Python
statistics
bioinformatics
lifelines

Pythonを使ってKaplan-Meier曲線を描こう

生存時間分析とは

 臨床研究を行う上で避けて通れないのが予後の検討ですが、本稿ではその代表的な手法であるKaplan-Meier曲線とlog-rank検定をPythonライブラリを用いて行いたいと思います。
 生存時間分析は、ある集団において死亡などのイベントが生じるまでの期間を短い順に並べて、その長短の度合いを評価するものです。人間には寿命がありますから、観察期間を無限にすると、最終的な生存率は0%になります。そのため、一定の観察期間を設けて、その間の経時的な生存率の推移(生存時間関数)を比較・評価することが一般的です。
 イベントが生じる前に意図的に観察を終了することをcensored(打ち切り)と表現し、イベントが発生した個体と区別することが、この分析手法の特徴です。
 今回はlifelinesというPythonライブラリを使用して、Group_A・Group_B 2群間の生存時間関数の統計学的比較を行います。

データの準備

 生存時間分析を行う際、Group_AとGroup_Bの各群において、生存期間とcensoredの有無をデータセットとして準備します。注意点として、lifelinesではイベント発生の有無を引数として関数に渡す仕様になっているので、censoredとは真偽値が逆になります。今回は下図のように、死亡イベントがあった場合を1として入力データを準備しました。

data_set.jpg

 xlrdモジュールを使用し、life_A・life_B、ev_A・ev_Bをそれぞれ生存期間とイベント発生の有無を格納したリストとして定義しました。その後の処理がしやすいように、生存期間のリストはNumpyのndarrayに変換しています。

#!/user/bin/env python
# coding:utf-8

import xlrd
import numpy as np
from lifelines.statistics import logrank_test
from lifelines import KaplanMeierFitter
import matplotlib.pyplot as plt

book = xlrd.open_workbook("kap_data.xlsx")
sheet = book.sheet_by_index(0)
N = sheet.nrows -1

life_A, life_B = [], []
ev_A, ev_B = [], []

for v in range(1,N+1):
    if sheet.cell(v,0).value != sheet.cell(1,0).value and sheet.cell(v,0).value != "":
        B = sheet.cell(v,0).value
        life_B.append(sheet.cell(v,1).value)
        ev_B.append(bool(sheet.cell(v,2).value))
    elif sheet.cell(v,0).value != "":
        A = sheet.cell(v,0).value
        life_A.append(sheet.cell(v,1).value)
        ev_A.append(bool(sheet.cell(v,2).value))

life_A = np.array(life_A)
life_B = np.array(life_B)

 次に観察期間を自由に変更できるよう、ユーザー入力を受け付ける処理を記述します。最終的にイベントが発生した個体でも、設定した観察期間内に生存していればイベントを0に書き換える処理を行っています。

obs_dur = input("Enter the Observational Period:")
if obs_dur == "":
    obs_dur = max(max(life_A), max(life_B)) #入力無ければ最大期間を設定
else:
    obs_dur = float(obs_dur)

for v in range(len(life_A)):
    if life_A[v] > obs_dur:
        ev_A[v] = 0
for v in range(len(life_B)):
    if life_B[v] > obs_dur:
        ev_B[v] = 0

life_A = np.minimum(life_A, obs_dur * np.ones(len(life_A)))
life_B = np.minimum(life_B, obs_dur * np.ones(len(life_B)))

 ここまでで必要なデータの成形が整いました。では早速Kaplan-Meier曲線を表示してみましょう。グラフを見やすくするために信頼区間は非表示とし、打ち切り症例(censored)はグラフ上にマークするようにしています。

axis = plt.subplot(111) #おまじない
kmf = KaplanMeierFitter()
kmf.fit(life_A, event_observed=ev_A, label=str(A))
ax = kmf.plot(ax=axis, ci_show=False, show_censors=True)
kmf.fit(life_B, event_observed=ev_B, label=str(B))
ax = kmf.plot(ax=axis, ci_show=False, show_censors=True)

plt.show()
plt.ion() #インタラクティブモードで起動

 上記コードをシェルから実行すると下のようにグラフが表示されます。
figure_1.png

 うまくいきましたね。どうやらGroup_Aのほうが各観察時点における生存率は高そうです(ただし今回の作例はGroup_Aに打ち切り症例が多いだけなので、実はどこかで野垂れ死んでいる可能性はあります)。
 視覚的に差が分かったので、今度はこの2群の生存率全体について統計学的な差異があるかどうか、log-rank testを用いて解析してみましょう。

ログランク検定

 ログランク検定はカイ二乗検定の一種で、ノンパラメトリックに2群の生存時間を比較することができます。lifelinesにはこの検定を行う関数が用意されているので、生存時間とイベント発生の各配列を引数として渡すだけで結果を得ることができます。 

result = logrank_test(life_A, life_B, event_observed_A=ev_A, event_observed_B=ev_B)
resut.print_summary() #print(result)でも出力は同様

 コードをシェルから実行すると、結果は下記のようになります。

Results
   t 0: -1
   test: logrank
   alpha: 0.95
   null distribution: chi squared
   df: 1

   __ p-value ___|__ test statistic __|____ test result ____|__ is significant __
         0.03822 |              4.295 |      Reject Null    |        True

 フムなるほど簡単でした。p値 = 0.038ということで、2群間には有意な生存時間の差があるようです。ちなみに3群以上の比較をしたい場合には pairwise_logrank_test() という関数も用意されているので、同様の手順で解析が可能です(Bonferroni補正により有意水準を低く抑える方式のようです)。

終わりに

 自分の知識の整理を兼ねて、生存時間分析の過程を書いてみました。
 SPSSやJMPのような高価なソフトを使わなくても、私たちにはPythonがあります。
 私はプログラミングも統計学も素人ですが、この記事が少しでも研究者の皆様の懐を温める一助になれば幸甚です。