生存分析の基礎

生存分析タブで使われている統計理論の背景です。操作方法は生存分析のページを参照してください。

生存時間データと打ち切り

生存分析は「イベントが発生するまでの時間」を分析する手法です。「生存」という名前は医学に由来しますが、対象は死亡に限りません。機械の故障までの時間、顧客の解約までの期間、再犯までの日数など、何らかのイベントが起きるまでの時間を扱う問題全般に適用できます。

生存時間データの特徴は 打ち切り（censoring） の存在です。観察期間内にイベントが発生しなかった対象（たとえば臨床試験の終了時にまだ生存していた患者、追跡から脱落した患者）は、「少なくともここまではイベントが起きなかった」という不完全な情報しか持ちません。

打ち切りを単純に除外すると、イベントが起きやすい対象だけが残り、生存時間を過小推定します。打ち切りのある対象を「イベントなし」として扱うと、真の生存時間が分からないまま過大推定します。生存分析の手法はこの打ち切りを適切に扱うために設計されています。ただし、打ち切りがイベント発生と独立であること（独立打ち切り）が前提です。独立打ち切りのもとでは打ち切りが推定に追加の情報をもたらさず（非情報的）、KM 推定量や Cox モデルの推定が偏りません。副作用の悪化で追跡から脱落するなど、打ち切りの発生がイベントの起きやすさと関連している場合、この前提が崩れて推定にバイアスが生じます。MIDAS が扱う打ち切りは右打ち切り（観察期間終了時やフォローアップ喪失による打ち切り）のみです。左打ち切り（観察開始時点で既にイベントが起きていたが、正確な発生時刻が分からない場合）や区間打ち切り（イベントの発生時刻が2つの観察時点の間にしか特定できない場合）には対応していません。

通常の回帰では扱えない理由

打ち切りがなければ、生存時間を応答変数とする通常の回帰分析が使えます。しかし打ち切りデータは「真の値は観測値以上である」という不等式情報であり、通常の残差（ $y_i - \hat{y}_i$ ）を定義できません。生存分析は尤度関数の構成にこの不等式情報を組み込むことで、打ち切りを正しく扱います。

生存関数とハザード関数

生存時間 $T$ の分布は2つの関数で特徴づけられます。

生存関数 $S(t) = P(T > t)$ は、時点 $t$ までイベントが発生しない確率です。 $S(0) = 1$ から始まり、時間とともに単調に減少します。

ハザード関数 $h(t)$ は、時点 $t$ まで生存していた条件のもとで、その直後にイベントが発生する瞬間的な強度です:

h(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}

ハザードは確率ではなく単位時間あたりの強度なので、1を超えることがあります。生存関数とハザード関数は $S(t) = \exp\left(-\int_0^t h(u)\,du\right)$ の関係で結ばれており、一方が分かれば他方が決まります。

Kaplan-Meier 推定量

Kaplan-Meier 推定量はノンパラメトリックな生存関数の推定法です。分布の形を仮定せず、観測されたイベント時刻から直接 $S(t)$ を推定します。

イベントが発生した時刻を $t_1 < t_2 < \cdots < t_k$ とし、各時刻 $t_i$ でのリスク集合（まだイベントを経験していない対象の数）を $n_i$ 、イベント発生数を $d_i$ とすると:

\hat{S}(t) = \prod_{t_i \le t} \left(1 - \frac{d_i}{n_i}\right)

各イベント時刻での「生き残り率」を累積的に掛け合わせています。打ち切りはリスク集合の減少を通じて反映されます。打ち切りの発生した時点でその対象はリスク集合から離脱しますが、イベント数 $d_i$ には数えません。

同一時刻に複数のイベントがある場合（同着）、 $d_i$ にその時刻のイベント数を、 $n_i$ にはその時刻の直前のリスク集合サイズを用います。

非情報的打ち切りのもとで、KM 推定量は $S(t)$ の一致推定量です。ただし追跡終了間際では打ち切りによりリスク集合が小さくなるため、推定の分散が大きくなり不安定になります。推定量の分散はデルタ法を適用した Greenwood の公式で求めます:

