« [統計][R]Rによる混合効果ロジスティック回帰-数値的な積分 | トップページ | [その他]NBAプレーオフあと1つずつ »

2008.04.13

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

 ここ何回か、Rによる混合効果ロジスティック回帰について書いてきた(といっても出てくるのは、glm()とglmmML()だが)。この2つについても、AICや対数尤度、あるいはdevianceの計算という、いままで書いてきたことと似た点についていくつか補足しておく。
 ロジスティック回帰では、ソフトウェアへのデータの与え方にいくつかのやり方がある。ここのところ使っていた、各容器の中にいる個体たちの生存/死亡なら、容器ごとに、生存個体数と死亡個体数を与えるやり方もあれば、1個体ごとに生存/死亡(そして容器の番号なども)を与えるやり方もある。(その1)から(その4)で使っていたのは前者である。

 Rのglm()という関数はかなりよく使われているので、当然の知識かもしれないが、上記の2つのやり方で、内容的には同じデータを与えたとき、glm()の返す結果は表面上は違っているところがある。とても簡単な例でやってみる。
まず使うデータだが、以下のmydata1gである。これは”前者”にあたる。
> mydata1g
alive dead x1
1 4 6 1
2 7 3 2
3 8 2 4

もう1つはmydata1gを個体ごとに直した以下のmydata1である。glm()で二項分布という仮定で行うなら(つまりoverdispersionもanderdispersionもないなら)、内容的には同じである。
> mydata1
alivedead x1
1 1 1
2 1 1
3 1 1
4 1 1
5 0 1
6 0 1
7 0 1
8 0 1
(長いので途中は省略)
27 1 4
28 1 4
29 0 4
30 0 4
以下のようにそれぞれロジスティック回帰する
> res1g.glm<-glm(cbind(alive,dead)~x1,family=binomial,data=mydata1g)
> res1.glm<-glm(alivedead~x1,family=binomial,data=mydata1)

> summary(res1g.glm)

Call:
glm(formula = cbind(alive, dead) ~ x1, family = binomial, data = mydata1g)

Deviance Residuals:
1 2 3
-0.3984 0.6178 -0.2580

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7318 0.8246 -0.887 0.3748
x1 0.5816 0.3518 1.653 0.0983 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 3.74390 on 2 degrees of freedom
Residual deviance: 0.60697 on 1 degrees of freedom
AIC: 12.41

Number of Fisher Scoring iterations: 4

> summary(res1.glm)

Call:
glm(formula = alivedead ~ x1, family = binomial, data = mydata1)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.886 -1.114 0.608 1.001 1.242

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7318 0.8245 -0.888 0.3748
x1 0.5816 0.3518 1.653 0.0983 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 36.293 on 28 degrees of freedom
AIC: 40.293

Number of Fisher Scoring iterations: 3

推定値や検定結果などは当然のことだが同じなのだが、AICやdevianceなどが異なっている。devianceもNullとResidualの差はどちらでも同じである。
 まず、AICのちがいが計算されている対数尤度の違いのみによることは、
> logLik(res1g.glm)
'log Lik.' -4.205005 (df=2)
> logLik(res1.glm)
'log Lik.' -18.14627 (df=2)
としてみると、対数尤度の差がAICのちがいのちょうど半分であることから確かめられる。

 それぞれの対数尤度は以下のように計算されている。”容器”ごとの場合の、ロジスティック回帰による生存確率の予測値は
> res1g.glm$fitted
1 2 3
0.4625256 0.6062115 0.8312628
である。そこで、対数尤度のうち、生存数×log(生存確率)+死亡数×log(1-生存確率)の部分は
> sum(mydata1g$alive*log(res1g.glm$fitted)+mydata1g$dead*log(1-res1g.glm$fitted))
[1] -18.14627
となる。これに二項係数の部分、
> sum(lchoose(mydata1g$alive+mydata1g$dead,mydata1g$alive))
[1] 13.94126
を加えたのが、 -4.205である。個体ごとというときには、二項係数の部分が入っていないわけである。

 ”容器”ごとの場合の飽和モデルは、各”容器”にそれぞれ尤度がもっとも大きくなる生存確率を割り当てたもので、alive/(alive +dead)である。p0がその尤度がもっとも大きくなる生存確率で、以下のようになる。

alive dead x1 total p0
1 4 6 1 10 0.4
2 7 3 2 10 0.7
3 8 2 4 10 0.8

生存数×log(生存確率)+死亡数×log(1-生存確率)の部分は-17.84278となり、二項係数の部分(上に書いたのと同じく 13.94126)を加えると、-3.9015となる。この飽和モデルの対数尤度とロジスティック回帰の対数尤度(上記の-4.205)のちがいの2倍が、Residual devianceになっている。

 ”個体”ごとのときには、それぞれの個体の結果を完全に予測してしまうのが飽和モデルになるから、尤度は1で対数尤度は0になる。そこで、-18.14627という対数尤度の2倍がResidual devianceになっている。これは、glmmML()のときの計算の仕方そっくりである(glmmML()のデータの与え方が拡張されてきた経過を考えると、このように計算しているのがわかる気がする)。


 

|

« [統計][R]Rによる混合効果ロジスティック回帰-数値的な積分 | トップページ | [その他]NBAプレーオフあと1つずつ »

コメント

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

トラックバック


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

« [統計][R]Rによる混合効果ロジスティック回帰-数値的な積分 | トップページ | [その他]NBAプレーオフあと1つずつ »