Survival Analysis with Python

備忘録としてまとめます。

生存分析は、生物の死や機械システムの故障など、あるイベントが発生するまでの期間を分析するための統計学的な手法の1つです。
ある時期を過ぎて生存する人の割合はどのくらい？あるイベントが起きた人のうち、その人たちはどのくらいの割合で天寿をまっとうする？ある特性は、生存の確率をどのように増加または減少させる？このような疑問に答える手法を模索するための学問ですね。

#Kaplan-Meier法を用いた予測

ある時期に発生するイベントの割合から発生確率を予測します。
データはこちらを使います。
「Tongue」
https://drive.google.com/file/d/1fpiIaKFMc9QPbOamKdxO1sWROzf1DP1Z/view?usp=sharing
元データ

sksurvを使う

#sksurvを使う方法
!pip install scikit-survival
#してから
#CSVファイルをこのノートへロードします。
#ファイルを開くダイアログを表示
from google.colab import files
file = files.upload()

#ファイルをCSVとして読み込みます
import pandas as pd
import io
data = pd.read_csv(io.StringIO(file['tongue.csv'].decode('utf-8')),header=0)
#中身の大きさと先頭5行のデータの内容を確認してみます。
print('データの形式：{}'.format(data.shape))
print(data.head())

#イベントのデータ型がintだと駄目なようなので、真偽値に変換しておきます。
#一先ず、グループ１のみを計算対象にします。
f = data.type==1
T = data[f]['time']
event = data[f]['delta']
eventBool = []
for i in event:
  if(event[i]==1):
    eventBool.append(True)
  else :
    eventBool.append(False)

#予測します。
from sksurv.nonparametric import kaplan_meier_estimator
x, y = kaplan_meier_estimator(eventBool, time)

#グラフにします
import matplotlib.pyplot as plt
plt.step(x, y, where="post")
plt.ylim(0, 1)
plt.show()

lifelinesを使う

!pip install lifelines # このライブラリも便利です。https://lifelines.readthedocs.io/en/latest/Quickstart.html

#からの予測器を初期化して
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()

#真偽値変換無しでそのまま予測器の引数に渡します
f = data.type==1
T = data[f]['time']
C = data[f]['delta']

kmf.fit(T, event_observed=C)

#グラフにします
kmf.plot(title='Tumor DNA Profile 1')

比較

sksurvでも出来ますが、lifelinesライブラリを使って、グループ2（つまりtype:2）と比較をしてみます。

f2 = data.type==2
T2 = data[f2]['time']
C2 = data[f2]['delta']

ax = plt.subplot(111)

kmf.fit(T, event_observed=C, label=['Type 1 DNA'])
kmf.survival_function_.plot(ax=ax)
kmf.fit(T2, event_observed=C2, label=['Type 2 DNA'])
kmf.survival_function_.plot(ax=ax)

plt.title('Lifespans of different tumor DNA profile')

logrank-test

生存時間解析において、2群の生存時間に差があるかを検定するノンパラメトリックな検定手法のこと。
一般化Wilcoxon検定に対して相対的に後期に起きた死亡を重く評価するという特徴を持つ。時間経過と共に生存率曲線の差が開いてくるような場合、ログランク検定は一般化Wilcoxon検定に比べて検出力が高くなる。（https://bellcurve.jp/statistics/glossary/719.html）

from lifelines.statistics import logrank_test
summary_= logrank_test(T, T2, C, C2, alpha=99)
summary_.print_summary()

出力

               t_0 = -1
 null_distribution = chi squared
degrees_of_freedom = 1
             alpha = 99
 test_statistic    p  -log2(p)
           2.79 0.09      3.40

#ハザード比の算出(Nelson-Aalen kernelを利用)

