1.緑本におけるofset項の説明
いわゆる緑本におけるポアソン回帰の説明において、「割り算を使わないためのoffset項について」（6.6 割算値の統計モデリングはやめよう)、次の説明があります。
（前提条件：私はあまり数学は得意ではありません。大学は理系でしたが数学はほとんど使いませんでした）
ここでは、Oは調査地iにおける観測値、Aは調査地iにおける面積とします。
密度Diは、
$$
\frac{O_i}{A_i}=D_i
$$
で表せます。
密度は正の値をとるので、調査地iに係る説明変数xを使ったモデリングをすると、上記を変形し、指数関数と組み合わせると
$$
O_i=A_i×D_i=A_iexp(\beta_1+\beta x_i)
$$
となり、これをさらに変形して
$$
O_i=exp(\beta_1+\beta_2 x_i + logA_i)
$$
となり、このlogAiをオフセット項と呼び、割り算の代わりに使用すると本の中では説明されています。しかし、いきなりlogがでてきますし、この変形の意味が、今ひとつ自分の頭では腑に落ちませんでした。
2.offset項の意味を考える
まずポアソン回帰のモデルを示します
$$
Y_{i} \sim Poisson(\lambda)
$$
ここでλはポアソン分布の期待値を表しています。このλを上記の定義の密度で表します。
$$
\lambda=\frac{O_i}{A_i}
$$
ポアソン回帰における期待値は正の値をとるので、この関係をgの関数で表すと、
$$
exp(\lambda)=g
$$
両辺を対数化して
$$
\lambda=log(g)
$$
これをLink関数と呼んでいます。数学的に常に正の値を与えるために指数関数を使っています。確率のように0-1の範囲に入る場合はロジスティック関数を使うようです。
たいがいLINK関数の説明はここで終わっているのだけど、なぜ指数関数を選ぶのかというところは、単に常に正になるということを示すのに数学的に扱いやすいというだけの理由だと思っています。
ここでLink関数を使って、先ほどの調査地iに係る説明変数xを使った線形モデリングを考えます。
$$
log(g_i)=log(\frac{O_i}{A_i})=\beta_1+\beta_2 x_i
$$
対数の割り算は引き算になるから、次のように変形できます。
$$
logO_i-logA_i=\beta_1+\beta_2 x_i
$$
右辺にlogAiを移動させると
$$
logO_i=logA_i+\beta_1+\beta_2 x_i
$$
ここで、初めて緑本で説明されたオフセット項がでてきました。
本当にこれで割り算の変わりになっているのでしょうか？
簡単に検算して確かめてみます。同じモデリングだとxの値は面積の値と反比例するはずです。
ここで、係数をそれぞれ、b1=1 b2=1とします。
Oが10で面積Aがそれぞれ1,2とします。
log(10)=log(1)+1+x
x=log(10)-log(1)-1
x=1.3
log(10)=log(2)+1+x
x=log(10)-log(2)-1
x=0.61
A=1の時、x=1.3, A=2の時、x=0.61となりました。
ただし、LINK関数を使ったときに対数化したので、指数化してスケールを元に戻す必要があります。
exp(1.3) =3.669297
exp(0.61)=1.840431
ほぼ面積に反比例しました。確かに割り算になっているようです。
3.シミュレーションでの解析事例
例
A地区とB地区では、ある種の雑草の変異が異なるか調べたい。
調査はそれぞれ50回実施したが、10a当たりの調査となったので、
調査本数が異なる。（調査数の単位は100とします。変異種の単位は1です）
A地区とB地区では、変異は異なっているのだろうか？
#A地区の単位当たり変異種の密度は8
denA=rpois(50,8)
area_1=floor(runif(50,10,100))
obs1=denA*area_1
#B地区の単位当たり変異種の密度は5
denB=rpois(50,5)
area_2=floor(runif(50,30,300))
obs2=denB*area_2
factor = rep(c("A","B"),each=50)
dat = data.frame("変異種"=c(obs1,obs2),"調査数"=c(area_1,area_2),"場所"=factor)
dat$場所=as.factor(dat$場所)
調査結果は一つのデータフレームに収納する。
A地区とB地区の結果の一部を示す。
head(dat,5)
tail(dat,5)
ポアソン回帰(GLM)で分析する。調査数はlog()でオフセット項とする
fit = glm( 変異種~ 場所, offset = log(調査数), family=poisson(link = "log"),data=dat)
summary(fit)
>Call:
>glm(formula = 変異種 ~ 場所, family = poisson(link = "log"),
> data = dat, offset = log(調査数))
>
>Deviance Residuals:
> Min 1Q Median 3Q Max
>-19.0265 -6.7443 -0.4323 4.5847 29.4774
>
>Coefficients:
> Estimate Std. Error z value Pr(>|z|)
>(Intercept) 2.098947 0.006611 317.50 <2e-16 ***
>場所B -0.507206 0.008523 -59.51 <2e-16 ***
>---
>Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>(Dispersion parameter for poisson family taken to be 1)
>
> Null deviance: 12261.4 on 99 degrees of freedom
>Residual deviance: 8871.5 on 98 degrees of freedom
>AIC: 9663.6
>
>Number of Fisher Scoring iterations: 4
結果は、有意にA地区とB地区は異なりAの方が変異割合が高い
それぞれの単位密度を推定すると
exp(fit$coefficients[1])
exp(fit$coefficients[1]+fit$coefficients[2])
>(Intercept)
> 8.157576
>(Intercept)
> 4.912296
上がA,下がBの単位当たりの変異数の推定量になる
真の値はAが8、Bが5なので、ほぼ正確に推定できている
（調査数の単位は100なので、100本当たりの変異出現密度）
4.まとめ
ポアソン回帰における割り算の代りに用いるofset項は、link関数を用いることで対数化し、割り算を引き算の関係にすることの意味でした。
（まあ、数学が得意な人には何を当たり前のことを言っているのかと言われそうですが...）
なにか数学的に素人なので間違いがあればご指摘いただけたらありがたいです。
5.参考
久保拓哉(2012).データ解析のための統計モデリング入門-一般化線形モデル・階層ベイズモデル・MCMC.岩波書店
Richard McElreath(2020).Statistical Rethinking 2nd Edition