前回の記事はこちら
Pythonで学ぶ生存時間解析1 -生存時間データとは
https://qiita.com/Goriwaku/items/8d00696d853da73505bd

今回使用しているデータについては前回の記事を参考にしていただくか、私のgit(https://github.com/goriwaku/survival_analysis )からダウンロード(whas500.csv)してください。

生存関数とは

　前回の記事で取り上げたWHAS500を使用して生存関数について解説していきます。連続変数としての時間tに対して、確率変数としての生存時間をTと定義します。この時、Tの累積分布関数はランダムに抽出された被験者がある時間t以下であることを表し、以下のように定義されます。

F(t)=Pr(T \leq t)

ここで、生存関数はある時間tよりも大きな生存時間Tを得る確率で表すことが可能で、

S(t)=Pr(T > t)

となります。このとき、任意の時点tにおいて、生存を死亡の余事象として考えることによって確率の公理により以下の等式が成立します。

S(t)=1-F(t)

ここからは、この生存関数の推定、特に右側打ち切りが存在しているときの生存時間関数の推定量であるKaplan-Meier推定量について見ていきます。

生存関数の推定

　実際に生存関数を推定していきます。例として、以下のような5人分の生存時間と生存状態から生存関数を推定することを考えます。

時点0においては、5人とも生存しているので$\hat{S}(t)=1.0$となります。6日目に被験者番号1が死亡するため、6日目の少し前から始まり、6日目で終了する微小区間$(6-\delta, 6]$を考えると、この微小区間で死亡する条件付確率の推定値が$1/5$で、生存確率は$4/5$となる。任意の時点において生存している被験者は死亡のリスクにさらされているリスク集合であり、この人数をnumber at riskと表します。生存時間の推定量は、(その時点まで生きている確率)×(上記微小区間で生存する条件付確率)によってあらわすことができます。以下、区間$I_i$を上の表の時間を昇順にして区切ったもの(ex: $I_0 = [0, 6)$)、$n_i$をその時点でのnumber at risk、$d_i$をその時点で死亡した人数とします。このとき、何らかのイベントが起こっている日に近い微小区間において、推定された死亡確率は$d_1/n_1$、生存確率は$(n_1 - d_1)/n_1$と表すことができます。よって、6日目、および14日目の生存確率を推定すると、

\hat{S}(6)=1 \times \frac{4}{5}=\frac{4}{5}\\
\hat{S}(14)=1 \times \frac{4}{5} \times \frac{3}{4} = \frac{3}{5}

となります。
次のイベントでは、21日目において被験者番号３が死亡に依らない打ち切りでat risk集合から脱落します。その微小区間において、$n_3=3$、$d_3=0$となることから、生存関数の推定値は、

\hat{S}(21)=1 \times \frac{4}{5} \times \frac{3}{4} \times \frac{3}{3} = \frac{3}{5}

当然ですが、区間$I_3$では死亡が発生していないため、生存関数の推定値は変化しません。以下の区間も同様に生存関数の推定値を求めることができます。この手法によって求めた生存関数の推定値をKaplan-Meier推定値といいます。

ここでKaplan-Meier推定値を階段関数としてプロットしてみましょう。

test_data = [[6, 1],[44, 1],[21, 0],[14, 1],[62, 1]]
test_data = sorted(test_data, key=lambda x: x[0], reverse=False)

s = 1
n = 5
pre = 0
for t, censor in test_data:
    plt.plot((pre, t), (s, s), color='blue')
    if censor == 1:
        plt.plot((t, t), (s, s * (n - 1) / n), color='blue')
        s = s * (n - 1) / n
    n -= 1
    pre = t
plt.show()

先ほどの表を2次元配列として格納し、扱いやすいようにtによって昇順でソートしました。2次元配列に対してソートを施したいときはkey=lambda x: x[0]などとキーワード引数によって何でソートするのかを明示します。もしcensorでソートしたいのであればx[1]を入れてあげれば大丈夫です。その後、for文で生存確率を計算し、それをプロットしたものが以下の図になります。

このように、Kaplan=Meier曲線は階段関数として出力されます。死亡が観測された時点で減少し、その間は一定です。$t=62$において生存者は0であるため、最後は$\hat{S}(62)=0$となります。今回の例では同一のtにおいて複数死亡することはありませんでしたが（このことをtieといいます)、tieが起こっても乱数を発生させてランダムに順序を割り振ることで同様の議論を展開できます。さらに、KM推定量においてtieデータを一律に扱ったとしても最終的な推定量は等しいため、調整を行う必要性は基本的にはありません。実務上は少ないですが、極端に多くのtieが発生している場合には、KM推定量ではなく離散時間モデルの使用を検討したほうが良いでしょう。また、最終観測時間が右側打ち切りの場合、それ以降の生存時間の推定量は定義できないことにも注意が必要です。

Kaplan-Meier推定量の一般化

　生存時間解析において、Kaplan-Meier推定量は非常によく使われるため、ここで一般的な定式化を与えておきます。tが有限の時、ソートすることは一般性を失わないため、tを昇順にソートして$t_i$時点における死亡によるセンサーかどうかの2値変数$c_i$、$t_i$時点のat risk人数$n_1$、観測される死亡数$d_i$のもとで、時点tでの生存関数のKaplan-Meier推定量は、

\hat{S}(t)=\prod_{t_i \leq t}\frac{n_i - d_i}{n_i}

ただし、$t \leq t_1$のとき$\hat{S}(t)=1$となります。

Kaplan-Meier曲線のplot -lifelines編

　先ほどは自力でKaplan-Meier曲線をプロットしましたが、わざわざ自分で実装しなくてもpythonにはlifelinesと呼ばれる生存時間解析用のライブラリが存在します。この節では、そのlifelinesを使用して前回のWHAS500のKaplan-Meier曲線を描画してみます。
lifelinesはデフォルトでは含まれていないので、一度も使用したことのない方はターミナルもしくはコマンドプロンプトで最初に以下を実行してください。

pip install lifelines

では、ここからpythonコードを示します。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lifelines as ll
from lifelines import KaplanMeierFitter


whas = pd.read_csv('whas500.csv')
kmf = KaplanMeierFitter()
kmf.fit(whas.LENFOL, event_observed=whas.FSTAT)
kmf.plot_survival_function()
print(kmf.survival_function_)
plt.show()

lifelinesのKaplanMeierFitterのfitの第一引数に時間tを、キーワード引数event_observedに観測したいイベント（死亡）が発生したかどうかの2値変数を入れて、フィッティングさせてあげます。そうすると、モデルの.plot_survival_fuction()関数によって簡単にKM曲線をプロットすることができます。また、推定された生存関数を得るにはsurvival_function_を使用します。pandasのDataFrame型で時点と推定値のペアを返してくれます。また、survival_functionをcumulative_densityに置き換えることで前述の累積分布関数$F(t)$を得ることも可能です。実際にKM曲線をプロットしたものが以下になります。

図の水色の幅は信頼区間です。

以上でKaplan-Meier推定量とその実装の解説を終えます。次回は、Kaplan-Meier推定量の解釈や2群比較を行いたいと思いますので、引き続き読んでいただけると幸いです。

参考文献・参照リンク

生存時間解析入門　Hosmer DW, Lemeshow S, May S
LIFELINES　https://lifelines.readthedocs.io/en/latest/index.html
京都大学OCW　京都大学大学院医学研究科 聴講コース 臨床研究者のための生物統計学「生存時間解析の基礎」　https://www.youtube.com/watch?v=NmZaY2tDKSA&feature=emb_title