共分散構造分析[R編]につまづき、CATDAP1に逃避

SEMパッケージは収束条件が厳しい感じで、完全情報最尤推定法でつまずいたのをきっかけに、このブログもまたまた随分放置していました()。

　全く話が変わるのですが、随分以前、統計数理研究所のセミナーで坂本先生のcatdapの話を聞きました。先生は、データが集まったら、とりあえずcatdapで変数選択してみるみたいな紹介をされたように記憶しています。便利かもと思って、作ったのが下の関数。実は、香川大学経済学部教授でいらした故堀啓三先生がSPSSマクロで作成されていたcatdapを下敷きに、ほぼ盗作？
　いまや本家本元の統計数理研究所でcatdap1とcatdap2のスクリプトが公開されている(ここ)ので、下にあるような怪しげなスクリプトを使うメリットはあるはずないよね、でへっ。いや、絶対に使ってはいけないぞ。
　実は、人に見せる集計表を選択するとき、よく、このcatdapを使っているのだけれど、自分はAICを説明しきれないのでcatdapは日陰者あつかいです。変数選択は頭の中でひねり出したことになってます。

下記スクリプトはcatdap01というからには、質的変数にしか使えないのです。

catdap01(y.data, x.data, zero.cell = exp(-1), nvar = 2, printout = 10)

• x.data : 従属変数ベクトル
• y.data : 独立変数候補のデータフレーム
• nvar : 考慮する独立変数候補の数(4層まで)
• zero.cell : 0になったセルに入れる値
• printout : コンソールに出力するデータの数

最近のPCはパワフルなので、変数候補が10くらいあっても3層ぐらいなら、茶を飲んで待てるくらい(2014年現在)で終了すると思います。

catdap01(Titanic[,4], Titanic[,-4], nvar = 1)だと、統計数理研究所のcatdapのヘルプにあるcatdap1と同じになり、
catdap01( (Titanic[,4], Titanic[,-4], nvar = 3)だと、統計数理研究所のcatdap2を使って、
catdap2(Titanic[,4], 2, Titanic[,-4])と同じ値が出力されるはず。

nvar=1　のとき、表示に問題があるので修正(2014/03/20)

catdap01 <- function (y.data, x.data, zero.cell = exp(-1), nvar = 2, printout = 10)
{
    if (!is.factor(y.data)) {
        y.data <- as.factor(y.data)
    }
    printout <- trunc(printout)
    if (printout < 1) {
        printout <- 1
    }
    nvar <- min(nvar, 4)
    nexp <- ncol(x.data)
    aicbox <- c(0)
    varbox <- list("independent")
    ct <- list()
    varlist <- NULL
    cat("\n\ncombining variables...\n")
    varlist <- combn(1:nexp, 1, simplify = FALSE)
    if (nvar > 1) {
        varlist <- c(varlist, combn(1:nexp, 2, simplify = FALSE))
    }
    if (nvar > 2) {
        varlist <- c(varlist, combn(1:nexp, 3, simplify = FALSE))
    }
    if (nvar > 3) {
        varlist <- c(varlist, combn(1:nexp, 4, simplify = FALSE))
    }
    cat("tabulating ", length(varlist), "combinations....\n")
    fun1 <- function(x, zero.cell) {
        x[x == 0] <- zero.cell
        return(x)
    }
    calcaic <- function(x) {
        if (is.vector(x))
            x <- as.matrix(x)
        rn <- ncol(x)
        cn <- nrow(x)
        csum <- colSums(x)
        rsum <- rowSums(x)
        total <- sum(x)
        aic.temp <- sum(x * log((total * x)/(rsum %*% t(csum))))
        aic.temp <- (-2) * aic.temp + 2 * (rn - 1) * (cn - 1)
        return(aic.temp)
    }
    rvarlist <- lapply(varlist, function(x) rev(x) + 1)
    rvarlist <- lapply(rvarlist, function(x) c(1, x))
    all.var <- data.frame(y.data, x.data)
    temp.table <- NULL
    temp.table <- lapply(rvarlist, function(x) data.frame(xtabs(~.,
        data = all.var[, x])))
    temp.table <- lapply(temp.table, function(x) matrix(x$Freq,
        length(x$Freq)/length(levels(y.data)), length(levels(y.data)),
        byrow = TRUE))
    temp.table <- lapply(temp.table, function(x) x[rowSums(x) >  0, ])
    temp.table <- lapply(temp.table, function(x) fun1(x, zero.cell))
    cat("calculating AICs...\n\n")
    aicbox <- sapply(temp.table, function(x) calcaic(x))
    aicorder <- order(aicbox)
    aicbox <- aicbox[aicorder]
    varlist <- varlist[aicorder]
    vars <- sapply(varlist, function(x) c(colnames(x.data)[rev(x)],
        "", "", "")[1:nvar])
    if(is.vector(vars){ vars <- as.matrix(vars) }else{ vars <- t(vars) }
    output.mat <- data.frame(round(aicbox, 3), vars)
    rownames(output.mat) <- NULL
    colnames(output.mat) <- c("AIC", paste("explanatory variables",
        1:nvar))
    cat("\n<<< catdap result>>>\n\n")
    print(output.mat[1:min(printout, 100), ])
    result <- list(output = output.mat, aics = aicbox, explvar = as.list(varlist))
    return(invisible(result))
}