データ解析のための統計モデリング入門 GLMの応用範囲をひろげる 読書メモ3
2017年 04月 14日
このブログ記事は『データ解析のための統計モデリング入門』（久保拓弥 著、岩波書店）という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。
今回は第６章、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 の分布も大きな値を取るようになっています。
＃人気の記事
＃タグ
