Coxモデル診断の目的と残差の基本概念
通常の線形回帰では、残差は観測値と予測値の差として直接定義されます。しかし生存時間解析では、打ち切り観測について真のイベント時刻が確認されないため、同様の定義を適用できません。打ち切り観測が含まれる状況で残差をどのように定義するかが、Coxモデル診断における根本的な課題です。
Coxモデルの部分尤度対数 $\ell_p(\beta)$ は次の形をとります。
$$
\ell_p(\beta) = \sum_{i:\,\delta_i=1} \left[ X_i^\top \beta – \log \sum_{j \in R(t_i)} \exp(X_j^\top \beta) \right]
$$
ここで $R(t_i)$ は時点 $t_i$ でのリスク集合、$\delta_i$ はイベント発生指示変数です。スコア関数 $U(\beta) = \partial \ell_p / \partial \beta$ と情報行列 $I(\beta) = -\partial^2 \ell_p / \partial \beta \partial \beta^\top$ は、残差近似と影響点分析の理論的基盤を構成します。
マーチンゲール残差は、各観測のイベント指示変数から推定累積ハザードを差し引いた量として次のように定義されます。
$$
M_i = \delta_i – \hat{H}(t_i \mid X_i)
$$
ここで $\hat{H}(t_i \mid X_i)$ は観測時点 $t_i$ における推定累積ハザードです。マーチンゲール過程の性質から、非情報的打ち切り仮定のもとで $E[M_i] = 0$ が成立します。この零期待値性が残差診断の理論的根拠となります。
残差分析はCoxモデル構築サイクルにおいて3つの補完的な役割を担います。第一に比例ハザード仮定の検証、第二に共変量の関数形確認（非線形効果の検出）、第三に影響点の識別（推定値に過度な影響を与える観測の特定）です。
非情報的打ち切り仮定のもとでのみ残差解釈が有効であり、打ち切り機構がイベントリスクと関連する場合、残差の零期待値性は成立しません。また、マーチンゲール残差は正規分布に従わないため、正規性に基づく診断手法（正規Q-Qプロット、Shapiro-Wilk検定等）を直接適用できません。
マーチンゲール残差と逸脱度残差
マーチンゲール残差の計算では、Breslowの基準累積ハザード推定量 $\hat{H}_0(t)$ を用いて次のように展開されます。
$$
M_i = \delta_i – \hat{H}(t_i \mid X_i) = \delta_i – \hat{H}_0(t_i)\exp(X_i^\top \hat{\beta})
$$
マーチンゲール残差は $(-\infty, 1]$ の範囲をとります。イベントが発生した観測は正値、打ち切り観測は負値となる傾向があります。
逸脱度残差はマーチンゲール残差を対称化した量であり、外れ値や問題観測の一次スクリーニングに用いられます。
$$
D_i = \text{sign}(M_i)\sqrt{-2\left[M_i + \delta_i \log(\delta_i – M_i)\right]}
$$
逸脱度残差は概ね対称な分布をとり、$|D_i|$ が2〜3を超える観測が外れ値の候補として識別されます。なお、$\delta_i = 0$ の場合は $\delta_i \log(\delta_i – M_i) = 0$ として計算されます。
共変量の関数形を確認する際は、マーチンゲール残差と連続共変量をプロットし、LOWESS（局所重み付き散布図平滑化）曲線を重ねます。LOWESS曲線がゼロ水平線から体系的に乖離するパターン（U字型や逆U字型）は、その共変量がlog-ハザード上で非線形な効果を持つことを示唆します。
（Fig1. マーチンゲール残差 vs 連続共変量のプロット（LOWESS平滑化付き）による非線形性の検出）
スコア残差は各観測のスコア方程式への寄与量（ベクトル量）を表し、主に影響点分析の基礎として用いられます。Schoenfeld残差は各イベント時点での観測共変量と期待共変量の差として定義され、時間との相関を利用した比例ハザード仮定の検証に特化しています
推定済みモデルが局所的に正しいことを前提として残差を解釈する必要があります。打ち切り割合が高いデータセットでは、マーチンゲール残差が全体的に負に偏る傾向があり、閾値判断の際に注意が必要です。イベント数が少ない場合、逸脱度残差の分布は不安定になります。
| 残差の種類 | 定義の概要 | 主な用途 | 値の範囲・特性 |
|---|---|---|---|
| マーチンゲール残差 | $M_i = \delta_i – \hat{H}(t_i \mid X_i)$ | 共変量の関数形確認（LOWESSプロット）、外れ値検出の補助 | $(-\infty, 1]$、非対称、非情報的打ち切り仮定下で零期待値 |
| 逸脱度残差 | $D_i = \text{sign}(M_i)\sqrt{-2[M_i + \delta_i\log(\delta_i – M_i)]}$ | 外れ値・問題観測の一次スクリーニング | 概ね対称、$|D_i| > 2$〜3 で外れ値候補 |
| スコア残差 | 各観測のスコア方程式への寄与量（ベクトル量） | 影響点分析（dfbeta計算の基礎） | ベクトル、全観測の和が零ベクトル |
| Schoenfeld残差 | 各イベント時点での観測共変量と期待共変量の差 | 比例ハザード仮定の時間依存性検証 | イベント発生数と同数、時間との相関が零で仮定成立 |
影響点分析：dfbetaとCook距離の拡張
影響点分析では、個々の観測を除去した場合に回帰係数推定値がどの程度変化するかを定量化します。完全な再適合（1観測削除ごとにモデルを再推定）は計算コストが高いため、Coxモデルでは情報行列の逆行列とスコア残差ベクトルの積による近似が用いられます。
$$
\text{dfbeta}_i \approx I(\hat{\beta})^{-1} L_i
$$
ここで $L_i$ は観測 $i$ のスコア残差ベクトル、$I(\hat{\beta})$ は推定情報行列です。dfbeta$_i$ は $\hat{\beta}$ と同じ単位を持ち、観測 $i$ を除外した場合の係数変化の近似量を表します。
dfbetaと同じ単位で共変量間を比較することはできないため、推定標準誤差による標準化を施したdfbetasが用いられます。
$$
\text{dfbetas}_i = \frac{\text{dfbeta}_i}{\widehat{\text{se}}(\hat{\beta})}
$$
dfbetasは無次元量であり、共変量の測定単位が異なる場合でも影響度の横断比較が可能です。一般的な閾値として $|\text{dfbetas}_i| > 2/\sqrt{n}$ が用いられます。観測番号を横軸、dfbetasを縦軸にとった指数プロットにより、閾値を超える観測を視覚的に識別できます。
（Fig2. dfbeta指数プロット：各観測のCox回帰係数推定値への影響度（水平参照線 ±2/√n））
影響点と外れ値は概念的に区別されます。外れ値は応答変数（生存時間）の観点から異常な観測であり、影響点は係数推定値への寄与が大きい観測です。外れ値であっても回帰係数への影響が小さい場合があり、逆に外れ値でなくても共変量の極端な値によって影響点となる場合があります。複数の共変量を持つモデルでは、dfbetaはベクトル量となるため、共変量ごとの指数プロットを個別に作成して解釈します。
複数の影響点が相互に隠蔽し合うマスキング効果が生じる場合、単純な指数プロットでは全影響点を検出できません。また、dfbetaは近似式であるため、小サンプルでは実際の観測削除法による推定値との乖離が生じます。
AIC・BICによる変数選択と逐次選択法
Coxモデルでは、部分尤度に基づくAIC（赤池情報量規準）を用いてモデルを比較します。
$$
\text{AIC} = -2\log\hat{L}_p(\hat{\beta}) + 2k
$$
ここで $k$ はモデルに含まれる共変量の数です。AICは予測精度とモデルの複雑さのトレードオフを捉える指標として機能します。
BIC（ベイズ情報量規準）はイベント発生数 $d$ を用いたペナルティを導入します。全観測数ではなくイベント数を基準とする定式化が、部分尤度に基づくCoxモデルに対して推奨されます。
$$
\text{BIC} = -2\log\hat{L}_p(\hat{\beta}) + k\log d
$$
ペナルティ $k\log d$ はイベント数が多いほど大きくなり、AICの $2k$ より強くモデルの複雑さを抑制する傾向があります。
逐次変数選択法には、前向き選択（変数を一つずつ追加）、後向き除去（変数を一つずつ削除）、ステップワイズ選択（追加と削除を繰り返す）の3種類があります。p値基準による変数選択は多重検定問題を引き起こすため、AIC基準での選択が推奨されます。変数選択の目的が予測モデル構築か交絡制御かによって選択基準が異なり、交絡制御目的では先験的知識に基づく変数の強制投入が必要です。比較するモデル間で欠損値の扱いが異なると観測数が変動するため、必ず同一のデータセット（同一観測数・イベント数）に適合されたモデル間で比較を行う必要があります。
AICおよびBICによる比較は、同一データセットに適合されたモデル間でのみ有効です。逐次選択によって選ばれたモデルの係数推定値には過大推定バイアスが生じ、信頼区間は過小評価される傾向があります。また、共変量間の高い相関（多重共線性）が存在する場合、ステップワイズ選択は反復によって不安定な結果をもたらします。
予測性能の評価：C-indexとBrier score
Coxモデルの識別能を評価する指標として、HarrellのC-indexが広く用いられます。C-indexは全比較可能ペアのうち、予測リスクの高低順がイベント発生の時間順序と一致するペアの割合として定義されます。
$$
C = \frac{\text{一致ペア数}}{\text{比較可能ペア総数}}
$$
ペア $(i, j)$ が比較可能とみなされる条件は、$\delta_i = 1$（観測 $i$ がイベントを経験）かつ $t_i < t_j$、または両観測でイベントが確認され $t_i \neq t_j$ が成立する場合です。一致ペアは、イベントが先に発生した観測の予測ハザードが他方より高い場合を指します。
C-indexは二値結果に対するAUCの生存分析への拡張です。$C = 0.5$ はランダム予測に相当し、$C = 1.0$ は完全識別を意味します。臨床的に有用とみなされる目安は疾患領域によって異なり、固定された閾値の適用には注意が必要です。打ち切り割合の変化によって比較可能ペアの構成が変わり、C-index値が変動する打ち切り依存性の問題があります。この問題に対し、逆確率打ち切り重み付けを用いたUno C統計量は打ち切りの影響を補正した識別能の推定量として提案されています。
統合Brierスコアは、各時点での予測二乗誤差の時間平均として定義され、予測確率の較正精度を評価します。
$$
\text{IBS} = \frac{1}{\tau} \int_0^{\tau} \frac{1}{n}\sum_{i=1}^n \left(\mathbf{1}[T_i > t] – \hat{S}(t \mid X_i)\right)^2 dt
$$
打ち切りデータへの適用では逆確率打ち切り重み付けによる補正が必要です。C-indexは観測間の相対的リスク順序（識別能）を評価しますが、予測確率の絶対的な精度（較正）を評価しません。識別能と較正は独立した性能指標であり、C-indexが高くても較正が不良なモデルは個人レベルの予後予測に誤った確率値を与える可能性があります。
交差検証とブートストラップによる過適合評価
モデルの見かけの予測性能は過適合によって過大評価されます。内部検証はこの楽観主義を定量化し、新規データに対して期待される性能の推定値を提供します。$k$ 分割交差検証では、データを $k$ 個のフォールドに分割し、$k-1$ フォールドでモデルを適合して残り1フォールドで性能を評価するプロセスを繰り返します。生存分析での適用では、打ち切り観測の割り当てが各フォールドで均等になるよう層別化が推奨されます。
ブートストラップ内部検証の手順は次の通りです。まず元データから復元抽出でブートストラップ標本を生成し、そのブートストラップ標本上でモデルを適合して見かけのC-index $C_\text{apparent}^{(b)}$ を計算します。次に、同じモデルを元データに適用してテストC-index $C_\text{test}^{(b)}$ を計算します。この操作を $B$ 回（通常200回以上）繰り返します。
楽観主義の推定量は各反復での差の平均として計算されます。
$$
\widehat{\text{optimism}} = \frac{1}{B}\sum_{b=1}^B \left( C_\text{apparent}^{(b)} – C_\text{test}^{(b)} \right)
$$
楽観主義補正後のC-indexは次のように算出されます。
$$
C_\text{corrected} = C_\text{apparent} – \widehat{\text{optimism}}
$$
（Fig3. ブートストラップ内部検証：見かけのC-indexと楽観主義補正後C-indexの分布比較）
楽観主義の大きさは過適合の程度を反映します。TRIPOD/PROBAST等の予後モデル報告基準では外部検証の実施が推奨されており、内部検証のみでは研究母集団外への予測性能の一般化は保証されません。小サンプルでは $k$ 分割交差検証の推定量の分散が大きくなるため、ブートストラップ法がより安定した推定を与えます。
疫学コホート研究への応用と競合モデルとの比較
前向き疫学コホート研究でCoxモデルによる予後モデルを構築する場合、次の分析フローが標準的です。まず研究目的と先験的知識に基づいて候補共変量を設定し、欠損値の処理方針（完全ケース解析または多重代入法）を決定します。次にAICまたはBIC基準による変数選択とともに、交絡因子については専門的判断に基づく強制投入を行います。その後、マーチンゲール残差プロットによる関数形確認、dfbeta解析による影響点検出、比例ハザード仮定の検証を実施します。最後にブートストラップ内部検証で楽観主義を定量化し、補正C-indexを報告します。
交絡因子と効果修飾因子の扱いは分析目的によって異なります。予測目的のモデル構築ではC-indexへの寄与を基準とした変数選択が中心となりますが、ハザード比の解釈を目的とする場合は未計測交絡の影響を明示する必要があります。観察研究では因果関係と相関関係の区別が不可欠であり、観察されたハザード比の推定値は未計測交絡因子の存在によって真の効果量と乖離している可能性があります。
競合リスクが存在するコホート研究では、通常のCox解析が競合イベントを単純打ち切りとして扱うため、特定の原因によるイベントの累積発症率が過大推定される問題が生じます。競合リスクが想定される場合、Fine-Grayサブ分布ハザードモデルへの移行が検討されます。Fine-Grayモデルは累積発症率関数を直接モデル化することで、競合リスクを考慮した予後予測を提供します。
共変量の数がイベント数に対して多い場合（イベント数/共変量数が10未満を目安として）、ペナルティ付きCoxモデルであるLasso-Coxへの移行が考慮されます。Lasso-Coxは部分尤度に $L_1$ ペナルティを付加することで自動的に変数選択とパラメータ縮小を実施しますが、選択された共変量のハザード比には縮小バイアスが残る点と、ペナルティパラメータの選択によって結果が変化する点に留意が必要です。Schoenfeld残差診断は比例ハザード仮定の検証に特化した役割を担い、本記事で扱う残差診断手順と相互補完的に機能します。
コホート研究での非情報的打ち切り仮定の現実妥当性を検討することが重要です。受診中断や転居による追跡脱落がイベントリスクと関連する場合、この仮定は成立せず、残差診断結果の解釈に制約が生じます。観察研究での未計測交絡因子の存在はC-indexを過大評価させ、内部検証で得られた予測性能は他の研究集団への外挿が制限されます。
