「ある製品が故障するまでの時間はどれくらいか」「がん患者が治療を受けてから再発するまでの期間はどの程度か」「顧客がサービスを契約してから解約するまでの期間はどう分布しているか」— このような「イベントが起こるまでの時間」を分析する統計手法が生存時間分析（survival analysis）です。英語では”time-to-event analysis”とも呼ばれ、医学・工学・社会科学で広く使われています。
生存時間分析の最大の特徴は、打ち切り（censoring）と呼ばれる不完全な観測を適切に扱えることです。たとえば、5年間の臨床試験で追跡中の患者が研究終了時にまだ生存している場合、その患者の正確な生存時間はわかりませんが、「少なくとも5年以上は生存している」という部分的な情報は得られます。打ち切りを無視して分析すると、生存時間を過小評価するなどの深刻なバイアスが生じます。
生存時間分析を理解すると、以下のような応用が可能です。
- 臨床医学: 新薬の生存期間改善効果の評価や、治療法の比較に不可欠です。FDAの新薬承認審査でも生存時間解析が標準的に要求されます
- 信頼性工学: 製品の寿命分布の推定、故障率の解析、予防保全のスケジュール設計に使われます
- マーケティング: 顧客の生涯価値（LTV）の予測、解約率（チャーン率）の分析に利用されます
- 社会科学: 失業期間、婚姻の持続期間、再犯までの期間などの分析に適用されます
本記事の内容
- 生存時間分析の基本概念（生存関数、ハザード関数、累積ハザード関数）
- 打ち切りデータの種類とその重要性
- パラメトリック生存モデルの概観
- カプラン・マイヤー推定量の直感的理解と最尤推定としての導出
- グリーンウッドの公式による分散推定の導出
- 信頼区間の構成法（線形スケール・対数変換・対数-対数変換）
- ネルソン・アーレン推定量との比較
- 中央生存時間と制限付き平均生存時間（RMST）
- Pythonによるカプラン・マイヤー曲線の実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 確率分布と累積分布関数 — CDF・PDFの定義
- 最尤推定の理論 — カプラン・マイヤーの導出に使用
- 条件付き確率 — 条件付き生存確率の概念
生存時間分析の基本概念
生存関数
生存時間 $T$ は、ある起点（診断日、製造日、契約日など）からイベント（死亡、故障、解約など）が発生するまでの時間を表す非負の連続確率変数です。
生存関数（survival function）$S(t)$ は、時刻 $t$ 以降もイベントが発生していない確率です。
$$ \begin{equation} S(t) = P(T > t) = 1 – F(t) \end{equation} $$
ここで $F(t) = P(T \leq t)$ は累積分布関数（CDF）です。生存関数は $S(0) = 1$（最初は全員が「生存」）から始まり、$t \to \infty$ で $S(t) \to 0$ に単調減少します。
日常の言葉で言えば、$S(5) = 0.7$ は「5年後に生存している確率が $70\%$」を意味します。$S(t)$ は医学では「5年生存率」、工学では「信頼度」（reliability）と呼ばれることもあります。
生存関数の基本的な性質をまとめておきましょう。
- $S(0) = 1$（時刻0では全員が生存）
- $\lim_{t \to \infty} S(t) = 0$（十分長い時間が経てば全員がイベントを経験）
- $S(t)$ は単調非増加（生存確率は時間とともに下がるか、一定）
- $S(t)$ は右連続（離散的なイベント時刻で階段的に下がる場合も含む）
性質2について補足すると、治癒モデル（cure model）のように「一定割合の個体は永久にイベントを経験しない」場合は $\lim_{t \to \infty} S(t) = p_{\text{cure}} > 0$ となります。しかし標準的な生存時間分析では性質2を仮定します。
生存関数から導かれる重要な量として、分位点があります。$p$ 分位生存時間は $S(t_p) = 1 – p$ を満たす $t_p$ であり、特に $p = 0.5$ のとき中央生存時間（median survival time）と呼ばれます。中央生存時間は生存時間の代表値として最もよく報告される統計量です。
ハザード関数
ハザード関数（hazard function, hazard rate）$h(t)$ は、時刻 $t$ まで生存した個体が、次の微小時間 $\Delta t$ の間にイベントを経験する「瞬間的なリスク」を表します。
$$ \begin{equation} h(t) = \lim_{\Delta t \to 0} \frac{P(t \leq T < t + \Delta t \mid T \geq t)}{\Delta t} = \frac{f(t)}{S(t)} \end{equation} $$
ここで $f(t) = F'(t) = -S'(t)$ は確率密度関数です。ハザード関数は確率ではなく「率」（rate）であるため、1を超える値を取ることもあります。単位は「1/時間」であり、$h(t) = 0.05$ ならば「時刻 $t$ で生存している個体が次の1単位時間にイベントを経験する瞬間的な率が5%/年」のように解釈します。
ハザード関数は、密度関数 $f(t)$ と異なり「$t$ まで生存していること」を条件としている点が重要です。これにより、「現在まで生き延びた個体にとってのリスク」を直接表現できます。
ハザード関数と生存関数の間には以下の基本的な関係があります。
$$ S(t) = \exp\left(-\int_0^t h(u)\,du\right) = \exp(-H(t)) $$
ここで $H(t) = \int_0^t h(u)\,du$ は累積ハザード関数（cumulative hazard function）です。
この関係式の導出を行いましょう。$h(t) = f(t)/S(t)$ かつ $f(t) = -S'(t)$ なので、
$$ h(t) = \frac{-S'(t)}{S(t)} = -\frac{d}{dt}\ln S(t) $$
両辺を $[0, t]$ で積分すると、
$$ \int_0^t h(u)\,du = -[\ln S(u)]_0^t = -\ln S(t) + \ln S(0) $$
$S(0) = 1$ より $\ln S(0) = 0$ なので、
$$ H(t) = -\ln S(t) \quad \Longleftrightarrow \quad S(t) = e^{-H(t)} $$
この関係は生存関数・ハザード関数・累積ハザード関数の間の完全な等価性を示しています。いずれか1つが決まれば他の2つも決まります。
ハザード関数の形状は、「加齢」や「摩耗」の効果を反映します。
- 増加ハザード: $h(t)$ が $t$ とともに増加する場合。時間が経つほどイベントのリスクが高まります。老化、摩耗、疲労が進行する状況に対応します
- 減少ハザード: $h(t)$ が $t$ とともに減少する場合。初期のリスクが高く、生き残った個体のリスクは低くなります。初期故障（infant mortality）、手術後の早期死亡に対応します
- 一定ハザード: $h(t) = \lambda$（定数）の場合。リスクが時間に依存せず、指数分布に対応します。指数分布の無記憶性と直接結びついています
- バスタブ型ハザード: 初期に減少→一定→後期に増加するU字型。工業製品でよく見られるパターンで、初期故障・偶発故障・摩耗故障の3つのフェーズに対応します
ハザード関数の形状を理解することは、適切なパラメトリックモデルの選択や、Cox比例ハザードモデルの仮定の検証に不可欠です。
打ち切りデータ
生存時間分析で最も重要な概念が打ち切り（censoring）です。打ち切りとは、対象のイベント時間を正確に観測できない状況のことです。
右打ち切り（right censoring）が最も一般的で、以下のような状況で生じます。
- タイプI打ち切り（研究終了による打ち切り）: 研究期間が事前に決められており、期間終了時にイベントが起きていない個体が打ち切りになります。すべての打ち切り時刻が同じです
- タイプII打ち切り: イベントが一定数発生した時点で研究を終了します。残りの個体は打ち切りになります
- ランダム打ち切り: 各個体が固有の打ち切り時刻を持ちます。追跡不能（転居、連絡喪失）や競合リスク（関心のないイベント）がこれに当たります。最も一般的なタイプです
打ち切りのあるデータの $i$ 番目の観測は、$(t_i, \delta_i)$ のペアで表されます。 – $t_i = \min(T_i, C_i)$: 観測時間（真のイベント時間 $T_i$ と打ち切り時刻 $C_i$ の小さい方） – $\delta_i = \mathbf{1}(T_i \leq C_i)$: イベント指示変数（$\delta_i = 1$ でイベント発生、$\delta_i = 0$ で打ち切り）
打ち切りの重要な仮定は非情報的打ち切り（noninformative censoring, independent censoring）です。これは、打ち切り時刻 $C_i$ がイベント時刻 $T_i$ と条件付き独立であることを意味します。つまり、打ち切りの時刻が将来のイベント発生リスクに関する情報を含まないという仮定です。
たとえば、「状態が悪化して通院できなくなり追跡不能になった」場合、打ち切りは情報的です（その患者は追跡不能でなければ早くイベントを経験した可能性が高い）。情報的打ち切りが存在するとカプラン・マイヤー推定量にバイアスが生じます。
他のタイプの打ち切りと切断
右打ち切りが最も一般的ですが、他のタイプの打ち切りも存在します。
左打ち切り（left censoring）: イベントが観測開始前に既に発生していた場合です。例えば、HIV検査で陽性と判明した患者の感染時期が不明な場合、感染時間は左打ち切りです。$T_i < C_i$ のとき $C_i$ のみが観測され、$T_i$ は $C_i$ 以下であることのみわかります。
区間打ち切り（interval censoring）: イベントが2つの観測時点の間に発生したことは分かるが、正確な時刻が不明な場合です。例えば、6ヶ月ごとの検診でがんが発見された場合、発症は前回検診と今回検診の間のどこかですが、正確な発症日は不明です。$(L_i, R_i]$ の区間のみが観測されます。
切断（truncation）は打ち切りと混同されやすい概念ですが、明確に異なります。
左切断（left truncation, delayed entry）: ある時間まで生存した個体のみが観測される場合です。例えば、60歳以上の退職者を対象とした研究では、60歳までに死亡した人はサンプルに含まれません。
打ち切りでは「個体は観測されるが、イベント時間の情報が不完全」なのに対し、切断では「条件を満たさない個体がそもそもサンプルに入らない」という違いがあります。切断を無視すると生存関数の推定に深刻なバイアスが生じるため、リスク集合の定義を適切に修正する必要があります（左切断された個体は切断時刻以前はリスク集合に含めない）。
なぜ打ち切りを無視できないか
打ち切りの重要性を具体例で理解しましょう。10人の患者を5年間追跡し、以下のデータが得られたとします。
- イベント（死亡）: 1年, 2年, 3年, 4年に各1人（計4人）
- 打ち切り: 1.5年, 2.5年, 3.5年に各1人（計3人）
- 5年時点で生存: 3人
3つの「間違った」分析方法を考えてみましょう。
間違った分析1（打ち切りを除外）: 打ち切りの3人 + 5年生存の3人を除外して、イベント4人 / 観測7人 = 57%の5年死亡率と計算。これは打ち切りまでの生存情報が捨てられており、選択バイアスが生じます。
間違った分析2（打ち切り = イベント）: 打ち切りの3人を「死亡」として扱い、7人 / 10人 = 70%の5年死亡率。これは生存時間を過小評価します。
間違った分析3（打ち切り = 生存）: 打ち切りの3人を「5年生存」として扱い、4人 / 10人 = 40%の5年死亡率。これは打ち切り後にイベントが起きた可能性を無視します。
正しい分析: カプラン・マイヤー法を使い、各時点でのリスク集合を動的に更新しながら生存確率を推定する。
打ち切りデータを適切に扱う方法として、カプラン・マイヤー推定量を導出しましょう。
パラメトリック生存モデルの概観
カプラン・マイヤー推定量はノンパラメトリックな方法ですが、生存時間の分布に特定の形を仮定するパラメトリック生存モデルも広く使われています。代表的な分布とそのハザード関数の特徴を整理しましょう。
指数分布: $S(t) = e^{-\lambda t}$, $h(t) = \lambda$（定数ハザード）。最もシンプルなモデルで、ポアソン過程のイベント間隔に対応します。無記憶性 $P(T > s + t \mid T > s) = P(T > t)$ を持ちます。
ワイブル分布: $S(t) = \exp(-(\lambda t)^p)$, $h(t) = p\lambda^p t^{p-1}$。$p > 1$ で増加ハザード（老化）、$p < 1$ で減少ハザード（初期故障）、$p = 1$ で指数分布に帰着します。ハザードの単調な変化をモデル化できる柔軟なモデルです。
対数正規分布: $T$ の対数が正規分布に従うモデルです。ハザード関数が非単調（最初は増加し、後に減少）な形を取ります。一部のがんの生存時間データに適しています。
ゴンペルツ分布: $h(t) = ae^{bt}$。ハザードが指数的に増加するモデルで、人間の死亡率のモデリングに古くから使われています（ゴンペルツの法則）。
対数ロジスティック分布: $S(t) = 1/(1 + (\lambda t)^p)$, $h(t) = p\lambda^p t^{p-1}/(1 + (\lambda t)^p)$。ハザード関数が $p > 1$ のとき増加してから減少する形を取ります。加速故障時間モデル（AFT）でよく使われます。
パラメトリックモデルの利点は、少ないパラメータで生存関数全体を記述できることと、データの範囲を超えた外挿が可能なことです。欠点は、分布の仮定が間違っていると推定にバイアスが生じることです。実務ではまずカプラン・マイヤー曲線で生存関数の形状を確認し、適切なパラメトリックモデルを選択するのが一般的です。
カプラン・マイヤー推定量
直感的な理解
カプラン・マイヤー推定量（Kaplan-Meier estimator, product-limit estimator）は、打ち切りを含むデータから生存関数 $S(t)$ を推定するノンパラメトリックな方法です。1958年にカプラン（E. L. Kaplan）とマイヤー（P. Meier）によって提案されました。
アイデアは次の通りです。生存関数 $S(t) = P(T > t)$ を、各イベント時刻での条件付き生存確率の積として分解します。
時刻 $t$ まで生存するには、最初のイベント時刻を越え、2番目のイベント時刻を越え、…、$t$ 以前のすべてのイベント時刻を越える必要があります。確率の乗法定理を使うと、
$$ S(t) = P(T > t_1) \times P(T > t_2 \mid T > t_1) \times P(T > t_3 \mid T > t_2) \times \cdots $$
各因子「時刻 $t_j$ を越える確率（$t_{j-1}$ を越えたという条件のもとで）」は、「$t_j$ の直前にリスクにある人数のうち、$t_j$ でイベントを起こさなかった人数の割合」で推定できます。これが「積極限」（product-limit）という名前の由来です。
数学的な定義
イベントが発生した異なる時刻を昇順に $t_{(1)} < t_{(2)} < \cdots < t_{(m)}$ とします（$m$ は異なるイベント時刻の数）。
各時刻 $t_{(j)}$ について: – $d_j$: 時刻 $t_{(j)}$ でイベントが発生した個数 – $n_j$: 時刻 $t_{(j)}$ の直前にリスクにある（まだイベントも打ち切りも起きていない）個体数（リスク集合のサイズ）
カプラン・マイヤー推定量は次のように定義されます。
$$ \begin{equation} \hat{S}(t) = \prod_{j: t_{(j)} \leq t} \left(1 – \frac{d_j}{n_j}\right) \end{equation} $$
$t < t_{(1)}$ のときは $\hat{S}(t) = 1$ です。各因子 $(1 - d_j/n_j)$ は、「$t_{(j)}$ の直前にリスクにある $n_j$ 人のうち、$t_{(j)}$ を生き延びた割合」であり、条件付き生存確率 $P(T > t_{(j)} \mid T \geq t_{(j)})$ の推定量です。
条件付きハザード確率（離散ハザード）を $\hat{q}_j = d_j/n_j$ と書くと、$\hat{S}(t) = \prod_{j:t_{(j)} \leq t} (1 – \hat{q}_j)$ です。
打ち切りの扱い
カプラン・マイヤー推定量が打ち切りを扱うメカニズムは、リスク集合（risk set）の動的な更新にあります。
リスク集合 $R(t)$ は「時刻 $t$ の直前にイベントを経験する可能性がある個体の集合」です。具体的には、
$$ R(t) = \{i : t_i \geq t\} $$
すなわち、時刻 $t$ 以前にイベントも打ち切りも起きていない個体の集合です。リスク集合のサイズ $n_j = |R(t_{(j)})|$ は次のように更新されます。
- 最初は全員（$n_1 = n$）
- 時刻 $t_{(j-1)}$ と $t_{(j)}$ の間に打ち切られた個体と、$t_{(j-1)}$ でイベントが起きた個体を除く
打ち切りが起きた時刻ではカプラン・マイヤー曲線はジャンプしません（生存確率の推定値は変わらない）が、リスク集合が小さくなるため、以降の推定の不確実性が増加します。
打ち切りの扱いにおいて重要なのは、打ち切り時刻とイベント時刻が同じ場合の処理です。標準的な慣例では、タイ（同時刻）がある場合、イベントが先に発生したと仮定します。つまり、イベントを起こした個体は $n_j$ に含め、同時刻で打ち切られた個体は $n_j$ には含めるがイベント数 $d_j$ にはカウントしません。
具体的な計算例
小さなデータで手計算をしてみましょう。6人の患者のデータは以下の通りです。
| 患者 | 観測時間 | イベント |
|---|---|---|
| A | 2 | 1 (死亡) |
| B | 3 | 0 (打ち切り) |
| C | 5 | 1 (死亡) |
| D | 5 | 1 (死亡) |
| E | 8 | 0 (打ち切り) |
| F | 10 | 1 (死亡) |
イベント時刻は $t_{(1)} = 2$, $t_{(2)} = 5$, $t_{(3)} = 10$ です。
$t_{(1)} = 2$: $n_1 = 6$（全員がリスク）, $d_1 = 1$
$$ \hat{S}(2) = 1 – \frac{1}{6} = \frac{5}{6} \approx 0.833 $$
$t_{(2)} = 5$: 患者Aはイベント済み（$n$ から除外）、患者Bは時刻3で打ち切り（$n$ から除外） → $n_2 = 4$, $d_2 = 2$
$$ \hat{S}(5) = \frac{5}{6} \times \left(1 – \frac{2}{4}\right) = \frac{5}{6} \times \frac{1}{2} = \frac{5}{12} \approx 0.417 $$
$t_{(3)} = 10$: 患者Eは時刻8で打ち切り → $n_3 = 1$（患者Fのみ）, $d_3 = 1$
$$ \hat{S}(10) = \frac{5}{12} \times \left(1 – \frac{1}{1}\right) = 0 $$
患者Bの打ち切り情報がどのように使われているか注目しましょう。患者Bは時刻3まで生存していたため、$t_{(1)} = 2$ のリスク集合には含まれますが（$n_1 = 6$）、$t_{(2)} = 5$ のリスク集合からは除外されます（$n_2 = 4$）。もし打ち切りを無視して全員をリスク集合に入れれば、$n_2 = 5$ となり、$\hat{S}(5) = 5/6 \times 3/5 = 1/2 = 0.5$ と異なる結果になります。打ち切りを適切に扱うことで、より正確な推定が得られるのです。
最尤推定としての導出
カプラン・マイヤー推定量は、実はノンパラメトリックな最尤推定量（NPMLE）です。この事実を示しましょう。
生存関数 $S(t)$ がイベント時刻 $t_{(1)}, t_{(2)}, \dots, t_{(m)}$ でのみジャンプする階段関数であると仮定します。各イベント時刻での条件付き死亡確率を
$$ q_j = P(T = t_{(j)} \mid T \geq t_{(j)}), \quad j = 1, \dots, m $$
とおきます。$\bm{q} = (q_1, q_2, \dots, q_m)$ がパラメータです。
打ち切りを含むデータの尤度関数を構成しましょう。時刻 $t_{(j)}$ でイベントが発生する確率は $q_j$（条件付き）、$t_{(j)}$ を越えて生存する確率は $1 – q_j$（条件付き）です。尤度は各イベント時刻での二項尤度の積として、
$$ L(\bm{q}) = \prod_{j=1}^{m} q_j^{d_j}(1 – q_j)^{n_j – d_j} $$
ここで $n_j – d_j$ は「$t_{(j)}$ の直前にリスクにあった $n_j$ 人のうちイベントを起こさなかった人数」で、$t_{(j)}$ を超えて生存した人と $t_{(j-1)}$ と $t_{(j)}$ の間で打ち切られた人を含みます。
各 $q_j$ は他の $q_k$ とは独立に最大化できます（尤度が $q_j$ ごとに分解されるため）。$q_j$ について微分すると、
$$ \frac{\partial \ln L}{\partial q_j} = \frac{d_j}{q_j} – \frac{n_j – d_j}{1 – q_j} = 0 $$
$$ \hat{q}_j = \frac{d_j}{n_j} $$
これは二項分布の最尤推定量です。したがって、
$$ \hat{S}(t) = \prod_{j: t_{(j)} \leq t} (1 – \hat{q}_j) = \prod_{j: t_{(j)} \leq t}\left(1 – \frac{d_j}{n_j}\right) $$
はノンパラメトリック最尤推定量です。この導出は、カプラン・マイヤー推定量が「自然な」推定量であることを保証しています。
分散推定と信頼区間
グリーンウッドの公式の導出
カプラン・マイヤー推定量の分散は、グリーンウッドの公式（Greenwood’s formula）で推定されます。
$$ \begin{equation} \widehat{\text{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{j: t_{(j)} \leq t} \frac{d_j}{n_j(n_j – d_j)} \end{equation} $$
この公式はデルタ法（delta method）を用いて導出されます。$\hat{S}(t) = \prod_{j}(1 – \hat{q}_j)$ の対数をとると
$$ \ln\hat{S}(t) = \sum_{j: t_{(j)} \leq t} \ln(1 – \hat{q}_j) $$
各 $\hat{q}_j = d_j/n_j$ は二項分布に基づく推定量で、異なるイベント時刻の推定量は条件付き独立です。$\hat{q}_j$ の分散は、
$$ \text{Var}(\hat{q}_j) \approx \frac{q_j(1 – q_j)}{n_j} $$
デルタ法を適用します。$g(q) = \ln(1 – q)$ に対して $g'(q) = -1/(1-q)$ なので、
$$ \text{Var}(\ln(1 – \hat{q}_j)) \approx [g'(q_j)]^2 \cdot \text{Var}(\hat{q}_j) = \frac{1}{(1 – q_j)^2} \cdot \frac{q_j(1 – q_j)}{n_j} = \frac{q_j}{n_j(1 – q_j)} $$
$q_j$ を $\hat{q}_j = d_j/n_j$ で推定すると、
$$ \widehat{\text{Var}}(\ln(1 – \hat{q}_j)) = \frac{d_j/n_j}{n_j(1 – d_j/n_j)} = \frac{d_j}{n_j(n_j – d_j)} $$
$\ln\hat{S}(t)$ の分散は各項の分散の和です（独立性より）。
$$ \widehat{\text{Var}}(\ln\hat{S}(t)) = \sum_{j: t_{(j)} \leq t} \frac{d_j}{n_j(n_j – d_j)} $$
最後に、$\hat{S}(t) = \exp(\ln\hat{S}(t))$ に再びデルタ法を適用して、
$$ \widehat{\text{Var}}(\hat{S}(t)) = \hat{S}(t)^2 \cdot \widehat{\text{Var}}(\ln\hat{S}(t)) = \hat{S}(t)^2 \sum_{j: t_{(j)} \leq t} \frac{d_j}{n_j(n_j – d_j)} $$
これがグリーンウッドの公式です。標準誤差は $\text{SE}[\hat{S}(t)] = \sqrt{\widehat{\text{Var}}[\hat{S}(t)]}$ です。
信頼区間
点ごとの $95\%$ 信頼区間は複数の方法で構成できます。
線形スケール（最も単純）:
$$ \hat{S}(t) \pm 1.96 \times \text{SE}[\hat{S}(t)] $$
ただし、この方法では信頼区間が $[0, 1]$ の範囲を超えることがあります。特に $\hat{S}(t)$ が0や1に近い場合に問題が生じます。
対数変換:
$\ln\hat{S}(t)$ の信頼区間を構成し、指数変換で $[0, 1]$ の範囲に戻します。
$$ \hat{S}(t) \cdot \exp\left(\pm \frac{1.96 \times \text{SE}[\hat{S}(t)]}{\hat{S}(t)}\right) $$
上界が1を超える可能性はありますが、下界は自動的に正になります。
対数-対数変換（推奨、多くの統計ソフトのデフォルト）:
$\ln(-\ln\hat{S}(t))$ に対して信頼区間を構成します。この変換は $[0, 1]$ を $(-\infty, \infty)$ に写すため、信頼区間が $[0, 1]$ に収まることが保証されます。
$$ \hat{S}(t)^{\exp(\pm 1.96 \cdot w)}, \quad w = \frac{\text{SE}[\hat{S}(t)]}{\hat{S}(t) \cdot |\ln\hat{S}(t)|} $$
小標本での被覆確率（信頼区間が真の値を含む確率）の観点からは、対数-対数変換に基づく信頼区間が最も良い性能を示すことが知られています。
ネルソン・アーレン推定量
累積ハザード関数の推定
カプラン・マイヤー推定量は生存関数を直接推定しますが、累積ハザード関数 $H(t) = -\ln S(t)$ を推定する別のアプローチもあります。ネルソン・アーレン推定量（Nelson-Aalen estimator）は累積ハザード関数のノンパラメトリック推定量です。
$$ \begin{equation} \hat{H}(t) = \sum_{j: t_{(j)} \leq t} \frac{d_j}{n_j} \end{equation} $$
カプラン・マイヤー推定量が積で構成されるのに対し、ネルソン・アーレン推定量は和で構成されています。
直感的には、$d_j/n_j$ はイベント時刻 $t_{(j)}$ での瞬間的なハザード率の推定量であり、これを足し合わせることで累積ハザードを推定しています。$H(t) = \int_0^t h(u) \, du$ の離散近似です。
ネルソン・アーレン推定量の分散
ネルソン・アーレン推定量の分散は次のように推定されます。
$$ \widehat{\text{Var}}[\hat{H}(t)] = \sum_{j: t_{(j)} \leq t} \frac{d_j}{n_j^2} $$
グリーンウッドの公式と比較すると、分母が $n_j^2$ のみで $n_j(n_j – d_j)$ ではない点が異なります。$d_j \ll n_j$（リスク集合が十分大きい）の場合は $n_j(n_j – d_j) \approx n_j^2$ なので両者はほぼ一致しますが、$d_j$ が $n_j$ に近い場合（リスク集合が小さい場合）に差が生じます。
カプラン・マイヤー推定量との関係
ネルソン・アーレン推定量から生存関数を推定することもできます。
$$ \hat{S}_{\text{NA}}(t) = \exp(-\hat{H}(t)) = \exp\left(-\sum_{j: t_{(j)} \leq t} \frac{d_j}{n_j}\right) $$
カプラン・マイヤー推定量との関係は、$\ln(1-x) \approx -x$（$x$ が小さいとき）から理解できます。
$$ \ln\hat{S}_{\text{KM}}(t) = \sum_j \ln\left(1 – \frac{d_j}{n_j}\right) \approx -\sum_j \frac{d_j}{n_j} = -\hat{H}(t) = \ln\hat{S}_{\text{NA}}(t) $$
$d_j/n_j$ が小さい（リスク集合が大きい）場合、カプラン・マイヤー推定量とネルソン・アーレン推定量はほぼ一致します。ネルソン・アーレン推定量の利点は、マルチンゲール理論に基づく統計的性質が扱いやすいことです。
中央生存時間と制限付き平均生存時間
生存時間の要約統計量として、以下が重要です。
中央生存時間: $\hat{S}(t) = 0.5$ を満たす最小の $t$ です。カプラン・マイヤー曲線が0.5を横切る時刻として読み取れます。生存時間の分布は通常右に歪んでいるため、中央値は平均値よりも代表性が高い指標として好まれます。中央生存時間が定義できない場合（全観測期間を通じて $\hat{S}(t) > 0.5$ の場合）は「not reached」と報告されます。
制限付き平均生存時間（RMST, Restricted Mean Survival Time）: 特定の時間 $\tau$ までのカプラン・マイヤー曲線の下の面積です。
$$ \text{RMST}(\tau) = \int_0^\tau \hat{S}(t)\,dt $$
RMSTは「$\tau$ までの期間に期待される生存時間」と解釈できます。たとえば、$\tau = 5$ 年での2群のRMSTの差 $\Delta\text{RMST} = \text{RMST}_{\text{治療}} – \text{RMST}_{\text{対照}}$ は「5年間で治療群は対照群に比べて平均して何ヶ月長く生存するか」を表し、臨床的に直感的な解釈が可能です。
RMSTは近年注目を集めている指標で、ハザードの比例性を仮定しないため、Cox比例ハザードモデルの仮定が成り立たない場合の代替指標として使われます。
カプラン・マイヤー推定量の統計的性質
カプラン・マイヤー推定量の統計的性質をまとめておきましょう。
一致性: サンプルサイズ $n \to \infty$ のとき、$\hat{S}(t) \xrightarrow{p} S(t)$（確率収束）。打ち切りが非情報的であれば、推定値は真の生存関数に収束します。
漸近正規性: 適切に標準化すると $\sqrt{n}(\hat{S}(t) – S(t))$ は正規分布に収束します。これがグリーンウッドの公式に基づく信頼区間の理論的根拠です。
ノンパラメトリック最尤推定量（NPMLE）: カプラン・マイヤー推定量は、離散的なイベント時刻でのみジャンプする階段関数の中で、尤度を最大化するものです。
ジャンプのタイミング: カプラン・マイヤー曲線はイベント時刻でのみジャンプし、打ち切り時刻ではジャンプしません。これは打ち切りの情報がリスク集合を通じてのみ推定に反映されることを意味しています。
自己一致性（self-consistency）: カプラン・マイヤー推定量はTurnbull（1976）の自己一致アルゴリズムの解でもあり、EMアルゴリズムと関連しています。
Pythonによる実装
カプラン・マイヤー推定量のスクラッチ実装
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
def kaplan_meier(times, events):
"""
カプラン・マイヤー推定量の計算
Parameters
----------
times : array 観測時間
events : array イベント指示変数（1=イベント, 0=打ち切り）
Returns
-------
result : dict KM推定量の各値
"""
times = np.asarray(times, dtype=float)
events = np.asarray(events, dtype=int)
n = len(times)
# ユニークなイベント時刻（イベント発生時のみ）
event_times = np.unique(times[events == 1])
event_times.sort()
# 各イベント時刻でのリスク集合とイベント数
km_times = [0.0]
km_survival = [1.0]
km_var = [0.0]
survival = 1.0
var_sum = 0.0
for t in event_times:
# リスク集合: t以上の観測時間を持つ個体数
n_risk = np.sum(times >= t)
# イベント数: ちょうどtでイベントが発生した数
d = np.sum((times == t) & (events == 1))
if n_risk > 0 and n_risk > d:
survival *= (1 - d / n_risk)
var_sum += d / (n_risk * (n_risk - d))
elif n_risk == d:
survival = 0.0
km_times.append(t)
km_survival.append(survival)
km_var.append(survival**2 * var_sum if survival > 0 else 0.0)
km_times = np.array(km_times)
km_survival = np.array(km_survival)
km_se = np.sqrt(np.array(km_var))
# 95%信頼区間（対数-対数変換）
ci_lower = np.ones_like(km_survival)
ci_upper = np.ones_like(km_survival)
for i in range(len(km_survival)):
s = km_survival[i]
se = km_se[i]
if 0 < s < 1 and se > 0:
w = se / (s * abs(np.log(s)))
ci_lower[i] = s ** np.exp(1.96 * w)
ci_upper[i] = s ** np.exp(-1.96 * w)
elif s == 1:
ci_lower[i] = 1.0
ci_upper[i] = 1.0
else:
ci_lower[i] = 0.0
ci_upper[i] = 0.0
# 中央生存時間
median_idx = np.where(km_survival <= 0.5)[0]
median_survival = km_times[median_idx[0]] if len(median_idx) > 0 else None
return {
"times": km_times, "survival": km_survival,
"se": km_se, "ci_lower": ci_lower, "ci_upper": ci_upper,
"median": median_survival
}
# --- テストデータ: 臨床試験 ---
np.random.seed(42)
# 治療群: ワイブル分布（中央生存時間が長い）
n_treat = 50
t_treat = np.random.weibull(1.5, n_treat) * 30
censor_treat = np.random.uniform(20, 60, n_treat)
observed_treat = np.minimum(t_treat, censor_treat)
event_treat = (t_treat <= censor_treat).astype(int)
# 対照群: ワイブル分布（中央生存時間が短い）
n_control = 50
t_control = np.random.weibull(1.5, n_control) * 20
censor_control = np.random.uniform(20, 60, n_control)
observed_control = np.minimum(t_control, censor_control)
event_control = (t_control <= censor_control).astype(int)
# KM推定
km_treat = kaplan_meier(observed_treat, event_treat)
km_control = kaplan_meier(observed_control, event_control)
print("=" * 55)
print("カプラン・マイヤー推定の結果")
print("=" * 55)
for name, km, n_obs, events in [
("治療群", km_treat, n_treat, event_treat),
("対照群", km_control, n_control, event_control)
]:
print(f"\n{name}: n={n_obs}, イベント={np.sum(events)}, "
f"打ち切り={n_obs - np.sum(events)}")
if km["median"] is not None:
print(f" 中央生存時間: {km['median']:.1f}")
else:
print(f" 中央生存時間: not reached")
# --- 可視化 ---
fig, ax = plt.subplots(figsize=(10, 6))
# 治療群
ax.step(km_treat["times"], km_treat["survival"], where="post",
linewidth=2, color="blue", label="Treatment")
ax.fill_between(km_treat["times"], km_treat["ci_lower"],
km_treat["ci_upper"], step="post",
alpha=0.15, color="blue")
# 対照群
ax.step(km_control["times"], km_control["survival"], where="post",
linewidth=2, color="red", label="Control")
ax.fill_between(km_control["times"], km_control["ci_lower"],
km_control["ci_upper"], step="post",
alpha=0.15, color="red")
# 打ち切りマーク
for km, obs, ev, color in [(km_treat, observed_treat, event_treat, "blue"),
(km_control, observed_control, event_control, "red")]:
cens_times = obs[ev == 0]
for ct in cens_times:
idx = np.searchsorted(km["times"], ct, side="right") - 1
idx = min(idx, len(km["survival"]) - 1)
ax.plot(ct, km["survival"][idx], "|", color=color,
markersize=8, markeredgewidth=1.5, alpha=0.5)
ax.axhline(0.5, color="gray", linestyle="--", linewidth=0.8, alpha=0.5)
ax.set_xlabel("Time (months)", fontsize=12)
ax.set_ylabel("Survival probability", fontsize=12)
ax.set_title("Kaplan-Meier survival curves with 95% CI", fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, None)
ax.set_ylim(0, 1.05)
plt.tight_layout()
plt.show()
このカプラン・マイヤー曲線から、以下のことが読み取れます。
-
階段関数の形状が生存関数の推定を示している。各ステップ（下降）はイベント（死亡・故障など）が発生した時刻に対応しています。平坦な部分は、打ち切りのみが発生した区間であり、その間はイベントが起きていないことを反映しています
-
帯状の領域は95%信頼区間を示している。時間が経つにつれてリスク集合が小さくなるため、信頼区間は広がっていきます。これは、残っている観測数が少なくなるにつれて推定の不確実性が増すことを反映しています。特に曲線の末端では信頼区間が非常に広くなるため、解釈には注意が必要です
-
縦のチックマーク（|）は打ち切りの時点を示している。これらの個体はその時点までは生存していたが、その後の状態は不明です。打ち切りはカプラン・マイヤー曲線のジャンプには寄与しませんが、リスク集合を減少させることで以降の推定に影響します
-
2群の曲線の分離度合いが治療効果を視覚的に示している。治療群（青）の曲線が対照群（赤）よりも上にあれば、治療群の方が生存が良好であることを示唆します。この差が統計的に有意かどうかは、ログランク検定で検定します
ネルソン・アーレン推定量との比較
import numpy as np
import matplotlib.pyplot as plt
def nelson_aalen(times, events):
"""ネルソン・アーレン推定量"""
times = np.asarray(times, dtype=float)
events = np.asarray(events, dtype=int)
event_times = np.sort(np.unique(times[events == 1]))
na_times = [0.0]
na_cumhaz = [0.0]
na_var = [0.0]
cumhaz = 0.0
var_sum = 0.0
for t in event_times:
n_risk = np.sum(times >= t)
d = np.sum((times == t) & (events == 1))
if n_risk > 0:
cumhaz += d / n_risk
var_sum += d / n_risk**2
na_times.append(t)
na_cumhaz.append(cumhaz)
na_var.append(var_sum)
return {
"times": np.array(na_times),
"cumhaz": np.array(na_cumhaz),
"survival": np.exp(-np.array(na_cumhaz)),
"se": np.sqrt(np.array(na_var))
}
na_treat = nelson_aalen(observed_treat, event_treat)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# KM vs NA 生存関数
ax = axes[0]
ax.step(km_treat["times"], km_treat["survival"], where="post",
linewidth=2, color="blue", label="Kaplan-Meier")
ax.step(na_treat["times"], na_treat["survival"], where="post",
linewidth=2, color="red", linestyle="--", label="Nelson-Aalen")
ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Survival probability", fontsize=12)
ax.set_title("KM vs Nelson-Aalen (treatment)", fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
# 累積ハザード関数
ax = axes[1]
km_cumhaz = -np.log(np.clip(km_treat["survival"], 1e-10, 1))
ax.step(km_treat["times"], km_cumhaz, where="post",
linewidth=2, color="blue", label="$-\\ln\\hat{S}_{KM}$")
ax.step(na_treat["times"], na_treat["cumhaz"], where="post",
linewidth=2, color="red", linestyle="--", label="Nelson-Aalen")
ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Cumulative hazard", fontsize=12)
ax.set_title("Cumulative hazard estimates", fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この比較から、以下のことが確認できます。
-
カプラン・マイヤー推定量とネルソン・アーレン推定量はほぼ一致している。リスク集合が大きい場合（序盤）は2つの推定量はほとんど区別がつきません。差が生じるのはリスク集合が小さくなる終盤のみです
-
累積ハザード関数の比較（右図）: カプラン・マイヤー推定量から計算した $-\ln\hat{S}_{\text{KM}}(t)$ とネルソン・アーレン推定量 $\hat{H}_{\text{NA}}(t)$ がほぼ一致しています。これは $\ln(1-x) \approx -x$ の近似が良好であることを示しています
生存時間分析の発展的トピック
本記事ではカプラン・マイヤー推定量を中心に解説しましたが、生存時間分析にはさらに多くの重要なトピックがあります。
回帰モデル
生存時間に共変量（年齢、性別、治療法など）が影響する場合、回帰モデルが必要です。Cox比例ハザードモデルは、ハザード関数を $h(t \mid \bm{x}) = h_0(t) \exp(\bm{\beta}^{\top}\bm{x})$ と表現するセミパラメトリックモデルです。ベースラインハザード $h_0(t)$ の形状を仮定しないため、非常に柔軟です。
生存曲線の比較
2群以上の生存曲線が統計的に異なるかを検定するために、ログランク検定が最も広く用いられます。ログランク検定はハザード関数の比が時間を通じて一定（比例ハザード）であるという仮定のもとで最も検出力が高い検定です。
競合リスク
複数のタイプのイベントが存在する場合（例: がんによる死亡と他の原因による死亡）、競合リスク分析が必要です。通常のカプラン・マイヤー推定量は、関心のあるイベント以外を打ち切りとして扱いますが、これは累積発生率（cumulative incidence）を過大推定するバイアスを生じさせます。累積発生関数（CIF）とFine-Grayモデルがこの問題に対処します。
まとめ
本記事では、生存時間分析の基本概念とカプラン・マイヤー推定量の理論を解説しました。
- 生存関数 $S(t) = P(T > t)$ は時刻 $t$ 以降の生存確率を表し、ハザード関数 $h(t) = f(t)/S(t)$ は時刻 $t$ での瞬間的なリスクを表します。両者は $S(t) = \exp(-H(t))$ で結びついています
- 打ち切りデータは生存時間分析に固有の概念であり、正確なイベント時間が観測できない場合の部分的な情報を活用します。非情報的打ち切りの仮定が分析の妥当性の鍵です
- カプラン・マイヤー推定量は条件付き生存確率の積として $\hat{S}(t) = \prod(1 – d_j/n_j)$ と定義され、ノンパラメトリック最尤推定量です。打ち切りはリスク集合の動的更新を通じて適切に処理されます
- グリーンウッドの公式 $\widehat{\text{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum d_j/[n_j(n_j-d_j)]$ で分散を推定し、対数-対数変換による信頼区間が推奨されます
- ネルソン・アーレン推定量は累積ハザード関数 $\hat{H}(t) = \sum d_j/n_j$ の推定量で、マルチンゲール理論との親和性が高いです
- 制限付き平均生存時間（RMST）はハザードの比例性を仮定しない要約指標として近年注目されています
次のステップとして、以下の記事も参考にしてください。
- Cox比例ハザードモデル — 共変量を含む生存時間モデル
- ログランク検定 — 生存曲線の群間比較
- ワイブル分布 — パラメトリック生存モデルの代表格