9. 部分尤度：Cox モデルの推定理論

部分尤度の動機と統計的位置づけ

Cox比例ハザードモデルは、ハザード関数を乗法構造

$$
h(t|\mathbf{x}) = h_0(t)\exp(\boldsymbol{\beta}^\top\mathbf{x})
$$

として定式化します。ここで$h_0(t)$は基準ハザード関数（ベースラインハザード）、$\boldsymbol{\beta}$は回帰係数ベクトル、$\mathbf{x}$は共変量ベクトルです。このモデルはセミパラメトリック構造を持ち、関心パラメータ$\boldsymbol{\beta}$のほかに、無限次元のニュアンスパラメータである$h_0(t)$を含みます。

完全尤度を構成すれば$\boldsymbol{\beta}$と$h_0(t)$を同時に推定することも原理的には可能です。しかし$h_0(t)$は非常に自由度の高い関数であり、これを最大化すると各観測時刻で無限大のスパイクが生じて解が崩壊します。また、$h_0(t)$を有限次元のパラメトリック族に制限することは、そのモデル設定の誤りに対してロバストでなくなるという代償を伴います。完全尤度最尤推定では、$h_0(t)$の推定に要するコストが$\boldsymbol{\beta}$の推定精度を損なう可能性があります。

この問題に対し、完全尤度を概念的に

$$
L = L_{\text{partial}}(\boldsymbol{\beta}) \times L_{\text{baseline}}(h_0, \boldsymbol{\beta})
$$

と分解する考え方が有効です。$L_{\text{partial}}$は$h_0(t)$を含まない$\boldsymbol{\beta}$のみの関数であり、$L_{\text{baseline}}$はベースライン成分を担います。この分解により、$h_0(t)$を推定せずに$\boldsymbol{\beta}$の推定が可能になります。これが部分尤度の核心的発想です。

部分尤度は通常の周辺尤度とも条件付き尤度とも異なります。周辺尤度はニュアンスパラメータを積分消去しますが、部分尤度は各イベント時刻での条件付き確率の積によって$h_0(t)$を代数的に消去します。この構成は、イベントの発生順序に着目した条件付けに基づいており、セミパラメトリック推定の中核をなす方法論です。プロファイル尤度は$h_0(t)$を$\boldsymbol{\beta}$の関数として最大化消去する方法であり、部分尤度とは消去の手続きが異なりますが、実用上は近似的に同等な推定結果を与えます。

部分尤度関数の定義と構造

$n$個の個体について生存時間データを観測し、$D$件のイベント発生があったとします。イベント発生時刻の順序統計量を

$$
t_{(1)} \le t_{(2)} \le \cdots \le t_{(D)}
$$

と表記します。時刻$t_{(j)}$におけるリスク集合$R(t_{(j)})$は、その時刻直前まで生存（観察打ち切りを含む）している個体の集合であり、

$$
R(t_{(j)}) = \{i : t_i \ge t_{(j)}\}
$$

と定義されます。リスク集合の構成において、観察打ち切り個体は打ち切り時刻まではリスク集合に含まれ、その後除外されます。この取り扱いは無情報打ち切り仮定のもとで統計的に正当化されます。

時刻$t_{(j)}$でイベントが発生した個体を$j$番目とすると、リスク集合内のいずれかの個体がこの時刻にイベントを経験するという条件のもとで、実際にそれが個体$j$である条件付き確率は

$$
\frac{\exp(\boldsymbol{\beta}^\top\mathbf{x}_j)}{\displaystyle\sum_{i \in R(t_{(j)})}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)}
$$

です。分子・分母に$h_0(t_{(j)})$が共通して現れるため代数的に相殺され、この確率は$h_0(t)$に依存しません。部分尤度は、この条件付き確率の$D$個の積として定義されます。

$$
L_P(\boldsymbol{\beta}) = \prod_{j=1}^{D}\frac{\exp(\boldsymbol{\beta}^\top\mathbf{x}_j)}{\displaystyle\sum_{i \in R(t_{(j)})}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)}
$$

対数をとると対数部分尤度は

$$
\ell_P(\boldsymbol{\beta}) = \sum_{j=1}^{D}\left[\boldsymbol{\beta}^\top\mathbf{x}_j – \log\sum_{i \in R(t_{(j)})}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)\right]
$$

と表現されます。各因子を互いに独立と見なすことは厳密には近似ですが、イベント時刻での条件付けにより$h_0(t)$の情報が十分に消去されており、残余の依存性は実用上無視できることが理論的に正当化されます。計数過程の記法では、この構造はマーチンゲール理論の枠組みで厳密に扱われます。

