TIM Labs

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

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第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)

サンプルデータに含まれる変数は xAy です。 説明変数は xA です。 応答変数は y です。 ポアソン回帰のリンク関数は λ = exp(b1 + b2 * x + log(A)) という形です。 なので、log(A) 部分がオフセット項になります。 本の例で言うと A は面積のような、y との割算で密度を計算できるような値です。

本の中では1000打数300安打の打者と10打数3安打の打者を同じ三割打者と見なしてよいかという問題に触れられています。 サンプルデータの変数名で言うと、A = 1000y = 300 の場合と、A = 10y = 3 の場合とを同じと見なしてよいかという問題です。 感覚的には、10回中の3回打つ打者より、1000回中の300回打つ打者の方がより正確に三割打者であると言えます。 なので、これらの状況には何らかの差が出てきて欲しいと思うところです。

この問題の状況を見るために、ポアソン回帰を2つの場合について行いました。 ケース1においては A1 から 10 までのランダムな値を取ります。 ケース2においては A100 から 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 の大きさの違いが推定値のばらつきに影響を与えているのが分かります。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/615

コメントする

このブログ記事について

このページは、sato.yuichiroが2017年4月25日 12:00に書いたブログ記事です。

ひとつ前のブログ記事は「Chainer:GPUを使うプログラムはどう書くか」です。

次のブログ記事は「カラー画像のデータセットを探し求めて」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。

月別 アーカイブ