UPDATE: 2022-09-05 12:58:27
このノートではポアソン回帰のパラメタをスコア法で計算する方法についてまとめておく。
ポアソン回帰分析の詳しい説明はさておき、ポアソン分布は下記の通り定義される。
\[ f(Y_{i}=y_{i}) = \frac{exp[-\lambda_{i}]\lambda_{i}^{y_{i}}}{y_{i}!}, \lambda_{i}>0,y=0,1,2, \dots \]
まず、ここでは下記の通り定義する。
\[ \boldsymbol{X} = \begin{bmatrix} x_{11} & \cdots & x_{1p} \\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{np} \end{bmatrix}, \boldsymbol{x} = \begin{bmatrix} x_{i1} \\ \vdots \\ x_{ip} \end{bmatrix}, \boldsymbol{\beta} = \begin{bmatrix} \beta_{1} \\ \vdots \\ \beta_{p} \end{bmatrix}, \boldsymbol{ x_{i}^{T} \boldsymbol{\beta}} = x_{i1}\beta_{1} + \cdots + x_{ip}\beta_{p} \]
ポアソン回帰分析は説明変数と回帰係数パラメタの積に指数を取ったもので\(\lambda\)を計算する。
\[ f(Y_{i}=y_{i}) = \frac{exp[-\lambda_{i}]\lambda_{i}^{y_{i}}}{y_{i}!}, \lambda_{i}=exp[\boldsymbol{ x_{i}^{T} \boldsymbol{\beta}}] \]
ポアソン分布を最尤法で計算してパラメタを求めることを考える。
\[ L(\boldsymbol{\beta})=\prod \frac{exp[-\lambda_{i}]\lambda_{i}^{y_{i}}}{y_{i}!} \]
尤度関数を対数尤度関数に変換し、
\[ log L(\boldsymbol{\beta})=\sum (y_{i} \boldsymbol{ x_{i}^{T} \boldsymbol{\beta}} - exp[\boldsymbol{ x_{i}^{T} \boldsymbol{\beta}}] - log y_{i}!) \]
\(\beta\)について微分する。この結果を0とすることで、対数尤度関数\(logL(β)\)が\(β\)について最大化されるときの1階条件となる。
\[ \begin{eqnarray} \frac{ \partial log L(\boldsymbol{\beta}) }{ \partial \beta} &=& \sum (y_{i} - exp[\boldsymbol{ x_{i}^{T} \boldsymbol{\beta}}])\boldsymbol{ x_{i}} = 0 \\ &=& \boldsymbol{ X^{T}}[\boldsymbol{y} - exp[\boldsymbol{X} \boldsymbol{\beta}]]= \boldsymbol{0} \end{eqnarray} \]
各パラメタの要素を書き下すと、下記の通りである。\(\boldsymbol{g}\)は勾配のGradientを意味する。対数尤度関数の各パラメタの1階微分は勾配であるため。
\[ \boldsymbol{g(\beta)}= \begin{bmatrix} \frac{ \partial log L(\boldsymbol{\beta}) }{ \partial \beta_{1}} \\ \vdots \\ \frac{ \partial log L(\boldsymbol{\beta}) }{ \partial \beta_{p}} \end{bmatrix}= \begin{bmatrix} [y_{1} - exp[\boldsymbol{x}^{T}_{1} \boldsymbol{\beta}]x_{11} + \cdots + [y_{n} - exp[\boldsymbol{x}^{T}_{n}] \boldsymbol{\beta}]x_{1n} \\ \vdots \\ [y_{1} - exp[\boldsymbol{x}^{T}_{1} \boldsymbol{\beta}]x_{p1} + \cdots + [y_{n} - exp[\boldsymbol{x}^{T}_{n}] \boldsymbol{\beta}]x_{pn} \end{bmatrix}= \begin{bmatrix} 0 \\ \vdots \\ 0 \end{bmatrix} \]
これまでに扱ってきたニュートン法で計算することを考える。\(β=β_{k}\)として、\(k\)回目の更新で得られた値を意味する。
\[ \boldsymbol{g(\beta_{k})}= \begin{bmatrix} \frac{ \partial log L(\boldsymbol{\beta}_{k}) }{ \partial \beta_{1(k)}} \\ \vdots \\ \frac{ \partial log L(\boldsymbol{\beta}_{k}) }{ \partial \beta_{p(k)}} \end{bmatrix} \]
ニュートン法の計算に必要なヘシアン行列\(\boldsymbol{H(\beta)}\)は下記となる。
\[ \boldsymbol{H(\beta)}=\frac{ \partial \boldsymbol{g(\beta)} }{ \partial \boldsymbol{\beta}}= \begin{bmatrix} \frac{ \partial^{2} logL(\boldsymbol{\beta})}{ \partial \beta_{1} \partial \beta_{1}} & \cdots & \frac{ \partial^{2} logL(\boldsymbol{\beta})}{ \partial \beta_{1} \partial \beta_{p}} \\ \vdots & \ddots & \vdots \\ \frac{ \partial^{2} logL(\boldsymbol{\beta})}{ \partial \beta_{p} \partial \beta_{1}} & \cdots & \frac{ \partial^{2} logL(\boldsymbol{\beta})}{ \partial \beta_{p} \partial \beta_{p}} \end{bmatrix} \]
また、\(k\)は、\(k\)回目の更新で得られた値を意味する。
\[ \boldsymbol{H(\beta_{k})}= \begin{bmatrix} \frac{ \partial^{2} logL(\boldsymbol{\beta_{k}})}{ \partial \beta_{1(k)} \partial \beta_{1(k)}} & \cdots & \frac{ \partial^{2} logL(\boldsymbol{\beta_{k}})}{ \partial \beta_{1(k)} \partial \beta_{p(k)}} \\ \vdots & \ddots & \vdots \\ \frac{ \partial^{2} logL(\boldsymbol{\beta_{k}})}{ \partial \beta_{p(k)} \partial \beta_{1(k)}} & \cdots & \frac{ \partial^{2} logL(\boldsymbol{\beta_{k}})}{ \partial \beta_{p(k)} \partial \beta_{p(k)}} \end{bmatrix} \]
以上より、ニュートン法の反復計算は次のようになる。
\[ \boldsymbol{\beta}_{k+1} = \boldsymbol{\beta}_{k} - [\boldsymbol{H(\beta}_{k})]^{-1} \boldsymbol{g(\beta}_{k}) \]
ただ、ポアソン回帰分析では、ヘシアン行列を使わず、Fisherスコアリング法を用いる。詳しい話は一般化線形モデルの書籍を参照。
Fisherスコアリング法は、ヘシアン行列\(\boldsymbol{H(\beta)}\)の代わりに負の符号を持つフィッシャー情報行列\(\boldsymbol{I(\beta)}\)を用いて反復計算を行う。
\[ \boldsymbol{\beta}_{k+1} = \boldsymbol{\beta}_{k} + [\boldsymbol{I(\beta}_{k})]^{-1} \boldsymbol{g(\beta}_{k}) \]
フィッシャー情報行列\(\boldsymbol{I(\beta)}\)は下記の通り計算できる。
\[ \begin{eqnarray} \boldsymbol{I(\beta)} &=& - E[\boldsymbol{H(\beta)}] \\ &=& - E \left[ \frac{\partial^{2} logL(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}} \right] \\ &=& \sum \boldsymbol{ x_{i}^{T}exp[x_{i}^{T} \beta]x_{i}} \\ &=& \boldsymbol{ X^{T} W X} \\ &=& \begin{bmatrix} x_{11} & \cdots & x_{np} \\ \vdots & \ddots & \vdots \\ x_{1p} & \cdots & x_{np} \end{bmatrix} \begin{bmatrix} exp[x_{1}^{T} \beta] & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & exp[x_{n}^{T} \beta] \end{bmatrix} \begin{bmatrix} x_{11} & \cdots & x_{1p} \\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{np} \end{bmatrix} \end{eqnarray} \]
以上より、ポアソン回帰分析は、フィッシャー情報行列\(\boldsymbol{I(\beta)}\)を用いて繰り返してパラメタを更新することで、回帰係数を求めることができる。
ここからは回帰係数の分散に関する話。回帰係数の推定量\(\hat{ \beta }\)の漸近分布は、
\[ \boldsymbol{ \hat {\beta} } \sim N[\boldsymbol{ \beta }, [\boldsymbol{ I(\beta) }]^{-1}] \]
であり、回帰係数の推定量\(\hat{ \beta }\)̂の分散共分散行列フィッシャー情報行列より、下記となる。
\[ \boldsymbol{ V[\hat \beta] } = [\boldsymbol{ I(\beta) }]^{-1} = (\boldsymbol{ X^{T} W X} )^{-1} \]
分散共分散行列の対角要素が各パラメタの分散となり、平方根を計算すれば、標準誤差が得られる。
\[ \boldsymbol{ V[\hat \beta] } = \begin{bmatrix} V[\hat \beta_{1}] & \cdots & C[\hat \beta_{1}, \hat \beta_{p}] \\ \vdots & \ddots & \vdots \\ C[\hat \beta_{p}, \hat \beta_{1}] & \cdots & V[\hat \beta_{p}] \end{bmatrix} \]
以上より、回帰係数を標準化して、漸近分布が正規分布に従うことを利用して、回帰係数の検定を行う。
\[ z_{i} = \frac{\boldsymbol{\hat \beta}_{i}}{ \sqrt { V[\boldsymbol{\hat \beta}_{i}] } } \]
まずはglm関数を利用して、ポアソン回帰分析を実行する。サンプルデータは下記よりお借りする。
library(tidyverse)
d <- read.csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv") %>%
mutate(f = ifelse(f == 'C', 0, 1),
x0 = 1) %>%
select(y, x0, x1 = x, x2 = f)
head(d)## y x0 x1 x2
## 1 6 1 8.31 0
## 2 6 1 9.44 0
## 3 6 1 9.50 0
## 4 12 1 9.07 0
## 5 10 1 10.16 0
## 6 4 1 8.32 0計算したい回帰係数の値を予め確認しておく。この値が得られるであろうパラメタの数値である。
Y <- as.vector(d$y)
Xnms <- c('x0', 'x1','x2')
X <- as.matrix(d[, names(d) %in% Xnms])
summary(glm(Y ~ X - 1, family = poisson(link = 'log')))##
## Call:
## glm(formula = Y ~ X - 1, family = poisson(link = "log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3977 -0.7337 -0.2023 0.6795 2.4317
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Xx0 1.26311 0.36963 3.417 0.000633 ***
## Xx1 0.08007 0.03704 2.162 0.030620 *
## Xx2 -0.03200 0.07438 -0.430 0.667035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1946.276 on 100 degrees of freedom
## Residual deviance: 84.808 on 97 degrees of freedom
## AIC: 476.59
##
## Number of Fisher Scoring iterations: 4さきほどの数式の部分をRで実装して実行する。glm関数と同じようなパラメタが計算できている。
iter <- 50
tolerance <- 10^-5
B <- rep(1, ncol(X))
for(i in 1:iter){
# 更新前のパラメタを保存
pre_B <- B
# 偏回帰係数と説明変数の積で計算されるλハット
lam <- exp(X %*% B)
# 勾配の計算
grad <- t(X) %*% (Y - lam)
# フィッシャー情報行列の計算用λハットの対角行列
W <- diag(nrow(X))
diag(W) <- lam
# フィッシャー情報行列を計算
I <- t(X) %*% W %*% X
# フィッシャー情報行列の逆行列と勾配を使ってパラメタを更新
B <- pre_B + solve(I) %*% grad
# 計算結果を表示
cat(sprintf("%d times: (x0=%2.3f, x1=%2.3f, x2=%2.3f)\n", i, B[1], B[2], B[3]))
if(sqrt(crossprod(B - pre_B)) < tolerance) {
cat(sprintf("The algorithm converged.\n%d times: (x0=%2.3f, x1=%2.3f, x2=%2.3f)", i, B[1], B[2], B[3]))
break
}
}## 1 times: (x0=0.000, x1=1.000, x2=1.000)
## 2 times: (x0=-0.998, x1=1.000, x2=1.000)
## 3 times: (x0=-1.995, x1=1.000, x2=0.999)
## 4 times: (x0=-2.985, x1=0.999, x2=0.998)
## 5 times: (x0=-3.959, x1=0.997, x2=0.995)
## 6 times: (x0=-4.888, x1=0.992, x2=0.987)
## 7 times: (x0=-5.699, x1=0.978, x2=0.964)
## 8 times: (x0=-6.202, x1=0.941, x2=0.906)
## 9 times: (x0=-5.975, x1=0.849, x2=0.768)
## 10 times: (x0=-4.375, x1=0.653, x2=0.505)
## 11 times: (x0=-1.497, x1=0.361, x2=0.189)
## 12 times: (x0=0.655, x1=0.144, x2=0.011)
## 13 times: (x0=1.226, x1=0.084, x2=-0.030)
## 14 times: (x0=1.263, x1=0.080, x2=-0.032)
## 15 times: (x0=1.263, x1=0.080, x2=-0.032)
## 16 times: (x0=1.263, x1=0.080, x2=-0.032)
## The algorithm converged.
## 16 times: (x0=1.263, x1=0.080, x2=-0.032)これまでと同じく、optim関数でも同じ値が得られるかを確認しておく。最大化する対数尤度関数は下記のとおり。
\[ log L(\boldsymbol{\beta})=\sum (y_{i} \boldsymbol{ x_{i}^{T} \boldsymbol{\beta}} - exp[\boldsymbol{ x_{i}^{T} \boldsymbol{\beta}}] - log y_{i}!) \]
optimにこの対数尤度関数を渡し、フィッシャーのスコア法がないので、準ニュートン法で最適化を行う。
# lfactorial(n) = log(factorial(n))
loglik <- function(B) sum(Y * X %*% B - exp(X %*% B) - log(factorial(Y)))
B <- rep(1, ncol(X))
set.seed(1989)
res <-
optim(
par = B,
fn = loglik,
control = list(fnscale = -1), # 最大化
hessian = TRUE,
method = "BFGS" #準ニュートン法
)
# 標準誤差
se <- sqrt(diag(solve(-res$hessian)))
# summary風
cbind(coef = res$par,
se = se,
zval = res$par / se )## coef se zval
## [1,] 1.26499789 0.36961119 3.422510
## [2,] 0.07988111 0.03703513 2.156901
## [3,] -0.03189318 0.07437914 -0.428792