（旧サイトより移行）過分散データ：GLM負の二項分布、GLMMによる推定をAICでモデル選択することは可能か？

# 2014.02.20追記：下記のlmer()は現在はglmer()に相当します

かつて、数値実験をしてみて「できないようだ」と書いたのですが、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に固定してみるか…また追々試してみたいと思います。

（旧サイトより移行）連続量を単位あたりに直してGLMする場合の対処：gaussian("log")、Gamma("log")、offsetなどなど）

# 2014.02.20 旧サイトを閉じるため、移植しました。
# なお、下記のようなケースではゴンペルツ曲線などを用いた成長モデルを使った方がいいかなと最近思ってます。

最近、実験系の統計を引き受けていて気になったので、連続量を単位あたりでGLM推定する方法について検証してみました。

単位あたりの量をGLMで解析する時、個体数や頻度などの整数の場合の対処はoffset項の利用で解決するのが常套手段ですね（密度や○×率のような割り算した値ではなく、元の値はいじらずに単位量を係数1として説明変数に加える。当ブログでもかつて紹介：http://nhkuma269.blog77.fc2.com/blog-entry-9.html）。
例えば、同じく個体密度0.1でも、1/10と100/1000とでは意味合いが違うが、割り算すると両者は同一密度として扱われてしまう。前者では1個体の増減が大きな誤差を生むが、後者ではほとんど影響なし。

では、元々が連続量のものを単位あたりの量にする場合も同様にoffset項による対処がよいのだろうか？例えば、サンゴの枝あたりのクロロフィル量とか、成長量（実験後重量 / 実験前重量）のような量は、必ず単位量（この場合は枝）のバラツキによる誤差が出てしまう（なるべく条件は揃えるように努力はしているだろうが）。これを割り算してしまうと誤差が統計結果に悪影響するだろう。

CrawleyのR統計本によると、単位量は共変量として説明変数に加えるか、log(目的の量 / 単位量)というように、割り算した上で対数変換せよ。と書かれていた。前者は要因の1つとしてカウントするならばそれでいいとして、後者はすでに古典的な対処法でしょう…（本自体がすでに一時代前のものなので仕方ないですが）。

ついでに、論文などで時々見掛ける、family=gaussian(link="log")、正規分布でリンク関数が対数という一見すると対数正規分布っぽいけれど違うらしい（対数正規分布の分散は平均と共に増大するようだ→平均や分散の式をチェックされたし）、こいつの性質もチェックしてみます。もう一つおまけに、古典統計で対数変換する際によく使うlog(x + 1)変換もチェックしてみました。ただし、これはlog(0)が計算できない問題を回避する目的に限定します。


比較するモデルは下記の6バージョン。これを応答変数が対数正規分布とガンマ分布の場合のそれぞれについて一次回帰式のパラメータ推定（切片と回帰係数）をします：
（色は下記の図とおおよそ対応）
m11：正規分布モデル、割合にした上でlog(x + 1)の対数変換という古典統計の常套手段（グラフ上では黒ライン）
m12：正規分布モデル、割合にした上でlog(x + 1e-10)変換、1の代わりにごく小さい小数(1.0 x 10^(-10))でlog(0)を回避
m13：正規分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定（以下、log(0)回避には1e-10を使用）
m14：正規分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定
m23：ガンマ分布モデル、割合にしているが、変数変換はせず、連結関数を対数に指定（ガンマ分布は正の値しか取れないので、1と2の変数変換のモデルは作れない）
m24：ガンマ分布モデル、単位量をoffset(log(off))として係数1の説明変数に加え、連結関数を対数に指定


