第6章一般化線形モデルとモデル選択

この章の目標

• 指数型分布族とリンク関数の言葉で、線形回帰・ロジスティック回帰・ポアソン回帰を一つの枠組み（一般化線形モデル）として書ける
• 正準リンクの対数尤度が凹であることを証明し、ニュートン法が反復重み付き最小二乗法（IRLS）になることを導ける
• ロジスティック回帰の係数をオッズ比として正しく解釈し、逸脱度と尤度比検定でモデルを比べられる
• 完全分離のとき最尤推定量が存在しないことを証明し、実務での兆候と対処を説明できる
• 過学習の仕組みを式で説明し、リッジ・ラッソ・交差検証・AIC・BIC が何をしているかを説明できる

前提：第3章、第4章、第5章。6.3 節と 6.7 節では 01 微分積分学 第7章（テイラーの定理、コンパクト集合上の最大値）を使う。交差検証とデータ漏洩は 24 機械学習の数理 第1章でも詳しく扱う。

会員が来月解約するかどうか（0 か 1）、1 日に届く問い合わせの件数（0, 1, 2, …）を説明変数で予測したいとする。線形回帰をそのまま使うと、予測確率が $[0, 1]$ からはみ出したり、予測件数が負になったりする。しかも、確率 $\pi$ の 0/1 データの分散は $\pi(1 - \pi)$、ポアソン分布に従う件数の分散は平均に等しく、第5章の等分散の仮定 (A2) は初めから成り立たない。一般化線形モデル (generalized linear model, GLM) は、データの分布を指数型分布族から選び、平均を説明変数の一次式とリンク関数で結ぶことで、これらを線形回帰と同じ枠組みで扱う。推定は最尤法で行い、その計算は第5章の最小二乗法の繰り返しに帰着する。

本章の後半では「説明変数をいくつ使うか」を考える。第5章で見たように、変数を増やすほど手元のデータへの当てはまりは良くなるが、新しいデータの予測は悪くなりうる（過学習）。これを抑える正則化（リッジ・ラッソ）と、モデルを比べる道具（交差検証・AIC・BIC）が、それぞれ何を推定・近似しているのかを正確に述べる。

6.1 指数型分布族

定義 6.1（指数型分布族, exponential family）$\Theta \subset \mathbb{R}$ を開区間、$\phi > 0$ とする。確率関数または密度が

$$f(y; \theta, \phi) = \exp\left(\frac{y\theta - b(\theta)}{\phi} + c(y, \phi)\right) \qquad (\theta \in \Theta)$$

の形の分布の族を（分散パラメータつきの）指数型分布族という。ここで $b$ は $\Theta$ 上の $C^2$ 級関数、$c$ は $\theta$ によらない関数で、各 $\theta \in \Theta$ について $f(\cdot; \theta, \phi)$ は確率関数（密度）である。$\theta$ を自然パラメータ (natural parameter)、$\phi$ を分散パラメータ (dispersion parameter) という。

例 6.2

1. 正規分布 $N(\mu, \sigma^2)$：$\log f = \frac{y\mu - \mu^2/2}{\sigma^2} - \frac{y^2}{2\sigma^2} - \frac{1}{2}\log(2\pi\sigma^2)$ なので、$\theta = \mu$、$b(\theta) = \theta^2/2$、$\phi = \sigma^2$。
2. ベルヌーイ分布 $B(1, \pi)$：$f(y) = \pi^y(1 - \pi)^{1 - y} = \exp\bigl(y\log\frac{\pi}{1 - \pi} + \log(1 - \pi)\bigr)$（$y = 0, 1$）なので、$\theta = \log\frac{\pi}{1 - \pi}$、$b(\theta) = \log(1 + e^{\theta})$、$\phi = 1$。試行回数 $m$ が既知の二項分布 $B(m, \pi)$ も、$b(\theta) = m\log(1 + e^{\theta})$ として同じ形に書ける。
3. ポアソン分布 $\operatorname{Po}(\lambda)$：$f(y) = \exp(y\log\lambda - \lambda - \log y!)$ なので、$\theta = \log\lambda$、$b(\theta) = e^{\theta}$、$\phi = 1$。

命題 6.3（平均と分散）$Y$ が定義 6.1 の分布に従うとき、$E[Y] = b'(\theta)$、$\operatorname{Var}(Y) = \phi b''(\theta)$ である。

