Language

生存分析の基礎

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

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

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

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

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

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

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

生存関数とハザード関数

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

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

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

h(t)=lim⁡Δt→0P(t≤T<t+Δt∣T≥t)Δth(t) = \lim_{\Delta t \to 0} \frac{P(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}

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

Kaplan-Meier 推定量

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

イベントが発生した時刻を t1<t2<⋯<tkt_1 < t_2 < \cdots < t_k とし、各時刻 tit_i でのリスク集合（まだイベントを経験していない対象の数）を nin_i、イベント発生数を did_i とすると:

S^(t)=∏ti≤t(1−dini)\hat{S}(t) = \prod_{t_i \le t} \left(1 - \frac{d_i}{n_i}\right)

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

非情報的打ち切りのもとで KM 推定量は S(t)S(t) の一致推定量です。推定量の分散はデルタ法を適用した Greenwood の公式で求めます:

Var⁡^[S^(t)]=S^(t)2∑ti≤tdini(ni−di)\widehat{\operatorname{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{t_i \le t} \frac{d_i}{n_i(n_i - d_i)}

この分散から信頼区間を構築します。MIDAS はデフォルトで log 変換法を使用し、exp⁡(log⁡S^(t)±z⋅SE/S^(t))\exp(\log \hat{S}(t) \pm z \cdot \text{SE} / \hat{S}(t)) として計算します。log 変換により、信頼区間が [0,1][0, 1] の範囲を逸脱しにくくなります。

Log-rank 検定

Log-rank 検定は、2つ以上の群のハザード関数が等しいかどうかを検定するノンパラメトリック検定です。

各イベント時刻 tit_i で、群 jj の観測イベント数 dijd_{ij} と期待イベント数 eij=nij⋅di/nie_{ij} = n_{ij} \cdot d_i / n_i を計算します。nijn_{ij} は時刻 tit_i における群 jj のリスク集合の大きさです。帰無仮説のもとでは、各時点のイベントがリスク集合の群構成比に応じて配分されると期待されます。この配分は超幾何分布に従い、分散もそこから導かれます。2群の場合、O1+O2=∑diO_1 + O_2 = \sum d_i かつ E1+E2=∑diE_1 + E_2 = \sum d_i であるため O1−E1=−(O2−E2)O_1 - E_1 = -(O_2 - E_2) が成り立ち、一方の群の情報で検定統計量を構成できます:

χ2=(O1−E1)2V1\chi^2 = \frac{(O_1 - E_1)^2}{V_1}

O1=∑id1iO_1 = \sum_i d_{1i}（群1の総観測イベント数）、E1=∑ie1iE_1 = \sum_i e_{1i}（総期待イベント数）、V1=∑in1i⋅n2i⋅di(ni−di)ni2(ni−1)V_1 = \sum_i \frac{n_{1i} \cdot n_{2i} \cdot d_i (n_i - d_i)}{n_i^2(n_i - 1)}（分散）です。分母は期待値 E1E_1 ではなく、超幾何分布から導かれる分散 V1V_1 です。3群以上の場合は分散共分散行列を用いた二次形式 χ2=U′V−1U\chi^2 = \mathbf{U}' \mathbf{V}^{-1} \mathbf{U} に拡張します。帰無仮説のもとでこの統計量は自由度 g−1g - 1（gg は群数）のカイ二乗分布に近似的に従います。

Log-rank 検定は重み付き検定の一種であり、ハザード比が時間を通じて一定であるとき（比例ハザード仮定のもとで）、Wilcoxon 型などの他の重み付き検定と比較して最も検出力が高くなります。生存曲線が交差するような状況では検出力が下がります。

Cox 比例ハザードモデル

モデルの定式化

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

h(t∣X)=h0(t)exp⁡(β1X1+β2X2+⋯+βpXp)h(t \mid X) = h_0(t) \exp(\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p)

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

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

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

比例ハザード仮定

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

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

Schoenfeld 残差

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

rij=Xij−Xˉj(t(i))r_{ij} = X_{ij} - \bar{X}_j(t_{(i)})

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

Xˉj(t(i))=∑k∈R(t(i))Xkjexp⁡(Xk′β^)∑k∈R(t(i))exp⁡(Xk′β^)\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)}

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

Schoenfeld 残差の合計はスコア関数（対数部分尤度の β\beta に関する勾配）に一致します (Schoenfeld, 1982)。最尤推定値 β^\hat\beta ではスコア関数は理論的にゼロなので、残差の合計もゼロになります。実際には Newton-Raphson が有限回で打ち切られるため、収束許容誤差の範囲内でゼロに近い値をとります。

スケーリング済み Schoenfeld 残差は、生の残差を分散共分散行列でスケーリングしたものです:

ri∗=d⋅V^ ri+β^r^*_i = d \cdot \hat{V} \, r_i + \hat\beta

dd はイベント総数、V^=I^(β^)−1\hat{V} = \hat{I}(\hat\beta)^{-1} は推定された分散共分散行列、rir_i と ri∗r^*_i はそれぞれイベント時点 ii での生残差ベクトルとスケーリング済み残差ベクトルです。ri∗r^*_i の第 jj 成分は βj(t(i))\beta_j(t_{(i)}) の推定値として解釈できます。比例ハザード仮定が成り立つ場合、ri∗r^*_i は時間に対して系統的な傾向を示しません (Grambsch & Therneau, 1994)。

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

  • Grambsch-Therneau 検定: スケーリング済み Schoenfeld 残差と時間の関連を検定します。共変量ごとの検定と全体の検定を報告します
  • スケーリング済み Schoenfeld 残差プロット: rij∗r^*_{ij} を時間に対してプロットし、LOESS 回帰線を重ねます
  • log(-log(S(t))) プロット: 群別の Kaplan-Meier 推定値を log⁡(−log⁡(S(t)))\log(-\log(S(t))) 対 log⁡(t)\log(t) にプロットします。比例ハザード仮定の下では曲線は近似的に平行になります

部分尤度

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

L(β)=∏i:eventexp⁡(Xi′β)∑j∈R(t(i))exp⁡(Xj′β)L(\beta) = \prod_{i:\text{event}} \frac{\exp(X_i'\beta)}{\sum_{j \in \mathcal{R}(t_{(i)})} \exp(X_j'\beta)}

各因子は、時刻 t(i)t_{(i)} のリスク集合における条件付き確率 h(t(i)∣Xi)/∑jh(t(i)∣Xj)h(t_{(i)}|X_i) / \sum_j h(t_{(i)}|X_j) に対応します。h(t∣X)=h0(t)exp⁡(X′β)h(t|X) = h_0(t)\exp(X'\beta) を代入すると h0(t(i))h_0(t_{(i)}) が分子と分母で約分され、β\beta の推定に h0(t)h_0(t) を知る必要がなくなります。部分尤度は通常の尤度ではありませんが、漸近的に通常の最尤推定量と同じ性質（一致性、漸近正規性）を持つことが示されています (Cox, 1975)。

ハザード比の解釈

exp⁡(βj)\exp(\beta_j) はハザード比（HR）として解釈されます:

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

ハザード比の信頼区間が1を含まない場合、データは「ハザードに差がない」という帰無仮説と整合しません。信頼区間の幅は推定の精度を反映し、狭い区間はより精密な推定を、広い区間はデータから得られる情報が限られていることを示します。ハザード比は効果の方向と大きさを直接読み取れるため、p 値単体よりも解釈に有用です。

検定統計量

Cox モデルでは「全ての β=0\beta = 0」という帰無仮説に対する3種類の検定統計量が報告されます。いずれも帰無仮説の下で自由度 pp（共変量数）の χ2\chi^2 分布に漸近的に従い、大標本では同じ値に収束します。

尤度比検定

Λ=2[ℓ(β^)−ℓ(0)]\Lambda = 2\bigl[\ell(\hat\beta) - \ell(0)\bigr]

帰無モデル（β=0\beta = 0）と適合モデルの対数部分尤度の差に基づきます。尤度面の全体的な形状を反映するため、Wald 検定や Score 検定のような局所近似に依存しません。

Wald 検定

W=β^′ I^(β^) β^W = \hat\beta' \, \hat{I}(\hat\beta) \, \hat\beta

最尤推定値 β^\hat\beta と推定された情報行列 I^(β^)\hat{I}(\hat\beta) に基づきます。これは対数尤度の MLE での曲率を使った局所近似です。係数テーブルの各共変量の p 値は、この多変量 Wald 検定の1変量版 zj=β^j/SE(β^j)z_j = \hat\beta_j / SE(\hat\beta_j) に基づいています。モデル全体の Wald 統計量と個別の zz 値を見比べることで、「個別には有意だがモデル全体では有意でない」のような矛盾を検出でき、多重共線性や推定の不安定さに気づく手がかりになります。

Wald 検定の弱点は、β^\hat\beta が大きいときに尤度面の非対称性を捉えられないことです。

Score 検定

S=U(0)′ I(0)−1 U(0)S = U(0)' \, I(0)^{-1} \, U(0)

U(0)U(0) は β=0\beta = 0 でのスコア関数（対数部分尤度の勾配）、I(0)I(0) はその点での情報行列です。β^\hat\beta を計算する必要がないため、Newton-Raphson が収束しない場合でも計算できます。

収束しない典型的な状況は、共変量がイベントをほぼ完全に予測する場合です。たとえば小規模な臨床試験で治療群にイベントが集中すると、β^→∞\hat\beta \to \infty に向かって反復が発散します。ハザード比 exp⁡(β)\exp(\beta) をモデル化しているため、効果が極めて大きいと係数が発散するのです。このとき尤度比検定と Wald 検定は計算できませんが、Score 検定は「共変量と生存の間に関連がある」という証拠を示せます。ただし効果の大きさは推定できません。対処法として Firth 補正や正確法がありますが、MIDAS は現在これらに対応していません。

See also

参考文献