結論から言うと、ケースバイケースな複雑な結果になりました…。実際の利用では、m13とm23、またはm14とm24をAICで比較するのがよいでしょう。
・やはりデータの対数変換は止めた方がよい、連結関数：logを用いるべきです。対数変換したモデルでは推定が大きく乱れる場合があり信頼できない（図1_1, 1_2のm11黒、m12赤）。
・データが対数正規分布の時、割合の正規分布モデル（m13緑）は、offsetの正規分布モデル（m14青）よりも推定の分散が小さかった（図1_1、2_1）。ただし、データがガンマ分布の時は、m14の方がm13よりも推定の分散が小さかった（図1_2、2_2）。両者の違いはさほど大きくなかった。
・ガンマ分布モデルでは、割合（m23水色）と、offset（m24ピンク）とは推定値がまったく同じ！（そのため、m23の水色は完全にマスクされている）
・つまり、この数値実験の限りでは、予想に反して連続量を割合や比率にした量をGLMで解析する際、割り算をした値を用いても大した問題はないことになる。offsetの使用による論文中での説明のややこしさや、記述可能なモデルの可塑性などを考えると…割り算していい気がしてきました。しかし、なぜ割り算にしても大丈夫だったのか、背景にある数学的なロジックは朧気なままです…たぶん、元の値と割算値とで確率分布が変わらないからだと思います（よくある整数値の割り算の場合は、小数点を取るようになるので本来のポアソンや二項分布が使用不可になるが）。例外があるとすれば、単位量の方にも確率誤差（その値が属している確率分布からのズレ）がある場合でしょう。その場合はoffsetすべきということかと。
・ちなみに、常套手段なlog(x + 1)変換（m11黒）の推定は危険です。とくに平均値の小さい推定では危険（図1_1、1_2）。log(0)の回避には1の代わりにごく小さな値（1e-10くらい）を足すのがよいでしょう。ただし、古典統計で整数を対数変換する場合は、平均と分散の関係を調整する意図があるので + 1 のままで）。

なお、ここでやっているのは入り口と出口がちゃんと一致するかを確認するだけの数値実験に過ぎません。参考にする際には自己責任でお願いします。



# 以下、こんな関数でGLM推定を1000回繰り返し、パラメータ推定の精度をチェックしました。
# 平均 exp(alpha + beta*X)という一次回帰を考え、こちらで指定したalphaとbetaを用いて、対数正規分布とガンマ分布のデータを発生させ、各モデルによってalphaとbetaを逆推定し、初めに指定した値を再現できるかどうかチェックします。
# alpha, beta, distributionは1が対数正規分布、2がガンマ分布、n.itrは繰り返し数を表します。

cont.skew <- function(alpha, beta, distribution, n.itr) {
estim <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの
set.seed(1)
X <- rep(c(1:10)*0.1, each=10)
off <- runif(length(X), min=1, max=10) # 下記のoffsetで用います。数字に特に意味はなし
mean <- exp(alpha + beta*X)*off # 平均：log(Y/off) ~ alpha + beta*X、なので。
sd <- exp(0.5)
for (n in 1:n.itr) {
Y <- rbind(rlnorm(100, meanlog=log(mean), sdlog=log(sd)), # 対数正規分布
rgamma(100, shape = mean^2/sd^2, scale = sd^2/mean)) # ガンマ分布（cf. rgamma()のhelp）
Y <- Y[distribution, ]
D <- data.frame(X, Y, off)
m11 <- glm(log(Y/off + 1) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル（割合、対数変換1）
m12 <- glm(log(Y/off + 1e-10) ~ X, family=gaussian(link="identity"), D) # 正規分布モデル（割合、対数変換2）
m13 <- glm((Y/off + 1e-10) ~ X, family=gaussian(link="log"), D) # 正規分布モデル（割合、連結関数：対数）
m14 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=gaussian(link="log"), D) # 正規分布モデル（オフセット、連結関数：対数）
m23 <- glm((Y/off + 1e-10) ~ X, family=Gamma(link="log"), D) # ガンマ分布モデル（割合、連結関数：対数）
m24 <- glm((Y + 1e-10) ~ X + offset(log(off)), family=Gamma(link="log"), D) # ガンマ分布モデル（オフセット、連結関数：対数）
estim2 <- c(rbind(coef(m11), coef(m12), coef(m13), coef(m14), coef(m23), coef(m24)))
estim <- rbind(estim, estim2)
} # forループ、ここまで
par(mfcol=c(1,2))
plot(density(estim[,1]), lwd=4, xlim=c(alpha-0.5, alpha+2), ylim=c(0,3), main="alpha (intercept)")
for (i in 2:6) { lines(density(estim[,i]), lwd=4, col=i)
abline(v=alpha) } # alpha値の図示
plot(density(estim[,7]), lwd=4, xlim=c(beta-1, beta+1), ylim=c(0,3), main="beta (coefficient)")
for (j in 2:6) { lines(density(estim[,j+6]), lwd=4, col=j)
abline(v=beta) } # beta値の図示
estim.m <- matrix(apply(estim, 2, mean), ncol=2) # alpha, betaの推定値を計算
colnames(estim.m) <- c("alpha", "beta")
estim.m} # 推定値を表示（cont.skew関数、ここまで）

