7. GLMのパラメータ推定：IRLSと最尤法

GLMにおける最尤推定の枠組み

一般化線形モデルのパラメータ推定は、指数型分布族に属する確率分布の対数尤度関数を最大化する問題として定式化されます。観測値 $y_i$ が指数型分布族に従うとき、その確率密度関数（または確率質量関数）は次の標準形で表されます。

$$
f(y_i; \theta_i, \phi) = \exp\left\{ \frac{y_i \theta_i – b(\theta_i)}{\phi} + c(y_i, \phi) \right\}
$$

$\theta_i$ は自然パラメータ（正準パラメータ）、$b(\theta_i)$ は対数分配関数、$\phi$ は分散パラメータです。対数分配関数 $b(\theta)$ は平均と分散の生成関数としての役割を担い、$b'(\theta_i) = \mu_i$ によって正準パラメータと平均パラメータが結びつけられます。また $b”(\theta_i)$ は分散関数 $V(\mu_i)$ に対応します。観測値 $y_i$ は $\theta_i$ に関する十分統計量を構成します。

$n$ 個の独立な観測に基づく対数尤度関数は次のように書かれます。

$$
\ell(\beta) = \sum_{i=1}^{n} \left\{ \frac{y_i \theta_i – b(\theta_i)}{\phi} \right\} + \text{const}
$$

一般化線形モデルでは、平均 $\mu_i = E[y_i]$ とリンク関数 $g$ を通じて線形予測子 $\eta_i = \mathbf{x}_i^T \beta$ が $g(\mu_i) = \eta_i$ の関係で結びつけられます。$\ell(\beta)$ の $\beta_j$ に関する偏微分は連鎖則によって次のように展開されます。

$$
\frac{\partial \ell}{\partial \beta_j} = \sum_{i=1}^{n} \frac{y_i – \mu_i}{\phi \cdot V(\mu_i)} \cdot \frac{\partial \mu_i}{\partial \eta_i} \cdot x_{ij}
$$

この正規方程式は $\beta$ に関して非線形であり、正規分布と恒等リンクの組み合わせ以外では閉形式解が得られないため、反復的な数値解法が不可欠です。

前提条件として、全観測値が互いに独立であること、および指定した分布族がデータ生成過程と合致していること（モデル正定値）が必要です。分布族を誤特定した場合、最尤推定量は一致性を失う可能性があります。

スコア方程式とFisher情報行列

最尤推定の1階条件は、対数尤度 $\ell(\beta)$ の $\beta$ に関する勾配ベクトル、すなわちスコア関数 $U(\beta)$ をゼロと置くことで得られます。

$$
U(\beta) = \frac{\partial \ell(\beta)}{\partial \beta} = \mathbf{X}^T \mathbf{D} \boldsymbol{\Sigma}^{-1} (\mathbf{y} – \boldsymbol{\mu})
$$

ここで $\mathbf{D} = \text{diag}(\partial \mu_i / \partial \eta_i)$、$\boldsymbol{\Sigma} = \text{diag}(\phi V(\mu_i))$ です。最尤推定量 $\hat{\beta}$ はスコア方程式

$$
U(\hat{\beta}) = \mathbf{0}
$$

の解として定義されます。正則条件（積分と微分の交換可能性）の下では $E_\beta[U(\beta)] = \mathbf{0}$ が成立します。

期待Fisher情報行列は、スコアの外積の期待値として次のように定義されます。

$$
\mathcal{I}(\beta) = E\left[ U(\beta) U(\beta)^T \right] = -E\left[ \frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^T} \right]
$$

これは観測情報行列 $\hat{\mathcal{I}}(\beta) = -\partial^2 \ell(\beta) / \partial \beta \partial \beta^T$ と区別されます。前者は $\beta$ の真値における期待値であるのに対し、後者は実現した標本に基づく量です。Cramér-Rao下界の理論から、正則条件下の不偏推定量の共分散行列は $\mathcal{I}(\beta)^{-1}$ に対して半正定値の意味で下から抑えられます。Fisher情報行列は漸近信頼区間の構築においても中心的な役割を果たします。

前提条件として、正則条件（積分と微分の交換可能性）と対数尤度の二階微分可能性が必要です。完全多重共線性が存在する場合、$\mathcal{I}(\beta)$ が特異になり逆行列が存在しなくなります。有限標本では観測情報行列と期待情報行列は一般に一致せず、どちらを用いるかによって標準誤差の推定値が異なる点に注意が必要です。

Newton-Raphson法とFisher Scoring

スコア方程式を数値的に解く基本的な手法がNewton-Raphson法です。現在の推定値 $\beta_t$ に対する更新式は次の形をとります。

$$
\beta_{t+1} = \beta_t – \left[ H(\beta_t) \right]^{-1} U(\beta_t)
$$

