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が低くなるのが分かります。