# 図1_1： 対数正規分布のパラメータ推定、平均値が小さい場合
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=1)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# m12赤、m23&m24ピンク、m13緑、の順に推定がよかった。m11黒は大きくずれた。
# 推定値（1~6の順にm11～m24）、m23とm24（一番下の2つ）は完全に推定値が一致！
# alpha beta
# [1,] 0.3131626 0.3959152
# [2,] -1.0009638 0.9999089
# [3,] -0.8811454 1.0011294
# [4,] -0.8879904 1.0073653
# [5,] -0.8793282 0.9996438
# [6,] -0.8793282 0.9996438

# 図1_2：ガンマ分布のパラメータ推定、平均値が小さい場合
cont.skew(alpha=-1, beta=1, n.itr=1000, distribution=2)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# m14青、m13緑、m23&m24ピンク、の順に推定がよかった。m11黒とm12赤は大きくずれた。
# なぜかガンマ分布の推定なのにガンマ分布モデル（m23&m24）が最強にならない…

# 図2_1：対数正規分布のパラメータ推定、平均値が大きい場合
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=1)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# 平均値が大きくなると、モデルによる推定の違いは小さくなる（中心極限定理による正規分布への収束）
# m12赤、m23&m24ピンク、m11黒、の順に推定がよかった。

# 図2_2：ガンマ分布のパラメータ推定、平均値が大きい場合
cont.skew(alpha=2, beta=1, n.itr=1000, distribution=2)
# m11：黒、m12：赤、m13：緑、m14：青、m23：水色、m24：ピンク（水色はピンクと完全に重なって表示されない）
# 平均値が大きくなると、モデルによる推定の違いは小さくなる（中心極限定理による正規分布への収束）
# m11黒以外は真の値にほぼ収束

（旧サイトより移行）過分散データ：GLM負の二項分布、GLMMによる解析の比較

## 2014.02.20 追記：旧サイトを閉じるため、このページを移植しました。
なお、現在は lmer()は正規分布専用になっており、その他の分布ではglmer()関数を使用します。

##（2012.10.17 追記：GLM関数群とGLMM関数群との間でのモデル選択について、数値実験結果を追加・修正しました）

# 生態学会で話していて、過分散データの取り扱いが話題に上った。負の二項分布を使えという意見と、GLMMを使えという意見と両方があるが、実際どっちの方がよいというのはあるのだろうか？（cf. 2014.6.24追記： GLM負の二項分布とは実態はランダム効果がガンマ分布となったGLMM）

