ロジスティック回帰モデルと一般化線形モデル

数理統計　第9章 §9.4　★★ 指数型分布族・リンク関数・スコアと Fisher 情報・Newton–Raphson（IRLS）・オッズ比の推測

ロジスティック回帰

応答が成功か失敗かの2値 $Y_i\in\{0,1\}$ のとき、$Y_i$ を正規分布でモデル化するのは不自然である。平均 $\pi_i=P(Y_i=1)$ を $\bm{x}_i^\top\bm{\beta}$ と等しいとおくと、確率の範囲である0以上1以下からはみ出す。そこで $\pi_i$ を対数オッズに変換したものを線形にする。

$$\log\frac{\pi_i}{1-\pi_i}=\bm{x}_i^\top\bm{\beta},\qquad\text{すなわち}\quad \pi_i=\frac{\mathrm{e}^{\bm{x}_i^\top\bm{\beta}}}{1+\mathrm{e}^{\bm{x}_i^\top\bm{\beta}}} .\tag{9.13}$$

これをロジスティック回帰モデルと呼ぶ。$x_{ij}$ が1増えると、対数オッズが $\beta_j$ 増える。よって $\mathrm{e}^{\beta_j}$ はオッズ比である。

$Y_i$ が独立にベルヌーイ分布 $\mathrm{Bin}(1,\pi_i)$ に従うとき、対数尤度は

$$\ell(\bm{\beta})=\sum_{i=1}^n\{Y_i\log\pi_i+(1-Y_i)\log(1-\pi_i)\}=\sum_{i=1}^n\{Y_i\eta_i-\log(1+\mathrm{e}^{\eta_i})\},\qquad \eta_i=\bm{x}_i^\top\bm{\beta} .$$

次の2つの微分を使う。

$$\frac{\mathrm{d} }{\mathrm{d} \eta}\log(1+\mathrm{e}^{\eta})=\pi,\qquad \frac{\mathrm{d} \pi}{\mathrm{d} \eta}=\pi(1-\pi) .$$

すると、スコアとヘッセ行列は

$$\frac{\partial \ell}{\partial \bm{\beta}}=\sum_{i=1}^n(Y_i-\pi_i)\bm{x}_i=X^\top(\bm{Y}-\bm{\pi}),\qquad \frac{\partial^2\ell}{\partial\bm{\beta}\,\partial\bm{\beta}^\top}=-\sum_{i=1}^n\pi_i(1-\pi_i)\bm{x}_i\bm{x}_i^\top=-X^\top WX\tag{9.14}$$

となる。$W=\mathrm{diag}(\pi_i(1-\pi_i))$ である。ヘッセ行列は $Y$ を含まないので、Fisher 情報行列はそのまま $I_n(\bm{\beta})=X^\top WX$ である。これは正定値なので $\ell$ は狭義凹である。

最尤推定値は、存在すればただ1つで、尤度方程式 $X^\top(\bm{Y}-\bm{\pi})=\bm{0}$ の解である。ただし $\pi_i$ が $\bm{\beta}$ の非線形関数なので、一般には陽に解けない。

一般化線形モデル

正規線形モデルとロジスティック回帰は、次の枠組みの特別な場合である。

証明

