3次のクロス表のログリニア分析が、さっとスタイリッシュ(笑)にできたらいいなあと思いつつ色々試行錯誤するもうまくいかず。とりあえず眠いので作りかけの部分をメモしておく。
# 三つの変数 A, B, & C をもとにクロス表 ### x <- table(A , B , C) ### モデル総当たり式でフィットさせる ### # (これがもっと数行でできないものか...) M <- NULL M[[1]] <- loglin(x, list(c(1,2,3)), param=T) M[[2]] <- loglin(x, list(c(1, 2), c(1, 3), c(2, 3)), param=T) M[[3]] <- loglin(x, list(c(1, 2), c(1, 3)), param=T) M[[4]] <- loglin(x, list(c(1, 2), c(2, 3)), param=T) M[[5]] <- loglin(x, list(c(1, 3), c(2, 3)), param=T) M[[6]] <- loglin(x, list(c(1, 2), 3), param=T) M[[7]] <- loglin(x, list(c(2, 3), 1), param=T) M[[8]] <- loglin(x, list(c(1, 3), 2), param=T) M[[9]] <- loglin(x, list(1, 2, 3), param=T) # 各モデルを図示 lab.llm <- c("[123]" ,"[12][13][23]" ,"[12][13]","[12][23]","[13][23]" ,"[12][3]" ,"[23][1]","[13][2]","[1][2][3]") names(M) <- lab.llm ### 結果の出力 ### # p値とAICを算出します hoge.p <- NULL; hoge.aic <- NULL for(i in 1:9){ hoge.p[i] <- 1-pchisq(M[[i]]$lrt, M[[i]]$df) hoge.aic[i] <- M[[i]]$pearson-2*M[[i]]$df } res <- cbind( round(hoge.p,4), hoge.aic) rownames(res) <- lab.llm colnames(res) <- c("p","AIC") res
こちらのページを参考にさせていただきました。深謝!→Rで対数線形モデル - きわめて個人的なR言語のメモ - R言語ユーザーグループ
12月10日追記、無意味にガラガラポンの関数化
myloglin <- function(x){ M <- NULL M[[1]] <- loglin(x, list(c(1,2,3)), param=T) M[[2]] <- loglin(x, list(c(1, 2), c(1, 3), c(2, 3)), param=T) M[[3]] <- loglin(x, list(c(1, 2), c(1, 3)), param=T) M[[4]] <- loglin(x, list(c(1, 2), c(2, 3)), param=T) M[[5]] <- loglin(x, list(c(1, 3), c(2, 3)), param=T) M[[6]] <- loglin(x, list(c(1, 2), 3), param=T) M[[7]] <- loglin(x, list(c(2, 3), 1), param=T) M[[8]] <- loglin(x, list(c(1, 3), 2), param=T) M[[9]] <- loglin(x, list(1, 2, 3), param=T) hoge.G2 <- hoge.df <- hoge.p <- hoge.aic <- NULL for(i in 1:9){ hoge.G2[i] <- M[[i]]$pearson hoge.df[i] <- M[[i]]$df hoge.p[i] <- 1-pchisq(M[[i]]$lrt, M[[i]]$df) hoge.aic[i] <- M[[i]]$pearson-2*M[[i]]$df } hoge <- cbind(hoge.G2, hoge.df, hoge.p, hoge.aic) colnames(hoge) <- c("G^2","df","P","AIC") rownames(hoge) <- c("{123}","{12} {13} {23}" ,"{12} {13}","{12} {23}","{13} {23}" ,"{12} {3}","{23} {1}","{13} {2}" ,"{1} {2} {3}") return(hoge) }