生存時間分析のメモ(~Python ライブラリ lifelines~)

$$\newcommand{\Prob}[1]{\mathbb{P} \lbrack #1 \rbrack}$$

生存時間分析の基礎事項についてまとめてみた。pythonの生存時間分析ライブラリであるLifelinesを使った分析例も載せています.

生存時間分析

生存時間分析は，統計学の一つの分野で「ある事象(event)が起きるまでの期間(duration)を扱う」ことを目的とします．生存時間分析の適用例を挙げるならば

・(医学の例)生活習慣(飲酒，喫煙)と余命との関係を調べる

・(マーケティングの例)ある商品を購入した顧客が，再びその商品を購入するまでの期間を調べる

などの場面が挙げられます．

　生存時間分析の特徴として，打ち切りデータを扱うことが出来る点があります．打ち切りデータとは，観測したい事象の発生時刻が正確には分からない(少なくともある時刻以降である"等)データのことです．例えば，余命を調べるために実験に参加している被験者が5年間実験に参加したのちに実験の辞退を行った場合などが該当します．以下，生存時間の基礎となる事項について見ていきましょう．

生存関数，ハザード関数

記号

$\tau$ : $\lbrack 0, \infty )$ に値をとる確率変数
$f(\tau)$ : $\tau$ の確率密度

このとき， $\tau$ の生存関数(survival function)とは，次で定義される．

$\displaystyle{S(t) := \Prob{\tau \geq t} = \int_{\lbrack t, \infty )} f(\tau) d\tau \cdots (1.1.1).}$

また， $\tau$ のハザード関数(hazard function)とは，次で定義される関数 $h(t)$ である，

$\displaystyle{h(t) := \frac{f(t)}{S(t)} = - \frac{d}{dt} \log{S(t)}\cdots (1.1.2){}_.}$

確率密度，生存関数，ハザード関数の定義については，つぎのまとめを理解しておけば，問題ナシでしょう．

$f(t)$

$\overset{\text{微分}}{\underset{積分}{\leftrightharpoons}}$

$S(t)$

$\overset{\text{積分してexp}}{\underset{対数微分}{\leftrightharpoons}}$

$h(t)$

例1.1.1:指数分布(Exponential distribution)

確率変数 $\tau$ のハザード関数 $h(t)$ が $h(t) = \alpha \quad (\alpha$ > $\text{0:定数})$ で与えられるとき， $\tau$ の確率密度関数 $f(\tau)$ は次で与えられる，

$f(\tau) = \alpha \exp(-\alpha \tau).$

(導出) $S(t)$ を， $\tau$ の生存関数とする．ハザード関数の定義(1.1.2)より，

$\displaystyle{- \frac{d}{dt} \log{S(t)} = h(t) = \alpha .}$

両辺積分して，

$\displaystyle{- \log{S(t)} = \int h(\tau) d\tau = \alpha t + \text{積分定数}}$

$S(0)=1$ より積分定数は，0なので結果 $S(t) = \exp(-\alpha t)$ .確率密度関数は，

$\displaystyle{f(t) = - \frac{d \ S(t)}{dt} = \alpha \exp(-\alpha t)}.$

例1.1.2:ワイブル分布(Weibull distribution)

確率変数 $\tau | \lambda , \rho$ の生存関数 $S(t | \lambda , \rho)$ が，次で与えられるような分布をワイブル分布( $\it{Weibull\ distribution}$ )という,

$\displaystyle{ S(t | \lambda , \rho) = \exp (- (\lambda t){}^\rho)}$ .

また，このとき確率密度関数 $f(\tau | \lambda , \rho)$ ，およびハザード関数 $h(\tau | \lambda , \rho)$ は次で与えられる，

$f(\tau | \lambda , \rho) = \lambda (\lambda \tau) {}^{\rho -1}\exp (- (\lambda \tau){}^\rho)$

$h(\tau | \lambda , \rho) = \lambda \rho (\lambda \tau)^{\rho -1}$ .

パラメトリックモデル

この，サブトピックでは最尤法によるパラメトリックモデルのあてはめについて説明をします．まず，尤度関数を定義します． $\tau | \theta , (\theta:\text{パラメータ})$ を，非負の値をとる確率変数とすします．また， $\tau | \theta$ の確率密度関数，および生存関数をそれぞれ $f(\tau|\theta), S(\tau|\theta)$ とする．

$N$ 個体の観測データ $\underline{t} = \{t_1, \ldots, t_N\}$ があるとする．また，観測データにたいして $\underline{\delta} = \{\delta_i\quad | i= 1\ldots, N\}$ を

$\delta_i = \begin{cases} 1\quad (t_i \text{が打ち切りでない場合}) \\ 0　\quad (t_i \text{が打ち切りの場合})\end{cases}$

と定義しておく．以下， $N$ 個体の観測データは独立と仮定する．

このとき，尤度関数(likelihood function) $\mathcal{L}(\theta)$ は次のように求まる，

$\mathcal{L}(\theta)$

$= \Prob{\underline{t} | \theta}$

$=\displaystyle{\prod_{i:t_i \text{は打ち切りデータでない}}\Prob{t_i | \theta}\ \prod_{i:t_i \text{は打ち切りデータ}} \Prob{\tau \geq t_i | \theta} }$

$=\displaystyle{\prod_{i:t_i \text{は打ち切りデータでない}}f(t_i | \theta)\ \prod_{i:t_i \text{は打ち切りデータ}} S(t_i | \theta)}$

