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

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


2017年 04月 21日

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの五回目です。 この章では様々なタイプのGLMについて説明がされています。 その中で、ロジスティック回帰についてAICを使ってモデル選択をする様子が説明されています。 なので、サンプルデータを用意して、いろいろなモデルについてAICを計算してくれるコードを用意しました。 コードはRで書きました。
N <- 10
dataSize <- 100
createData <- function(createY) {
createY <- match.fun(createY)
x <- runif(dataSize, 0, 20)
c <- runif(dataSize, 0, 20)
y <- createY(x)
data.frame(x, c, y)
}

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

data <- dataGenerator()

calcAIC <- function(formula) {
fit <- glm(formula, data = data, family = binomial)
fit$aic
}

cat("model\t\taic\n")
cat("null\t\t", calcAIC(cbind(y, N - y) ~ 1), "\n")
cat("x\t\t", calcAIC(cbind(y, N - y) ~ x), "\n")
cat("c\t\t", calcAIC(cbind(y, N - y) ~ c), "\n")
cat("x + c\t\t", calcAIC(cbind(y, N - y) ~ x + c), "\n")
cat("xc\t\t", calcAIC(cbind(y, N - y) ~ x : c), "\n")
cat("x + xc\t\t", calcAIC(cbind(y, N - y) ~ x + x : c), "\n")
cat("c + xc\t\t", calcAIC(cbind(y, N - y) ~ c + x : c), "\n")
cat("x + c + xc\t", calcAIC(cbind(y, N - y) ~ x * c), "\n")
サンプルデータは `x`、`c`、`y` を変数として含んでいます。 説明変数は `x` で応答変数は `y` です。 そこに `y` と相関のない変数として `c` を用意しました。 `c` は必要ない変数なので、モデルに組み込むとバイアスが大きくなるはずです。 コードを実行すると以下のような結果が出力されます。
model		aic
null		 885.6348
x			 264.6336
c			 887.2804
x + c		 265.7918
xc			 667.2655
x + xc		 265.8728
c + xc		 381.9791
x + c + xc	 267.7918
確かに `x` だけがモデルに使われている時に一番AICが低くなるのが分かります。