glm()
)
な単純なハナシなんだけど,
すごく難しかった.
また多くの教訓をえてしまった.
stepAIC()
は信用するな,
というもの.
理由は
stepAIC()
が発見してくれたやつより,
さらに良い (多くの場合「さらに簡単な」)
モデルがすぐに見つけられる場合があるから.
stepAIC()
はちょー強力・便利きわまりない道具なんだけど,
自動「モデル簡単化」の手順に制約を設けていて,
「あまり簡単ではないモデル」
が最良モデルとして選ばれる場合がある.
具体的な状況としては,
たとえば二種の植物種 (s) とすみ場所 (h)
における開花について調べてるときに
stepAIC(glm((開花) ~ s * h, ...))
で選ばれる推定結果は
s == B
は開花しやすい (推定された係数の大きさ
+b とする)
h == Y
は何か効果ある? (係数の大きさがゼロにちかい)
s:h == sB:hY
は開花しにくくなる
(大きさ -b
…… これは s == B
の効果を打ち消すだけのために出現した量)
logit(開花) ~ (Intercept) + s + h + s:h
という 4 パラメーターモデルが最良になるだろう.
h
(すみ場所)
の係数なんかはほぼゼロの値であっても残ってしまう.
stepAIC()
では s:h
が残ると h
の効果ありなしについては検討しないからだ.
stepAIC()
なやつより)
よい線形モデルの例としては
logit(開花) ~ (Intercept) + sBhY
という 2 パラメーターモデル,
この
sBhY
は
s == B
かつ
h == Y
のときのみ TRUE
になるような boolean 変数である.
stepAIC()
が悪いわけではなく,
線形モデル・一般化線形モデルとはそういったものなのである.
しかしながら,
上の例で示している現象には強い非線形性
があるので
(変数そのまま列挙しただけの素朴な)
線形モデルのたぐいではうまく表現できていない
……
そういうことだ.
しかし我々はアタマがかたいのでこれを一般化線形モデルでなんとかしよう,
とするわけだ.
glm()
使うマシな方法としては
(stepAIC()
ではなく)
「s == B
かつ h == Y
」
の場合だけちがう,
という「場合わけ」を R に自動的に探索させる
(計算生態学自由集会 2004
の久保発表参照)
というのも可能だ.
しかしこの方法はちょいめんどくさい.
glm()
連動の z 検定の結果を参考にしつつ,
手さぐりの探索を試みてみよう
(z 検定そのものは変数選択には使えない
……
理由: これは「検定」だから).
stepAIC()
を素朴に用いてみる
logit(開花) ~ (Intercept) + s + h + s:h
,
h
がほぼゼロ,
s
と s:h
が打ち消されるような,
解釈不可能な複雑なモデルが選択される
summary()
してみる
……
すると係数の推定値とその推定のばらつき (SE),
この両者から導出される z 値と z 検定の p-value
が表示される
h
の |z| 値が小さい)
summary()
結果みながら,
(上でのべた)
sBhY
のような新しい説明変数が作れないか検討してみる
(このあたり,
一種の非線形モデリングを試みているのである)
glm((開花) ~ sBhY, ...)
を推定させてみる
stepAIC()
最良モデルより
さらに小さければ,
この logit(開花) ~ (Intercept) + sBhY
が (AIC で評価される世界において)
よりよいモデルである
(stepAIC()
の能力の限界を
手動で補正している,
ということ)
stemAIC(glm())
するといろいろ面倒な
……