# 過分散 overdispersion：ポアソン分布や二項分布のモデルに当てはめる場合、これらの分布型が平均と分散が表裏一体なのだが（ポアソンでは平均=分散)、実際のデータは大抵はそうでないので分散が過剰となることが多い、という問題。
# GLMのsummary出力結果で、下の方に、こんなのがあって…
# Residual deviance: 5094.7 on 71 degrees of freedom
# …ここで、5094.7/71の値が、1.5を超えているならば、何らかの対処をした方がよいだろう（cf. モデル選択ペンギン本、Zuur et al (2009) Mixed effects models and extensions in Ecology with R. Springer。これのp. 225-6には"1.5"が目安と書かれている）

# 過分散しているデータとは下記のように裾野が片方に広く延びている（上：ふつうのポアソン分布、中：負の二項分布、下：正規ノイズを加えたポアソン分布）。負の二項分布はゼロのあたりに集中、正規ノイズは正負にばらつく。

par(mfcol=c(3,1)) # グラフ領域を縦に3分割
hist(rpois(10000, lambda=exp(1)), freq=F, xlim=c(0,100), breaks=seq(0,20,by=2), col=1) # ポアソン分布（平均 exp(1)）
hist(rnegbin(10000, mu=exp(1), theta=1/4), freq=F, xlim=c(0,100), breaks=seq(0,10000,by=2), col=2) # 負の二項分布（平均 exp(1)、θ 1/4。θはガンマ分布のshapeパラメータに相当、shape=平均^2 / 分散）
hist(rpois(10000, lambda=exp(1 + rnorm(10000, mean=0, sd=4^0.5))), freq=F, xlim=c(0,100), breaks=seq(0,10000,by=2), col=3) # ポアソン分布（平均 exp(1)）＋正規ノイズ（平均0、分散4）

# では、ポアソン分布のデータと、負の二項分布をしたデータ、ポアソン分布に正規乱数ノイズを加えたデータを人工的に作成し、それぞれGLM (poisson)、GLM（負の二項分布）、GLMM（データ数分の個体差があると見なす）で推定した結果を比べてみよう（実行するためのプログラムは末尾に）。


OD(alpha=2, beta=0.5, siml=100, distribution=1) # データがポアソン分布の時
#（マシンパワーによってはけっこう時間が掛かる、アラートも沢山出るが無視）
# ランダム効果がN数と同じ数だけあるよ！？というアラートが大量に出てくるが、推定はちゃんとできているので無視
# Number of levels of a grouping factor for the random effects is *equal* to n, the number of observations

# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果
#（３本のラインがあるはずですが完全に重なっていて、視覚的には識別不能）

# models: glm.pois glm.nb glmm.pois
# aic.d 710.5256 712.1975 712.2007 # devianceに基づくAIC
# prob.d 0.9500 0.0400 0.0100 # 各モデルが選択される確率
# aic.rd 102.8072 101.6366 104.4823 # residual devianceに基づく AIC
# prob.rd 0.6000 0.4000 0.0000 # 各モデルが選択される確率

# GLM (poisson)、GLM（負の二項分布）、GLMMともalpha、betaの推定値はほぼ真の値となった（ということは、単なるGLM（poisson）を使う理由はないのではないか！？）。AIC（deviance）によるモデル選択は確率0.95で真のモデルであるGLM（ポアソン分布）を選択したが、AIC（residual deviance）ではAICの推定値が僅差で負の二項分布モデルが最小だし、誤って負の二項分布モデルを選択する確率が0.4もあった。とはいえ、いずれのモデルもほぼ真の値を推定できているので問題なしか！？



OD(alpha=2, beta=0.5, siml=100, distribution=2) # データが負の二項分布の時
# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果

# models: glm.pois glm.nb glmm.pois
# aic.d 36960.67 902.9074 922.3840 # devianceに基づく AIC
# prob.d 0.00 1.0000 0.0000 # 各モデルが選択される確率
# aic.rd 36572.47 121.0299 534.1787 # residual devianceに基づく AIC
# prob.rd 0.00 1.0000 0.0000 # 各モデルが選択される確率

