かつて、数値実験をしてみて「できないようだ」と書いたのですが、GLMM側（lmerやglmmML）のAIC計算方法をGLM側（glmやglm.nb）へと合わせるように修正すれば可能であることを確かめました。じつは、両者でAICの計算方法が異なっていたことが分かりました。道理でいつもGLMMのAICが小さ過ぎるわけです…。
★以前アップしていた数値実験を大幅に修正・加筆しておきました（cf.「過分散データ：GLM負の二項分布、GLMMによる解析の比較」）
なお、glmやglm.nb と、lmerやglmmMLのAICの計算方法の違いを検証するところまでは、検索すると複数のブログで見つかります。
ところがさらに検証してみたところ、逆にglmやglm.nbのAICの方をlmerやglmmMLのAIC計算方法へと修正すると、不当にglm、とくにglm.nbのAICが小さくなり、誤ったモデル選択となる確率が大幅に増大することも分かりました（ほぼ逆転してしまう！）。もっとも、ちゃんと数値実験をやっているとはいえ、なぜそうなるのか理論的裏付けがよく分からないので、なにか思いついたら検証したいと思います。
ちなみに、まったくのポアソン分布データのパラメータ推定でも、glm、glm.nb、lmerはいずれも真の値をほぼ推定できていました。過分散問題がけっこう厄介なことを思うと、もはや glm(, family=poisson) を紹介する意義はほとんど無いのかもしれません…。
数値実験では、乱数発生させた、ポアソン分布データと、負の二項分布データ、ポアソン分布に正規乱数ノイズを加えたデータの三者で、GLM（ポアソン分布：glm）とGLM（負の二項分布：glm.nb）、GLMM（ポアソン分布：lmer）が真の値を推定できているかをクロスチェックし、同時にAICが正しい推定をちゃんと反映できているかをもチェックしました。lmerの方は本来のGLMMの使い方とは異なり、各データを異なるid（1データ1id）としているので過分散の対処に用途は限定しています（こういう使い方は、help(lmer)や、Warton & Hui (2011) Ecology 92:3–10、和文では粕谷先生のGLM本「一般化線形モデル」で紹介されています）。
混乱しやすいのでメモ的に整理します。
まず、AIC()を用いた時の、glmやglm.nbと、lmerやglmmMLの計算方法の違い：
glmやglm.nbでは、
AIC = deviance + 2*推定するパラメータ数
lmerやglmmMLでは、
AIC = residual deviance + 2*推定するパラメータ数
また、deviance()は、いずれの場合もresidual devianceを計算する。
一方、logLik()では、glmやglm.nbではdeviance/(-2)、lmerやglmmMLではresidual deviance/(-2)
lmerやglmmMLのAICをglmの方へ合わせるには、下記のようにする（上述の理由により、逆はやめておきましょう）。
# 例えば、こんな２つのモデルがあるとき（id=1~N数）
pois <- glm(Y ~ X, family=poisson)
pois.m <- lmer(Y ~ X + (1|id), family=poisson)
glmのAICはふつうに、AIC(pois)、で求められる。
GLMMのAICは、上述の違いに基づくと、こうすればglmのと同じように求められる：
AIC(pois.m) - deviance(pois) -2*logLik(pois)
# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数
#（ランダム効果を抜く以外は同等であるglmモデルを用意してやる必要がある）
# ちなみに、AICがそのように変換できることの検証として、GLMMでidを全部 1 にしてやると（1グループのみ）、グループ間の分散は 0 と推定されます。切片と回帰係数の推定も、単なるGLMのそれと同一の推定結果となります。ただし、分散値は推定するパラメータ数としてカウントされるので、変換後のGLMMのAICはGLMのAICよりもキッカリ2だけ増加した値となります。これは、AICの計算式にある、2*パラメータ数（ここでは、2*1）、の部分に相当します。
しかし、なぜGLMMのパッケージがresidual devianceに基づいたAIC計算にしているのか、その意図が何なのかわかりません。こういう計算方法の修正をしないでくれと言っているようにも見えます。さしあたり、過分散データへの対処という限定的な使用法について数値実験をしてみた限りでは問題は無いようです。
こういう煩わしさを思うと、もはや個人的にはBUGSに全面移行してしまおうかなどと考えたりします…Stanの発展も待ち遠しいです。
# （当時いただいたコメント）
2012/11/15(木) 23:39:11 | URL | 高橋
興味深く読ませていただきました。
ちょっと気になったのですが、
＞# (residual deviance + 2*パラメータ数) - residual deviance + deviance = deviance + 2*パラメータ数
この式の２つのresidual devianceは別物だと思います。
residual deviance = (当該モデルのdeviance) - (飽和モデルのdeviance)
なので、
＞- deviance(pois) -2*logLik(pois)
は、
- ((ポアソンモデルのdeviance) - (ポアソンモデルの飽和モデルのdeviance)) + (ポアソンモデルのdeviance)
となります。
したがって、
＞AIC(pois.m) - deviance(pois) -2*logLik(pois)
は、
(混合モデルのdeviance) - (混合モデルの飽和モデルのdeviance) + 2*パラメータ数 + (ポアソンモデルの飽和モデルのdeviance)
となり、「混合モデルの飽和モデル」が「ポアソンモデルの飽和モデル」と等しい場合（多分そうなのだと思いますが）、混合モデルのdevianceに基づいてAICを求めているということになると思います。
glm.nb()の関数定義を見ると、推定されたThetaのもとで飽和モデルのdevianceを計算しているようです。この値は、ポアソンモデルにおける飽和モデルのdevianceよりずっと大きくなるので、residual devianceを用いてAICを計算すると、負の二項分布の常勝となるのでしょう。
ちなみに飽和モデルの対数尤度は、応答変数をyとすれば、
sum(dpois(y, y, log = TRUE))
で求められます。負の二項分布の場合は、glm.nb()が返すThetaパラメータを使って、
sum(dnbinom(y, mu = y, size = Theta, log = TRUE))
> 高橋さま
コメントありがとうございます、おかげで理解が深まりました。この検証に興味を持っていただき嬉しいです。
ポアソンモデルと、ポアソン混合モデルの飽和モデルが等しい時にのみ成り立ちうるというご指摘、たしかにその通りでした。
これはもう数値実験を見る限りそうなるとしか答えようがないのですが、1グループのみの混合モデルの変換後AICから2を引いた値とポアソンモデルのAICが一致する、というのが一応の証拠になっていると考えています（こちらの検証コード掲載は省いてしまっていますが）。
glm.nbの飽和モデルの求め方がそうなっているとは気づいていませんでした。言われてみて改めて関数定義をチェックするとたしかにそれらしきコードを見つけることができました。
混合モデルならば、グループを1つだけにするなどにより、ポアソンモデルとのすり合わせを試みることができましたが、負の二項モデルで同じようなすり合わせをするとなると…glm関数でThetaを1に固定してみるか…また追々試してみたいと思います。
0 件のコメント:
コメントを投稿