生存分析の基礎

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

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

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

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

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

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

打ち切りがなければ、生存時間を応答変数とする通常の回帰分析が使えます。しかし打ち切りデータは「真の値は観測値以上である」という不等式情報であり、通常の残差（ $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$ には数えません。

非情報的打ち切りのもとで 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)}

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

Log-rank 検定

Log-rank 検定は、2つ以上の群の生存曲線が同じかどうかを検定するノンパラメトリック検定です。

各イベント時刻 $t_i$ で、群 $j$ の観測イベント数 $d_{ij}$ と期待イベント数 $e_{ij} = n_{ij} \cdot d_i / n_i$ を計算します。 $n_{ij}$ は時刻 $t_i$ における群 $j$ のリスク集合の大きさです。帰無仮説のもとでは、各時点のイベントがリスク集合の群構成比に応じて配分されると期待されます。この配分はハイパー幾何分布に従い、分散もそこから導かれます。2群の場合、 $O_1 + O_2 = \sum d_i$ かつ $E_1 + E_2 = \sum d_i$ であるため $O_1 - E_1 = -(O_2 - E_2)$ が成り立ち、一方の群の情報で検定統計量を構成できます:

\chi^2 = \frac{(O_1 - E_1)^2}{V_1}

$O_1 = \sum_i d_{1i}$ （群1の総観測イベント数）、 $E_1 = \sum_i e_{1i}$ （総期待イベント数）、 $V_1 = \sum_i \frac{n_{1i} \cdot n_{2i} \cdot d_i (n_i - d_i)}{n_i^2(n_i - 1)}$ （分散）です。分母は期待値 $E_1$ ではなく、ハイパー幾何分布から導かれる分散 $V_1$ です。3群以上の場合は分散共分散行列を用いた二次形式 $\chi^2 = \mathbf{U}' \mathbf{V}^{-1} \mathbf{U}$ に拡張します。帰無仮説のもとでこの統計量は自由度 $g - 1$ （ $g$ は群数）のカイ二乗分布に近似的に従います。

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

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$ の推定後、Breslow 推定量によって $h_0(t)$ をノンパラメトリックに推定し、特定の共変量値に対する生存関数 $S(t|X)$ を求めることもできます。ただし MIDAS は現在、ベースラインハザードの推定を出力していません。

このモデルの共変量 $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 残差や log(-log) プロットなどの方法がありますが、MIDAS は現在これらの診断機能を備えていません。

部分尤度

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)$ を知る必要がなくなります。部分尤度は通常の尤度ではありませんが、漸近的に通常の最尤推定量と同じ性質（一致性、漸近正規性）を持つことが示されています (Cox, 1975)。

ハザード比の解釈

$\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$ はハザードに影響しない

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

同時イベント（タイ）の処理

複数の対象が同じ時刻にイベントを経験した場合、部分尤度の正確な計算は組み合わせ的に困難になります。MIDAS は Efron (1977) の近似を使用しています。Efron 法はタイのある対象をリスク集合から段階的に除外するため、タイのない場合の正確な部分尤度に対して Breslow 法よりも正確な近似を与えます。

参考文献

Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), 187-220.
Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481.
Cox, D. R. (1975). Partial likelihood. Biometrika, 62(2), 269-276.
Struthers, C. A., & Kalbfleisch, J. D. (1986). Misspecified proportional hazard models. Biometrika, 73(2), 363-369.
Efron, B. (1977). The efficiency of Cox's likelihood function for censored data. Journal of the American Statistical Association, 72(359), 557-565.