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

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

N <- 10
dataSize <- 100

createData <- function(createY) {
createY <- match.fun(createY)
x <- runif(dataSize, 0, 20)
y <- createY(x)
data.frame(x, y)
}

dataGenerator <- function()
createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.5 * x))))

data <- dataGenerator()

fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

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

xs <- seq(5, 20, by = 3)
ys <- 0:N
title <- "Logistic regression"
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
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) {
q <- linkFunction(x)
y <- dbinom(ys, N, q)
lines(ys, y, type="l")
points(ys, y, pch = x)
}
legend("topleft", legend = legends, pch = xs)

作ったサンプルデータにロジスティック回帰をやって、リンク関数のパラメータを推定しているのは以下の部分です。

fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

こうして得られた `b1` と `b2` を使って以下のリンク関数を作れば、説明変数 `x` に対して、事象の発生確率 `q` が得られます。

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

これで事象を `N` 回繰り返した時の応答変数 `y` の分布が求まるので、あとはそれをプロットします。 コードを実行すると、下図のようなグラフがプロットされます。 推定される `b2` は正の数なので、`x` が大きくなるにしたがって `y` の分布も大きな値を取るようになっています。