仮定として、イベント時刻は理論上連続な確率変数であり同着（ties）は生じないことを前提とします。また、打ち切りの発生機構はその後の生存時間の分布に依存しないという無情報打ち切り仮定が必要です。

（Fig1. イベント時刻ごとのリスク集合の縮小と部分尤度の積構造（模擬臨床試験データ））

同着処理：Breslow 近似と Efron 近似

実データでは、観測時刻の離散的な測定精度（日単位・月単位など）や丸め処理により、複数の個体が同一時刻にイベントを経験する同着が頻発します。連続時間の仮定に基づく基本的な部分尤度はそのままでは適用できないため、同着に対処する近似が必要となります。

時刻$t_{(j)}$で$d_j$件の同着イベントが発生したとし、その個体の集合を$D_j$、共変量の和を$\mathbf{s}_j = \sum_{k \in D_j}\mathbf{x}_k$とします。Breslow 近似の対数部分尤度は

$$
\ell_P^{\text{Breslow}}(\boldsymbol{\beta}) = \sum_{j=1}^{D^*}\left[\boldsymbol{\beta}^\top\mathbf{s}_j – d_j\log\sum_{i \in R(t_{(j)})}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)\right]
$$

です。ここで$D^*$は同着をグループ化した後のイベント時刻数です。この式では分母のリスク集合の和を$d_j$回同一の値で使用しており、同着内での順序を無視した粗い近似となっています。

Efron 近似はより精緻で、分母を重み付き調整します。

$$
\ell_P^{\text{Efron}}(\boldsymbol{\beta}) = \sum_{j=1}^{D^*}\left[\boldsymbol{\beta}^\top\mathbf{s}_j – \sum_{k=1}^{d_j}\log\!\left(\sum_{i \in R(t_{(j)})}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i) – \frac{k-1}{d_j}\sum_{i \in D_j}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)\right)\right]
$$

Efron 近似は、$d_j$件の同着イベントが任意の順序で発生したと仮定し、リスク集合から同着イベント者を逐次的に除去する重み付けを導入します。これにより同着内の順序情報の欠損を部分的に補正できます。Exact（離散）法は同着内の全順列にわたる組み合わせ和を正確に計算するものですが、$d_j$が大きくなると計算量が指数的に増大し、大規模データでは非現実的となります。

ties 割合が高い場合の Breslow 近似では、リスク集合の過大評価により係数$\boldsymbol{\beta}$の推定に系統的バイアスが生じます。Efron 近似はこのバイアスを大幅に低減し、Exact 法との差異も小さいことが確認されているため、標準的な実装で広く推奨されています。Exact 法の計算量爆発のため、同着割合が高いデータ（医療ビッグデータ等）では Efron 近似が実践的な第一選択となります。

（Fig2. 同着割合と Breslow/Efron 近似の係数推定バイアス比較（モンテカルロシミュレーション））

スコア方程式と情報行列の導出

$\boldsymbol{\beta}$の推定は$\ell_P(\boldsymbol{\beta})$を最大化する問題に帰着します。スコア関数$U(\boldsymbol{\beta})$は対数部分尤度を$\boldsymbol{\beta}$で微分して得られます。

$$
U(\boldsymbol{\beta}) = \frac{\partial \ell_P}{\partial \boldsymbol{\beta}} = \sum_{j=1}^{D}\left[\mathbf{x}_j – \bar{\mathbf{x}}(t_{(j)},\boldsymbol{\beta})\right]
$$

ここで$\bar{\mathbf{x}}(t_{(j)},\boldsymbol{\beta})$はリスク集合内の共変量の重み付き期待値

$$
\bar{\mathbf{x}}(t_{(j)},\boldsymbol{\beta}) = \frac{\displaystyle\sum_{i \in R(t_{(j)})}\mathbf{x}_i\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)}{\displaystyle\sum_{i \in R(t_{(j)})}\exp(\boldsymbol{\beta}^\top\mathbf{x}_i)}
$$

です。スコア関数はイベント発生者の実観測共変量$\mathbf{x}_j$と、リスク集合内の重み付き期待共変量$\bar{\mathbf{x}}$の差の総和として解釈できます。イベントが実際に発生した個体の共変量が、リスク集合内の期待値と一致する方向にパラメータが調整されます。スコア方程式$U(\boldsymbol{\beta})=\mathbf{0}$の解が部分尤度推定量$\hat{\boldsymbol{\beta}}$です。

