はじめに
生存時間解析について、ハザード関数の導出からCox比例ハザードモデルまでを掻い摘んで説明します。
加えて、PythonライブラリであるLifelinesを使って実際のデータで挙動を見ていきます。
細かい理論的な部分にはあまり興味がない方が多いかもしれませんが、モデルを使うにあたって注意するべき部分がいくつかあるので、その理解のためにも必要と思われる部分はまとめて説明しています。
医学統計の専門家ではないので見当違いな説明や不足部分があるかもしれません。その際は補足いただければと思います。
ハザード関数
生存時間$T$の分布関数を$F(t)$、密度関数を$f(t)$とし、生存時間関数を$S(t)$とすると、これらは下記のように表せます。
\begin{align}
S(t) &= Pr\{ T \geq t \}\\
F(t) &= 1 - S(t) \\
f(t) &= \frac{dF(t)}{dt}\\
\end{align}
次に瞬間死亡率であるハザード関数を以下のように定義すると、$S(t)$を用いて次のように表せます。
\begin{align}
\lambda(t) &= \lim_{\Delta t \rightarrow 0} \frac{Pr \{ t \leq T \leq t + \Delta t \mid T \geq t \}}{\Delta t}
\\
&= \lim_{\Delta t \rightarrow 0} \frac{Pr \{ t \leq T \leq t + \Delta t \ \ and \ \ T \geq t \}}
{P\{ T \geq t\}\Delta t}
\\
&= \frac{1}{S(t)} \lim_{\Delta t \rightarrow 0} \frac{F(t+\Delta t) - F(t)}
{\Delta t}
\\
&= - \frac{S(t)^{'}}{S(t)}
\end{align}
この微分方程式を初期条件$S(t=0) = 1$で解けば、ハザード関数と生存関数の関係は
\begin{align}
S(t) = \exp \left( -\int_{0}^{t} \lambda(u) du \right)
\end{align}
として表す事ができます。
- 重要ではないかもしれない豆知識（統計検定）
この導出は2015年の統計検定1級で出題されていたりします。
ここでハザード$\lambda(t)$関数を時間依存しない定数$\lambda(t) = \lambda$として与えると下記の指数分布となります。
$$
\begin{align}
S(t) &= e^{-\lambda t} , \
f(t) = \lambda e^{-\lambda t}
\end{align}
$$
また、ハザード$\lambda(t)$関数を時間依存するように$\lambda(t) = \lambda p (\lambda t)^{p-1}$として与えると、下記のワイブル分布となります。
\begin{align}
S(t) = \exp \{ -(\lambda t)^p \} ,\ \
f(t) = \lambda p (\lambda t)^{p-1} \exp \{ -(\lambda t)^p \}
\end{align}
イメージとしては指数分布 $\subset$ ワイブル分布 $\subset$ 生存時間関数です。
指数分布やワイブル分布のように明示的にハザード関数を与える生存関数をパラメトリックモデルといいます。
ノンパラメトリック生存解析（Kaplan-Meier法）
- 重要ではないかもしれない豆知識（名前の由来と変化）
Kaplan-Meier推定量はKaplanとMeierの名前に由来しますが、彼らはこの推定量をproduct limit estimatorと呼んでいたそうです(1958)。
しかし現在ではKaplan-Meier推定量のほうが通りがよいのでそう呼ばれています。
この60年間に何があってそうなってしまったのでしょうか…
明示的なハザード関数を与えずにデータから生存関数を推定するノンパラメトリック方法としてKaplan-Meier法があります。
KM法はただ単に分布の仮定が必要でなくなったことが特徴ではありません。それを説明する為にセンサーデータについて述べます。
- センサーデータ
研究終了までに生存、途中で意図的に観察を打ち切られた、意図せず追跡不能となったデータのことです。
生存時間の解析には各個体のエンドポイントが必要ですが、一般にすべてのエンドポイントを観測することは難しく、打ち切りデータが発生してしまいます。これらセンサーデータが含まれる場合はパラメトリックモデルで推定することが難しくなります。
Kaplan法はそれに適応できるという点がもう一つの特徴です。
- KM法による生存時間分布推定
$0 < t_1 < t_2 <$…$ < t_j <$… とすると、$t_j < t \leq t_{j+1}$となる$t$では、生存時間分布を下記のように算出できます。
\begin{align}
\hat S(t) &= \hat S(t_j) \left( 1 - \frac{d_j}{n_j} \right) \\
SE(t) &= \hat S(t_j) \left\{ \sum_i \frac{d_i}{n_i(n_i - d_i)}\right\}^{1/2}
\end{align}
（二つ目の式はGreenwoodの公式による推定量の標準誤差です）
ただし、$\sum$は$t_j$以前に発生した死亡$t_i \leq t_j$についての和、$d_j$、$n_j$は$t_j$での死亡発生数と観察対象数を示します（観察対象数はセンサーデータを含まない）。また、$\hat S(t)$は漸近的に$N(S(t), \ \ SE(t)^2)$に従うので、95%信頼区間を求めることができます。
このようにKM法では、データ上で死亡が発生した時点に対して各$S(t_i)$を計算することで生存関数を推定しています。
PythonによるKaplan-Meier法の実行（Lifelines）
生存時間解析にはPythonパッケージであるLifelinesがあります。
サンプルデータを用いてLifelinesを実行します。
- サンプルデータ
サンプルデータはSurvival Times after Mastectomy of Breast Cancer Patientsを用います。
乳癌患者の乳房切除後の生存期間です。
data = pd.read_csv("mastectomy.csv", index_col=0)
data["metastized"] = [1 if k == "yes" else 0 for k in data["metastized"].values]
data.head()
eventは観測されたか（センサーでないか）どうかであり、metastized(正式はmetastasized)は転移の有無です。
Stringだったのでbool値に変換しています。
for k, col in enumerate(data.columns[1:]):
plt.subplot(220 + k + 1)
plt.hist(data[col])
plt.ylabel(col)
ヒストグラムをみると転移している患者が多い事がわかります。
plt.hist(data[data.metastized == 1].time, normed=True, alpha=0.7, label="metastized")
plt.hist(data[data.metastized == 0].time, normed=True, alpha=0.7, label="normal")
plt.legend()
また転移の有無による生存時間は差があるようです。
当然ですが転移が発生している患者の方が生存時間が短い傾向があります。
- KMの実行
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(data["time"], data["event"], timeline=np.linspace(0,200,200))
kmf.plot()
第一引数に生存時間、第二引数にセンサー有無を渡せば推定してくれます（青い区間は信頼区間です）。
これでセンサーデータを含むデータの生存時間解析ができました。めでたしめでたし！
とはならないですね。
- 問題 共変量はどこへいったのか？
KM法では共変量をモデルへ組み込むことができません。
それでは実用性がないじゃないか、という事では決してなく、後に共変量を含めたモデルを構成する際に活躍します。
また、共変量であるかどうかを検定する方法はあります。これを用いればある変数が生存時間に影響を与えているのかどうかを検定できるので、モデルそのものに共変量を組み込まなくても意義のある分析は可能です。
ログランク検定
- 重要ではないかもしれない豆知識（名前の由来）
$j$番目に死ぬ個体はrank $j$の個体といいます。そしてその影響度は$\sum_{k=0}^j \frac{1}{n-k}$となります。
この値は総観測個体$n$が大きくなれば$\log (n)$に近くなります。
rank $j$の個体に$\log$に近いscoreを与える ⇒ ログランク検定 です。
Mantelさんが1966に用いました。
二群の生存時間分布の差が有意かどうかを検定する手法です。
生存時間解析に用いるデータにはセンサーデータが含まれているため、t検定や分散分析を用いることができません。そのため、順位のみ用いるMantel-Haenxel検定法の考え方に基き、死亡時間の順位のみ比較することで検定します。イメージとしてはマン・ホイットニー検定のようなものです。
- ログランク検定の実行
from lifelines.statistics import logrank_test
data_A = data[data["metastized"] == 1]
data_B = data[data["metastized"] == 0]
results = logrank_test(data_A["time"], data_B["time"], data_A["event"], data_B["event"])
results.print_summary()
有意な差があると見ていいです。
このように、KM法では共変量を考慮しないモデルを作成するか、ログランク検定を用いて共変量が存在するかどうか検定することまではできます。
- 注意点
ログランク検定に関わらず、検定をする二群に対して別の共変量がある場合、それが対象群に均等に配分されていることを確認する必要があります。専門的に言えば、ランダム化比較試験（RCT）であるかどうかという事です。
そうでない場合は傾向スコア等を用いて層別化するという手法がとられることがあります。
これついては紙面を改めて記事にできればと思います。
では共変量をモデルに組み込むにはどうすれば良いのかという事ですが、それについては後述するCox比例ハザードモデルが存在します。
Cox比例ハザードモデル
- 重要ではないかもしれない豆知識（考案者）
Coxさんが1972年に考案しました。
CoxさんといえばBoc-Cox変換が有名でしょうか。また、ロジスティック回帰の考案者でもあります。
彼はこのモデル推定の正当性について厳密な数学的証明を与えませんでした。それは1980年に確率過程を用いて正当化されることになります。
（デルタ関数が超関数を用いて数学的に正当化されたのと似ていますね。）
Cox比例ハザードモデルを説明する前に比例ハザード性について説明します。
- 比例ハザード性
すべての可能な$t>0$に対し、二つのハザード関数$\lambda_1(t)$、$\lambda_2(t)$の間に次の関係式が成立する場合、比例ハザード性が成立するという。（ただし$c$は時間$t$に依存しない定数）
$$
\begin{align}
\lambda_1(t) = c\lambda_2(t)
\end{align}
$$
一言でいえば、ある危険度（瞬間死亡率）が別のある危険度の定数倍で表せることです。
これはかなり強い仮定ですが、モデルに組み込む場合はさらに対数線形性を仮定します。
- Cox比例ハザードモデル
\begin{align}
\lambda(t \mid \boldsymbol{z}) &= \lambda_0(t) r(\boldsymbol{z}) \\
\log r(\boldsymbol{z}) &= \sum_i^k \beta_i z_i
\end{align}
ここで、$\lambda_0(t)$はベースラインハザード、$r(\boldsymbol{z})$は相対危険度関数と呼ばれるものです。
一つ目の式が比例ハザード性であり、二つ目の式が対数線形性です。
（対数線形性はロジスティック回帰で見覚えがある方も多いのではないでしょうか。）
なぜこのような仮定が必要なのかというと、それはモデル推定に必要な計算を可能にするためです。
モデル推定方法
Cox比例ハザードモデルは部分尤度による計算とKM法を組み合わせてモデルを推定します。
- 部分尤度最大化
観察された死亡数を$D$、観察された死亡時間を$t_i \ \ (1 < i < D)$、$t_i$の直前まで観察されていた集合を$R_i$とすると、部分尤度$L$は次のように表せます。
\begin{align}
L &= \prod_i^D L_i \\
L_i &= \frac{\lambda (t_i \mid \boldsymbol{z}_{(i)})}{\sum_{j \in R_i} \lambda(t_i \mid \boldsymbol{z}_{(j)})}
\end{align}
イメージとしては各死亡時刻において確率を計算している感じです。KM法の推定方法と似ています。
ここで比例ハザード性と対数線形性により、
\begin{align}
L_i &= \frac{\lambda (t_i \mid \boldsymbol{z}_{(i)})}{\sum_{j \in R_i} \lambda(t_i \mid \boldsymbol{z}_{(j)})}
= \frac{r(\boldsymbol{z}_{(i)}) }{\sum_{j \in R_i} r( \boldsymbol{z}_{(j)})} \\
&= \frac{\exp (\boldsymbol{\beta}^T\boldsymbol{z}_{(i)}) }{\sum_{j \in R_i} \exp (\boldsymbol{\beta}^T\boldsymbol{z}_{(j)})}
\end{align}
となり、$L_i$は$t$や$\lambda_0$と無関係になります。
それにより各死亡時刻における共変量の値で$L_i$を表すことができます。
そして、この積を取った次の式をCoxの部分尤度といいます。
\begin{align}
L(\beta) = \prod_i \left( \frac{\exp (\boldsymbol{\beta}^T\boldsymbol{z}_{(i)}) }{\sum_{j \in R_i} \exp (\boldsymbol{\beta}^T\boldsymbol{z}_{(j)})} \right), \ \ i \in D.
\end{align}
あとはこれを最大化するという最尤推定をすれば$\boldsymbol{\beta}$を推定することができます。
最尤推定自体は一般的な方法と同様なので詳細は省きます。
- $S_0(t)$の推定
回帰係数の推定値$\boldsymbol{\hat \beta}$を真値として$S_0(t)$を求めます。
$ w_j = \exp (\boldsymbol{\beta}^T\boldsymbol{z}_{(j)})$、$W = \sum w$とすると、
\begin{align}
q_i &= \left( 1- \frac{w_{(i)}}{W_i} \right) ^{\frac{1}{w_{(i)}}}, \ \ for \ \ i=1, \dots ,D, \\
\hat S_0(t) &= \hat S_0(t_i) q_i , \ \ for \ \ t_i < t \leq t_{i+1}.
\end{align}
この$\hat S_0(t)$はKM法の$\hat S(t) = \hat S(t_j) \left( 1 - \frac{d_j}{n_j} \right)$に対応しています。
- タイがある場合
同時点で複数のイベントが起きることをタイといいます。
多くの統計ソフトでは部分尤度はBreslowによる尤度を用いて計算していますが、これはタイが全体の死亡数に比べて少ない場合には計算が簡単でよい近似であるからです。
タイが多い場合ズレが生じてしまいます。そのようなデータに対してこの近似法を用いるのは注意が必要です。
PythonによるCox比例ハザードモデルの実行（Lifelines）
cph = CoxPHFitter()
cph.fit(data, duration_col="time", event_col="event")
cph.print_summary()
cph.plot_covariate_groups('metastized', [0, 1])
係数をみると影響度は高いと言えます。ベースラインハザードと比較すれば、転移による影響がどのように表れてくるのか把握することができます。
このConcordanceというのはC値と呼ばれるもので適合度を示しています。
詳細は省きますが、順序のみを考慮した適合度であるので一般的な精度のような数値と比較すると高いと感じると思います。
さて、共変量も考慮してモデルを推定して生存関数を書くことができました。適合度もランダムより精度が高いです。
くぅ疲、これにて完了、と言いたいところですが、まだ検証するべきことが残っています。
比例ハザード性の検証
差異を引き起こす共変量をモデルへ組み込み、生存関数の推定を行った所までは良いのですが、そもそも比例ハザード性が成り立っているかどうか、というチェックをしていませんでした。
このチェックをしていないと、間違ったモデルへデータを適合させてしまう恐れがあります。
- 層別log-logプロットの平行性
$S(t) = \exp \left( -\int_{0}^{t} \lambda(u) du \right)$と、$\lambda(t \mid \boldsymbol{z}) = \lambda_0(t) r(\boldsymbol{z})$を用いれば、
\begin{align}
S(t \mid \boldsymbol{z}) = S_0(t)^{r(\boldsymbol{z})}
\end{align}
となるので、両辺の対数を取り、$-$をかけてさらに対数を取ると、
\begin{align}
\log \{ - \log S(t \mid \boldsymbol{z}) \} = \log \{ - \log S_0(t) \} + \log r(\boldsymbol{z})
\end{align}
となります。ここで対数線形性により、さらに、
\begin{align}
\log \{ - \log S(t \mid \boldsymbol{z}) \} = \log \{ - \log S_0(t) \} + \boldsymbol{\beta}^T \boldsymbol{z}
\end{align}
となります。
この式より、共変量に対して層別にKMモデルを当てはめてlog-logプロットをすれば、各モデルが平行になることがわかります。
つまり、共変量については層別log-logプロットの平行性を確認できれば、それは比例ハザード性を満たしていると言えます。
- 層別log-logプロットの実行
def verify_proportional_hazard(pre_data, column, col_type="bool"):
if col_type == "bool":
pre_data_0 = pre_data[pre_data[column] == 0]
pre_data_1 = pre_data[pre_data[column] == 1]
elif col_type == "numerical":
pre_data_0 = pre_data[pre_data[column] <= pre_data[column].median()]
pre_data_1 = pre_data[pre_data[column] > pre_data[column].median()]
else:
print("col_type Error.")
kmf0 = KaplanMeierFitter()
kmf0.fit(pre_data_0['time'], event_observed=pre_data_0['event'])
kmf1 = KaplanMeierFitter()
kmf1.fit(pre_data_1['time'], event_observed=pre_data_1['event'])
fig, axes = plt.subplots()
kmf0.plot_loglogs(ax=axes)
kmf1.plot_loglogs(ax=axes)
axes.legend(['False', 'True'])
plt.title("Proportional Hazard of " + column)
plt.show()
verify_proportional_hazard(data, "metastized")
グラフから見れば対数線形性は成立していると判断できます。
ちなみに対数線形性が成立しない場合は、層別プロットが重なってしまったり、傾きがずれていたりします。
NG例1
NG例2
※厳密に行うのであればTime関数を利用した方法等がありますが、それらについては今回は触れないことにします。
以上によりKM法に対する共変量の組み込みとその仮定の検証までを完了しました。
先ほど推定したCox比例ハザードモデルは比例ハザード性も満たしているので妥当なモデル化と言っていいと思います。
このデータセットに対しては本当の本当におわりと言っても良いでしょう。
比例ハザード性が成立しない場合
先に示した比例ハザードの仮定はとても強い仮定です。成立しない場合は多々あります。
この場合はいくつか手段があります。
- 層別にCox比例ハザードモデルを構築する
- 時間依存を取り込む
- 線形項に非線形変換を適用する
これらは発展途中な部分もありますが応用化は進んでいっているようです。
（lifeliesの比例ハザードモデルには層別化を渡す引数があります。）
変数選択
共変量が一つという事は実際の現場でははまれです。いくつかの変数の中から選択しなければいけないと思います。この変数選択部分は分析において最も難しい部分の一つです。
なぜなら変数が複数になれば交互作用等の検証を行う必要が出てくるからです。
ただし、変数選択は一般的な戦略が存在するので、Cox比例ハザードの場合もそういった手法を利用しています。
具体的に言えばステップワイズという手法で、最初にすべての変数を入れておく後ろ向き法(backward)と、最初はすべての変数を入れておかない前向き法(forward)があります。
これについては機会があれば紙面を改めて記載したいと思います。
おわりに
生存時間解析の記事が少なく、基礎的な部分の解説があまりなかったので今回の記事を作成しました。
内容としては下記専門書の前半半分程度の内容なので、言及できていない部分はまだまだあります。
生存時間解析は医療や製品故障に関わりのない方にとっては身近ではないので、データサイエンス全体としてみればマイナーな分野かもしれませんが、様々な分野で活躍できるポテンシャルがある手法だと思います。
下記に参考書籍とリンクを張ります。
Cox比例ハザードモデル
Lifelines
David Cox
以上です。