ポアソン回帰(カウントデータ)

ポアソン回帰はカウントデータあるいはイベントの発生率をモデル化する際に用いられる。このページではカウントデータをモデル化する例を示す。

サンプルデータはある感染症の新規感染者数を記録したデータである。感染者が確認された年を 1 期間とし、以降 12 期間(年間)についてその新規感染を記録している。

期間を x、各期間に対応する新規感染者数を y に代入する。

x <- 1:12
y <- c(6, 5, 14, 14, 21, 31, 38, 51, 86, 136, 169, 234)

誤差構造をポアソン分布とした時、期待値と分散はそれぞれ次のように求められる(計算略)。

\[P(y) = \frac{\lambda ^{y}}{y!}e^{-\lambda}\] \[E[Y] = Var[Y] = \lambda\]

また、デザイン行列を X、パラメーターを β、リンク関数を g とすると、モデルは以下のように構築される。

\[E[Y] = \lambda\] \[ g(\lambda) = \mathbf{x}\mathbf{\beta}= \begin{pmatrix} 1 & x\end{pmatrix}\begin{pmatrix}\beta_{1}\\ \beta_{2}\end{pmatrix} \]

ここで、感染者数(Y)はカウントデータであるため、マイナスの値を取ることはない。しかし、デザイン行列の値次第で上式の右辺はマイナスとなることもある。そのため、リンク関数の逆関数を、上の式の右辺を常にプラスにするような関数を選ぶ必要がある。ポアソン回帰では一般的に対数関数をリンク関数として用いる。

\[ g(\lambda) = \log(\lambda) = \mathbf{x}\mathbf{\beta}= \begin{pmatrix} 1 & x\end{pmatrix}\begin{pmatrix}\beta_{1}\\ \beta_{2}\end{pmatrix} \]

R でポアソン回帰を行ってみる。関数は glm 関数を用いる。familypoisson を指定する。

m <- glm(y ~ x, family = poisson(link = "log"))

summary 関数で解析結果の要約情報を見ることができる。

summary(m)
## Call:
## glm(formula = y ~ x, family = poisson(link = "log"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.14929  -0.48356  -0.02623   0.26379   1.44696  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.29894    0.14658   8.861   <2e-16 ***
## x            0.34870    0.01454  23.990   <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: 821.4922  on 11  degrees of freedom
## Residual deviance:   6.2131  on 10  degrees of freedom
## AIC: 75.058
## 
## Number of Fisher Scoring iterations: 4

Coefficients の項目をみると (Intercept) と x の値がそれぞれ 1.29894 と 0.34870 に対応している。また、構築されたモデルの式から、期間 x と λ の間は以下の関係を持つ。

\[ E[Y] = \lambda = \exp\left( \beta_{1} + \beta_{2}x \right) \]

これを利用して、ポアソン回帰線を散布図に書き入れると以下のようになる。

plot(x, y, xlim = c(0, 15), ylim = c(0, 400),
     xlab = "period", ylab = "number of infected people")
par(new = T)
curve(exp(1.29894 + 0.34870 * x), 
      xlim = c(0, 15), ylim = c(0, 400), xlab = "", ylab = "")

※ サンプルに利用したデータは参照 [4] で公開されている「HIV患者とAIDS患者報告数の年次推移」から取得したもので、ポアソン回帰の使い方を説明するために特定の期間(1985 年~ 1996 年)のデータだけを利用した。全期間のデータを用いるとポアソン回帰から大きく外れる。

References

  1. 久保 拓弥. 確率と情報の科学 データ解析のための統計モデリング入門 ― 一般化線形モデル・階層ベイズモデル・MCMC. 2012.
  2. Dobson AJ. An Introduction to Generalized Linear Models. Second Edition. 2002.
  3. Carolyn AJ. Poisson Regression or Regression of Counts (& Rates). PDF
  4. AIDS/STI -related database Japan. HIV患者とAIDS患者報告数の年次推移. AIDS/STI