# GLM (poisson)は両方向に大きく広がった。GLM（負の二項分布）はalphaもbetaも真の値に近かった。GLMMはbetaの値は真の値に近かったが分布の裾が広かった、alphaは盛大にズレた。
モデル選択はすべてGLM（負の二項分布）となった。



# ポアソン分布に正規乱数ノイズを加えたデータを同様に推定してみる
OD(alpha=2, beta=0.5, siml=100, distribution=2)

# 黒：GLM (poisson)、赤：GLM（負の二項分布）、緑：GLMM。100回分の推定結果

# models: glm.pois glm.nb glmm.pois
# aic.d 845220.9 1309.4785 1278.7899 # devianceに基づく AIC
# prob.d 0.0 0.0100 0.9900 # 各モデルが選択される確率
# aic.rd 844616.2 136.0312 674.1263 # residual devianceに基づく AIC
# prob.rd 0.0 1.0000 0.0000 # 各モデルが選択される確率

# GLM (poisson)は真の値から大きく外れた。GLM（負の二項分布）はalphaは盛大にズレた、betaは真の値に近かったが裾野が広かった。GLMMはalphaもbetaも真の値に近かった。AIC（deviance）によるモデル選択は0.99の確率で真の母集団のモデルGLMMを選択したが、AIC（residual deviance）によるモデル選択はすべてのケースで誤ってGLM（負の二項分布）を選択してしまった…。



# 結論として、過分散データと言っても一辺倒ではなく、負の二項分布、正規ノイズ、どちらに近いかによって結果が変わってくるという、ごく当然の結果となった。glm.nb、GLMMのどちらがよいかは、結局データ依存ということになりました。

# 選択基準としては…glm.nb、GLMMでAICを基にモデル選択は、一手間加えれば可能であるようだ。また、glmとglm.nbの方をresidual devianceに基づくAICに修正する計算も試してみたところ（2*logLik(pois) + deviance(pois) + AIC(pois)）、誤ってglm.nbモデルが選択される確率が大きく増大することがわかった。なので、GLMとGLMMのモデルをAICで比較する場合には、GLMMのAICをdevianceベースのAICに修正して用いる必要がある（-2*logLik(pois) - deviance(pois) + AIC(pois.mm)）ということですね。（注：poisはglm(,family=poisson)のモデル、pois.mmはlmer(,family=poisson)のモデル、両者の違いはランダム効果の有無のみ → 下記のプログラム中にモデル式あり）


# しかし、lmerやglmmMLのデフォルトがresidual devianceに基づくAICを計算していること自体が、そういう比較を推奨しないというメッセージにも思える。

# もっとも、ポアソン分布の場合は今回のような選択肢があるけれど、二項分布の場合はGLMMを使うくらいしか対処の方法がない（怪しげなQAICの使用を回避するならば）。

# 2014.06.24 追記：実際のデータのモデル化を考えた時、固定効果のみで推定値を求めるときGLMM (poisson)とglm.nbは異なる推定値になってしまうのは気になるが、もはやモデルの違いと割り切るしかない気がしている。双方ともGLMMなので、誤差はランダム効果でオーバーフィットしてしまうため、AICなどで比較して差異があっても実質的なモデルの良し悪しを反映しているのだろうか？現実のデータの場合は真の答えは分からないし、モデルを使用している以上は仮定を捨てられないので、悩んでも仕方がないと振り返って思います。

# 参照：
# lmerはresidual devianceに基づいたAICを返しているので、直接比較をするにはdevianceに基づいた計算にし直す必要がある
# deviance: -2*logLik()
# residual deviance: deviance()
# AIC(glm())の場合は、deviance + 2*推定するパラメータ数
# AIC(lmer())では、residual deviance + 2*推定するパラメータ数





# 上記の数値実験を実行するためのプログラムはこちら：
# 平均 exp(alpha + beta*x)とし、真の値はalpha = 2、beta = 0.5、あと負の二項分布の歪みのパラメータθ= 1/4、正規乱数ノイズの分散4も設定
# （ただし、θや分散の値については今回はチェックはしない）