証明. $\Theta$ は開区間なので、$\lvert t \rvert$ が小さければ $\theta + t\phi \in \Theta$ である。このとき（離散型なら積分を和に読み替えて）

である。モーメント母関数が $0$ の近くで有限なので、第1章 定理 1.15 より $E[Y] = M_Y'(0) = b'(\theta)$、$E[Y^2] = M_Y''(0) = \phi b''(\theta) + b'(\theta)^2$。$\square$

例 6.2 では、正規分布で $b'(\theta) = \mu$、$\phi b'' = \sigma^2$、ベルヌーイ分布で $b'(\theta) = e^{\theta}/(1 + e^{\theta}) = \pi$、$b''(\theta) = \pi(1 - \pi)$、ポアソン分布で $b' = b'' = \lambda$ となり、よく知られた値と一致する。分布が 1 点に集中しなければ $b'' > 0$ なので、$\mu = b'(\theta)$ は $\theta$ の狭義単調増加関数で、$\theta$ と平均 $\mu$ は 1 対 1 に対応する。分散を平均の関数として $\operatorname{Var}(Y) = \phi V(\mu)$ と書いたとき、$V$ を分散関数という（正規 $V = 1$、ベルヌーイ $V(\mu) = \mu(1 - \mu)$、ポアソン $V(\mu) = \mu$）。

6.2 一般化線形モデル

定義 6.4（一般化線形モデル）$Y_1, \dots, Y_n$ は独立で、$Y_i$ は自然パラメータ $\theta_i$、共通の分散パラメータ $\phi$ の指数型分布族に従うとする。平均 $\mu_i = E[Y_i]$ と説明変数 $x_i \in \mathbb{R}^p$ が、狭義単調で微分可能な関数 $g$ によって

$$g(\mu_i) = \eta_i, \qquad \eta_i = x_i^{\top}\beta$$

