データ解析のための統計モデリング入門 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]

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

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

これで事象を N 回繰り返した時の応答変数 y の分布が求まるので、あとはそれをプロットします。
コードを実行すると、下図のようなグラフがプロットされます。

logistic-regression.png

推定される b2 は正の数なので、x が大きくなるにしたがって y の分布も大きな値を取るようになっています。