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

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


2017年 04月 25日

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの六回目です。 この章では様々なタイプのGLMについて説明がされています。 その中でオフセット項わざという手法について説明がされているので、サンプルデータに対してこの手法を試してみるコードを用意しました。 コードはRで書きました。
dataSize <- 100
createData <- function(n, createY) {
createY <- match.fun(createY)
x <- runif(dataSize, 0, 5)
A <- runif(dataSize, n, 10 * n)
y <- createY(x, A)
data.frame(x, A, y)
}

linkFunction <- function(x, A) rpois(dataSize, exp(0.1 + 0.3 * x + log(A)))
dataGenerator1 <- function() createData(1, linkFunction)
dataGenerator2 <- function() createData(100, linkFunction)

doGlm <- function(dataGen) {
dataGen <- match.fun(dataGen)
data <- dataGen()
fit <- glm(y ~ x, offset = log(A), data = data, family = poisson)
c(fit$coefficients[1], fit$coefficients[2])
}

showResult <- function(dataGen) {
dataGen <- match.fun(dataGen)
m <- replicate(100, doGlm(dataGen))
meanm <- apply(m, 1, mean)
sdm <- apply(m, 1, sd)
cat("mean b1", meanm[1], "\n")
cat("mean b2", meanm[2], "\n")
cat("sd b1", sdm[1], "\n")
cat("sd b2", sdm[2], "\n")
}

cat("case1\n")
showResult(dataGenerator1)
cat("\ncase2\n")
showResult(dataGenerator2)
サンプルデータに含まれる変数は `x`、`A`、`y` です。 説明変数は `x` と `A` です。 応答変数は `y` です。 ポアソン回帰のリンク関数は `λ = exp(b1 + b2 * x + log(A))` という形です。 なので、`log(A)` 部分がオフセット項になります。 本の例で言うと `A` は面積のような、`y` との割算で密度を計算できるような値です。 本の中では1000打数300安打の打者と10打数3安打の打者を同じ三割打者と見なしてよいかという問題に触れられています。 サンプルデータの変数名で言うと、`A = 1000` で `y = 300` の場合と、`A = 10` で `y = 3` の場合とを同じと見なしてよいかという問題です。 感覚的には、10回中の3回打つ打者より、1000回中の300回打つ打者の方がより正確に三割打者であると言えます。 なので、これらの状況には何らかの差が出てきて欲しいと思うところです。 この問題の状況を見るために、ポアソン回帰を2つの場合について行いました。 ケース1においては `A` は `1` から `10` までのランダムな値を取ります。 ケース2においては `A` は `100` から `1000` までのランダムな値を取ります。 `x` はどちらの場合も `0` から `5` までのランダムな値です。 `y` は `λ = exp(0.1 + 0.3 * x + log(A))` つまり `λ = A * exp(0.1 + 0.3 * x)` というリンク関数で得られる `λ` を平均に持つポアソン分布からサンプリングされます。 なので、`A` が100倍になれば、`y` の平均も100倍になります。 コードを実行すると次のような結果が得られます。
case1
mean b1 0.09576121
mean b2 0.3013721
sd b1 0.05874369
sd b2 0.01892569

case2
mean b1 0.1002007
mean b2 0.299965
sd b1 0.006423389
sd b2 0.00181942
推定されるパラメータの平均値はどちらも大差がありませんが、その標準偏差は10倍近い差があります。 確かに `A` の大きさの違いが推定値のばらつきに影響を与えているのが分かります。