データ解析のための統計モデリング入門 GLMの応用範囲をひろげる 読書メモ6
2017年 04月 25日
このブログ記事は『データ解析のための統計モデリング入門』（久保拓弥 著、岩波書店）という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。
今回は第６章、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回打つ打者の方がより正確に三割打者であると言えます。
なので、これらの状況には何らかの差が出てきて欲しいと思うところです。
この問題の状況を見るために、ポアソン回帰を２つの場合について行いました。
ケース１においては A は 1 から 10 までのランダムな値を取ります。
ケース２においては 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 の大きさの違いが推定値のばらつきに影響を与えているのが分かります。
＃人気の記事
＃タグ
