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

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


2017年 04月 14日

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、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` の分布が求まるので、あとはそれをプロットします。 コードを実行すると、下図のようなグラフがプロットされます。 logistic-regression.png 推定される `b2` は正の数なので、`x` が大きくなるにしたがって `y` の分布も大きな値を取るようになっています。