なお、現在は 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) }
0 件のコメント:
コメントを投稿