前回は生存関数とハザード関数の関係性について紹介しました。
生存関数の推定には、確率分布を仮定しない（ノンパラメトリックな）推定法を用います。
これらについて定義を見直し、lifelinesを使ってカプラン・マイヤー推定量を算出・プロットする方法についてまとめてみました。
経験分布による推定量 / 経験生存関数
打ち切りのない生存時間データがあるとき、このデータの生存時間関数 $S(t)$ は、下記の式で推定することができます。
$$ \hat{S}(t) = \frac{（生存時間 \geq t） である被験者数}{データセットに含まれる被験者数} $$
経験生存関数では、時間 $t$ における生存関数の推定に $t$ 以前に生存時間が打ち切りとなった被験者の情報を使用することはできないため、打ち切りのない生存データにしか用いることができません。
データの読み込み
今回はlifelinesパッケージで提供されているテストデータ（世界中の政治指導者の生涯データ）を用いて生存時間関数の推定について見ていきます。
なお、このデータでは、誕生に対応するイベントは在職期間の始まり、死亡に対応するイベントは個人の退職となっています。
from lifelines.datasets import load_dd data = load_dd() data.head() # ctryname cowcode2 politycode un_region_name ... regime start_year duration observed # 0 Afghanistan 700 700.0 Southern Asia ... Monarchy 1946 7 1 # 1 Afghanistan 700 700.0 Southern Asia ... Civilian Dict 1953 10 1 # 2 Afghanistan 700 700.0 Southern Asia ... Monarchy 1963 10 1 # 3 Afghanistan 700 700.0 Southern Asia ... Civilian Dict 1973 5 0 # 4 Afghanistan 700 700.0 Southern Asia ... Civilian Dict 1978 1 0 # # [5 rows x 12 columns] T = data["duration"] # 在職期間 E = data["observed"] # 退職が観測されたか否か
生命表法/保険数理法推定量
死亡がカレンダー日付のある単位で収集されるような大規模試験では、生命表法/保険数理法推定量が用いられることがあります。この推定量では、観測の期間を連続した区間に分割し、個々の区間での生存率を掛け合わせせたものとして定義します。
区間を $m$ 個に分割し（$j$ 番目を $t'_j$ から $t'_{j+1}$までの区間とする）、この区間における死亡数を $d_j$、打ち切り数を$c_j$であらわします。
更に$j$ 番目の区間の始点で生存しており、死亡のリスクがある被験者数を $n_j$ とします。
このとき、区間内で打ち切りが一様に発生するという仮定を置くと、区間内で死亡のリスクがある被験者数は、
$$
n_j' = n_j - c_j / 2
$$
となります。
また、$j$ 番目の区間における生存確率は、$(n_j' - d_j) / n_j'$ となります。
被験者が $t_k$ を超えて生存した確率を考えると、これは $k$ 番目の区間の始点を超えて生存する確率と、 $k-1$ 個それぞれの区間を生存した確率の積となります。
したがって、$t_k \leq t \lt t_{k+1}'$ 番目の区間における生命表法による生存関数の推定量は、
$$ S^*(t) = \prod_{j=1}^k\frac{n_j'-d_j}{n_j'} $$
とあらわせます。
この定義では、生存時間を任意にグループ化することで生存関数を導出するため、情報の損失が生じてしまいます。
特に症例数が少ない場合、この影響は大きくなります。
カプラン・マイヤー推定量
グループ化していない生存時間データを分析する手法として、頻繁に用いられるのがカプラン・マイヤー推定になります。
カプラン・マイヤー推定量を算出する際には、1つ1つの死亡時刻がそれぞれの区間の始点となるようにデータを分割します。
一般に $n$ 人の被験者について、観測された生存時間を $t_1, t_2, ..., t_n$ とします（$t_1 \lt t_2 \lt ... \lt t_n$）。これらの被験者について、 $r \leq n$ を満たす $r$ 個の死亡時刻があるとしたとき、 $j$ 番目を $t_j$と表します。また、 $t_j$ の直前に生存していた被験者の数を $n_j$、この時間で死亡する $d_j$ と表します。
このとき、区間 $t_j$ を通じて生存する確率は $\frac{n_j-d_j}{n_j}$となり、生存関数 $\hat{S}(t)$ は、
$$ \hat{S}(t) = \prod_{j=1}^k\frac{n_j-d_j}{n_j} $$ と表せます。
データに打ち切りがない場合のカプラン・マイヤー推定量を考えると、
$$ \hat{S}(t) = \frac{n_2}{n_1} \times \frac{n_3}{n_2} \times \cdots \times \frac{n_{k+1}}{n_k} $$
$k = 1, 2, \cdots, r-1$ について $n_{k+1}/n_1$ となります。
即ちカプラン・マイヤー推定量は打ち切りを考慮できるように経験生存関数を一般化したものといえます。
カプラン・マイヤー推定量の算出
カプラン・マイヤー推定量を算出、プロットするには、KaplanMeierFitter関数を使用します。
from lifelines import KaplanMeierFitter kmf = KaplanMeierFitter() kmf.fit(durations=T, event_observed=E) kmf.plot()
なお、生存関数の推定値は、kmf.survival_function内に格納されています。
kmf.survival_function_.head() # KM_estimate # timeline # 0.0 1.000000 # 1.0 0.721792 # 2.0 0.601973 # 3.0 0.510929 # 4.0 0.418835
ネルソン・アーレン推定量
生存関数は累積ハザード関数と1対1に対応しているため、累積ハザード関数の推定値から生存関数を推定することができます。 この推定値は Nelson (1986) が最初に提案したことから，ネルソン（Nelson）推定量，ネルソン‐アーレン（Nelson-Aalen）推定量などと呼ばれます。 前の記事でかいたように、累積ハザード関数$H(t)=\int_0^t{h(u)du}$ と生存時間関数 $S(t)$ の間には下記のような関係が成り立ちます。 $$ S(t) = \exp\{-H(t)\} $$ ところで、各時点でのハザード関数 $h(t)$ は下記であらわされるため、
$$ h(t)=\lim_{\delta t \to 0}\frac{P(t\leq T<t+\delta t|T\geq t)}{\delta t} $$
$\delta t$ を各イベント間の期間とみなすと、累積ハザード関数 $H(t)$ は $$ \hat{H}(t) = \sum_{j=1}^k{(d_j/n_j)} $$ で近似的に表せます。 この累積ハザード関数から算出した生存関数がネルソン・アーレン推定量となります。 $$ \tilde{S}(t) = \prod_{j=1}^k\exp{(-d_j/n_j)} $$ カプラン・マイヤー推定量はネルソン・アーレン推定量の近似ということができます。 $$ e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots $$ から、xの値が小さいとき、 $$ e^x \simeq 1 + x $$ となります。 終点に近い時点を除き、$d_j/n_j$は十分に小さくなるため、 $$ \exp{(-d_j/n_j)} \simeq 1 + (-d_j/n_j) = (n_j -d_j)/n_j $$ が成り立ち、この右辺はカプラン・マイヤー推定量の足し合わせ項に等しくなっています。 任意の $x$ について、 $e^x \geq 1+x$ なので、ネルソン・アーレン推定量はカプラン・マイヤー推定量よりも大きい値となります。
ネルソン・アーレン推定量の算出
lifelinesでは、累積ハザード関数の推定量 $\hat{H}(t)$ がNelsonAalenFitterとして実装されています。
from lifelines import NelsonAalenFitter naf = NelsonAalenFitter() naf.fit(durations=T, event_observed=E) naf.plot()
累積ハザード関数の推定値は、naf.cumulative_hazard_内に格納されています。
naf.cumulative_hazard_.head() # NA_estimate # timeline # 0.0 0.000000 # 1.0 0.325912 # 2.0 0.507356 # 3.0 0.671251 # 4.0 0.869867
カプラン・マイヤー推定量とネルソン・アーレン推定量の比較
カプラン・マイヤー推定量とネルソン・アーレン推定量から求めた生存関数について、両者をプロットすると実質的にほとんど差がないことが分かります。
plt.figure() ax = kmf.survival_function_.plot(drawstyle="steps-post", style=['r--']) np.exp(-naf.cumulative_hazard_).plot(drawstyle="steps-post", linewidth=1, ax=ax, style=['b--']) plt.show()
