データ解析のための統計モデリング入門 一般化線形混合モデル(GLMM) 読書メモ7


2017年 06月 06日

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。
今回は第7章、一般化線形混合モデル(GLMM)についてのまとめの七回目です。

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。
その中でポアソン分布と正規分布を混ぜた分布について説明されています。
なので、ポアソン分布と正規分布を掛けた値をグラフにしてくれるコードを用意しました。
コードはRで書きました。

linkFunction <- function(b1, b2, x) {
  function(r) {
    exp(b1 + b2 * x + r)
  }
}

distribution <- function(y, s, linkFunction) {
  l <- match.fun(linkFunction)
  function(r) {
    dpois(y, l(r)) * dnorm(r, mean = 0, sd = s)
  }
}

xl <- "r"
yl <- "Probability"
ys <- seq(0, 8, by = 2)
rs <- seq(-5, 5, by = 0.5)
b1 <- 0.1
legends = paste0("y = ", ys)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    l <- linkFunction(b1, b2, 2)
    title <- paste0("lambda  = exp(", b1, " + ", b2, " * 2 +  r), r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.2), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      d <- distribution(ys[i], s, l)
      lines(rs, d(rs), type = "l")
      points(rs, d(rs), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。
r は平均 0、標準偏差 s の正規分布から生成される変数です。
あとは、正規分布から r が生成される確率と、平均 λ のポアソン分布を掛け合わせれば、欲しい値になります。

コードを実行すると20枚のグラフがプロットされます。
例えば次のようなグラフが得られます。

mixed_poission.png

y によってグラフの山の大きさが変わってくるのが分かります。
混合分布はこれを r で積分したものです。
なので、グラフ上で面積の大きな y が確率の高い分布になります。