2010年8月21日土曜日

自分用反復測定データにおけるGLMMまとめ

現時点ではlme4パッケージのlmer関数を用いるのが一般的と思われる。
測定日が一定間隔でなかったり、測定日によって繰り返し数が異なる(missing dataがある)ようなデータセットの場合、repeated measures ANOVAはあまり有用ではない。
そこでGLMMの登場。
繰り返し(rep)をランダム効果にすることで、missing dataがある場合でも手元にある全データを活用して統計できる。
固定効果として、測定日(day)、生物種(sp)の2つあったとする。
どんな固定構造のモデルが最も当てはまりが良いか探索することにする。

> attach(dataset)
> sp=factor(sp)
> rep=factor(rep)
> library(lme4)

(1)最適モデルの探索
> model0 <- lmer(Y ~ 1 + (1|rep), data=dataset, REML=FALSE)
#固定効果なしモデル
> model1 <- lmer(Y ~ sp*day+(1|rep), data=dataset, REML=FALSE)
#フルモデル
> model2 <- lmer(Y ~ sp+day+(1|rep), data=dataset, REML=FALSE)
#交互作用なしモデル
> model3 <- lmer(Y ~ sp+(1|rep), data=dataset, REML=FALSE)
#種のみモデル
> model4 <- lmer(Y ~ day+(1|rep), data=dataset, REML=FALSE)
#調査日のみモデル

REML=FALSEにしたのは、固定構造(fixed structure)の異なるモデルを比較する場合、REMLでなくMLで推定したモデルでなければならないためである。
lmerではFALSEにするとML、TRUEにするとREMLを使った推定となる。
(参考文献) Zuur et al. 2009 "Mixed Effects Models and Extensions in Ecology with R"
Pinheiro & Bates 2000 "Mixed Effects Models in S and S-Plus"

(1-1)パターン1 手動でAIC比較
> AIC(model0, model1, model2, model3, model4)

(1-2)パターン2 自動でAIC比較
フルモデル(model1)のみあればよい。
> library(MuMIn)
> dredge(model1)
これで固定効果の全組み合わせを自動的に計算して、AICやらAICcを一覧にしてくれる。

(1-3)パターン3 尤度比検定
> anova(model0, model1, model2, model3, model4)
ただし3つ以上のモデルを比較する場合、多重比較の問題が出てくるので、この後bonferroni補正などが必要。
具体的にどうやるかは、勉強中。
AICによるモデル選択に惹かれつつあるのでこのまま放置かも。

(2) 固定効果の多重比較: どれとどれが異なるのか?
種(sp)のみが固定効果として残り(model3)、どの種とどの種が異なるのか検定したい場合について。

(2-1)多重比較法
> library(multcomp)
> summary(glht(model3, linfct=mcp(sp="Tukey")))
これはTukey検定を行う場合。他の方法も使える。多分。

(2-2)多重比較的モデル選択(勝手に命名)
A, B, Cという3種の比較を行うとする。
データセットに列を追加し、ダミー変数を入れて、モデル選択を行う。
例えばAだけがB, Cと異なり、BとCは差がない、というモデルなら、種Aの行には「A」というダミー変数を、種BとCの行には「BC」というダミー変数を入れる。
同様に、A=B≠C、A≠B≠Cなど、あらゆるパターンの分だけ列を追加してダミー変数を入れていく。
rep sp day Y ABC A_BC AB_C AC_B A_B_C
1 A 1 Y ABC A AB AC A
2 A 1 Y ABC A AB AC A
~
1 A X Y ABC A AB AC A
2 A X Y ABC A AB AC A
~
1 B X Y ABC BC AB B B
2 B X Y ABC BC AB B B
~
1 C X Y ABC BC C AC C
2 C X Y ABC BC C AC C
~

このあと(1)と同様にモデル選択する。
> modelABC <- lmer(Y ~ ABC+(1|rep), data=dataset, REML=FALSE)
#3種で差がないモデル
> modelA_BC <- lmer(Y ~ A_BC+(1|rep), data=dataset, REML=FALSE)
#A≠B=Cのモデル
> modelAB_C <- lmer(Y ~ AB_C+(1|rep), data=dataset, REML=FALSE)
#A=B≠Cのモデル
> modelAC_B <- lmer(Y ~ AC_B+(1|rep), data=dataset, REML=FALSE)
#A=C≠Bのモデル
> modelA_B_C <- lmer(Y ~ A_B_C+(1|rep), data=dataset, REML=FALSE)
#A≠C≠Bのモデル

> AIC(modelABC, modelA_BC, modelAB_C, modelAC_B, modelA_B_C)

2010年8月20日金曜日

ML vs. REML

Mixed Model勉強中。
今トレンド(?)のlmerですが、デフォルトはREMLになってます。
しかしいろいろ調べてみると、固定効果部分が異なるモデルを比較する場合、REMLを使っちゃいかんようです。
REMLは固定構造を考慮せずランダム構造のみの尤度を最大化するのがその理由らしい。
というか固定構造が変わるとREMLの評価基準も変わってしまうようだ。
固定構造を考慮しない推定法で、固定構造の異なるモデルの比較を行うこと自体ナンセンス、ということだろうか。
ここの解釈がまだアヤシイ。

いろいろ見てたらTREEにこんな論文を見つけた。
Bolker et al. 2008 "Generalized linear mixed models: a practical guide for ecology and evolution"

この本もいい。ほすぃ。
Zuur et al. 2009 "Mixed Effects Models and Extensions in Ecology with R"