ページ

2016年3月10日木曜日

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

これまでのあらすじ:行を統合する関数、統合前後のAIC変化を計算する関数を定義し、うまく動作することを確認した。

今回は、列の統合関数を作成(といっても、コピペで済むけれど)、全体の自動カテゴリー統合プログラムoptrecodを完成させようと思う。このプログラムの挙動はこんな感じ。
  1. 分割表(飽和モデルそのもの)を引数として受け取る。また、オプションとして繰り返し計算をするか、一段階のみの計算とするかの論理値も受け取る。今回は試作品として、オプションはFALSE固定。
  2. 行の組み合わせのうち、AICのマイナスが最小となる組み合わせを探索する。もしもAICがマイナスにならない、または組み合わせ候補が存在しない場合は、スキップする。
  3. 列について同様。
  4. 行、列のいずれがAICが小さくなるか判定。最小の方で分割表をカテゴリー統合する。行、列のいずれもスキップされていた場合は何もしない。
  5. (繰り返しオプションがTRUEの場合、1に戻る。)FALSEの場合、結果を出力して、終了する。

列統合関数
まずは列統合についての関数定義。

# -----------------------------
# mcatCols:分割表の列を統合する
# -----------------------------
mcatCols<-function(x,k,l){
 # チェックする
 nc <- ncol(x)
 if(nc > 2 && k <= nc && l <= nc && k != l){
  x[,k] <- x[,k] + x[,l]
  z <- x[,-l]
  return(z)
 } else {
  return(NULL)
 }
}

# -----------------------------
# dmcatCols:分割表統合によるAIC変化
# -----------------------------
dmcatCols<-function(x,k,l){
 nc <- ncol(x)
 if(nc > 2 && k <= nc && l <= nc && 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 * nrow(x)
  return(retval)
 } else {
  return(NULL)
 }
}
動作は以下のとおり。
> mcatCols(t11_14D,3,4)
   
      1   2  3 5 6 7  8 9
  1 119 299 66 9 1 7 15 7
  2  23  96  2 2 1 2  2 2
  3  17  76  8 2 0 2  4 0
  4   3  20  2 1 0 0  1 2
  5   9  31 11 1 0 1  4 2
  6   1   0  0 0 1 0  0 0
> dmcatCols(t11_14D,3,4)
[1] -26.76791

カテゴリー自動統合関数
次に自動統合のための関数を定義してみた。あくまでも試作品として。 いろいろ計算してみたが、次のような問題があるようだ。

  1. 何も考えずに統合してゆくと、どんどん統合されてしまう。クラスター分析で最近隣法での挙動に似ており、一番度数の大きいカテゴリーがどんどん膨らんでゆく。これを避けるには、一番度数の少ないカテゴリー同士が統合されるようにしなければならない。
  2. どこに統合するかについては、「その他」カテゴリーがあるならば、そちらに優先して統合したい。
  3. 目的について、最小AICのモデルを探索することに主眼を置いていたが、それは分析目的からして適当ではないのではないか。あくまでも、独立化によってAICが悪化する度合いが一番大きいものが良いのではないか。

今後の予定
もう一度、分析目的の点から考え直してみたい。