ページ

2016年3月7日月曜日

【R】分割表の自動カテゴリー統合プログラム : optrecod (その5)

これまでのあらすじ:比較するべきものが2択ではなく、4択であることをはっきりさせた。カテゴリー統合に伴ってAICがどう変化するかを計算した。

\[
\begin{array}{c c c} {\cal M}_{1} & {\tiny ind} \longrightarrow & {\cal M}_{0} \\ \downarrow {\tiny mcat} & & \downarrow {\tiny mcat} \\ {\cal M}_{1}^{\prime} & {\tiny ind} \longrightarrow & {\cal M}_{0}^{\prime} \end{array}
\]
今回はカテゴリー統合を計算するR関数を二種類、開発する。

今回は、k行とl行を統合した分割表を返す関数mcatRows(x,k,l)と、統合前後のAIC変化を計算するdmcatRows(x,k,l)を開発してみる。列についても同様に出来る。

Rソースコード

まずはmcatRows関数から。

# -----------------------------
# mcatRows:分割表の行を統合する
# -----------------------------
mcatRows<-function(x,k,l){
 # チェックする
 nr <- nrow(x)
 if(nr > 2 && k <= nr && l <= nr && k != l){
  x[k,] <- x[k,] + x[l,]
  z <- x[-l,]
  return(z)
 } else {
  return(NULL)
 }
}
例題6.1を実行するとこうなる。
> ex6.1
     [,1] [,2]
[1,]  749   83
[2,]  445  636
> mcatRows(ex6.1,1,2)
NULL
2行以下なので、計算しないでNULLを返している。

実習データについても実行してみる。
> t11_14D
   
      1   2   3   4   5   6   7   8   9
  1 119 299  31  35   9   1   7  15   7
  2  23  96   2   0   2   1   2   2   2
  3  17  76   2   6   2   0   2   4   0
  4   3  20   2   0   1   0   0   1   2
  5   9  31   8   3   1   0   1   4   2
  6   1   0   0   0   0   1   0   0   0
> mcatRows(t11_14D,3,4)
   
      1   2  3  4 5 6 7  8 9
  1 119 299 31 35 9 1 7 15 7
  2  23  96  2  0 2 1 2  2 2
  3  20  96  4  6 3 0 2  5 2
  5   9  31  8  3 1 0 1  4 2
  6   1   0  0  0 0 1 0  0 0
ちゃんと動作しているようだ。

次にdmcatRows関数。
# -----------------------------
# dmcatRows:分割表統合によるAIC変化
# -----------------------------
dmcatRows<-function(x,k,l){
 nr <- nrow(x)
 if(nr > 2 && k <= nr && l <= nr && k != l){
  z <- x[c(k,l),]
  a <- sum(lchoose(z[1,] + z[2,],z[1,]))
  z[z==0] <- exp(-1)
  b <- sum(z * log(z)) - sum((z[1,] + z[2,])*log(z[1,] + z[2,]))
  retval <- 2 *(a + b) - 2 * ncol(x)
  return(retval)
 } else {
  return(NULL)
 }
}
肝はlchoose関数という組み込み関数を使っていることと、対数計算の直前にゼロ度数対策を施していること。

実行するとこうなる。
> dmcatRows(ex6.1,1,2)
NULL
> -catdap(t11_14D)
[1] -3.886714
> dmcatRows(t11_14D,3,4)
[1] -38.73089
> -catdap(mcatRows(t11_14D,3,4))
[1] 0.8423838
実習データについて結果を解釈すると、 \[ \begin{array}{c c c} {\cal M}_{1} & {\tiny ind} \longrightarrow & {\cal M}_{0} \\ & -3.89 & \\ \downarrow {\tiny mcat} & & \downarrow {\tiny mcat} \\ -38.73 & & \\ {\cal M}_{1}^{\prime} & {\tiny ind} \longrightarrow & {\cal M}_{0}^{\prime} \\ & +0.84 & \end{array} \] つまり、飽和モデルM1からそのまま独立性の検定を行うと(M1 → M0)帰無仮説が棄却されない可能性が高いが、行を統合後のモデルM1'では独立性の帰無仮説はぎりぎり棄却できるかどうかといったところ(少しマシになった)。AICの水準で言えばM1'が最も最適。

次回の予告

次は、列についても同様の関数を完成させた上で、自動カテゴリー統合の基本方針を確認する。