« [その他]NBAプレーオフ間近 | トップページ | [その他]NBA »

2008.04.08

[統計][R]Rによる混合効果ロジスティック回帰(その1)

 Rで混合効果ロジスティック回帰というのは、いろいろなところで取り上げられていて、今さら感もある話題なのだが、細かいところではときどき「これは何」と思うこともある。以下、シンプルなものを使って考えてみる。
 まず、サンプルデータだが、以下のようにmydata11という名前をつけたデータフレームに入れておく。容器の中に9個体か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は、対数尤度の-2倍にパラメーター数を加えたものだから、6.781ということは対数尤度は-0.39くらいとなる。各容器ごとの生存確率をパラメーターだと考えても、普通に考えると一番対数尤度が大きいのは、それぞれの容器の生存率が生存確率に等しいときのはずである。そのときの対数尤度として-0.39みたいなこんな大きな値になるのだろうかと感じるわけだが、計算してみるとそのときの対数尤度は-6.51くらいである(二項係数の部分も入れた値、入れなければもっと小さくなる)。-0.39はそれよりもだいぶ大きいので、そんな値をとるのは変だ、というのが、問題の1つの発端である。

|

« [その他]NBAプレーオフ間近 | トップページ | [その他]NBA »

コメント

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: [統計][R]Rによる混合効果ロジスティック回帰(その1):

« [その他]NBAプレーオフ間近 | トップページ | [その他]NBA »