ニュートン・ラフソン法によるポアソン回帰

以下のサンプルデータを利用して、ポアソン回帰を行い、そのパラメーターをニュートン・ラフソン法により推測する。

サンプルデータはポアソン分布に従うように乱数生成したものである。ある実験区画に入り込むある絶滅危惧種の個体数を計数したデータである。計数は 1 年間 4 回行い、これを 5 年間続けた。このデータを利用して、年が経つ毎に個体数がどのように増えたのかをポアソン回帰によりモデル化する。

計数    1回目 2回目 3回目 4回目
1 年目    6     9    11    11
2 年目    21   14    13    20
3 年目    38   36    23    29
4 年目    30   33    46    40
5 年目    64   48    54    44

右図は、横軸を経過年数、縦軸を観測値としてプロットした散布図である。観測値が小さい時、そのばらつきも小さい。逆に、観測値が大きくなると、そのばらつきも大きくなっていることが確認できる。

ポアソン分布

ポアソン回帰を行うとき、リンク関数(g)およびパラメーター(β1、β2)を以下のように設定する。

\[ g(\lambda) = \log\lambda = \beta_{1} + \beta_{2} \cdot year \]

観測値ベクトル y 、パラメーターベクトル β およびデザイン行列 X は以下のように書ける。

\[ \mathbf{y} = \begin{pmatrix} 6 \\ 9 \\ 11 \\ 11 \\ 21 \\ \vdots \\ 47 \end{pmatrix}, \mathbf{X} = \begin{pmatrix} 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ 1 & 2 \\ \vdots & \vdots \\ 1 & 5 \end{pmatrix}, \mathbf{\beta} = \begin{pmatrix} \beta_{1} \\ \beta_{2} \end{pmatrix} \]

ニュートン・ラフソン法を利用する際に上で挙げた観測値ベクトル、デザイン行列、パラメーターベクトルのほかに WU なども必要とする(参照 1, 参照 2)。

\[ \begin{eqnarray} \mathbf{U} &=& \sum _{i=1}^{n} \frac{(y_{i} - \mu_{i})\mathbf{x}^{T}_{i \cdot}}{g'(\lambda_{i})Var[y_{i}]} \\ &=& \sum _{i=1}^{n} (y_{i} - \mu_{i})\mathbf{x}^{T}_{i \cdot} \\ \mathbf{W} &=& \begin{pmatrix} \frac{1}{(g'(\mu_{1}))^{2} Var[y_{1}]} & 0 & \cdots & 0 \\ 0 & \frac{1}{(g'(\mu_{2}))^{2} Var[y_{2}]} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{1}{(g'(\mu_{n}))^{2} Var[y_{n}]} \\ \end{pmatrix} \\ &=& \begin{pmatrix} \lambda_{1} & 0 & \cdots & 0 \\ 0 & \lambda_{2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_{n} \\ \end{pmatrix} \end{eqnarray} \]

行列の準備はこれで終わり。次に観測値、デザイン行列およびパラメーターを R の変数に保存する。

Y <- c( 6,  9, 11, 11,
       21, 14, 13, 20,
       38, 36, 23, 29,
       30, 33, 46, 40,
       61, 48, 54, 44)
X <- cbind(
       rep(1, length(Y)),
       c(1, 1, 1, 1,
         2, 2, 2, 2,
         3, 3, 3, 3,
         4, 4, 4, 4, 
         5, 5, 5, 5))

次に、for 文を利用して 99 回ニュートン・ラフソン法を繰り返して実行する。パラメーターの初期値はそれぞれ 1 と 10 に設定する。全 100 セットのパラメーターは行列型の beta 変数に保存される。

beta <- matrix(0, ncol = 2, nrow = 100)
beta[1, ] <- c(1, 10)
for (m in 2:100) {
  lambda <- exp(X %*% beta[m - 1, ])
  W <- diag(lambda[, 1])
  XtWX <- t(X) %*% W %*% X
  U <- t(Y - lambda) %*% X
  beta[m, ] <- beta[m - 1, ] + solve(XtWX) %*% t(U)
}

推定されたパラメーターをプロットすると、約 50 回繰り返したあとに、パラメーターが収束したと確認できる。

ポアソン分布

具体的にはパラメーターは以下の値に収束した。

tail(beta)
##           [,1]      [,2]
## [95,]  2.09693 0.3806021
## [96,]  2.09693 0.3806021
## [97,]  2.09693 0.3806021
## [98,]  2.09693 0.3806021
## [99,]  2.09693 0.3806021
## [100,] 2.09693 0.3806021

次に、同じことを R の glm 関数で行うと、結果は以下のようになる。(Intercept) の予測値が 2.09693、X2 の予測値が 0.38060 となり、上で求めた値と同じになったことが確認できる。

poi.model <- glm(Y ~ X, family = poisson(link = "log"))
summary(poi.model)
## Call:
## glm(formula = Y ~ X, family = poisson(link = "log"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8959  -0.8895  -0.2677   0.7155   2.3056  
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.09693    0.12565   16.69   <2e-16 ***
## X1                NA         NA      NA       NA    
## X2           0.38060    0.03193   11.92   <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: 180.724  on 19  degrees of freedom
## Residual deviance:  25.271  on 18  degrees of freedom
## AIC: 130.12
## 
## Number of Fisher Scoring iterations: 4

References

  1. P.M.E.Altham, Statistical Laboratory, University of Cambridge. Introduction to Generalized Linear Modelling. 2011. PDF
  2. Dobson AJ. An Introduction to Generalized Linear Models. Second Edition. 2002.
  3. ニュートン法. Wikipedia