4次のクロス表用ログリニア解析のためのR関数 - こにしき（言葉・日本社会・教育）

誰がこんな関数使うんだよのコーナー！

つくったのでアップします。ただし、モデル設定がかなり個人的な事情を反映してますので、ご使用の際はご注意下さい。

とりあえず以下をコピー＆ペースト

### 注意！！！
### 以下は、4次のクロス表で2番目、3番目、4番目の相互連関を最初から仮定してます

myloglin4 <- function(x){
 M <- NULL
 M[[1]] <- loglin(x, list(c(1,2,3,4)), param=T)
 M[[2]] <- loglin(x, list(c(1,2),c(1,3),c(1,4),c(2,3),c(3,4),c(2,4)), param=T)
 M[[3]] <- loglin(x, list(c(1,2),c(1,3)      ,c(2,3),c(3,4),c(2,4)), param=T)
 M[[4]] <- loglin(x, list(c(1,2),c(1,4)      ,c(2,3),c(3,4),c(2,4)), param=T)
 M[[5]] <- loglin(x, list(c(1,3),c(1,4)      ,c(2,3),c(3,4),c(2,4)), param=T)
 M[[6]] <- loglin(x, list(c(1,2)             ,c(2,3),c(3,4),c(2,4)), param=T)
 M[[7]] <- loglin(x, list(c(1,3)             ,c(2,3),c(3,4),c(2,4)), param=T)
 M[[8]] <- loglin(x, list(c(1,4)             ,c(2,3),c(3,4),c(2,4)), param=T)
 M[[9]] <- loglin(x, list(c(1)               ,c(2,3),c(3,4),c(2,4)), param=T)
hoge.p <- NULL; hoge.aic <- NULL
 for(i in 1:length(M)){
  hoge.p[i] <- 1-pchisq(M[[i]]$lrt, M[[i]]$df) 
  hoge.aic[i] <- M[[i]]$pearson-2*M[[i]]$df 
 }
 hoge <- cbind(hoge.p, hoge.aic)
 colnames(hoge) <- c("P","AIC")
 rownames(hoge) <- c("{1 2 3 4}","{1 2}{1 3}{1 4}..."
    ,"{1 2}{1 3}...","{1 2}{1 4}...","{1 3}{1 4}..."
    ,"{1 2}...","{1 3}...","{1 4}...","...")
return(hoge)
}

分析事例

たとえば、ある自治体の住民調査（サンプリング調査）をやったとします。
その結果、得られた人口分布表が以下のようなものだったとします。

> ### 架空のデータ
> a <- c(15, 63, 36, 17, 76, 36, 38, 65, 23, 43, 88, 31, 21
+     ,99, 52, 31, 132, 67, 46, 78, 35, 52, 97, 38, 46, 164
+     ,84, 67, 197, 109, 29, 67, 21, 15, 36, 13)
> 
> X <- array(a,dim=c(3,2,2,3))
> dimnames(X)[[1]] <- c("都市部","住宅地","農村部")
> dimnames(X)[[2]] <- c("男性","女性")
> dimnames(X)[[3]] <- c("非・大卒","大卒") 
> dimnames(X)[[4]] <- c("若年","中年","高齢") 
> 
> # 見づらいので、 ftable() で4次のクロス表を2次につぶす
> (X0 <- ftable(X,row.vars=c(2,3,4)) )   
                    都市部 住宅地 農村部
                                        
男性 非・大卒 若年      15     63     36
              中年      21     99     52
              高齢      46    164     84
     大卒     若年      38     65     23
              中年      46     78     35
              高齢      29     67     21
女性 非・大卒 若年      17     76     36
              中年      31    132     67
              高齢      67    197    109
     大卒     若年      43     88     31
              中年      52     97     38
              高齢      15     36     13

このとき、(1) どういう地域に住んでいるか、(2) 性別はなにか、(3) 最終学歴は大卒か否か、(4) 年齢は若いか中くらいかお年寄りか、という4つの変数をもとに集計されていたとします。

理論上、うえの4つの変数にはそれぞれ相互の連関があるかもしれません。

こんなかんじ

しかし、これだとすべての矢印に「はてな」（？）がついていて煩雑です。
しかも、たとえば、常識的に考えて、世代と学歴や、ジェンダーと学歴は関係があると考えられます。大学進学率は戦後全体でみれば徐々に増加しているので、若いほうが大卒が多いですし、（最近やや緩和したとはいえ）男女の高等教育への進学格差は大きいものがあります。

ということで、ジェンダー、世代、学歴には最初から連関を想定しちゃいます。

とうことで、「居住地域」とそれ以外の変数同士の関連があるかないかだけを見ていけばいいことになります。
こうするとかなりすっきり。
以下のように、
連関パタン（モデル）も8つくらいに減ります。（なお、準拠モデルとして「飽和モデル」というものも用意するので、上の自作関数では９個のモデルが比較されています）

さっきのデータを再確認

分析の前に、居住者の分布が、ジェンダー、学歴、世代でどう変動があるかを簡単にチェックしておきます。

### 居住地域の分布（パーセンテージ表示で基本的な分布を確認）
P0 <- X0/apply(X0,1,sum)            # 各属性ごとの居住地シェア
rownames(P0) <-  paste(
           rep(c(dimnames(X)[[2]]),times=1,each=6)
          ,rep(c(dimnames(X)[[3]]),times=2,each=3)
          ,rep(c(dimnames(X)[[4]]),times=4,each=1)
          ,sep="／")

par(mar=c(10,5,2,2))
 barplot(100*t(P0),las=2,ylab="シェア（％）",col=2:4
   ,legend=dimnames(X)[[1]])

こんな感じ。学歴と居住地の間に差がありそうです。ただ、ジェンダーによっても多少の差があるようにも見えます。

というわけで、先ほどのログリニア関数

> res <- myloglin4(X)
2 iterations: deviation 0 
5 iterations: deviation 0.09124073 
5 iterations: deviation 0.01192215 
4 iterations: deviation 0.07844824 
5 iterations: deviation 0.07874301 
4 iterations: deviation 0.07796009 
4 iterations: deviation 0.07576277 
4 iterations: deviation 0.07796009 
4 iterations: deviation 0.07796009 
>
> round(res,3)
                       P     AIC
{1 2 3 4}          0.000   0.000
{1 2}{1 3}{1 4}... 0.076  -7.469
{1 2}{1 3}...      0.175 -14.313
{1 2}{1 4}...      0.000  43.725
{1 3}{1 4}...      0.131 -11.278
{1 2}...           0.000  38.752
{1 3}...           0.260 -18.162
{1 4}...           0.000  39.754
...                0.000  34.771

一般的に、P値（危険率）が極度に小さくなくて（棄却されなくて）、かつ、AICが十分小さいものが良いモデルだと言われています。

なので、上記の結果では、(1) 居住地と(3) 学歴だけに連関を想定するモデル