ログリニア解析、試行錯誤

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言語ユーザーグループ

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)
}