と結ばれるモデルを一般化線形モデルという。$\eta_i$ を線形予測子 (linear predictor)、$g$ をリンク関数 (link function) という。$g = (b')^{-1}$、すなわち $\theta_i = \eta_i$ となるリンク関数を正準リンク (canonical link) という。

ロジスティック回帰では、$\Lambda(t) = 1/(1 + e^{-t})$（ロジスティック関数）として $\pi_i = P(Y_i = 1) = \Lambda(x_i^{\top}\beta)$ であり、予測確率は必ず $(0, 1)$ に入る。ポアソン回帰では $\mu_i = e^{x_i^{\top}\beta}$ は必ず正である。件数を観測期間や人数 $t_i$ あたりの率でモデル化したいときは、$\log\mu_i = \log t_i + x_i^{\top}\beta$ と係数 $1$ の項（オフセット）を加える。一般化線形モデルは「$y$ を変換して線形回帰する」こととは違う。$\log y$ を線形回帰すると $E[\log Y]$ をモデル化することになり、$\log E[Y]$ とは一致しないうえ、$y = 0$ を扱えない。

6.3 対数尤度とその凹性

正準リンクでは $\theta_i = x_i^{\top}\beta$ なので、対数尤度は

である。特に

がロジスティック回帰とポアソン回帰の対数尤度である。$m_i$ 回の試行のうち $y_i$ 回成功した集計データ（$Y_i \sim B(m_i, \pi_i)$）では、$\sum_i (y_i\eta_i - m_i\log(1 + e^{\eta_i}))$ に定数を加えたものになる。以下、集計データでは $b$ を $m_i\log(1 + e^{\theta})$ に読み替えれば、議論はそのまま成り立つ。

命題 6.5（スコアとヘッセ行列）正準リンクのとき、$\mu_i = b'(\eta_i)$、$W = \operatorname{diag}\bigl(b''(\eta_1), \dots, b''(\eta_n)\bigr)$ とおくと

$$\nabla\ell(\beta) = \frac{1}{\phi}X^{\top}(y - \mu), \qquad \nabla^2\ell(\beta) = -\frac{1}{\phi}X^{\top}WX$$

証明. 連鎖律より、$\beta \mapsto b(x_i^{\top}\beta)$ の勾配は $b'(\eta_i)x_i$、ヘッセ行列は $b''(\eta_i)x_ix_i^{\top}$ である。これらを足し合わせればよい。$\square$

命題 6.3 より $\phi W = \operatorname{diag}(\operatorname{Var}(Y_i))$ である。ヘッセ行列が $y$ を含まないので、観測されたフィッシャー情報量と期待値は一致する。尤度方程式 $X^{\top}(y - \hat{\mu}) = 0$ は、第5章の正規方程式 $X^{\top}e = 0$ の類似であり、定数項があれば $\sum_i \hat{\mu}_i = \sum_i y_i$ となる（問題 6.3）。

定理 6.6（正準リンクの対数尤度の凹性）正準リンクの一般化線形モデルで、$\Theta = \mathbb{R}$（正規・ベルヌーイ・二項・ポアソン分布はこれをみたす）とし、各 $Y_i$ の分布は 1 点に集中しないとする。このとき $\ell$ は $\mathbb{R}^p$ 上の凹関数である。さらに $\operatorname{rank} X = p$ ならば $\ell$ は狭義凹で、最尤推定量は存在すればただ一つであり、それは尤度方程式 $\nabla\ell(\beta) = 0$ の解と一致する。

証明. 命題 6.5 と $b'' > 0$ より、任意の $\beta, h$ について $h^{\top}\nabla^2\ell(\beta)h = -\frac{1}{\phi}\sum_i b''(x_i^{\top}\beta)(x_i^{\top}h)^2 \leq 0$ で、$Xh \neq 0$ なら $< 0$。テイラーの定理（01 第7章 定理 7.21）より、ある $\tau \in (0, 1)$ について

であり、$Xh \neq 0$ なら不等号は狭義である。$\beta_0, \beta_1 \in \mathbb{R}^p$、$0 < t < 1$、$\beta_t = (1 - t)\beta_0 + t\beta_1$ とする。(1) を $\beta = \beta_t$ と $h = \beta_0 - \beta_t$、$h = \beta_1 - \beta_t$ に使い、それぞれ $1 - t$ 倍、$t$ 倍して足すと、$(1 - t)(\beta_0 - \beta_t) + t(\beta_1 - \beta_t) = 0$ より

を得る。$\operatorname{rank} X = p$ かつ $\beta_0 \neq \beta_1$ なら $X(\beta_0 - \beta_t) = tX(\beta_0 - \beta_1) \neq 0$ なので不等号は狭義である。最大点が 2 つ $\beta_0 \neq \beta_1$ あれば、中点での値が最大値より大きくなって矛盾する。最後に、$\nabla\ell(\hat{\beta}) = 0$ なら (1) より任意の $h$ で $\ell(\hat{\beta} + h) \leq \ell(\hat{\beta})$ なので $\hat{\beta}$ は最大点であり、逆に最大点では勾配が $0$ である（01 の命題 7.23）。$\square$

正準リンク以外では、ヘッセ行列に $y_i - \mu_i$ を含む項が加わり、対数尤度は一般には凹とは限らない（プロビットリンク $g = \Phi^{-1}$ のように、別の理由で凹になる例はある）。定理 6.6 は最大点の存在は保証しない。存在しない典型例が 6.7 節の完全分離である。

6.4 ニュートン法と反復重み付き最小二乗法

$\ell$ を最大化するにはニュートン法を使う。現在の点 $\beta$ のまわりで $\ell$ を 2 次のテイラー多項式 $\ell(\beta) + \nabla\ell(\beta)^{\top}h + \frac{1}{2}h^{\top}\nabla^2\ell(\beta)h$ で近似し、$\nabla^2\ell(\beta)$ が負定値ならその最大点 $h = -\nabla^2\ell(\beta)^{-1}\nabla\ell(\beta)$ へ進む：

定理 6.7（ニュートン法は反復重み付き最小二乗法）正準リンクで $\operatorname{rank} X = p$ とする。$\beta^{(t)}$ での線形予測子・平均・重みを $\eta^{(t)} = X\beta^{(t)}$、$\mu_i^{(t)} = b'(\eta_i^{(t)})$、$W^{(t)} = \operatorname{diag}(b''(\eta_i^{(t)}))$ とし、作業応答 (working response)

$$z^{(t)} = \eta^{(t)} + (W^{(t)})^{-1}(y - \mu^{(t)})$$

を定めると、ニュートン法の更新は

$$\beta^{(t+1)} = (X^{\top}W^{(t)}X)^{-1}X^{\top}W^{(t)}z^{(t)}$$

であり、これは $\sum_i w_i^{(t)}(z_i^{(t)} - x_i^{\top}\beta)^2$（$w_i^{(t)}$ は $W^{(t)}$ の対角成分）を最小にする重み付き最小二乗法の解である。

証明. 命題 6.5 より（$\phi$ は約分されて）$\beta^{(t+1)} = \beta^{(t)} + (X^{\top}WX)^{-1}X^{\top}(y - \mu)$（添字 $(t)$ を省く）。$X^{\top}WX\beta^{(t)} = X^{\top}W\eta^{(t)}$ なので

一方 $\sum_i w_i(z_i - x_i^{\top}\beta)^2 = \lVert W^{1/2}z - W^{1/2}X\beta \rVert^2$ で、$W^{1/2}X$ の階数は $p$ なので、第5章の定理 5.4（正規方程式）より最小点は $(X^{\top}WX)^{-1}X^{\top}Wz$ である。$\square$

重みと作業応答を更新しながら最小二乗法を解き直すので、この計算法を反復重み付き最小二乗法 (iteratively reweighted least squares, IRLS) という。作業応答は $g(y_i)$ を $\mu_i$ のまわりで 1 次近似したもの $\eta_i + g'(\mu_i)(y_i - \mu_i)$（正準リンクでは $g'(\mu) = 1/b''(\theta)$）で、重み $w_i$ はその分散の逆数に比例する。ロジスティック回帰では $w_i = \pi_i(1 - \pi_i)$（集計データでは $m_i\pi_i(1 - \pi_i)$）、ポアソン回帰では $w_i = \mu_i$ である。正準リンク以外でも、ヘッセ行列をその期待値で置き換えたフィッシャーのスコア法は、$w_i = 1/(V(\mu_i)g'(\mu_i)^2)$、$z_i = \eta_i + g'(\mu_i)(y_i - \mu_i)$ とする同じ形の反復になる（計算は省略する）。

ニュートン法は最大点の近くでは局所的に 2 次の速さで収束するが（$-\ell$ の最小化とみて 23 最適化 第6章 定理 6.4）、遠くから始めると収束は保証されないので、実装では対数尤度が増えるまで更新幅を半分にするなどの工夫をする（同章 6.3 節の減衰ニュートン法）。推定値の精度については、適当な正則条件のもとで $n \to \infty$ のとき $\hat{\beta}$ は近似的に $N_p\bigl(\beta, \phi(X^{\top}WX)^{-1}\bigr)$ に従う（第3章 定理 3.32 の最尤推定量の漸近正規性を、多次元の、独立だが同一分布でない観測に拡張したもの。証明は省略する）。そこで $W$ を $\hat{\beta}$ で評価した $\operatorname{SE}(\hat{\beta}_j) = \sqrt{\phi[(X^{\top}\hat{W}X)^{-1}]_{jj}}$ を標準誤差とし、$z_j = \hat{\beta}_j/\operatorname{SE}(\hat{\beta}_j)$ を標準正規分布と比べる（ワルド検定）。

例 6.8（クーポンと購入）ある通販サイトで、メールに付けるクーポンの額 $x$（100 円単位で $0, 1, \dots, 5$）を変え、既存顧客（$s = 0$）と新規顧客（$s = 1$）にそれぞれ各額 200 通ずつ送って、購入した人数 $y$ を数えた（説明用の架空のデータ）。

12 群の購入数を $Y_i \sim B(200, \pi_i)$、$\log\frac{\pi_i}{1 - \pi_i} = \beta_0 + \beta_1x_i + \beta_2s_i$ とモデル化し、IRLS で当てはめる。

import numpy as np

x = np.tile(np.arange(6.0), 2)          # クーポン額（100 円単位）
s = np.repeat([0.0, 1.0], 6)            # 0 = 既存顧客, 1 = 新規顧客
m = np.full(12, 200.0)                  # 各群の送付数
y = np.array([6, 10, 19, 20, 18, 27, 16, 15, 26, 23, 27, 41], dtype=float)  # 購入数

X = np.column_stack([np.ones(12), x, s])
beta = np.zeros(3)
for it in range(1, 7):
    eta = X @ beta
    pi = 1 / (1 + np.exp(-eta))
    w = m * pi * (1 - pi)               # 重み（各群の購入数の分散）
    z = eta + (y - m * pi) / w          # 作業応答
    beta = np.linalg.solve(X.T @ (w[:, None] * X), X.T @ (w * z))  # 重み付き最小二乗
    print(it, beta.round(6))

pi = 1 / (1 + np.exp(-X @ beta))
w = m * pi * (1 - pi)
se = np.sqrt(np.diag(np.linalg.inv(X.T @ (w[:, None] * X))))
print("標準誤差", se.round(4), " z 値", (beta / se).round(2))

出力：

1 [-1.872381  0.082286  0.16    ]
2 [-2.683023  0.172323  0.334791]
3 [-2.998649  0.222944  0.431219]
4 [-3.033762  0.2292    0.442531]
5 [-3.034123  0.229266  0.442641]
6 [-3.034123  0.229266  0.442641]
標準誤差 [0.1632 0.041  0.1374]  z 値 [-18.6    5.59   3.22]

係数の変化の最大値は 4〜6 回目で $0.035$、$0.00036$、$4 \times 10^{-8}$ とほぼ 2 乗のオーダーで小さくなり、5 回でほぼ収束している。推定値は $\hat{\beta} = (-3.034, 0.2293, 0.4426)$、標準誤差は $(0.163, 0.0410, 0.137)$ である。尤度方程式から、購入数の当てはめ値の合計は既存顧客・新規顧客それぞれで実際の合計（100 人、148 人）に一致する。

6.5 オッズ比の解釈

確率 $\pi$ に対し $\pi/(1 - \pi)$ をオッズ (odds) という。ロジスティック回帰は対数オッズ $\log\frac{\pi}{1 - \pi}$ を $x^{\top}\beta$ とおくモデルなので、ほかの説明変数を固定して $x_j$ を $1$ 増やすと、対数オッズは $\beta_j$ 増え、オッズは $e^{\beta_j}$ 倍になる。$e^{\beta_j}$ をオッズ比 (odds ratio) という。この倍率がほかの説明変数の値によらないのはモデルの仮定であり、交互作用を入れれば変わりうる。

例 6.8 では、クーポン 100 円あたりのオッズ比は $e^{0.2293} = 1.258$（95% 信頼区間は $e^{0.2293 \pm 1.96 \times 0.0410}$ より $[1.161, 1.363]$）、新規顧客の既存顧客に対するオッズ比は $e^{0.4426} = 1.557$（$[1.189, 2.038]$）である。購入確率でいえば、既存顧客ではクーポンなしの 4.6% から 500 円で 13.1% に、新規顧客では 7.0% から 19.1% になる。確率の変化 $\partial\pi/\partial x_j = \beta_j\pi(1 - \pi)$ は $\pi$ によって変わり、$\pi(1 - \pi) \leq 1/4$ より $\lvert \beta_j \rvert/4$ を超えない。

6.6 逸脱度と尤度比検定

各観測に別々のパラメータを与え $\hat{\mu}_i = y_i$ とするモデルを飽和モデル (saturated model) といい、その対数尤度を $\ell_{\text{sat}}$ と書く。

定義 6.9（逸脱度, deviance）$D = 2\phi\bigl(\ell_{\text{sat}} - \ell(\hat{\beta})\bigr)$ をモデルの逸脱度という。

定義から次のように計算できる（$0\log 0 = 0$ とする）。

• 正規分布：$D = \sum_i (y_i - \hat{\mu}_i)^2 = \mathrm{RSS}$。逸脱度は残差平方和の一般化である。
• ポアソン分布：$D = 2\sum_i \bigl(y_i\log(y_i/\hat{\mu}_i) - (y_i - \hat{\mu}_i)\bigr)$。
• 二項分布（集計データ）：$D = 2\sum_i \Bigl(y_i\log\frac{y_i}{m_i\hat{\pi}_i} + (m_i - y_i)\log\frac{m_i - y_i}{m_i(1 - \hat{\pi}_i)}\Bigr)$。ベルヌーイ分布（$m_i = 1$）では $\ell_{\text{sat}} = 0$ なので $D = -2\ell(\hat{\beta})$。

計画行列 $X_0$ のモデル $M_0$ が $X_1$ のモデル $M_1$ に含まれる（$\mathcal{C}(X_0) \subset \mathcal{C}(X_1)$、パラメータ数 $p_0 < p_1$）とき、$\ell_{\text{sat}}$ が打ち消し合って

となり、逸脱度の差は尤度比検定（第4章 定義 4.13）の統計量そのものである。$\phi$ が既知で $M_0$ が正しければ、ウィルクスの定理（第4章 定理 4.14。ここでも適当な正則条件のもとで、独立で同一分布でない観測への拡張を認める）により、$n \to \infty$ で近似的に $\chi^2(p_1 - p_0)$ に従う。正規分布で $\sigma^2$ が未知なら、第5章 定理 5.16 の $F$ 検定が正確な検定を与える。

例 6.8 では、逸脱度はモデル全体で $D = 7.47$、新規顧客の項を除いたモデルで $18.03$、定数項だけのモデルで $50.45$ である。「新規か既存かで購入確率は変わらない」（$\beta_2 = 0$）の尤度比統計量は $18.03 - 7.47 = 10.56$（自由度 1、$p = 0.0012$）で、ワルド統計量 $z^2 = 3.22^2 = 10.38$ とほぼ同じ結論になる。

逸脱度の使い方には注意がいる。

• 値そのものはデータの持ち方に依存する。 例 6.8 を 2400 人分の 0/1 データとして当てはめると、推定値は同じだが飽和モデルが変わり、逸脱度は $1552.29$ になる。入れ子のモデルの差（$10.56$）は変わらない。
• 適合度の検定：集計データで各 $m_i$ が大きい（ポアソンなら各 $\mu_i$ が大きい）とき、モデルが正しければ $D$ は近似的に $\chi^2(n - p)$ に従う。例 6.8 では $D = 7.47$（自由度 9、$p = 0.59$）で、当てはまりの悪さを示す証拠はない。しかし 0/1 データ（$m_i = 1$）では $n$ が大きくてもこの近似は成り立たず（飽和モデルのパラメータが $n$ とともに増え、ウィルクスの定理の前提が崩れる）、$D$ を適合度の検定に使ってはいけない。
• 過分散：ポアソン回帰・二項回帰は分散が平均で決まる（$\phi = 1$）と仮定している。件数が大きいのに $D/(n - p)$ が 1 よりずっと大きければ、実際の分散がモデルより大きい（過分散）疑いがあり、標準誤差は過小になる。$\phi$ を推定する方法（擬似尤度）や負の二項分布のモデルを検討する。

6.7 完全分離と最尤推定量の非存在

ロジスティック回帰の当てはめで、係数と標準誤差が異常に大きくなり、計算が収束しないことがある。原因は説明変数が 0/1 を「分けすぎる」ことにある。以下、0/1 データ（ベルヌーイ分布）で、$s_i = 2y_i - 1 \in \lbrace -1, 1 \rbrace$ とおく。

定義 6.10（分離）$s_ix_i^{\top}v > 0$ がすべての $i$ で成り立つ $v \in \mathbb{R}^p$ が存在するとき、データは完全分離 (complete separation) しているという（$y_i = 1$ の点と $y_i = 0$ の点が超平面 $x^{\top}v = 0$ で完全に分かれる）。完全分離ではないが、$s_ix_i^{\top}v \geq 0$ がすべての $i$ で成り立つ $v \neq 0$ が存在するとき、準完全分離 (quasi-complete separation) しているという。

定理 6.11（最尤推定量の存在）ロジスティック回帰で $\operatorname{rank} X = p$ とする。

1. 完全分離していれば、すべての $\beta$ で $\ell(\beta) < 0$ かつ $\sup_{\beta}\ell(\beta) = 0$ であり、最尤推定量は存在しない。
2. 最尤推定量が存在するための必要十分条件は、完全分離も準完全分離もしていない（任意の $v \neq 0$ について $s_ix_i^{\top}v < 0$ となる $i$ がある）ことである。存在すればただ一つである。

証明. $1 - \Lambda(t) = \Lambda(-t)$ より、$y_i = 1$ の項 $\log\Lambda(x_i^{\top}\beta)$ と $y_i = 0$ の項 $\log(1 - \Lambda(x_i^{\top}\beta))$ はどちらも $\log\Lambda(s_ix_i^{\top}\beta)$ と書けて、

となる。各項は負で、$s_ix_i^{\top}\beta$ について単調増加である。

(1) $\ell < 0$ は明らか。完全分離を与える $v$ について、$u \to \infty$ のとき各項 $\log\Lambda(u \cdot s_ix_i^{\top}v) \to 0$ なので $\ell(uv) \to 0$。よって上限 $0$ は達成されない。

(2) 必要性：ある $v \neq 0$ ですべての $s_ix_i^{\top}v \geq 0$ なのに最大点 $\hat{\beta}$ があるとする。各項の単調性より $\ell(\hat{\beta} + v) \geq \ell(\hat{\beta})$ なので $\hat{\beta} + v$ も最大点となり、定理 6.6 の一意性に反する。十分性：単位球面 $S = \lbrace v \mid \lVert v \rVert = 1 \rbrace$ 上の連続関数 $F(v) = \max_i(-s_ix_i^{\top}v)$ は仮定より正の値をとる。$S$ はコンパクトなので $F$ は正の最小値 $c$ をもつ（01 第7章 定理 7.7）。$\beta \neq 0$ なら $v = \beta/\lVert \beta \rVert$ について $-s_ix_i^{\top}\beta \geq c\lVert \beta \rVert$ となる $i$ があり、その項だけを残して $\ell(\beta) \leq \log\Lambda(-c\lVert \beta \rVert)$。$\Lambda(-cR) = 2^{-n}$ となる $R \geq 0$ をとると、$\lVert \beta \rVert > R$ なら $\ell(\beta) < -n\log 2 = \ell(0)$。よって $\ell$ の最大値は閉球 $\lbrace \lVert \beta \rVert \leq R \rbrace$ 上での最大値に等しく、これはコンパクト集合上の連続関数なので存在する（01 の定理 7.7）。一意性は定理 6.6 による。$\square$

完全分離のとき IRLS を回すと、係数の比はほぼ一定のまま大きさだけが増え続け、対数尤度は $0$ に近づく（問題 6.5）。ソフトウェアは収束しない、あるいは「当てはめ確率が 0 か 1 になった」と警告し、巨大な係数と標準誤差を出力する。これは推定の失敗であって、「その変数の効果が非常に大きい」という発見ではない。

6.8 過学習と正則化

第5章の線形回帰で、説明変数を増やしたときに何が起こるかを式で見る。真の平均 $\mu = E[Y]$ は $\mathcal{C}(X)$ に入っていなくてもよいとし、誤差は (A1)(A2) をみたすとする。

命題 6.12（訓練誤差の楽観性）$Y = \mu + \varepsilon$、$\hat{y} = HY$（$H$ は $p$ 列の計画行列 $X$ のハット行列）とし、$Y' = \mu + \varepsilon'$ を同じ説明変数での新しい観測（$\varepsilon'$ は $\varepsilon$ と独立で同じ分布）とすると、

$$E\lVert Y - \hat{y} \rVert^2 = \lVert (I_n - H)\mu \rVert^2 + (n - p)\sigma^2, \qquad E\lVert Y' - \hat{y} \rVert^2 = \lVert (I_n - H)\mu \rVert^2 + (n + p)\sigma^2$$

である。特に、新しい観測に対する誤差は訓練誤差より平均して $2p\sigma^2$ 大きい。

証明. $Y - \hat{y} = (I_n - H)\mu + (I_n - H)\varepsilon$ で、交差項の期待値は $0$、$E\lVert (I_n - H)\varepsilon \rVert^2 = \sigma^2\operatorname{tr}(I_n - H) = (n - p)\sigma^2$（第5章 定理 5.12 と同じ計算）。また $Y' - \hat{y} = (I_n - H)\mu + \varepsilon' - H\varepsilon$ で、$\varepsilon'$ と $\varepsilon$ は独立で平均 $0$ なので交差項の期待値は $0$、$E\lVert \varepsilon' \rVert^2 = n\sigma^2$、$E\lVert H\varepsilon \rVert^2 = \sigma^2\operatorname{tr}H = p\sigma^2$。$\square$

$\lVert (I_n - H)\mu \rVert^2$（偏り）は列を増やすと減るが、$p\sigma^2$（分散）は増える。予測誤差は両者の和なので、変数を増やしすぎると悪化する。一方、訓練誤差はこの増加を $-p\sigma^2$ として逆向きに数えるので、変数を増やすほど良く見える。これが過学習 (overfitting) である。命題 6.12 から $\mathrm{RSS} + 2p\sigma^2$ は $E\lVert Y' - \hat{y} \rVert^2$ の不偏推定量であり、これを $\sigma^2$ で割って定数を調整したものをマローズの $C_p$ という。同じ考え方を尤度に一般化したのが 6.10 節の AIC である（同じ結果は 24 第1章 命題 1.18 でも、バイアス–バリアンス分解の例として扱う）。

過学習を抑える一つの方法は、係数が大きくなることに罰則を課す正則化 (regularization) である。以下、説明変数は標準化（平均 $0$・分散 $1$）し、$y$ も中心化して定数項を除いておく（定数項には罰則を課さない）。

命題 6.13（リッジ回帰, ridge regression）$\lambda > 0$ について、$\lVert y - X\beta \rVert^2 + \lambda\lVert \beta \rVert^2$ の最小点はただ一つで、$\hat{\beta}_{\lambda} = (X^{\top}X + \lambda I_p)^{-1}X^{\top}y$ である（$X$ の階数は問わない）。$X = UDV^{\top}$（$U$ は列が正規直交な $n \times r$ 行列、$D = \operatorname{diag}(d_1, \dots, d_r)$、$d_j > 0$、$V$ は列が正規直交な $p \times r$ 行列、$r = \operatorname{rank} X$）を特異値分解とすると

$$X\hat{\beta}_{\lambda} = \sum_{j=1}^r \frac{d_j^2}{d_j^2 + \lambda}u_ju_j^{\top}y$$

証明. $X$ の下に $\sqrt{\lambda}I_p$ を、$y$ の下に $0 \in \mathbb{R}^p$ を付け加えて

とおくと、目的関数は $\lVert \tilde{y} - \tilde{X}\beta \rVert^2$ で、$\tilde{X}$ の階数は $p$ である。第5章の定理 5.4 より最小点はただ一つで、$(\tilde{X}^{\top}\tilde{X})^{-1}\tilde{X}^{\top}\tilde{y} = (X^{\top}X + \lambda I_p)^{-1}X^{\top}y$。後半は、$V$ を $\mathbb{R}^p$ の正規直交基底に延長して計算すると $(X^{\top}X + \lambda I_p)^{-1}X^{\top} = V(D^2 + \lambda I_r)^{-1}DU^{\top}$ となるので、$X\hat{\beta}_{\lambda} = UD^2(D^2 + \lambda I_r)^{-1}U^{\top}y$。$\square$

最小二乗法の当てはめ値は $\sum_j u_ju_j^{\top}y$（$\mathcal{C}(X)$ への射影）なので、リッジ回帰は各特異ベクトル方向の成分を $d_j^2/(d_j^2 + \lambda)$ 倍に縮める。縮み方は特異値の小さい方向ほど大きい。特異値の小さい方向とは説明変数どうしがほぼ一次従属になる方向で、第5章で見た多重共線性で最小二乗推定量の分散が大きくなる方向である。偏りを入れる代わりに分散を減らすことで、平均二乗誤差では最小二乗法に勝つことができる。

定理 6.14（リッジ推定量の平均二乗誤差）(A1)(A2) を仮定し、$\operatorname{rank} X = p$ とする。$X = UDV^{\top}$（$V$ は $p$ 次直交行列）、$\alpha = V^{\top}\beta$ とおくと、$0 < \lambda \leq \sigma^2/\max_j\alpha_j^2$ をみたすすべての $\lambda$ について（$\beta = 0$ ならすべての $\lambda > 0$ について）

$$E\lVert \hat{\beta}_{\lambda} - \beta \rVert^2 < E\lVert \hat{\beta} - \beta \rVert^2$$

である（$\hat{\beta}$ は最小二乗推定量）。

証明. $V$ は直交行列なので $\lVert \hat{\beta}_{\lambda} - \beta \rVert = \lVert V^{\top}\hat{\beta}_{\lambda} - \alpha \rVert$ で、命題 6.13 の証明より $V^{\top}\hat{\beta}_{\lambda}$ の第 $j$ 成分は $d_ju_j^{\top}Y/(d_j^2 + \lambda)$。$E[u_j^{\top}Y] = u_j^{\top}UDV^{\top}\beta = d_j\alpha_j$、$\operatorname{Var}(u_j^{\top}Y) = \sigma^2$ より、この成分の平均二乗誤差（分散と偏りの 2 乗の和）は