対数密度の $\theta$ 微分は $(y-b'(\theta))/\phi$ である。スコアの期待値は0である（$\int f\,\mathrm{d} y=1$ を $\theta$ で微分して得られる、第6章 §6.3）。よって $E[Y]-b'(\theta)=0$。さらに Fisher 情報量の2通りの表し方 $E[(\partial_\theta\log f)^2]=-E[\partial_\theta^2\log f]$ を使うと

$$\frac{E[(Y-b'(\theta))^2]}{\phi^2}=\frac{b''(\theta)}{\phi},\qquad\text{すなわち}\quad V[Y]=\phi b''(\theta) .$$

主な例を並べる。

• 正規分布 $\mathrm{N}(\mu,\sigma^2)$：正準リンクは恒等写像。

$$\theta=\mu,\qquad b(\theta)=\theta^2/2,\qquad \phi=\sigma^2 .$$

• ベルヌーイ分布：正準リンクはロジット（確認テスト 9.4-1）。

$$\theta=\log\{\pi/(1-\pi)\},\qquad b(\theta)=\log(1+\mathrm{e}^\theta),\qquad \phi=1 .$$

• ポアソン分布 $\mathrm{Po}(\mu)$：$f(y)=\exp\{y\log\mu-\mu-\log y!\}$ より、正準リンクは対数。

$$\theta=\log\mu,\qquad b(\theta)=\mathrm{e}^\theta,\qquad \phi=1 .$$

正準リンクのときは $\theta_i=\bm{x}_i^\top\bm{\beta}$ なので、対数尤度は

$$\ell(\bm{\beta})=\frac1\phi\sum\{Y_i\bm{x}_i^\top\bm{\beta}-b(\bm{x}_i^\top\bm{\beta})\}+(\text{定数})$$

となる。よって

$$\frac{\partial \ell}{\partial \bm{\beta}}=\frac1\phi X^\top(\bm{Y}-\bm{\mu}),\qquad -\frac{\partial^2\ell}{\partial\bm{\beta}\,\partial\bm{\beta}^\top}=\frac1\phi X^\top WX,\quad W=\mathrm{diag}(b''(\theta_i)) .\tag{9.16}$$

ロジスティック回帰の式 (9.14) は、この $\phi=1$、$b''=\pi(1-\pi)$ の場合である。正規線形モデルでは $W=I$ で、尤度方程式は正規方程式に一致する。

Newton–Raphson 法と反復重みつき最小二乗

尤度方程式を Newton–Raphson 法で解く。現在の値 $\bm{\beta}^{(t)}$ のまわりでスコアを1次近似して $\bm{0}$ とおくと

$$\bm{\beta}^{(t+1)}=\bm{\beta}^{(t)}+(X^\top WX)^{-1}X^\top(\bm{Y}-\bm{\mu})$$

を得る（$W,\bm{\mu}$ は $\bm{\beta}^{(t)}$ で評価する）。これは次のように書き直せる。

$$\bm{\beta}^{(t+1)}=(X^\top WX)^{-1}X^\top W\bm{z},\qquad \bm{z}=X\bm{\beta}^{(t)}+W^{-1}(\bm{Y}-\bm{\mu}) .$$

すなわち、作業変数 $\bm{z}$ を重み $W$ で重みつき最小二乗回帰する操作の繰り返しになっている。この見方から反復重みつき最小二乗法（IRLS）と呼ばれる。正準リンクでは観測された情報行列と Fisher 情報行列が一致するので、Newton–Raphson 法と Fisher のスコア法は同じ反復になる。

推測：Wald 検定・尤度比検定・逸脱度

第6章で示した最尤推定量の漸近正規性により、$n$ が大きいとき $\hat{\bm{\beta}}$ は近似的に次に従う。

$$\hat{\bm{\beta}}\approx\mathrm{N}_p\bigl(\bm{\beta},\ \phi(X^\top WX)^{-1}\bigr).$$

$W$ を $\hat{\bm{\beta}}$ で評価する。この行列の対角成分から、標準誤差 $\mathrm{se}(\hat\beta_j)$ を作る。$z=\hat\beta_j/\mathrm{se}(\hat\beta_j)$ を $\mathrm{N}(0,1)$ と比べるのが Wald 検定である。

入れ子のモデルの比較には、第7章 §7.3 の尤度比検定を使う。データに完全に合わせた飽和モデル（$\hat\mu_i=Y_i$）の対数尤度を $\ell_{\mathrm{sat}}$ として、

$$D=2\{\ell_{\mathrm{sat}}-\ell(\hat{\bm{\beta}})\}$$

を逸脱度（deviance）と呼ぶ。正規線形モデル（$\phi=\sigma^2$）では $D=S_e/\sigma^2$ で、逸脱度は残差平方和の一般化である。入れ子のモデルで、母数が $q$ 個少ない $H_0$ を考える。そのもとでは、逸脱度の差 $D_0-D$（尤度比検定統計量）が漸近的に $\chi^2_q$ に従う。

確認テスト

オリジナル

1. 9.4-1 ベルヌーイ分布 $\mathrm{Bin}(1,\pi)$ の確率関数を (9.15) の形に書き、$\theta,\ b(\theta),\ \phi$ を求めよ。命題 9.11 から $E[Y]=\pi$、$V[Y]=\pi(1-\pi)$ を確かめよ。
2. 9.4-2 ポアソン回帰 $\log\mu_i=\bm{x}_i^\top\bm{\beta}$（$Y_i\sim\mathrm{Po}(\mu_i)$ 独立）を考える。$X$ が切片の列を含むなら、最尤推定値で $\sum\hat\mu_i=\sum Y_i$ が成り立つことを示せ。
3. 9.4-3 2群のポアソン回帰 $\log\mu=\beta_0+\beta_1x$ を考える。$x=0$ の10単位の合計件数は30、$x=1$ の10単位の合計件数は48であった。$\hat\beta_0,\ \hat\beta_1,\ \mathrm{se}(\hat\beta_1)$ を求め、$H_0:\beta_1=0$ を Wald 検定（有意水準 5%）で判定せよ。
4. 9.4-4 切片だけのロジスティック回帰（$n$ 個の観測、標本平均 $\bar y$）を考える。初期値 $\beta^{(0)}=0$ から Newton–Raphson 法を1回行った値を $\beta^{(1)}$ とする。これを $\bar y$ で表せ。$\bar y=0.7$ のとき $\beta^{(1)}$ と最尤推定値を比べよ。