« [香椎・千早]香椎浜団地入口交差点ー洋服の青山 | トップページ | [本]『日本のソフトウェア産業がいつまでもダメな理由』 »

2008.04.16

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

 他のことに(私が)使われているあいだに、East_Scrofaさんや足軽日記さんのところで、だいぶ進んでいる。

 また、余談であるが、さきに書いた2つのデータの与え方を、glm()だけでなくglmmML()にも試してみる。結果はここまでこの話題で書いてきたことから容易に想像されるので、「glm()のAICとglmmML()のAICはこんなときだったら比べられるんだろ」とすぐ分かってしまう方は読み飛ばしてください。使ったのはこんなデータである。
まず、1行が1cluster
> mydata1g
alive dead x1 total ID
1 4 6 1 10 1
2 7 3 2 10 2
3 8 2 4 10 3
4 3 4 0 7 4

次に内容的には同じデータの1行1個体版である。
> mydata1
alivedead x1 ID
1 1 1 1
2 1 1 1
3 1 1 1
4 1 1 1
5 0 1 1
6 0 1 1
7 0 1 1
8 0 1 1
9 0 1 1
10 0 1 1
11 1 2 2
(中略)
30 0 4 3
31 1 0 4
32 1 0 4
33 1 0 4
34 0 0 4
35 0 0 4
36 0 0 4
37 0 0 4
長いので少し略した。
> res.glm.i1<-glm(alivedead~x1,family=binomial,data=mydata1)
> res.glm.g1<-glm(cbind(alive,dead)~x1,family=binomial,data=mydata1g)
とまず、glm()で、次に
> res.glmmML.g1<-glmmML(cbind(alive,dead)~x1,family=binomial,cluster=ID,data=mydata1g)
> res.glmmML.i1<-glmmML(alivedead~x1,family=binomial,cluster=ID,data=mydata1)
とglmmML()でやってみる。

 以下の2つはglm()の結果で、先に書いたとおりで当たり前だが、点推定値や検定結果などまったく同じである。
> summary(res.glm.i1)

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

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8363 -1.1767 0.6401 0.9766 1.3939

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.88 0.3791
x1 0.4943 0.2657 1.86 0.0629 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 49.961 on 36 degrees of freedom
Residual deviance: 46.009 on 35 degrees of freedom
AIC: 50.009

Number of Fisher Scoring iterations: 4

> summary(res.glm.g1)

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

Deviance Residuals:
1 2 3 4
-0.6320 0.5252 -0.1188 0.2711

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.88 0.3791
x1 0.4943 0.2657 1.86 0.0629 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 4.71427 on 3 degrees of freedom
Residual deviance: 0.76285 on 2 degrees of freedom
AIC: 15.016

Number of Fisher Scoring iterations: 4

 以下の2つはglmmML()の結果で、こちらも、点推定値や検定結果などまったく同じである。
> summary(res.glmmML.i1)

Call: glmmML(formula = alivedead ~ x1, family = binomial, data = mydata1, cluster = ID)


coef se(coef) z Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.8795 0.3790
x1 0.4943 0.2657 1.8601 0.0629

Standard deviation in mixing distribution: 1.041e-08
Std. Error: 0.4001

Residual deviance: 46.01 on 34 degrees of freedom AIC: 52.01
> summary(res.glmmML.g1)

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


coef se(coef) z Pr(>|z|)
(Intercept) -0.4959 0.5638 -0.8795 0.3790
x1 0.4943 0.2657 1.8601 0.0629

Standard deviation in mixing distribution: 1.611e-09
Std. Error: 0.4001

Residual deviance: 0.7628 on 1 degrees of freedom AIC: 6.763

 さて、この例では、ほとんどランダム効果の項の効果はないのだが(したがってあまりわかりやすい例ではなくてすみません)、1行1個体という与え方のときのglm()とglmmML()が返すAICを見ると、これが同じ基準で計算されているらしいことがわかる。この例ではほとんどランダム効果項の効果がないので、1パラメーター多い分だけ、つまりAICなら2ちがってくるはずだが、実際にもそうなっている。
 glm()とglmmML()が返すAICは、1行1clusterというデータの与え方のときは基準値が違うのでそもそもそのままでは比較可能でなく、一方、1行1個体という与え方のときは同じ基準値で計算されているので比較しうるものになっている、ということのようだ。
 1行1clusterというデータの与え方のときは、glmmML()にclusterがないときのdevianceを計算させ(いわばglm()を代行させ)れば、同じ基準値によるAICはすぐ求められるわけである。

|

« [香椎・千早]香椎浜団地入口交差点ー洋服の青山 | トップページ | [本]『日本のソフトウェア産業がいつまでもダメな理由』 »

コメント

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

トラックバック


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

« [香椎・千早]香椎浜団地入口交差点ー洋服の青山 | トップページ | [本]『日本のソフトウェア産業がいつまでもダメな理由』 »