カプランマイヤー曲線について
はじめに¶
書こうと思ったきっかけ¶
実習中論文を頂いて読む機会が多いが、医学論文は薬は新しい治療方法の成績を比較しているのが多い。 だいたい Kaplan-Meier曲線 が使われているので(偏見)、論文を読むときにはもちろん、将来自分が論文を執筆するときに、何が行われているのか少しでも理解できていればいいかなと思いまとめてみた。
Kaplan-Meier 曲線とは¶
Wikipediaによると「治療を行った後の患者の生存率をグラフにしたもの」ということらしい。ただ、実際には問題がいくつかあって
- 観察期間は限られている(患者さんを何年もフォローすることはできない)
- 観察不能になることがある(病院移った、他の要因で死亡したなど)
- 患者さんによって観察開始時間は様々(みんな一緒に罹患したり治療を開始するわけではない) そういった問題をうまく対処してくれるのがKaplan-Meier曲線。
Kaplan-Meier曲線の縦軸は累積生存率で、
$$ S(t)=\prod_j^t \left(1-\frac{d_j}{n_j}\right) $$で表されます。ここで$d_j$は時間$j$における死亡者数で、$n_j$は時間$j$における生存者数です。
今回はpythonのコードでKaplan-Meier曲線を書いていこうと思います。
pythonのライブラリにlifelinesという便利なライブラリがあるのでそれを使っていこうと思います。installはpipでできます：pip install lifelines
データの準備¶
ある病気に対する治療法の予後のKaplan-Meier曲線を描きたいとします。 ある期間内に観察された患者の経過をみることになると思いますが、患者さんによって発症する時間も様々なので、イメージとして下の図のようになると思います。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from lifelines import KaplanMeierFitter
from lifelines.plotting import plot_lifetimes
from numpy.random import uniform, exponential
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
N = 50
current_time = 48
birth = np.random.randint(0,45,50)
actual_lifetimes = np.array([[int(exponential(20))+1, int(exponential(10))+1][uniform() < 0.5] for i in range(N)])
observed_lifetimes = np.minimum((actual_lifetimes+birth), current_time)
observed = (birth + actual_lifetimes) < current_time
plt.xlim(0, 50)
plt.vlines(48, 0, 50, lw=2, linestyles='--',colors='red')
plt.xlabel("time")
plt.title("Births and deaths of our population, at $t=48$")
plot_lifetimes(observed_lifetimes-birth,
event_observed=observed,
birthtimes=birth,
sort_by_duration=False)
#plot_lifetimes(observed_lifetimes-birth, event_observed=observed,birthtimes=birth)
print("Observed lifetimes at time %d:\n" % (current_time), observed_lifetimes)
赤線では期間内(t=48)に患者さんが死亡したケース、青線は観測不能になったケースです(観察期間終わったからそれ以上追えない)。
簡単のため、観測不能になったケースは観測期間が終わった場合のみとします。観察期間の長さでソートすると、下の図のようになります（defaultはソートされて出力されます）。
s = pd.Series(observed_lifetimes-birth,observed)
s = s.sort_values()
s.name = 'time'
observed, result_lifetimes = s.index, s.values
plt.xlim(0,50)
#plt.vlines(22, 0, 30, lw=2, linestyles='--',colors='red')
plt.xlabel("time")
plt.title("Births and deaths of our population, at $t=22$")
plot_lifetimes(result_lifetimes, event_observed=observed)
実際の計算¶
さっき作ったデータを使って、累積生存率を計算していきます。たとえばt=5における累積生存率を計算するとき、
pd.DataFrame(s).head(12)
rueが患者さんが死亡したのを観測できた事例、Falseが観察期間終了により観測できなかった事例です。観察期間の短い順に見ていくと、
- t=1のとき、50人のうち2人が死亡,1人が打ち切り。
- t=2のとき、死亡者、打ち切りなし。
- t=3の時は４７人のうち１人が死亡。
- t=4の時は46人のうち２人が死亡、２人が打ち切り。
- t=5の時は42人中１人死亡、３人が打ち切り。
よって$ S(5)$は
$$ S(5)=\left(1-\frac{1}{50}\right)\left(1-\frac{0}{47}\right)\left(1-\frac{1}{47}\right)\left(1-\frac{2}{46}\right) \left(1-\frac{1}{42}\right) \simeq 0.896 $$となります。
kmf = KaplanMeierFitter()
kmf.fit(result_lifetimes,event_observed=observed)
kmf.plot()
plt.title('Survival function of political regimes');
hormone therapyの有無で2群に分けKaplan-Meier曲線を描いた例¶
from lifelines.datasets import load_gbsg2
df = load_gbsg2()
ax = None
for name, group in df.groupby('horTh'):
kmf = KaplanMeierFitter()
kmf.fit(group['time'], event_observed=group['cens'],
label = 'hormone therapy =' + str(name))
# 描画する Axes を指定。None を渡すとエラーになるので場合分け
if ax is None:
ax = kmf.plot()
else:
ax = kmf.plot(ax=ax)
plt.title('Kaplan-Meier Curve')
plt.show()
Cox比例ハザードモデル(CoxPH model)の適応¶
from lifelines import CoxPHFitter
cph = CoxPHFitter()
df = pd.get_dummies(df,columns=['horTh','menostat'],drop_first=True)
df.drop('tgrade',axis = 1,inplace = True)
cph.fit(df, duration_col='time', event_col='cens', show_progress=True)
CoxPH回帰の結果のsummary
cph.print_summary()
CoxPH回帰の結果のplot
cph.plot()
- 前の記事 : pandasのpivot_tableを用いた高速データ処理
- 次の記事 : 2018年阪大理系数学をPythonで解いてみた
- 関連記事 :
