glmmとかとか。 - 知行合一日記

というか，正規分布だと個体ごとにデータが１つしかないと，ランダム効果を推測できないのか！？

まずポワソン分布の場合。

> a<-0.25
> b<-rnorm(50,0,2)　#ランダム誤差
> y<-rpois(50,exp(a+b))
> id<-c(1:50)
> dat<-data.frame(y,id)

> fglm<-glm(y~1,poisson,dat) #普通にglmでやると
> fglm

Call:  glm(formula = y ~ 1, family = poisson, data = dat) 

Coefficients:
(Intercept)  
      1.834  

Degrees of Freedom: 49 Total (i.e. Null);  49 Residual
Null Deviance:	    787.8 
Residual Deviance: 787.8 	AIC: 893.2

推定がずれている。
では，次にglmmMLでやってみる。

> library(glmmML)
> fglmmML<-glmmML(y~1,cluster=id,data=dat,poisson)
> fglmmML

Call:  glmmML(formula = y ~ 1, family = poisson, data = dat, cluster = id) 


              coef se(coef)      z Pr(>|z|)
(Intercept) 0.1722   0.3587 0.4799    0.631

Scale parameter in mixing distribution:  1.954 gaussian 
Std. Error:                              0.3087 

Residual deviance: 150.4 on 48 degrees of freedom 	AIC: 154.4

結構いい感じに推定されている。　aが0.17（ホントは0.25），bが1.954（ホントは2）。
で，lmer()だと

> flmer<-lmer(y~1+(1|id),family=poisson,dat,REML=F)　　#REMLは正規分布にしかきかないみたいです
> flmer
Generalized linear mixed model fit by the Laplace approximation 
Formula: y ~ 1 + (1 | id) 
   Data: dat 
   AIC   BIC logLik deviance
 154.4 158.2 -75.21    150.4
Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 3.8191   1.9543  
Number of obs: 50, groups: id, 50

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.1721     0.3145  0.5472    0.584

んー，~~若干ずれています。とくにbが。~~
ウソです。glmmMLと同じような推定値です。
lmerはσがstd.dev.みたいです。
と，familyをポワソンにした時はランダム効果をidにして，各idは１つしかないけど，σ^2を推定できてるわけです。

ところが，正規分布データの場合は，idが１つだとできない。

> y2<-rnorm(50,a+b,3)
> dat2<-data.frame(y2,id)
> flmer2<-lmer(y2~1+(1|id),dat2,family=gaussian,REML=F)
 エラー： length(levels(dm$flist1)) < length(Y) is not TRUE

と，推定できない。

そもそもσ^2は，そのランダム効果内のバラツキから推定してるから，正規分布の場合に推定できないのはわかるけど，ポワソンや二項分布の場合はなぜできるのだろう。やはり，分布のパラメータが１つ（ポワソンλや二項p）と２つ（正規μとσ）で決まっていることの違いなのでしょうか。わからん。

ちなみに，idごとにいくつかデータがあれば，（前もやったけど）

> b<-rnorm(10,0,1) #ランダム効果
> id2<-rep(1:10,each=5) #idごとに同じランダム効果を持つ。
> y3<-rnorm(50,a+b[id2],3)
> dat3<-data.frame(y3,id2)


> fglm3<-glm(y3~1,gaussian,dat3)
> fglm3

Call:  glm(formula = y3 ~ 1, family = gaussian, data = dat3) 

Coefficients:
(Intercept)  
    0.04821  

Degrees of Freedom: 49 Total (i.e. Null);  49 Residual
Null Deviance:	    432.2 
Residual Deviance: 432.2 	AIC: 253.7

次にlmer()。

> flmer3<-lmer(y3~1+(1|id2),dat3,REML=F)
> flmer3
Linear mixed model fit by maximum likelihood 
Formula: y3 ~ 1 + (1 | id2) 
   Data: dat3 
   AIC   BIC logLik deviance REMLdev
 253.7 259.4 -123.9    247.7   247.1
Random effects:
 Groups   Name        Variance Std.Dev.
 id2      (Intercept) 1.3854   1.1770  
 Residual             7.2586   2.6942  
Number of obs: 50, groups: id2, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.04821    0.53262 0.09051

そりゃできますよ。できますけど，切片がよくないなぁ。

> flmer4<-lmer(y3~1+(1|id2),dat3,family=gaussian)
> flmer4
Linear mixed model fit by REML 
Formula: y3 ~ 1 + (1 | id2) 
   Data: dat3 
   AIC   BIC logLik deviance REMLdev
 253.1 258.8 -123.5    247.7   247.1
Random effects:
 Groups   Name        Variance Std.Dev.
 id2      (Intercept) 1.7006   1.3041  
 Residual             7.2586   2.6942  
Number of obs: 50, groups: id2, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.04821    0.56142 0.08586

はて，バラツキεの方を大きくしすぎたかな。