ハザード比も算出できます。
ある臨床試験で検討したい新薬Aと比較対象の薬剤Bとを比べたとき、ハザード比が1であれば2つの治療法に差はなく、ハザード比が1より小さい場合には治療Aの方が有効と判定され、その数値が小さいほど有効であるとされます。　例えばA薬と対象のB薬を比較するというある臨床試験でハザード比が0.94という結果であれば、A薬はB薬よりリスクを6%減少させたという意味になります。(https://oncolo.jp/dictionary/hr)

from lifelines import NelsonAalenFitter

naf = NelsonAalenFitter()
naf.fit(T, event_observed=C)

naf.plot(title='Nelson-Aalen Estimate')

#最後に「Time Dependent AUC」

よくバイナリの予測器を作成する際にROC曲線を使いますが、このときのROCはあくまでもその時点で収集できたデータで作成した予測器の精度を示すROCです。
このROCを一定の期間ごとに計算し、それぞれの時点でデータを整えて分類モデルでROC（その時点の生存の正解を利用）を描き、AUCを計算し、経時的なAUCの変化を可視化したのがTime dependent AUCです。

これはsksurvライブラリを使って試してみます。
免疫グロブリン遊離L鎖κ/λ比のサンプルデータをロードします。
これは、以下のようなデータです。
7874サンプルで9つの項目があります。
age：年数
sex：F =女性、M =男性
sample.yr：血液サンプルが取得された暦年
kappa：FLCカッパ
lambda：FLCラムダ
flc.grp：被験者のFLCグループ
creatinine：血清クレアチニン
mgus：被験者がmonoclonal gammapothy（MGUS）と診断されたかどうか
chapter：国際疾病コードICD-9の章の見出しによる主な死因のグループ化

#データをロードします
from sksurv.datasets import load_flchain
x, y = load_flchain()

#今回はテスト用にデータを分けることにします
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=0)

#'age', 'creatinine', 'kappa', 'lambda'のみ学習に利用します
num_columns = ['age', 'creatinine', 'kappa', 'lambda']

#欠測を簡易的に単一補完します
from sklearn.impute import SimpleImputer
imputer = SimpleImputer().fit(x_train.loc[:, num_columns])
x_train = imputer.transform(x_train.loc[:, num_columns])
x_test = imputer.transform(x_test.loc[:, num_columns])

#正解ラベルからトレーニングとテストのそれぞれのデータセットの生存期間の最大最小を取得します。
y_events = y_train[y_train['death']]
train_min, train_max = y_events["futime"].min(), y_events["futime"].max()

y_events = y_test[y_test['death']]
test_min, test_max = y_events["futime"].min(), y_events["futime"].max()

#データの整合性を確認します。テストデータに用いるデータの生存期間はトレーニングに用いるデータの生存期間内になければなりません。
assert train_min <= test_min < test_max < train_max, \
    "time range or test data is not within time range of training data."

#今回は生存期間が、5から81までの間で15間隔で計算します。
import numpy as np
times = np.percentile(y["futime"], np.linspace(5, 81, 15))
print(times)

#time dependent aucを描きます。
from sksurv.metrics import cumulative_dynamic_auc
from sksurv.metrics import concordance_index_ipcw
import matplotlib.pyplot as plt

def plot_cumulative_dynamic_auc(risk_score, label, color=None):
    auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, risk_score, times)

    plt.plot(times, auc, marker="o", color=color, label=label)
    plt.xlabel("days from enrollment")
    plt.ylabel("time-dependent AUC")
    plt.axhline(mean_auc, color=color, linestyle="--")
    plt.legend()
    
for i, col in enumerate(num_columns):
    plot_cumulative_dynamic_auc(x_test[:, i], col, color="C{}".format(i))
    ret = concordance_index_ipcw(y_train, y_test, x_test[:, i], tau=times[-1])

上記のグラフから、予測力高い特徴量は「Age, kappa/lambda, creatinine」の順であることがわかります。

以上です。

p.s

さらに、sksurvの開発者が深層学習を利用した生存分析のアプローチも紹介しています。
勉強しなくては。。
SurvivalAnalysisForDeepLearning

#Reference

• https://en.wikipedia.org/wiki/Survival_analysis
• https://plot.ly/python/v3/ipython-notebooks/survival-analysis-r-vs-python/
• https://k-d-w.org/slides/pyconuk-2017/
• http://sinhrks.hatenablog.com/entry/2014/11/01/223504