※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。

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


2017年 05月 16日

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第7章、一般化線形混合モデル(GLMM)についてのまとめの一回目です。 この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 まずはロジスティック回帰にパラメータを追加したモデルについて説明されています。 なので、いろいろなパラメータについてそのモデルをプロットしてくれるコードを用意しました。 コードはRで書きました。
N <- 20

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

distribution <- function(x, r, linkFunction) {
l <- match.fun(linkFunction)
function(y) {
dbinom(y, N, l(x, r))
}
}

rs <- -2:2
ys <- 0:N
pchs <- 1:5
xl <- "y"
yl <- "Probability"
legends <- paste0("r = ", rs)
for (b1 in c(-0.5, 0.5, by = 0.2)) {
for (b2 in c(-0.5, 0.5, by = 0.2)) {
l <- linkFunction(b1, b2)
title <- paste0("logit(q) = ", b1, " + ",b2," * 2"," + r")
plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
for (i in 1:5) {
d <- distribution(2, rs[i], l)
lines(ys, d(ys), type = "l")
points(ys, d(ys), pch = pchs[i])
}
legend("topright", legend = legends, pch = pchs)
}
}
このモデルでは、二項分布のパラメータである事象の発生確率 `q` が以下のコードで与えられます。
linkFunction <- function(b1, b2) {
function(x, r) {
1 / (1 + exp(-b1 - b2 * x - r))
}
}
GLMのモデルと比べると、新しい変数として `r` が追加されています。 この `r` が説明変数 `x` とは違った、何とは特定されていない、余分な事象の発生確率への影響を表します。 コードを実行すると様々な `b1` と `b2` と `r` の値に対してモデルをプロットします。 `x` は `2` で固定にしました。 事象の試行回数は `N = 20` で、事象の発生回数は `y` です。 例えば下図のような図がプロットされます。 mixed.png 図を見ると、確かに `r` の増減に従って `y` の分布が影響を受けるのが分かります。