Python lifelines で生存分析
サイトのアクセス履歴をみていたら "Python", "生存分析", "Kaplan-Meier" なんかがちらほらあったので、知っている方法を書いてみる。
生存分析とは
いくつかのサンプルについて、何らかのイベントが起きるまでの時間とイベント発生率との関係をモデル化するもの。よく使われるのは
- 何かの病気になってからの死亡率
- 部品の故障率
など。
この説明がわかりやすい。
Python でやる方法
どちらか。lifelines のほうが手間はすくない。
- lifelines パッケージを使う
- rpy2 でデータを R に渡して survival パッケージ を使う
1. lifelines パッケージを使う
これを使えば Pure-Python で生存分析できる。使えるのは、以下のノンパラ、セミパラ系の手法。Weibull などのパラメトリック系はない (Scipy 使えばできる)。
- Kaplan-Meier
- Nelson-Aalen
- Cox 比例ハザード
インストール
pip install lifelines
ドキュメント
バージョン古いようなのでアップデートされてないかも。
Lifelines — lifelines 0.8.0.1 documentation
やり方
データは pandas の DataFrame で扱うのがよい。ここでは、lifelinesにもサンプルデータとして入っている lung (肺ガン) データを使う。
lungデータ
https://github.com/CamDavidsonPilon/lifelines/blob/master/datasets/lung.csv
データの意味はこちら。
ここでは、性別でグループ分けしたデータについて、イベントの発生時間 ("time" カラムの値) とそのときの状態 ("death" カラムの値, Trueなら死亡, Falseなら censored = 観測打ち切り) で Kaplan-Meierモデルを作る。
これは lifelines.KaplanMeierFitter クラスでできる。引数の詳細はドキュメント参照。
Introduction to using lifelines
import pandas as pd # 表示する行数を指定 pd.options.display.max_rows = 5 # データの読み込み df = pd.read_csv("lung.csv") # 元データには、死亡なら2, censored なら1が入っている # lifelines に渡すため、死亡なら True, censored なら False となるよう変換する df['death'] = df['status'] == 2 df # Unnamed: 0 inst time status age sex ph.ecog ph.karno pat.karno \ # 0 1 3 306 2 74 1 1 90 100 # 1 2 3 455 2 68 1 0 90 90 # .. ... ... ... ... ... ... ... ... ... # 226 227 6 174 1 66 1 1 90 100 # 227 228 22 177 1 58 2 1 80 90 # # meal.cal wt.loss death # 0 1175 NaN True # 1 1225 15 True # .. ... ... ... # 226 1075 1 False # 227 1060 0 False # # [228 rows x 12 columns] from lifelines import KaplanMeierFitter # グラフを描画する Axes ax = None # 性別でグループ分け for name, group in df.groupby('sex'): kmf = KaplanMeierFitter() kmf.fit(group['time'], event_observed=group['death'], label = 'sex=' + str(name)) # 描画する Axes を指定。None を渡すとエラーになるので場合分け if ax is None: ax = kmf.plot() else: ax = kmf.plot(ax=ax) # Kaplan-Meier曲線を描画 import matplotlib.pyplot as plt plt.show()
Kaplan-Meier曲線 (横軸が経過時間、縦軸が生存率 = イベントの未発生率) となる。性別 = 2 (女性) のほうが 同じ経過時間時点での生存率が高い傾向がある。
2. rpy2 でデータを R に渡して survival パッケージ を使う
pandas と rpy2 との連携はこちらを参照。
やり方
やりたいことにもよるが、結果を Python に戻してからの処理がちょっとめんどい。
import numpy as np import pandas as pd # 表示する行数を指定 pd.options.display.max_rows = 5 # rpy2 使って survival::lung データセットを呼び出す # 手元のデータを使うなら不要 import rpy2.robjects as robjects import pandas.rpy.common as com robjects.packages.importr('survival') lung = com.load_data('lung') # rpy2 形式に変換し、 Rへ渡す lungr = com.convert_to_r_dataframe(lung) robjects.r.assign('lungr', lungr); # R で生存分析実施 robjects.r('fitted <- survfit(Surv(time, status) ~ sex, data = lungr)'); # 結果を Python に戻し、各プロパティを numpy.arrayとして取得 fitted = robjects.r['fitted'] time = np.array(fitted.rx2('time')) surv = np.array(fitted.rx2('surv')) lower = np.array(fitted.rx2('lower')) upper = np.array(fitted.rx2('upper')) # strata のキーを順序どおりに取得するには一手間必要 robjects.r('strata_keys <- names(fitted$strata)') strata_keys = com.convert_robj(robjects.r['strata_keys']) strata = com.convert_robj(fitted.rx2('strata')) strata = reduce(lambda x, y: x + y, [ [k] * strata[k] for k in strata_keys]) # DataFrame作成 result = pd.DataFrame({'time': time, 'surv': surv, 'lower': lower, 'upper': upper, 'strata': strata}) result # lower strata surv time upper # 0 0.954230 sex=1 0.978261 11 1.000000 # 1 0.943423 sex=1 0.971014 12 0.999412 # .. ... ... ... ... ... # 204 0.001582 sex=2 0.011111 821 0.078022 # 205 NaN sex=2 0.000000 965 NaN # プロット import matplotlib.pyplot as plt # グラフを描画する Axes ax = None # 生存曲線と 信頼区間を別々に描画する必要がある # strata ごとに同じ色を使うために グループ名 : 色の辞書を作成 colors = {'sex=1': 'green', 'sex=2': 'blue'} for name, group in result.groupby('strata'): # pandas は index を x軸として plot するので、 # time の値を indexに設定 group = group.set_index('time') color = colors[name] # pandas.Series.plot で生存曲線自体を描画 ax = group['surv'].plot(ax=ax, color=colors[name]) # matplotlib.Axes.fill_between で信頼区間を描画 ax.fill_between(group.index, group['lower'], group['upper'], color=colors[name], alpha=0.3) plt.show()
