［統計］［R］Rによる混合効果ロジスティック回帰（その１）

　Rで混合効果ロジスティック回帰というのは、いろいろなところで取り上げられていて、今さら感もある話題なのだが、細かいところではときどき「これは何」と思うこともある。以下、シンプルなものを使って考えてみる。
　まず、サンプルデータだが、以下のようにmydata11という名前をつけたデータフレームに入れておく。容器の中に９個体か10個体をそれぞれ入れておいて、x1で表される条件に置いておいたら、何個体が生存（alive）で何個体が死亡（dead）なのかというデータである。よくあるシンプルな状況だと思うが、データ自体は計算例としての仮想的なものである（生存か死亡のどちらかだけで半死半生とかはないことにする）。x1は連続量である。

> mydata11
ID alive dead x1
1 1 2 8 1
2 2 4 6 2
3 3 6 4 3
4 4 8 2 4
5 5 4 5 3

（CRANのパッケージだけを考えても）混合効果ロジスティック回帰はたとえば、glmmMLパッケージの関数glmmMLでできる。少し前にデータの取り扱いに変更があった（久保さんのサイトなど参照）。以下のように分析してみる、
> glmm.res11<-glmmML(cbind(alive,dead)~1+x1,family=binomial,cluster=ID,data=mydata11)
> glmm.res11

Call: glmmML(formula = cbind(alive, dead) ~ 1 + x1, family = binomial, data = mydata11, cluster = ID)

coef se(coef) z Pr(>|z|)
(Intercept) -2.2550 0.9271 -2.432 0.0150
x1 0.8466 0.3296 2.569 0.0102

Standard deviation in mixing distribution: 2.947e-08
Std. Error: 0.3394

Residual deviance: 0.7813 on 2 degrees of freedom AIC: 6.781

こんな風に書いてもできる、
> glmm.res11x<-glmmML(cbind(alive,dead)~x1,family=binomial,cluster=ID,data=mydata11)
> glmm.res11x

Call: glmmML(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata11, cluster = ID)

coef se(coef) z Pr(>|z|)
(Intercept) -2.2550 0.9271 -2.432 0.0150
x1 0.8466 0.3296 2.569 0.0102

Standard deviation in mixing distribution: 2.947e-08
Std. Error: 0.3394

Residual deviance: 0.7813 on 2 degrees of freedom AIC: 6.781

もちろん、結果はまったく同じである。さて、ぱっと見て、あれれと思うのは（私が思ったということだが）AICの値である。AICは、対数尤度の-２倍にパラメーター数を加えたものだから、6.781ということは対数尤度は-0.39くらいとなる。各容器ごとの生存確率をパラメーターだと考えても、普通に考えると一番対数尤度が大きいのは、それぞれの容器の生存率が生存確率に等しいときのはずである。そのときの対数尤度として-0.39みたいなこんな大きな値になるのだろうかと感じるわけだが、計算してみるとそのときの対数尤度は-6.51くらいである（二項係数の部分も入れた値、入れなければもっと小さくなる）。-0.39はそれよりもだいぶ大きいので、そんな値をとるのは変だ、というのが、問題の１つの発端である。