$H(\beta_t) = \partial^2 \ell / \partial \beta \partial \beta^T |_{\beta_t}$ はHessian行列です。指数型分布族では対数尤度が $\beta$ に関して凹であるため、Hessianは負定値となります。この性質は、対数尤度の $\beta_t$ 周りの二次テイラー展開

$$
\ell(\beta) \approx \ell(\beta_t) + U(\beta_t)^T (\beta – \beta_t) + \frac{1}{2}(\beta – \beta_t)^T H(\beta_t)(\beta – \beta_t)
$$

を最大化することに対応しており、更新方向は対数尤度曲面の局所的な曲率情報を用いた二次近似に基づきます。

Fisher Scoringは、HessianをFisher情報行列の期待値 $-\mathcal{I}(\beta_t)$ で置き換えた手法です。

$$
\beta_{t+1} = \beta_t + \mathcal{I}(\beta_t)^{-1} U(\beta_t)
$$

正準リンク関数を使用する場合は $H(\beta_t) = -\mathcal{I}(\beta_t)$ が代数的に成立するため、Newton-Raphson法とFisher Scoringは厳密に同一の更新式となります。非正準リンクでは両者は一致しませんが、期待値をとることで各反復における数値的安定性が向上します。この安定化の効果は、Hessianが確率的に変動する標本依存の量である点に起因します。

前提条件として、対数尤度関数の局所的凹性（指数型分布族と正準リンクの組み合わせで保証）と、初期値が収束域内に存在することが必要です。$p \times p$ 行列の反転を要するため1反復あたり $O(p^3)$ の計算コストが発生し、高次元では計算負荷が大きくなります。初期値が大幅に離れている場合、局所解への収束や発散のリスクがあります。

IRLS（反復重み付き最小二乗法）の導出

IRLS（反復重み付き最小二乗法）は、Fisher Scoringを加重最小二乗問題として等価に書き換えた実装形式です。各反復において、調整従属変数（作業応答）$z_i$ を次のように定義します。

$$
z_i = \hat{\eta}_i + (y_i – \hat{\mu}_i)\left.\frac{d\eta}{d\mu}\right|_{\hat{\mu}_i}
$$

ここで $\hat{\eta}_i = \mathbf{x}_i^T \beta_t$ は現在の線形予測子の値です。$z_i$ は観測値 $y_i$ と現在の平均推定値 $\hat{\mu}_i$ の差分をリンク関数の微分で線形化したものであり、現在の推定点における接線近似として解釈されます。対応する作業重みは

$$
w_i = \left[ V(\hat{\mu}_i) \cdot \left(\frac{d\hat{\eta}_i}{d\hat{\mu}_i}\right)^2 \right]^{-1}
$$

と定義されます。対角重み行列 $W_t = \text{diag}(w_1, \ldots, w_n)$ を用いると、IRLS更新式は次のWLS標準形に帰着します。

$$
\beta_{t+1} = \left(\mathbf{X}^T W_t \mathbf{X}\right)^{-1} \mathbf{X}^T W_t \mathbf{z}_t
$$

この式がFisher Scoringの更新式と代数的に等価であることは、$\mathcal{I}(\beta_t) = \phi^{-1} \mathbf{X}^T W_t \mathbf{X}$ の関係から確認できます。正準リンクを使用する場合は $d\eta/d\mu$ の形が簡略化され、作業重みの計算式が $w_i = \phi^{-1} \cdot [b”(\theta_i)]^{-1}$ の形になります。作業応答 $z_i$ は各反復における「近似的な目標値」、作業重み $w_i$ は「各観測の局所的な信頼度」として機能しており、これらが対数尤度曲面の局所線形化近似における勾配と曲率情報をそれぞれ反映しています。

（Fig1. IRLSアルゴリズムの収束過程：反復回数と対数尤度の単調増加（ロジスティック回帰））

（Fig2. ロジスティック回帰における作業重みの構造：推定確率 μ と重み w = μ(1−μ) の関係）

前提条件として、分散関数 $V(\mu)$ が正値であること（指数型分布族の性質から保証される）、およびリンク関数が微分可能であることが必要です。非正準リンクを使用する場合は作業重みの計算式が複雑化し、デバッグが困難になります。重みの値域が極端に広い場合は $W_t$ の数値条件が悪化し、計算が不安定になる可能性があります。

収束診断と数値的安定性

IRLSの収束判定には、逸脱度の変化量を基準とする方法が広く用いられます。逸脱度 $D$ は飽和モデルの対数尤度 $\ell_s$ との差分として $D = -2(\ell(\hat{\beta}) – \ell_s)$ で定義され、収束基準は

$$
|D_{t+1} – D_t| < \varepsilon $$

と設定されます。$\varepsilon$ の選択は問題の精度要件に依存し、固定の閾値を全問題に一律適用することは適切ではありません。最大反復回数とトレランス $\varepsilon$ はトレードオフの関係にあり、トレランスを厳しく設定するほど収束に要する反復数が増加します。実際には最大反復回数に達した場合でも収束が成立していない可能性があり、逸脱度の推移を直接確認することが重要です。