library(MASS) # glm.nbの呼び出し
library(lme4) # lmerの呼び出し

OD <- function(alpha, beta, siml, distribution) {
# 切片、傾き、シミュレーション回数、データ分布（1：ポアソン分布、２：負の二項分布、３：ポアソン分布＋正規ノイズ）
estim <- aic.d <- aic.rd <- prob.d <- prob.rd <- numeric(0) # 後でデータをくっつけるためのイントロンみたいなもの
set.seed(36) # 乱数の種を指定（他の値にすると推定エラーが出るかもですglm.nbやlmerの問題）
for (s in 1:siml) {
Noise <- rnorm(100, 0, 4^0.5) # ノイズの分散を4とした
x <- rep(c(0:9), each=10)
y <- rbind(rpois(100, lambda=exp(alpha + beta*x)), # ポアソン分布
rnegbin(100, mu=exp(alpha + beta*x), theta=1/4), # 負の二項分布
rpois(100, lambda=exp(alpha + beta*x + Noise))) # 普通のポアソン分布にノイズを追加
y <- y[distribution,] # データ分布の選択
ID <- 1:length(y) # データの数だけ個体差（ID）を設定
D <- data.frame(x, y, ID)
pois <- glm(y ~ x, family=poisson, data=D) # GLM (poisson)
negb <- glm.nb(y ~ x, data=D) # GLM (負の二項分布)
pois.mm <- lmer(y ~ x + (1|ID), family=poisson, data=D) # GLMM
estim2 <- c(cbind(coef(pois), coef(negb), fixef(pois.mm)))
estim <- rbind(estim, estim2) # simlの数だけ推定値を重ねていく
aic.p <- AIC(pois) + 2*logLik(pois) + deviance(pois) # residual devianceに基づくGLMのAIC
aic.n <- AIC(negb) + 2*logLik(negb) + deviance(negb) # residual devianceに基づくGLM.nbのAIC
aic.m <- -2*logLik(pois)[1] - deviance(pois) + AIC(pois.mm) # devianceに基づくGLMMのAIC
Aic <- c(AIC(pois), AIC(negb), aic.m) # deviance に基づくAIC
Aic2 <- c(aic.p, aic.n, AIC(pois.mm)) # residual deviance に基づく AIC
aic.d <- rbind(aic.d, Aic)
aic.rd <- rbind(aic.rd, Aic2)
prob.d <- c(prob.d, which.min(Aic))
prob.rd <- c(prob.rd, which.min(Aic2))
}
aic.d <- apply(aic.d, 2, mean)
aic.rd <- apply(aic.rd, 2, mean)
prob.d <- factor(prob.d, levels=factor(c(1:3))) # モデルを示す番号をカテゴリー化
prob.d <- table(prob.d)/siml # AICモデル選択の結果を集計
prob.rd <- factor(prob.rd, levels=factor(c(1:3))) # モデルを示す番号をカテゴリー化
prob.rd <- table(prob.rd)/siml # AICモデル選択の結果を集計

select <- rbind(aic.d, prob.d, aic.rd, prob.rd)
colnames(select) <- c("glm.pois", "glm.nb", "glmm.pois")
par(mfcol=c(1,2))
plot(density(estim[,1]), lwd=4, xlim=c(-2,6), ylim=c(0,1), main="alpha (intercept)")
for (i in 2:3) { lines(density(estim[,(2*i-1)]), lwd=4, col=i)
abline(v=alpha) } # alpha値の図示
plot(density(estim[,2]), lwd=4, xlim=c(0,1), ylim=c(0,6), main="beta (coefficient)")
for (j in 2:3) { lines(density(estim[,(2*j)]), lwd=4, col=j)
abline(v=beta) } # beta値の図示
return(select) }

負の二項分布のパラメータ推定：R と WinBUGSの比較

