データ解析のための統計モデリング入門 GLMの応用範囲をひろげる 読書メモ2

このブログ記事は『データ解析のための統計モデリング入門』（久保拓弥　著、岩波書店）という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第６章、GLMの応用範囲をひろげるについてのまとめの二回目です。 この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、いろいろなパラメータについてロジスティック回帰のモデルをプロットしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

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

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

xs <- seq(0, 10, by = 2)
ys <- 0:N
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
for (b1 in seq(-0.5, 0.5, by = 0.2)) {
for (b2 in seq(-0.5, 0.5, by = 0.2)) {
l <- linkFunction(b1, b2)
title <- paste0("q = 1 / (1 + exp(-(", b1, " + ", b2, "x)))")
plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
for (x in xs) {
d <- distribution(x, l)
lines(ys, d(ys), type = "l")
points(ys, d(ys), pch = x)
}
legend("topright", legend = legends, pch = xs)
}
}

リンク関数はロジスティック関数です。 コードでは以下の部分です。

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

これが二項分布のパラメータの一つである事象の発生確率 `q` を与えます。 もう一つのパラメータは試行回数 `N` ですが、これは `20` にしました。 ということで、リンク関数のパラメータ `b1` と `b2` さえ与えれば、説明変数 `x` に対して応答変数 `y` の分布がどのように変化するかが分かります。 コードを実行すると様々な `b1` と `b2` の値に対してモデルをプロットします。 例えば下図のような図がプロットされます。 これは `b1 = -0.3`、`b2 = 0.3` の場合です。 `x` が `0` なら`q = 1 / (1 + exp(0.3)) = 0.425...` です。 `x` が `10` なら`q = 1 / (1 + exp(-2.7)) = 0.937...` です。 `b2` が正なので `x` の増加に従って `exp` の引数部分がどんどん小さくなって、`q` 全体としては `1` に近づくわけです。 事象の発生確率が高くなれば、`N` 回の試行でその事象が起きる回数は大きくなります。 図を見ると、確かに `x` の増加に従って `y` の分布が大きな値を取るようになっているのが分かります。