lmer研究（５）&multcomp研究（２）lmerで多重比較

まず、TukeyのHSDなどの多重比較は、正規分布のモデルしか使えないはず。それについてmultcompに関しての論文のp7にgeneralized linear modelとmixed modelについて書かれてある。

読む限りだと、他の分布でも使っていいようなことが書いてある。

また、mixed modelの場合は、REMLならばよい。と書いてあるように見えるけれど、上の論文にlmerで二項分布を使った例がある（検定はしていない）ので、使えないことはないのかもしれない。

一応、Lhaplace近似して尤度の近似値を求めているので、使ってもいいんではないのかなと思う。これは、lmerでAIC水準のモデル選択ができることと同じ根拠。

実際に動かしてみる。

ポアソン分布のデータを作る。

set.seed(813)
count<- c(rpois(30,10) + floor(rnorm(30,0,3)), rpois(30,10) + floor(rnorm(30,0,3))
         ,rpois(30,14) + floor(rnorm(30,0,3)) )
treat<- factor(c(rep("A",30),rep("B",30),rep("C",30)))
site<- rep(c(rep("alpha",10),rep("beta",10),rep("gamma",10)),3)
dat<- data.frame(count,treat,site)

昨日のデータとほとんど変わりません。
AとBが平均１０で、Cが平均14
siteはブロックを想定したランダム変数（平均０分散３）。ブロックは３つ。

正規分布なデータではないので、ほんとはできないのだが、試しにmultcompBoxplot()を使ってみる。

library(multcompView)
multcompBoxplot(count~treat,data=dat)

Aが重なるデータができました。

ここからが本題

library(lme4)
mod.lmer<- lmer(count~ treat+(1|site), family=poisson)
summary(mod.lmer)

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  2.27899    0.05842   39.01   <2e-16 ***
treatB      -0.05267    0.08373   -0.63   0.5293   
treatC       0.18938    0.07897    2.40   0.0165 *

これを、Tukeyの方法で多重比較する。

library(multcomp)
tukey1<- glht(mod.lmer,linfct=mcp(treat = "Tukey"))
summary(tukey1)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts
Fit: NULL

Linear Hypotheses:
           Estimate Std. Error z value p value   
B - A == 0 -0.05267    0.08373  -0.629  0.8040   
C - A == 0  0.18938    0.07897   2.398  0.0435 *
C - B == 0  0.24205    0.08013   3.021  0.0073 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Warning message:
In x$model$call : $ operator not defined for this S4 class, returning NULL

ポアソンモデルでやることにより、きれいにトリートメントA, BグループとCに分けることができた。

上で紹介したmultcompの論文の１８ページにlmer()の二項分布の例があり、glht()を用いて、予測値と９５％信頼区間を推定している。モデル選択の場合、やっていることは推定だと思うので、こっちの方が求めたいもの次第によってはいいのかもしれない。
今回の例でやってみる。

mod.lmer2<- lmer(count~ treat-1+(1|site), family=poisson) #切片なしモデルをつくる
K <- diag(length(fixef(mod.lmer2)))
ci2 <- confint(glht(mod.lmer2, linfct = K))
ci2$confint<- poisson()$linkinv(ci2$confint)
ci2


         Simultaneous Confidence Intervals

Fit: NULL

Estimated Quantile = 2.3878
95% family-wise confidence level

Linear Hypotheses:
       Estimate lwr     upr   
1 == 0  9.7524   8.4817 11.2134
2 == 0  9.2499   8.0147 10.6755
3 == 0 11.7992  10.3929 13.3959

こうして、予測値と信頼区間がもとまりました。
ここで、ひとつ注意。
lmerの場合、fitted()を使って予測値を求めることができる。しかし、glm()の場合fitted()で、リンク関数を外した状態で得られるのに対し、lmer()の場合、出てくる予測値はリンク関数が外されていない。そのため、リンク関数の逆関数（この場合はexp()）をしてあげる必要がある。その行程が、上のpoisson()$linkinv()の関数だ。これで、リンク関数を外した値を得ることができる。

これを図にする。

cci<- ci2$confint
plot(1:3, cci[,1], xaxt="n",ylab="",xlab="",ylim=c(min(cci[,2]),max(cci[,3])))
axis(1,1:3,c("A","B","C"))
arrows(1:3,cci[,2],1:3,cci[,3],length=.05,angle=90,code=3)
text(1:3,13,c("a","a","b"),col=c("blue","blue","red"),cex=2)

　

棒グラフと情報量は変わらないけれど、論文に棒グラフで平均と標準誤差を載せるよりましな気がする。

最後にひとつ注意を。

Tukeyの方法など、p値を求める検定がRのバージョン2.7.0でできなくなりました（今までの結果は2.6.2で求めた）。これは、lmerでやってはいけないからではなく、推測ですが2.7.0では、S4メソッドの何かしらの仕様が変わったからだと思います。
S4のlmer()にS3のglht()を無理やりなんとかしていた感が2.6.2からあった（warningがよく出る）ので、仕方ないのかと思います。
そこのところは、将来的にR側でなんとかするのか、multcompパッケージ側で対応するのだろうと思います。

過去のlmer研究
（１）自動akaike weightプログラム
（２）自動akaike weightプログラム２
（３）gaussian modelのp値
（４）AICについて

過去のmultcomp研究
（１）とりあえず正規分布でやってみる