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

サンプルデータは xcy を変数として含んでいます。
説明変数は 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が低くなるのが分かります。