$=\displaystyle{\prod_{i = 1, \ldots, N} f(t_i | \theta)^{\delta_i} S(t_i | \theta)^{1-\delta_i}}\cdots (1.2.1).$

例1.2.1:指数分布の場合

今までの記号を利用する． $\tau | \alpha \sim \text{Exp}(\alpha)$ (パラメータ $\alpha$ )の指数分布とする．このとき最尤推定量 $\alpha^{\text{MLE}}$ は，次で与えられる．

$\displaystyle{ \alpha {}^{\text{MLE}} = \frac{\sum_{i=1}^N t_i}{|\underline{\delta}|}}$

,ただし $|\underline{\delta}| = \sum_{i= 1}^N \delta_i$ (=打ち切りデータ以外のデータの個数)とする．

(導出)観測データ $\underline{t} = \{t_1, \ldots, t_N\}, \underline{\delta} = \{\delta_i | i= 1\ldots, N\}$ に対する，尤度は式(1.2.1)より，

$\displaystyle{\mathcal{L}(\alpha)=\prod_{i = 1, \ldots, N}\{\alpha \exp(-\alpha t_i)\}{}^{\delta_i} \{ \exp(-\alpha t_i) \}{}^{1-\delta_i}\}\} }$

$\displaystyle{ \qquad = \alpha {}^{|\underline{\delta}|} \exp \Bigl( - \alpha \sum_{i=1}^N t_i \Bigr) .}$

よって，最尤推定量 $\alpha^{\text{MLE}}$ は

$\displaystyle{0 = \frac{\partial \mathcal{L}}{\partial \alpha}(\alpha^{\text{MLE}})}$

$\displaystyle{ \Leftrightarrow \alpha {}^{\text{MLE}} = \frac{\sum_{i=1}^N t_i}{|\underline{\delta}|}}$

と求まる．

PythonパッケージLifelinesでパラメータ当てはめ

Lifelinesは，生存時間分析ようのライブラリです．以下，このパッケージを用いてワイブル分布による推定を行います．

(i)分析データ

肺がんデータ

https://github.com/CamDavidsonPilon/lifelines/blob/master/lifelines/datasets/lung.csv

・・・肺がん治療患者のデータ．10因子(性別や，年齢etc...)ごとに死亡までの時間が

保存されている．

(ii)分析の目的

肺がんデータを性別毎に分け，そのハザード関数をWeibull分布により推定する．

(iii)Nelson-Aalen推定量

Nelson-Aalen推定量とは，累積ハザード関数 $H(t) := \int_{\lbrack 0, t \rbrack} \lambda(\tau) d\tau$ の推定量の一つで次の式で定義されます．イベント発生までの時間 $\{t_1 \$ < $\ \ldots \$ < $\ t_N \}$ に対して，

$\displaystyle{\hat{H}(t) = \sum_{t_i \leq t} \frac{d_i}{n_i}}$

ただし，

$d_i := \text{時刻}t = t_i \text{におけるイベント発生数}$

$n_i := \text{時刻}t = t_i\text{までの生存数}$

です．今回は比較用にNelson-Aalen推定量もプロットを行います．

(iv)分析



#2018/12/28
from lifelines import WeibullFitter
from lifelines import datasets as ds
from lifelines import NelsonAalenFitter
import matplotlib.pyplot as plt
import numpy as np
#肺がんデータ(10因子)を読み込み#性別:男性(sex = 1), 女性(sex = 2)経過日数:time
#@Reference:https://github.com/CamDavidsonPilon/lifelines/blob/master/lifelines/datasets/lung.csv
df = ds.load_lung()
#@memo 男性
durations_male = (df[df['sex']==1])['time']
status_male = (df[df['sex']==1])['status']event_observed_male = []
#@todo statusから，イベント観測の有無をTrue Falseに変換

for status in status_male:
	if status == 1:
		event_observed_male.append(False)
	else :
		event_observed_male.append(True)

#@memo 女性
durations_female = (df[df['sex']==2])['time']
status_female = (df[df['sex']==2])['status']
event_observed_female = []
for status in status_female:
if status == 1:
		event_observed_female.append(True)
	else :
		event_observed_female.append(False)

 #weibull分布で当てはめ
 wf_male =WeibullFitter().fit(durations_male, event_observed_male, label='Weibull_filter_male')
 wf_female = WeibullFitter().fit(durations_female, event_observed_female, label='Weibull_filter_female')


#Nelson-Aalen estimator
naf_male = NelsonAalenFitter().fit(durations_male, event_observed_male, label='Nelson-Aalen_filter_male')
naf_female = NelsonAalenFitter().fit(durations_female, event_observed_female, label='Nelson-Aalen_filter_female')
#@memo param
print(wf_male.lambda_, wf_male.rho_)
print(wf_female.lambda_, wf_female.rho_)
 

#@todo ハザード関数の描画(性別毎)
ax = wf_male.plot()
naf_male.plot(ax = ax)
plt.title("Weibull_filter_versus_Nelson_Aalen_filter_male")
plt.xlim(0, 1000)
plt.ylim(0, 6.0)
ax = wf_female.plot()
naf_female.plot(ax = ax)
plt.title("Weibull_filter_versus_Nelson_Aalen_filter_female")
plt.xlim(0, 1000)
plt.ylim(0, 6.0)
plt.show()

推定結果のプロットは以下の通りです。(注:実線の周りの色付き部分は，95%Confidential Intervalです．)

f:id:kimigayoseishou:20181230164912p:plain