データ解析のための統計モデリング入門 GLMのモデル選択 読書メモ３

このブログ記事は『データ解析のための統計モデリング入門』（久保拓弥　著、岩波書店）という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第４章、GLMのモデル選択についてのまとめの三回目です。 この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 それを理解するために、ポアソン回帰のバイアスを見るためのコードを用意しました。 コードはRで書きました。

createData <- function() {
y <- rpois(10, lambda = 4)
data.frame(y)
}

logLikelihood <- function(data, b1) {
sum(log(dpois(data, lambda = exp(b1))))
}

printLogLikelihood <- function(data, i) {
cat(paste0("Input data is data", i, "\n"))
print(data[i]$y)

fit <- glm(formula = y ~ 1, family = poisson, data = data[i])
b1 <- fit$coefficients[1]

cat("Estimated b1 =", b1, "\n")
cat("Log likelihood for\n")
cat("data1\t\tdata2\t\tdata3\t\n")
cat(logLikelihood(data[1]$y, b1), "\t",
logLikelihood(data[2]$y, b1), "\t",
logLikelihood(data[3]$y, b1), "\n\n")
}

printData <- function(data, i) {
cat(paste0("data", i))
print(data[i]$y)
}

data <- c(createData(), createData(), createData())

for(i in 1:3) {
printData(data, i)
}
cat("\n")
for(i in 1:3) {
printLogLikelihood(data, i)
}

コードを実行すると、まずポアソン分布からサンプリングした `10` 個のデータからなるデータセットを3個作ります。 それぞれのデータに対してポアソン回帰をして、推定されたパラメータを使って他の2つのデータセットに対して対数尤度を計算します。 対数尤度を計算しているコードは以下のコードです。

logLikelihood <- function(data, b1) {
sum(log(dpois(data, lambda = exp(b1))))
}

`b1` が最尤推定されたパラメータです。 `data` が対数尤度を計算する対象のデータです。 平均 `exp(b1)` を持つポアソン分布からデータが得られる確率の対数を計算しています。 実行結果は以下のようになります。

data1 [1] 5 1 2 4 9 3 3 2 3 4
data2 [1] 7 4 3 3 3 8 4 6 3 1
data3 [1] 9 1 5 6 1 3 3 8 3 2

Input data is data1
[1] 5 1 2 4 9 3 3 2 3 4
Estimated b1 = 1.280934
Log likelihood for
data1		data2		data3
-20.59338 	 -21.43294 	 -24.32331

Input data is data2
[1] 7 4 3 3 3 8 4 6 3 1
Estimated b1 = 1.435085
Log likelihood for
data1		data2		data3
-21.04396 	 -20.95861 	 -24.00313

Input data is data3
[1] 9 1 5 6 1 3 3 8 3 2
Estimated b1 = 1.410987
Log likelihood for
data1		data2		data3
-20.91147 	 -20.97071 	 -23.99113

傾向として、ポアソン回帰に使ったデータセットへの尤度が一番大きくて、他のデータセットへの尤度は小さくなります。 時々は、他のデータセットへの尤度が大きくなります。 最尤推定は与えられたデータに最もよく当てはまるパラメータを返すので、まだ見ぬデータにもよく当てはまるかどうかは分からないということです。