「統計学関連なんでもあり」の過去ログ--- 046

No.21042　dredgeとmodel.selのAICが異なる　　【gibaco 】　2014/05/19(Mon) 12:15

済みませんが，教えて下さい。
以下のような，オフセット項，ランダム効果を加えたGLMM（glmmadmb）をつくり，AIC基準でモデル選択をしようとしています。

ana <- glmmadmb(monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)) + (1 | site), data =d, family= "nbinom")

モデル選択をするためにdredgeとmodel.selをそれぞれ使って，モデル間の比較を行いました。ところが，この二つのやり方ではlogLik， dfの値が異なり，当然ながらAICも違っています。また，各モデルのAIC順も異なっていました。個々のモデルについて，summaryを出すと， model.selの結果とは一致し，dredgeとは一致しません。どうしてこんなことが起きるのか，どちらを使うべきか，教えていただければと思います。

　色々試してみたので以下に結果の例をのせておきます。うまくファイル添付できなかったので，長くなってしまいました。申し訳ありません。ここではmodel.selにはフルモデルしか入れていないので，dredgeの結果についてもフルモデル（8番のモデル）のみ注目して頂ければと思います。
どうも，dredgeを使うと，ランダム効果がないものとしてAIC等を算出しているように見えます。ただし，dredgeの入っているMuMInパッケージ説明ではglmmadmbがつかえるようなことが書いてあります（http://cran.r-project.org/web/packages/MuMIn/MuMIn.pdf）。どう考えればいいのでしょうか。dredgeは使うべきではないのでしょうか？よろしくお願いします。

結果の例
その1　glmmadmbを使い，dredgeした場合

> ana <- glmmadmb(monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)) + (1 | site), data =d, family= "nbinom")
> dredge(ana, fixed = "offset(log(Rdcamt))", rank="AIC")
Global model call: glmmadmb(formula = monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)) +
(1 | site), data = d, family = "nbinom")
---
Model selection table
(Int) b15 mnD p15 off(log(Rdc)) df logLik AIC delta weight
4 -2.894 0.2754 -3.760 + 4 -850.338 1708.7 0.00 0.597
8 -2.965 0.2838 -4.018 0.030100 + 5 -850.138 1710.3 1.60 0.268
2 -2.991 0.2669 + 3 -853.137 1712.3 3.60 0.099
6 -2.995 0.2672 0.001188 + 4 -853.136 1714.3 5.60 0.036
5 -1.546 -0.155200 + 3 -878.674 1763.3 54.67 0.000
7 -1.488 -2.011 -0.144700 + 4 -878.002 1764.0 55.33 0.000
3 -1.644 -3.030 + 3 -882.215 1770.4 61.75 0.000
1 -1.752 + 2 -883.755 1771.5 62.83 0.000
Random terms (all models):
‘1 | site’

その2　glmmadmbを使い（その1と同じモデル），model.selを使って，フルモデルの統計量を算出。その1と異なる値になる。なお，summaryで算出されるAICや係数はこれ（その2）と同じ値でした。

> ana <- glmmadmb(monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)) + (1 | site), data =d, family= "nbinom")
> model.sel(ana, rank = AIC)
Model selection table
(Int) b15 mnD p15 off(log(Rdc)) df logLik AIC delta weight
ana　　 -3.114 0.2455 -3 0.08234 + 6 -834.01 1680 0 1
Random terms (all models):
‘1 | site’

その3　glmmadmb（その1）からランダム効果項を削除して計算し，dredgeした場合＝＞その1の結果と一致

> ana <- glmmadmb(monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)), data =d, family= "nbinom")
> dredge(ana, fixed = "offset(log(Rdcamt))", rank="AIC")
Global model call: glmmadmb(formula = monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)),
data = d, family = "nbinom")
---
Model selection table
(Int) b15 mnD p15 off(log(Rdc)) df logLik AIC delta weight
4 -2.894 0.2754 -3.760 + 4 -850.338 1708.7 0.00 0.597
8 -2.965 0.2838 -4.018 0.030100 + 5 -850.138 1710.3 1.60 0.268
2 -2.991 0.2669 + 3 -853.137 1712.3 3.60 0.099
6 -2.995 0.2672 0.001188 + 4 -853.136 1714.3 5.60 0.036
5 -1.546 -0.155200 + 3 -878.674 1763.3 54.67 0.000
7 -1.488 -2.011 -0.144700 + 4 -878.002 1764.0 55.33 0.000
3 -1.644 -3.030 + 3 -882.215 1770.4 61.75 0.000
1 -1.752 + 2 -883.755 1771.5 62.83 0.000
> model.sel(ana_deer, rank = AIC)

その4　glmmadmb（その1）からランダム効果項をむりやり削除して計算し，model.selを使って，フルモデルの統計量を確認＝＞その1とその3と結果一致

> ana <- glmmadmb(monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)) + (1 | site), data =d, family= "nbinom")
> model.sel(ana, rank = AIC)
Model selection table
(Int) b15 mnD p15 off(log(Rdc)) df logLik AIC delta weight
ana　　 -2.965 0.2838 -4.018 0.0301 + 5 -850.138 1710.3 0 1

混合モデルではなく，GLMで試してみました。
その5　glm.nbを使い（当然ランダム効果は入っていません），dredgeした場合＝＞その1，3，4と一致

> ana <- glm.nb(monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)), data =d)
> dredge(ana, fixed = "offset(log(Rdcamt))", rank="AIC") #
Global model call: glm.nb(formula = monkN ~ b150 + p150 + meanmnkhuntD + offset(log(Rdcamt)),
data = d, init.theta = 0.6566762117, link = log)
---
Model selection table
(Int) b15 mnD p15 off(log(Rdc)) df logLik AIC delta weight
4 -2.894 0.2754 -3.760 + 4 -850.338 1708.7 0.00 0.597
8 -2.965 0.2838 -4.018 0.030100 + 5 -850.138 1710.3 1.60 0.268
2 -2.991 0.2669 + 3 -853.137 1712.3 3.60 0.099
6 -2.995 0.2672 0.001188 + 4 -853.136 1714.3 5.60 0.036
5 -1.546 -0.155200 + 3 -878.674 1763.3 54.67 0.000
7 -1.488 -2.011 -0.144700 + 4 -878.002 1764.0 55.33 0.000
3 -1.644 -3.030 + 3 -882.215 1770.4 61.75 0.000
1 -1.752 + 2 -883.755 1771.5 62.83 0.000

その6　glm.nbを使い（その5と同じ），model.selを使って，フルモデルの統計量を算出＝＞その1，3，4，5と一致

> model.sel(ana, rank=AIC)
Model selection table
(Int) b15 mnD p15 off(log(Rdc)) df logLik AIC delta weight
ana　　 -2.965 0.2838 -4.018 0.0301 + 5 -850.138 1710.3 0 1