観測情報行列$I(\boldsymbol{\beta})$は対数部分尤度の負の2階微分として定義されます。

$$
I(\boldsymbol{\beta}) = -\frac{\partial^2 \ell_P}{\partial \boldsymbol{\beta}\partial\boldsymbol{\beta}^\top} = \sum_{j=1}^{D}\left[\overline{\mathbf{x}^{\otimes 2}}(t_{(j)},\boldsymbol{\beta}) – \bar{\mathbf{x}}(t_{(j)},\boldsymbol{\beta})\bar{\mathbf{x}}(t_{(j)},\boldsymbol{\beta})^\top\right]
$$

この行列は各イベント時刻での共変量の重み付き分散共分散行列の和として表現され、半正定値となります。対数部分尤度の凹性により$I(\boldsymbol{\beta})$は正定値となり、スコア方程式の解は一意に存在します。Newton-Raphson 反復アルゴリズムによる更新式は

$$
\boldsymbol{\beta}^{\text{new}} = \boldsymbol{\beta}^{\text{old}} + I(\boldsymbol{\beta}^{\text{old}})^{-1}U(\boldsymbol{\beta}^{\text{old}})
$$

です。実用上、観測情報行列と期待情報行列（Fisher情報行列）はほぼ同等な性能を示しますが、観測情報行列を用いたほうが数値的安定性に優れることが多く、多くの実装で採用されています。

推定量の漸近理論

部分尤度推定量$\hat{\boldsymbol{\beta}}$の漸近的性質は、通常の最尤推定理論と類比した形で確立されています。標本サイズ$n$の増加に伴いイベント数が正比例的に増加する規則的な増大条件のもとで、

$$
\sqrt{n}(\hat{\boldsymbol{\beta}} – \boldsymbol{\beta}_0) \xrightarrow{d} \mathcal{N}\!\left(\mathbf{0},\,\Sigma\right)
$$

が成立します。漸近分散$\Sigma$の推定量としては観測情報行列の逆行列$I(\hat{\boldsymbol{\beta}})^{-1}$が用いられます。

この一致性・漸近正規性の証明には、マーチンゲール中心極限定理が根拠として機能します。スコア関数$U(\boldsymbol{\beta}_0)$を計数過程に関するマーチンゲールの増分として表現し、中心極限定理を適用することで漸近正規性が導かれます。これは通常の独立同分布の仮定によらない証明であり、打ち切りや時変共変量を含む一般的な設定に適用可能です。

セミパラメトリック効率下界との関連について、Andersen-Gill の計数過程理論は、比例ハザード仮定が成立するもとで部分尤度推定量が$\boldsymbol{\beta}$に関する情報をほぼ完全に利用することを示します。すなわち、$h_0(t)$を特定のパラメトリック族に制限した完全尤度最尤推定量と比較したとき、部分尤度推定量の漸近的情報損失は理論上限定的です。

漸近理論が適用されるための条件として、比例ハザード仮定の成立と、イベント数の規則的な増大が必要です。有限標本、特に共変量1つあたりのイベント数（Events Per Variable, EPV）が5〜10未満の場合には、係数推定が過学習や数値不安定性を示す可能性があります。また、共変量次元が標本サイズを超える高次元設定（$p \gg n$）では通常の部分尤度推定量は発散し、罰則付き部分尤度（Lasso-Cox や Ridge-Cox）が必要となります。これらは対数部分尤度に正則化項を付加したものであり、$p \gg n$の設定でも安定した推定が可能です。

（Fig3. 部分尤度推定量の漸近正規性：標本サイズによる分布収束）

Wald 検定・尤度比検定・スコア検定

部分尤度に基づく仮説検定は主に3種類の統計量を通じて実施されます。帰無仮説$H_0: \boldsymbol{\beta} = \mathbf{0}$に対して、各統計量は漸近的に自由度$p$のカイ二乗分布に従います。

Wald 統計量は推定量$\hat{\boldsymbol{\beta}}$と推定情報行列$I(\hat{\boldsymbol{\beta}})$を用いて

$$
W = \hat{\boldsymbol{\beta}}^\top I(\hat{\boldsymbol{\beta}})\hat{\boldsymbol{\beta}} \xrightarrow{d} \chi^2_p
$$

と定義されます。対数部分尤度比統計量は帰無仮説下と対立仮説下の対数部分尤度の差の2倍として

