確認３

統計 R

こちらの続き。
random effectをどのように評価するか？

> leaf<-rep(1:4,each=4)
> disc<-rep(1:4,4)
> conc<-runif(16)*runif(16)*sqrt(leaf)
> data<-data.frame(leaf=factor(leaf),disc=factor(disc),conc)
> summary(data)
 leaf  disc       conc         
 1:4   1:4   Min.   :0.007103  
 2:4   2:4   1st Qu.:0.084563  
 3:4   3:4   Median :0.280326  
 4:4   4:4   Mean   :0.420890  
             3rd Qu.:0.580985  
             Max.   :1.523141

まずは尤度比検定。χ^2分布で近似する方法。

> mixmodel<-lmer(conc~1+(1|leaf),data)
> mixmodel
Linear mixed model fit by REML 
Formula: conc ~ 1 + (1 | leaf) 
   Data: data 
   AIC   BIC logLik deviance REMLdev
 15.14 17.46 -4.569    7.897   9.138
Random effects:
 Groups   Name        Variance Std.Dev.
 leaf     (Intercept) 0.176978 0.42069 
 Residual             0.052425 0.22897 
Number of obs: 16, groups: leaf, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)   0.4209     0.2180   1.931

バラツキが葉で大きいそうなことがつかめる。ML法でやってみると，

> mixMLmodel<-lmer(conc~1+(1|leaf),data,REML=F)
> mixMLmodel
Linear mixed model fit by maximum likelihood 
Formula: conc ~ 1 + (1 | leaf) 
   Data: data 
   AIC   BIC logLik deviance REMLdev
 13.78 16.10 -3.889    7.779   9.243
Random effects:
 Groups   Name        Variance Std.Dev.
 leaf     (Intercept) 0.129457 0.35980 
 Residual             0.052425 0.22897 
Number of obs: 16, groups: leaf, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)   0.4209     0.1888   2.229

となり，分散成分が多少低くなっている。
そして，null modelを作る。

> nullmodel<-lm(conc~1,data)
> as.numeric(2*(logLik(mixMLmodel)-logLik(nullmodel)))
[1] 10.35676 　　　#対数尤度比統計量。Devianceの差でもある。
> pchisq(10.35676,1,lower=F)　#それが自由度1のχ2分布に従うことを利用してp値を算出。
[1] 0.001290014

で，パラメトリックブートストラップだと

> pbst<-numeric(3000)
> for(i in 1:3000){
　　　　y<-unlist(simulate(nullmodel))
　　　　bnull<-lm(y~1)
　　　　bmodel<-lmer(y~1+(1|leaf),data,REML=F)
　　　　pbst[i]<-as.numeric(2*(logLik(bmodel)-logLik(bnull)))
　　　　　　　　}
> hist(pbst)

で尤度比統計量をシミュレーションしてくれる。
p値は，さきほどの尤度比統計量がどの程度の確率で出てくる値かで求めることができる。

> mean(pbst>10.3567)
[1] 0.0003333333

とかなり低いp値になる。そして標準誤差は

> sqrt(0.0003333*(1-0.0003333)/3000)
[1] 0.0003332611

とかなり低い。