$\mathbf{X}^T W \mathbf{X}$ の数値的安定性は条件数によって評価されます。

$$
\kappa(\mathbf{X}^T W \mathbf{X}) = \frac{\lambda_{\max}(\mathbf{X}^T W \mathbf{X})}{\lambda_{\min}(\mathbf{X}^T W \mathbf{X})}
$$

条件数が大きいほど浮動小数点演算における誤差が増幅され、係数推定の信頼性が低下します。説明変数間の多重共線性が強い場合や、作業重みの値域が極端に広い場合に条件数は大きくなります。

完全分離とは、あるベクトル $\mathbf{c}$ が存在して $\mathbf{c}^T \mathbf{x}_i > 0 \Leftrightarrow y_i = 1$（二値応答の場合）が成立する状況であり、この条件下ではロジスティック回帰の最尤推定量が存在しません。反復のたびに係数の絶対値が拡大し、IRLSが収束しない形として現れます。準完全分離はその緩和版であり、一部の観測のみで線形分離が成立する状況です。どちらの場合も係数の急激な拡大として診断できます。

（Fig3. 完全分離問題の図示：係数の反復発散と通常収束との対比）

凸問題（指数型分布族と正準リンクの組み合わせ）では初期値によらず大域最適解への収束が保証されます。収束判定閾値は問題の精度要件に応じて設定することが前提です。完全分離の場合は最尤推定量が存在せず、Firth法などの正則化推定への移行が必要です。多重共線性が強い場合は $\mathbf{X}^T W \mathbf{X}$ が悪条件化し、係数標準誤差が過大になります。

推定アルゴリズムの比較と実用上の指針

IRLSはFisher Scoringと代数的に等価であり、二次収束特性を持ちます。二次収束では各反復ごとに誤差がおおよそ二乗のオーダーで減少するため、通常5〜10反復程度で実用的な精度に収束します。これに対して一次法（勾配降下法）は線形収束であり、同等の精度を得るまでの反復数は大幅に多くなります。勾配降下法の更新式は次の形をとります。

$$
\beta_{t+1} = \beta_t + \alpha \, U(\beta_t)
$$

$\alpha$ はステップサイズです。二次収束は誤差のオーダーが $O(q^{2^t})$（$0 < q < 1$）であるのに対し、線形収束では $O(q^t)$ となります。精度の高い解を得るまでの反復数において、両者には指数的な差が生じます。

確率的勾配降下法（SGD）は各反復でミニバッチの勾配のみを使用するため大規模データに適用できますが、確率的収束のため最終解に基づく標準誤差や検定統計量の直接的な算出には追加の理論的枠組みが必要です。欠損データや潜在変数を含む一般化線形モデルでは、EMアルゴリズムが対数尤度の単調増加を保証しながら最尤推定量を求める枠組みとして機能します。正則化一般化線形モデル（$\ell_1$ 正則化など）に対しては、座標降下法が各パラメータを順次最適化する方式として有効であり、非微分点を含む目的関数にも適用できます。IRLSを $\ell_1$ 正則化に直接拡張しようとすると、WLS内部でのサブ問題が非平滑になるため実装が複雑になります。

IRLSは毎反復 $p \times p$ 行列演算を要するため $n \gg p$ の標準的な設定では効率的ですが、$p$ が大きい高次元データでは計算コストが制約となります。SGDの適用時には、最終解の統計的解釈（信頼区間・仮説検定）が通常の漸近理論と直接対応しない点に留意が必要です。

生物統計学への応用

臨床試験における二値転帰（治療成功・失敗）の治療効果推定では、ロジスティック回帰をIRLSで推定し、Fisher情報行列の逆行列から得られる標準誤差に基づいてWald信頼区間を構築した上で、治療群と対照群のオッズ比を報告します。Wald統計量は推定値 $\hat{\beta}_j$ をその推定標準誤差 $\sqrt{[\mathcal{I}(\hat{\beta})^{-1}]_{jj}}$ で除した量であり、漸近正規分布のもとで検定と信頼区間の構築に用いられます。

小標本（$n < 50$ 程度）ではFisher情報行列に基づく漸近正規近似の精度が低下し、Wald信頼区間が実際の不確実性に比べて過度に狭くなる場合があります。希少転帰や強い予測力を持つ共変量が存在する場合には完全分離が発生しやすく、通常の最尤推定では係数が発散するため、正則化推定（Firth法）への移行が必要になることがあります。また、推定されたオッズ比は観察データに基づく関連性の指標であり、治療と転帰の間の因果関係を直接示すものではありません。無作為化試験のデザインと交絡調整の評価を併せて行うことが、推定結果の解釈の信頼性を高める上で必要です。