$$
LR = 2\!\left[\ell_P(\hat{\boldsymbol{\beta}}) – \ell_P(\mathbf{0})\right] \xrightarrow{d} \chi^2_p
$$

と定義されます。スコア統計量は帰無仮説下でのスコア関数の値と情報行列を用いて

$$
S = U(\mathbf{0})^\top I(\mathbf{0})^{-1}U(\mathbf{0}) \xrightarrow{d} \chi^2_p
$$

と定義されます。3統計量は$n \to \infty$の極限で漸近的に等価であり、同一の漸近的検出力を持ちます。

単一係数$\hat{\beta}_k$に対する Wald 信頼区間は

$$
\hat{\beta}_k \pm z_{\alpha/2} \times \widehat{SE}(\hat{\beta}_k)
$$

として構成されます。ただし Wald 検定はパラメータ変換後の不変性を持たないという特性があります。例えば$\beta$ではなく$\exp(\beta)$（ハザード比）に対して Wald 統計量を適用すると、変換前と異なる結論を与えることがあります。

対数部分尤度比検定はプロファイル尤度信頼区間と対応関係にあります。プロファイル信頼区間は関心パラメータを固定して残りのパラメータを最大化したプロファイル対数尤度に基づいて構成され、Wald 信頼区間と異なりパラメータ変換後の不変性を持ちます。小標本では Wald 検定の被覆確率が名目水準を下回ることがあり、プロファイル信頼区間と Wald 信頼区間の乖離が問題となる場合があります。スコア検定は帰無仮説下での推定のみを必要とするため、制約付きモデルの検定において計算上の優位性を持ちます。

仮定・限界と生物統計学への応用

部分尤度推定が統計的に有効であるためには、いくつかの前提条件の成立が必要です。

最も基礎的な仮定は無情報打ち切り仮定です。打ち切りが発生する機構がその後の生存時間の分布に依存しない場合に無情報打ち切りが成立します。臨床試験では、試験終了による行政的打ち切りや追跡不能は通常無情報打ち切りとして扱えます。しかし毒性副作用による試験中断・脱落では、患者の健康状態（すなわち生存時間の予後因子）と打ち切り機構が連動するため、情報的打ち切りとなります。情報的打ち切りが生じると、部分尤度推定量に系統的バイアスが生じます。この場合、打ち切り機構を明示的にモデル化した競合リスクモデルや感度解析が必要となります。

比例ハザード仮定は、2群間のハザード比が観察期間を通じて一定であることを要求します。この仮定が成立しない場合、推定された$\hat{\boldsymbol{\beta}}$は時間平均的なハザード比の何らかの要約値として解釈されることになりますが、その解釈は明確ではなく、臨床的判断を誤らせる可能性があります。比例ハザード仮定の診断には Schoenfeld 残差のプロットや時間交互作用項の検定が用いられます。

EPV の観点では、共変量1つあたりのイベント数が10未満の場合に係数推定の不安定性が知られており、5未満では過学習のリスクが顕著となります。標本サイズ計画の段階で、所定の検出力と精度基準を満たすイベント数を確保することが求められます。

部分尤度推定では基準ハザード関数$h_0(t)$は直接推定されません。生存曲線や累積ハザードの推定には、別途 Breslow 推定量

$$
\hat{H}_0(t) = \sum_{j:\, t_{(j)} \le t}\frac{d_j}{\displaystyle\sum_{i \in R(t_{(j)})}\exp(\hat{\boldsymbol{\beta}}^\top\mathbf{x}_i)}
$$

を用いて累積ハザードを推定する必要があります。

生物統計学における実践的応用として、臨床試験における全生存時間（OS）および無増悪生存期間（PFS）エンドポイントの解析があります。化学療法群と標準療法群を比較する場合、部分尤度によりハザード比の点推定と95%信頼区間が構成され、対数部分尤度比検定により統計的検定が実施されます。ハザード比は時間依存性を持たない要約統計量として解釈されますが、これは比例ハザード仮定が成立する場合にのみ有効な解釈です。

完全尤度最尤推定との比較では、部分尤度は$h_0(t)$を推定しない分だけ情報を損失しますが、$h_0(t)$のノンパラメトリック推定に伴う推定誤差を回避できます。比例ハザード仮定が正しい場合、この情報損失は漸近的に無視できる程度であることが理論的に示されており、部分尤度はロバストな代替手段として有効です。高次元共変量設定における通常部分尤度の限界を克服するには、罰則付き部分尤度（Lasso-Cox や Ridge-Cox）が用いられます。