\widehat{\operatorname{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{t_i \le t} \frac{d_i}{n_i(n_i - d_i)}

この分散から信頼区間を構築します。Greenwood の公式から得た $\hat{S}(t)$ の標準誤差を $\text{SE}$ と書くと、 $\log \hat{S}(t)$ の標準誤差はデルタ法により $\text{SE}/\hat{S}(t)$ となります。MIDAS は log 変換法を使用し、 $\exp\!\bigl(\log \hat{S}(t) \pm z \cdot \text{SE}/\hat{S}(t)\bigr)$ として信頼区間を計算します。log 変換により、信頼区間が $[0, 1]$ の範囲を逸脱しにくくなります。

RMST（制限平均生存時間）

RMST (Restricted Mean Survival Time) は、制約時点 $\tau$ までの生存関数の面積です。

\text{RMST}(\tau) = \int_0^\tau S(t)\,dt

$E[\min(T, \tau)]$ 、すなわち追跡を $\tau$ で打ち切ったときの生存時間の期待値です。ハザード比と異なり、比例ハザード仮定を必要としません。

KM 推定量 $\hat{S}(t)$ は階段関数なので、積分は矩形の和になります。 $\tau$ 以下のイベント時刻を $t_1 < t_2 < \cdots < t_m$ とし、 $t_0 = 0$ , $\hat{S}(t_0) = 1$ , $t_{m+1} = \tau$ とすると:

\widehat{\text{RMST}}(\tau) = \sum_{i=0}^{m} \hat{S}(t_i)\bigl(\min(t_{i+1}, \tau) - t_i\bigr)

分散

KM 推定量の各時点での変動が、その後の面積を通じて RMST に伝播します。 $A(t_i) = \int_{t_i}^{\tau} \hat{S}(t)\,dt$ （時刻 $t_i$ から $\tau$ までの KM 曲線下面積）とすると:

\widehat{\text{Var}}\bigl[\widehat{\text{RMST}}(\tau)\bigr] = \sum_{t_i \le \tau} A(t_i)^2 \frac{d_i}{n_i(n_i - d_i)}

$d_i / \bigl(n_i(n_i - d_i)\bigr)$ は Greenwood の公式と同じ項で、時刻 $t_i$ での KM 推定量の変動に対応します。各時点の変動が面積 $A(t_i)$ の二乗で重み付けされるため、早期の時点ほど RMST の分散への寄与が大きくなります。

RMST は $n \to \infty$ で漸近正規性を持ちます。この性質に基づき、信頼区間は Wald 型で構築します: $\widehat{\text{RMST}} \pm z_{\alpha/2} \cdot \text{SE}$ 。

群間差

2 群の RMST 差 $\Delta = \text{RMST}_1 - \text{RMST}_2$ について、群が独立なので分散は各群の分散の和です。

\text{Var}(\hat\Delta) = \text{Var}(\widehat{\text{RMST}}_1) + \text{Var}(\widehat{\text{RMST}}_2)

信頼区間は同じく Wald 型です。3 群以上の場合はペアごとに差と信頼区間を算出しますが、多重性の調整は行いません。

Cox 比例ハザードモデル

モデルの定式化

Cox (1972) の比例ハザードモデルは、共変量がハザードに与える効果を推定するセミパラメトリックモデルです:

h(t \mid X) = h_0(t) \exp(\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p)

$h_0(t)$ はベースラインハザード（すべての共変量が0のときのハザード）、 $\exp(\beta_j)$ は共変量 $X_j$ が1単位増加したときのハザード比です。

「セミパラメトリック」と呼ばれるのは、 $\beta$ はパラメトリックに推定しますが、 $h_0(t)$ の関数形を指定しないためです。これにより、ベースラインハザードの分布を仮定する必要がなくなります。 $\beta$ の推定後、累積ベースラインハザード $H_0(t)$ をノンパラメトリックに推定し、特定の共変量値に対する生存関数 $S(t|X)$ を求められます。

S(t \mid X) = \exp\!\bigl(-H_0(t) \cdot \exp(\beta' X)\bigr)

ここで $H_0(t) = \int_0^t h_0(u)\,du$ です。 $H_0(t)$ の推定では、同着の処理に係数推定と同じ Efron 法を使います。MIDAS は $H_0(t)$ と、指定した共変量値での $S(t|X)$ を出力します。

このモデルの共変量 $X$ は各対象について観察期間を通じて固定された値です。時間とともに変化する共変量（時間依存共変量）を扱うには拡張が必要であり、MIDAS は現在この拡張に対応していません。

比例ハザード仮定

モデルの核心的な仮定は、共変量の効果が時間によらず一定であることです。つまり2つの対象のハザード比 $h(t \mid X_1) / h(t \mid X_2) = \exp(\beta'(X_1 - X_2))$ は $t$ に依存しません。

この仮定が成り立たない場合（たとえば治療効果が時間とともに薄れる場合）、 $\beta$ の推定値はリスク集合の構成とベースラインハザードに依存する重み付き平均としての意味しか持たず、解釈が困難になります。この平均の重みはベースラインハザードの形状に依存し事前に特定できないため、推定値が何を推定しているかを知ることができません (Struthers & Kalbfleisch, 1986)。

Schoenfeld 残差

比例ハザード仮定の診断に使われる残差です。イベント時点 $t_{(i)}$ ごとに、各共変量 $j$ について定義されます:

r_{ij} = X_{ij} - \bar{X}_j(t_{(i)})

$X_{ij}$ はイベントを経験した対象 $i$ の共変量 $j$ の値です。 $\bar{X}_j(t_{(i)})$ はその時点のリスク集合における共変量 $j$ の加重平均で、次のように定義されます:

\bar{X}_j(t_{(i)}) = \frac{\sum_{k \in \mathcal{R}(t_{(i)})} X_{kj} \exp(X_k'\hat\beta)}{\sum_{k \in \mathcal{R}(t_{(i)})} \exp(X_k'\hat\beta)}

$k$ はリスク集合 $\mathcal{R}(t_{(i)})$ の各対象、 $X_k$ は対象 $k$ の共変量ベクトル、 $X_k'\hat\beta$ はその線形予測子です。 $\exp(X_k'\hat\beta)$ が対象 $k$ の重みとなり、ハザードが高い対象ほど大きく寄与します。この重みは全ての共変量 $j$ に共通です。

同着の処理に Breslow 近似を用いる場合、Schoenfeld 残差の合計はスコア関数（対数部分尤度の $\beta$ に関する勾配）に一致します (Schoenfeld, 1982)。最尤推定値 $\hat\beta$ ではスコア関数はゼロなので、残差の合計もゼロになります。MIDAS は同着の処理に Efron 法を使用しているため厳密には成立しませんが、収束許容誤差の範囲内でゼロに近い値をとります。

スケーリング済み Schoenfeld 残差は、生残差を分散共分散行列でスケーリングし、漸近的に $\beta_j(t)$ の推定値として解釈できるように変換したものです:

r^*_{ij} = d \sum_{k=1}^{p} \hat{V}_{jk} \, r_{ik} + \hat\beta_j

$d$ はイベント総数、 $\hat{V}_{jk}$ は $\hat{V} = \hat{I}(\hat\beta)^{-1}$ （推定された $p \times p$ 分散共分散行列）の $(j,k)$ 成分です。 $d \cdot \hat{V}$ は漸近共分散行列のスケールを調整する役割を果たします。 $r^*_{ij}$ は $\beta_j(t_{(i)})$ の推定値として解釈できます。比例ハザード仮定が成り立つ場合、 $r^*_{ij}$ は時間に対して系統的な傾向を示しません。ただし個々の値は分散が大きいため、時間に対してプロットし LOESS 等で平滑化して傾向を読み取ります (Grambsch & Therneau, 1994)。

MIDAS は以下の診断を表示します（操作方法）:

Proportional Hazards Diagnostics: スケーリング済み Schoenfeld 残差と KM 時間変換の Pearson 相関係数 rho を共変量ごとに表示します
スケーリング済み Schoenfeld 残差プロット: $r^*_{ij}$ を時間に対してプロットし、LOESS 回帰線を重ねます
log(-log(S(t))) プロット: 群別の Kaplan-Meier 推定値を $\log(-\log(S(t)))$ 対 $\log(t)$ にプロットします。比例ハザード仮定の下では曲線は近似的に平行になります

rho の計算では時間軸に KM 時間変換 $g(t)=1-\hat{S}(t^-)$ を使いますが、残差プロットの横軸は生の時間です。 $g(t)$ は時間の単調増加変換なので傾向の向きは一致しますが、軸のスケールは異なります。

部分尤度

Cox モデルのパラメータ推定には部分尤度（partial likelihood）を使います。イベントが時刻 $t_{(i)}$ に発生した対象 $i$ について、その時点のリスク集合 $\mathcal{R}(t_{(i)})$ （まだイベントを経験していない全対象）の中から対象 $i$ がイベントを経験する条件付き確率を考えます:

L(\beta) = \prod_{i:\text{event}} \frac{\exp(X_i'\beta)}{\sum_{j \in \mathcal{R}(t_{(i)})} \exp(X_j'\beta)}

各因子は、時刻 $t_{(i)}$ のリスク集合における条件付き確率 $h(t_{(i)}|X_i) / \sum_j h(t_{(i)}|X_j)$ に対応します。 $h(t|X) = h_0(t)\exp(X'\beta)$ を代入すると $h_0(t_{(i)})$ が分子と分母で約分され、 $\beta$ の推定に $h_0(t)$ を知る必要がなくなります。上の式はイベント時刻が全て異なることを前提としています。同じ時刻に複数のイベントが発生する場合（同着, tied events）、条件付き確率の構成が一意でなくなるため近似が必要です。Breslow 法は同着イベントの各対象に同一のリスク集合を適用する近似、Efron 法はリスク集合を同着イベント間で段階的に減少させるより精度の高い近似です。MIDAS は Efron 法を使用しています。

部分尤度は通常の尤度ではありませんが、漸近的に通常の最尤推定量と同じ性質（一致性、漸近正規性）を持つことが示されています (Cox, 1975; Andersen & Gill, 1982)。

ハザード比の解釈

他の共変量を一定に保ったとき、 $\exp(\beta_j)$ はハザード比（HR）として解釈されます:

HR > 1: $X_j$ が1単位増えるとハザードが $(\text{HR} - 1) \times 100\%$ 増加
HR < 1: ハザードが $(1 - \text{HR}) \times 100\%$ 減少
HR = 1: $X_j$ はハザードに影響しない

信頼区間の幅は推定の精度を反映します。狭い区間はより精密な推定を、広い区間はデータから得られる情報が限られていることを示します。ハザード比は効果の方向と大きさを信頼区間と合わせて読み取れるため、p 値単体よりも情報量が多い指標です。

モデル適合度指標

Cox モデルでは、モデルの判別能力とモデル比較のための指標を報告します。

Concordance Index

Concordance index（Harrell's C）は、モデルが生存時間の順序をどの程度正しく予測できるかを測る指標です。

比較可能なペア（一方がイベントを経験し、他方がその時点でまだリスク集合にいる組み合わせ）について、リスクスコア $\hat\eta_i = X_i\hat\beta$ の大小関係がイベント発生順序と一致する割合として定義されます。0.5 は無情報（ランダムな予測と同等）、1.0 は完全な判別を意味します。

標準誤差は影響関数（infinitesimal jackknife）で推定します。各観測 $i$ に対して $\delta_i = (c_i - C \cdot n_i) / N$ を計算します。 $c_i$ は観測 $i$ を含む concordant ペア数、 $n_i$ は比較可能ペア数、 $C$ は concordance、 $N$ は比較可能ペアの総数です。 $SE = \sqrt{\sum \delta_i^2}$ です。

AIC

\text{AIC} = -2\ell(\hat\beta) + 2p

$\ell(\hat\beta)$ は部分対数尤度、 $p$ は共変量の数です。値が小さいほど予測と節約のバランスが良いモデルであることを示します。異なる共変量の組み合わせを比較する場合に使います。

参考文献

Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), 187-220. https://www.jstor.org/stable/2985181
Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481. https://www.jstor.org/stable/2281868
Cox, D. R. (1975). Partial likelihood. Biometrika, 62(2), 269-276. https://www.jstor.org/stable/2335362
Andersen, P. K., & Gill, R. D. (1982). Cox's regression model for counting processes: A large sample study. The Annals of Statistics, 10(4), 1100-1120. https://www.jstor.org/stable/2240714
Struthers, C. A., & Kalbfleisch, J. D. (1986). Misspecified proportional hazard models. Biometrika, 73(2), 363-369. https://www.jstor.org/stable/2336212
Schoenfeld, D. (1982). Partial residuals for the proportional hazards regression model. Biometrika, 69(1), 239-241. https://www.jstor.org/stable/2335876
Grambsch, P. M., & Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81(3), 515-526. https://www.jstor.org/stable/2337123

このページの Markdown 版もあります。