「データ解析のための統計モデリング入門」読書ノート 7章 GLMMとGLMを比較する

前回：「データ解析のための統計モデリング入門」読書ノート：http://yagays.github.com/blog/2012/10/14/review-glm-mcmc/

前回の読書ノートでは全体の流れを簡単に追ったが，今回は「データ解析のための統計モデリング入門」7章の一般化線形混合モデルのシミュレーションを動かしてみた．

7章の概要

6章までは，ある植物100個体から得られた生存種子数と葉数のデータから，その関係性をモデル化しGLMで推定してきた．7章では，個体差を取り入れたモデルを作成し，生存種子数の過分散にうまくあてはまるようなパラメータを推定する．ただし，個体差を表すパラメータを個体ごとに割り当てて最尤推定するのではなく，個体差がある分布に従うと仮定した上で，その分布のパラメータを推定する．

大雑把な目標としては，葉の数だけが生存種子数に効いているわけではないんだから，それぞれの個体で観測されていない/できないデータというのを考えようというもの．ただし，個体ごとに何らかの値を振って生存種子数を上手く表すことのできるモデルを組んだとしても予測には使えないので，じゃあ個体差自体を扱えるように分布を仮定してそこからランダムに値が選ばれるとする．そうすると，個体ごとにパラメータを持たせるんじゃなくて，個体のばらつきという一つのパラメータで表すことができるので，より現実的なモデルとなると同時に簡単にパラメータの推定ができるようになる．

問題設定

ここでは，ロジスティック回帰におけるGLMMを考える．種子生存数が二項分布

$\mathrm{Bin}(y_i|\beta_1,\beta_2,s) = {{N}\choose{y_i}} q_{i}^{y_i}(1-q_i)^{N-y_i}$

に従うとした上で，そのパラメータ $q_i$ を

$\mathrm{logit}(q_i) = \beta_1 + \beta_2 x_i + r_i$

とおく．ここで，個体差 $r_i$ は

$r_i \sim \mathcal{N}(0,s)$

平均0，標準誤差sの正規分布に従うとする．

実験

今回の実験は，なぜ応答変数が過分散の場合にGLMでは駄目でGLMMで推定すべきかを確かめるために，本書に記載されているGLMおよびGLMMの実験を1000回行なって，そのときの切片 $\beta_1$ と傾き $\beta_2$ の分布を見ることにする．目標は「「個体差」の統計モデリング」(pdf)の図11をプロットすること．

各実験は，葉数が2〜6枚の個体を20個体ずつ合計100個体のデータが得られたという条件で，それぞれの生存種子数をパラメータ $\beta_1 = -4$ ， $\beta_2 = 1$ ， $r_i \sim \mathcal{N}(0,3)$ の二項乱数から生成した．このように生成したデータを用いて，GLMMとGLMを実行した．

この実験を1000回繰り返し，その時の $\beta_1$ と $\beta_2$ の分布をプロットしたのが以下の図となっている．左図が切片 $\beta_1$ の分布，右図が傾き $\beta_2$ の分布で，それぞれ黒の線がGLMが推定した値の分布，赤の線がGLMMが推定した値の分布となっている．また，GLMとGLMMの推定値の平均を縦の点線で表している．

library(glmmML)

logistic <- function(z){ 1/(1+exp(-z))}
Niter <- 1000
sigma <- 3
glm_b <- matrix(0, nrow=Niter, ncol=2)
glmm_b <- matrix(0, nrow=Niter, ncol=2)

for(i in 1:Niter){
  N <- rep(8,100)
  x <- rep(2:6,each=20)
  re <- rnorm(100, 0, sigma)
  y <- rbinom(100, 8, prob=logistic(-4+x+re))
  id <- 1:100
  d <- data.frame(N,x,y,id)

  glm_b[i,] <- glm(cbind(y,N-y) ~ x, data=d, family=binomial)$coefficients
  glmm_b[i,] <- glmmML(cbind(y,N-y) ~ x, data=d, family=binomial, cluster=id)$coefficients
}

par(mfrow=c(1,2))
plot(density(glm_b[,1]), xlim=c(-8,0), ylim=c(0,0.8), lwd=2, main="Estimate of beta_1", xlab="")
par(new=T)
plot(density(glmm_b[,1]), xlim=c(-8,0), ylim=c(0,0.8), col="red", lwd=2, main="Estimate of beta_1", xlab="")
abline(v=mean(glm_b[,1]), lty=2)
abline(v=mean(glmm_b[,1]), lty=2, col="red")

plot(density(glm_b[,2]), xlim=c(0,2), ylim=c(0,3), lwd=2, main="Estimate of beta_2", xlab="")
par(new=T)
plot(density(glmm_b[,2]), xlim=c(0,2), ylim=c(0,3), col="red", lwd=2, main="Estimate of beta_2", xlab="")
abline(v=mean(glm_b[,2]), lty=2)
abline(v=mean(glmm_b[,2]), lty=2, col="red")
legend(1.1, 3.0, c("glm","glmmML"), col=c(1:2), lwd=2)

さて，この図の解釈だが，今回データセット作成に用いた実験条件は $\beta_1 = -4$ ， $\beta_2 = 1$ だったので，GLMMにおいて切片・傾きどちらもある程度正しい推定ができていることがわかる．逆にGLMでは，切片は過大推定，傾きは過小推定されている．また，推定値のばらつきである分布の山の裾の拡がりを見てみると，GLMMでは大きく，GLMでは小さくなっている．これはいわゆるバイアス・バリアンスのトレードオフによるものだと思われる．

7章の概要

問題設定

実験

参考