読者です 読者をやめる 読者になる 読者になる

StatsFragments

Python, R, Rust, 統計, 機械学習とか

Python lifelines で生存分析

pandas Python

サイトのアクセス履歴をみていたら "Python", "生存分析", "Kaplan-Meier" なんかがちらほらあったので、知っている方法を書いてみる。

生存分析とは

いくつかのサンプルについて、何らかのイベントが起きるまでの時間とイベント発生率との関係をモデル化するもの。よく使われるのは

  • 何かの病気になってからの死亡率
  • 部品の故障率

など。

この説明がわかりやすい。

Rと生存時間分析(1)

Python でやる方法

どちらか。lifelines のほうが手間はすくない。

  1. lifelines パッケージを使う
  2. 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

データの意味はこちら。

R: NCCTG Lung Cancer Data

ここでは、性別でグループ分けしたデータについて、イベントの発生時間 ("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 (女性) のほうが 同じ経過時間時点での生存率が高い傾向がある。

f:id:sinhrks:20141101215357p:plain

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()

f:id:sinhrks:20141101215936p:plain