過分散した整数データのモデリングによく使う負の二項分布、Rではglm.nb関数でできるが、WinBUGSの例をほとんど見ない。ちょうど使う用事があったので使い方をテストしてみた。

# サンプルデータは、平均 exp(2)、dispersion parameter 1.5 の負の二項分布から抽出した N = 1000 のデータ（後で対数リンク関数を使用するのでexp()で表記）
require(MASS)
Y <- rnegbin(1000, mu=exp(2), theta=1.5) 
# 注：乱数なので実行する度に結果は少しずつ変わる

# まずはRのglm.nb関数による推定
require(MASS)
summary(glm.nb(Y ~ 1))

Coefficients:

                  Estimate Std. Error z value Pr(>|z|)    

(Intercept)  1.97519    0.02771   71.28   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.5895) family taken to be 1)

# 平均 1.97519（expの中身に相当）、dispersion parameter 1.5895 と順当な推定

# 次にWinBUGS、負の二項分布だからdnegbin、と試してみるが…（thetaがdispersion parameterに相当しそう）

data <- list(N=1000, Y=Y)

parm <- list(theta=rlnorm(1), alpha=rnorm(1))

model <- function() {

theta ~ dgamma(1E-1, 1E-2)

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dnegbin(p[i], theta) # 負の二項分布では平均をそのまま使えない

p[i] <- theta/(theta+mu[i]) # pと平均muの関係（pは0~1の値を取る）

log(mu[i]) <- alpha # muは正の値

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000, debug=T) # 自作ラッパー関数wbugsを使用

# order of negative binomial *** must be an integer というエラー

# theta は負の二項分布の狭義の定義通り、整数でないといけないようだ

# 戸惑いつつも多項分布の事前分布で thetaが正の整数しか取らないようにし再計算

data <- list(N=1000, Y=Y, q=rep(1,20)/20)

parm <- list(theta=sample(c(1:5), 1), alpha=rnorm(1))

model <- function() {

theta ~ dcat(q[]) # 事前分布1~20の整数（多項分布使用）

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dnegbin(p[i], theta)

p[i] <- theta/(theta+mu[i]) 

log(mu[i]) <- alpha 

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)

               mean        sd          2.5%      25%        50%        75%        97.5%   Rhat   n.eff

theta        2.000       0.000    2.000    2.000       2.000       2.000       2.000          1        1

alpha       1.976       0.026    1.927    1.958       1.976       1.993       2.026          1  1500

deviance 6032.409 1.446 6031.000 6032.000 6032.000 6033.000 6036.000    1  1500

# 平均 alphaの推定は1.976と正しく推定できているのだが、thetaは当然のように整数値となった…

# Rの?rnbinomなどを読み返してみると、"dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer" とある。

# gamma mixing distributionのshapeパラメータに相当するということか。

# なので、ポアソン分布にガンマ分布を混ぜた階層ベイズとして負の二項分布の推定を試みる

data <- list(N=1000, Y=Y)

parm <- list(shape=rlnorm(1,log(1),log(2)), alpha=rnorm(1), lambda=rlnorm(1000,log(1),log(2)),F) 

model <- function() {

shape ~ dgamma(1E-1, 1E-2)

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dpois(lambda[i]) # 負の二項分布を階層ベイズで表現

lambda[i] ~ dgamma(shape, rate[i])

rate[i] <- shape / mu[i] # rate と shape、平均 mu の関係を表す

log(mu[i]) <- alpha 

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)

                    mean     sd     2.5%      25%      50%      75%    97.5%  Rhat   n.eff

shape           1.587  0.087    1.422    1.522    1.587    1.647    1.755 1.002  1200

alpha           1.975  0.028    1.917    1.957    1.975    1.993    2.029 1.003  1200

# lambdaは省略

# 平均 1.975、shape 1.587 と 元の値やglm.nbの推定とほぼ一致、これが正解のようだ

# それにしてもややこしい…負の二項分布には異なる定義が多すぎる！