ポアソン回帰で推定しているモノはλの式

某所の(1)ポアソン回帰モデルの説明が、(2)対数変換OLSと同じになっている気がします。違うものだと思うのですが、シミュレーションをして(1)と(2)の推定をして確認してみました。

1. モデル

ポアソン分布はパラメーター $\lambda$ で決定されるわけですが、 $\lambda$ を説明変数で説明するモデルになります。 $n$ 個のパラメータがあり、 $x$ を説明変数、 $\beta$ を係数として、以下のような式ですね。
$\lambda = e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n}$
$P(y=k) = \frac{\lambda^k e^{-\lambda}}{k!}$
被説明変数 $y$ の値が0以上の整数 $k$ のときの確率を、 $x$ は間接的に決定するわけですね。教科書的には最尤法を用いて求めることになるみたいですが、実用的にはリンク関数を用いて一般化線形回帰モデルで推定できるようです。

2. データ作成

まずはポアソン回帰モデル用のデータxとyを作成します。

set.seed(20130919)
x <- round(runif(100, max=3))
lambda <- exp(1.1 + 1*x)
y <- numeric(length(x))
for(i in 1:length(x)){
y[i] <- rpois(1, lambda[i])
}

3. Rで推定してみる

glm関数を使うだけですけど、推定できます。

r_pois <- glm(y ~ x, family=poisson)
summary(r_pois)

4. 対数変換OLSと比較してみる

比較のために、対数変換OLSをかけてみましょう*1。

r_ols <- lm(log(y) ~ x)
summary(r_ols)

当たり前ですが、残差平方和も異なります。

sum((y - fitted.values(r_pois))^2)
sum((y - exp(fitted.values(r_ols)))^2)

5. over-dispersionのテスト

tの切片項が0を棄却されたら、つまり有意に0以外ならばover-dispersionでポアソン分布でないことになります。

t <- fitted.values(r_pois)
z <- ((y - t)^2 - y)/(t*sqrt(2))
summary(lm(z ~ 0 + t))

人工的なデータなのですが、P値が0.323でした。

*1:実は対数化できるようにyはゼロ以上の数字になっています。

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。