はじめに
scouty AI LAB ではこれまで、生存時間解析に関して概要の記事とノンパラメトリックなモデルに関する解説の記事をそれぞれ公開してきました。
本記事では、その続編としてセミパラメトリックなモデル、特にその中で非常によく用いられるモデルであるCox比例ハザードモデルとその学習の過程について詳しく紹介します。
Cox比例ハザードモデル
Cox比例ハザードモデル (Cox Proportional Hazards model) は、何らかの対象に関する生存時間を分析する生存時間解析 (Survival Analysis) のモデルの一つです。ここで言う生存時間とは、患者の死亡、社員の離職、機械部品の故障といったイベント (event) が生起するまでの時間のことを指します。Cox比例ハザードモデルを用いると、各特徴量がイベントの生起にどれぐらいの影響力をもつのかを定量的に比較することや、将来イベントが発生する可能性を予測することなどが可能になります。
生存時間を解析したい対象（社員、患者、機械部品など）が \(p\) 個あるとき、 \(\{\boldsymbol{x}_i\}_{i = 1}^p\) が対象それぞれのもつ特徴量を表すことにします。また、 \(\boldsymbol{x}_i = (x_{i1}, x_{i2}, \dots, x_{in})^\top\) は \(i\) 番目の対象が持つ \(n\) 次元の特徴量ベクトルを指します。このとき、Cox比例ハザードモデルにおけるハザード関数 \(h(\boldsymbol{x}_i, t)\) は、以下のような式で定義されます：
$$
\begin{equation}
h(\boldsymbol{x}_i, t) = h_0(t) \exp \left({\sum_{j=1}^{n}{{\beta_j x_{ij}}}} \right).
\end{equation}
$$
ここで、上式内の各部分はそれぞれ以下のような意味を持ちます：
- \(h(\boldsymbol{x}_i, t)\): ハザード関数。対象 \(i\) に対し、ある時点 \(t\) の直前までイベントが生起していないという条件の下で、ちょうど時点 \(t\) でイベントが生起する可能性を表す条件付き確率。
- \(h_1(\boldsymbol{\beta}, \boldsymbol{x}_i) = \exp \left( {\sum_{j=1}^{n}{\beta_j x_{ij} }} \right)\): ハザード率。 特徴量 \(\boldsymbol{x}_i\) をもつ対象 \(i\) の瞬時イベント発生リスク。対象のハザード率は時間によって変化しないことを仮定して、右辺の式で対象ごとのハザード率の推定値を計算します。\(\boldsymbol{\beta} = (\beta_1, \beta_2, \dots, \beta_n)^\top\) は各特徴の重みで、学習によってこれを獲得します。文献によっては、\(h_2(\boldsymbol{\beta}, \boldsymbol{x}_i) = {\sum_{j=1}^{n}{\beta_j x_{ij} }}\)をリスク関数 (Risk Function) と呼ぶこともあります。
- \(h_0(t)\): 基準ハザード関数。時間 \(t\) にのみ依存する関数で、対象 \(i\) には依存しません。ハザード率に「時間による変化」を掛けることで、ハザード関数が時間によって変化するようにします。
データの形式
Cox比例ハザードモデルでは、以下のような形式の入力データを扱います：
$$
\begin{equation}
D = {(\boldsymbol{x}_i, t_i, c_i)}
\end{equation}\ \ (i = 1, 2, \dots , p )
$$
ここで、 \(D\) の各要素はそれぞれ以下のように説明されます：
- \(\boldsymbol{x}_i = (x_{i1}, x_{i2}, \dots, x_{in})^\top\): 対象 \(i\) がもつ \(n\) 次元の特徴量ベクトル 。
- \(t_i\): イベント発生時間または打ち切り時間。
- \(c_i\): イベントの発生を表す変数。イベントが発生している場合は \(c_i = 1\) 、そうでない場合は \(c_i = 0\) となる。
時間経過を伴う観察を行うとき、観察期間内に何らかの理由で対象が観察から脱落する場合や、観察終了時点までイベントが発生していない場合があります。これらの場合を打ち切りが生じた (censored) といい、脱落した時の時間や観察が終わった時間を打ち切り時間と言います。上記の \(c_i\) は、 \(i\) が打ち切られている場合は \(c_i = 0\), そうでない場合は \(c_i = 1\) となります。
また、データはタイデータ (tie data) の有無によって分類されます。タイデータとは、 \(i \neq j\) でありながら \(t_i = t_j\) となっているような \(i\) と \(j\) のことで、具体的には以下のいずれかの場合に該当します：
- 2人以上が同時にイベント発生している。
- 2人以上が同時に打切られている。
- 1つ以上のイベントと1つ以上の打切りが同時に生じている。
Cox比例ハザードモデルの学習
Cox比例ハザードモデルの学習では、学習データを用いて、ハザード率のパラメータ \(\boldsymbol{\beta}\) および基準ハザード関数 \(h_0(t)\) を推定します。
- 部分尤度最大化によるハザード率のパラメータ \(\boldsymbol{\beta}\) の推定
- \(h_0(t)\) の推定
以下ではこれらを順に説明します。
1. ハザード率のパラメータ推定
モデルのパラメータ推定のためには、しばしば最尤推定法が用いられます。詳細は省略しますが、データ \(\boldsymbol{x}_1, \boldsymbol{x}_2, \dots, \boldsymbol{x}_p\) が得られており、あるデータ \(\boldsymbol{x}_i\) があるモデル \(P\) から発生する確率 \(P(\boldsymbol{x}_i)\) が計算できるとき、すべての \(i\) に関する発生確率を掛け合わせることで尤度 \(L\) が以下のように定まります：
$$
L = \prod_{i=1}^p P(\boldsymbol{x}_i).
$$
データが手に入れば \(L\) はパラメータ \(\boldsymbol{\beta}\) にのみ依存する関数になります。そこで \(\boldsymbol{\beta}\) を変化させ、 \(L\) を最大にするようなパラメータを探せば、そのモデルがデータ \(\boldsymbol{x}_1, \boldsymbol{x}_2, \dots, \boldsymbol{x}_p\) を最もよく表現しているといえます。このような、尤度 \(L\) を最大にするようなパラメータを求める推定法を最尤推定と呼びます。
さて、Cox比例ハザードモデルでも最尤推定によってパラメータ \(\boldsymbol{\beta}\) の推定を行いたいところですが、残念ながら、基準ハザード関数 \(h_0(t)\) の形が定まっていないため、それによって以下の尤度 \(L\) を計算することができず、パラメータを推定することができません [Kalbfleisch and Prentice, 2011]。
$$
\begin{equation}
L(t, \boldsymbol{\beta}) = \prod_{i=1}^p P(\boldsymbol{x}_i) = \prod_{i=1}^p \left[ h_0(t_i) \exp({\boldsymbol{\beta}^\top \boldsymbol{x}_i}) \right]^{c_i}
\end{equation}
$$
もちろん、基準ハザード関数を先に計算する、あるいは一定の形に決めた後に最尤推定法を適応することは可能ですが、以下に示すCox比例ハザードモデルの基本的な考え方 [Cox, 1992] に反することになります：
- 基準ハザード関数を指定しないことで、時間 \(t\) の分布がどのような分布であってもモデルに含まれるようにする。
- 基準ハザード関数が分からない場合でも、特徴量がイベント発生に対する影響力の比較を可能にする。
これに対応するために、基準ハザード関数を用いずに部分尤度 (partial likelihood) を定義し、最大化することでパラメータ \(\boldsymbol{\beta}\) を近似的に推定します。パラメータの推定方法はタイデータが存在するか否かにより違います。以下では、タイデータが存在する場合、しない場合の順に推定方法を紹介します。
以後の説明で用いる表記を定義します:
- \(R(t_i)\): 時間 \(t_i\) でイベントが発生する可能性のある対象の集合。リスク集合と呼ばれる。
- \(D(t_i)\): 時間 \(t_i\) でイベントが発生した対象の集合。
- \(d_i\): 生存時間が \(t_i\) である人数。 \(d_i = \left| D(t_i) \right|\).
タイデータが存在しない場合の部分尤度の定義
タイデータが存在しないとき、部分尤度 \(L(\boldsymbol{\beta})\) は次のように定義されます:
$$
L(\boldsymbol{\beta}) = \prod_{i=1}^p\left[\frac{\exp({\boldsymbol{\beta}^\top\boldsymbol{x}_i})}{ \sum_{j \in R(t_i)}{\exp({\boldsymbol{\beta}^\top\boldsymbol{x}_j})}}\right]^{c_i}
$$
上式の括弧 \([]\) 内の分数は、特徴 \(\boldsymbol{x}_i\) をもつ被験者 \(i\) に時間 \(t_i\) でイベントが発生したという事実を表現した条件付き確率に該当します。この確率に対し、打ち切られているデータを分子で考慮しないようにするために \(c_i\) 乗を計算しています。 \(p\) 個全てのデータに関しての積を計算することで、部分尤度 \(L(\boldsymbol{\beta})\) を定義します。定義の詳細に関しては [Klein and Moeschberger, 2006] をご参照ください。
タイデータが存在する場合の部分尤度の定義
タイデータが存在する場合には部分尤度の定義が変化しますが、その定義にはいくつかの種類があります。以下では代表的な部分尤度の定義のみを紹介します。さらに詳しい内容については [Zhang, 2005] をご参照ください。
Exact法 (1972)
Cox 本人が提唱した正確な計算方法です。この方法は、時間 \(t_i\) においてタイである被験者 \(1, 2, \dots d_i\) の真の生存時間が微妙に異なると仮定しています。
つまり、これらの \(d_i\) 人の被験者は、測定に十分な精度がない、あるいは便宜的に元のデータが丸められた結果として情報が失われたなどの理由で、データの生存時間が同じであるかのように見えているけれども、実はイベント発生に順序がある、と解釈するということです。
例えば、データ上では時間 \(t_i\) で被験者 \(j\), \(k\) \(\in R(t_i)\) に同時にイベントが発生し、タイとして扱われているときを考えます。この場合、被験者 \(j\), \(k\) の真のイベント発生順序はわかりませんが、以下のいずれかになります：
- 事象 \(A_{j \rightarrow k}\): 被験者 \(j\) のイベントが先に発生し、その後被験者 \(k\) のイベントが発生
- 事象 \(A_{k \rightarrow j}\): 被験者 \(k\) のイベントが先に発生し、その後被験者 \(j\) のイベントが発生
そこで、2つの条件付き確率の和として、タイが生じている事象の発生確率は以下のように求められます。これを各イベント発生時刻について掛けあわせたものを部分尤度と定めます。
$$
P(A_{j \rightarrow k} \cup A_{k \rightarrow j}) = P(A_{j \rightarrow k}) + P(A_{k \rightarrow j}) \\
= \frac{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_j)}{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_1) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_j) + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_k) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_{d_i})} \times \frac{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_k)}{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_1) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_k) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_{d_i})} \\
+ \frac{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_k)}{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_1) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_j) + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_k) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_{d_i})} \times \frac{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_j)}{\exp(\boldsymbol{\beta}^\top\boldsymbol{x}_1) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_j) + \dots + \exp(\boldsymbol{\beta}^\top\boldsymbol{x}_{d_i})}
$$
自明ですが、タイが存在しないデータにExact法を適用して計算した部分尤度 \(L(\boldsymbol{\beta})\) は、上で定義したタイデータがない時の部分尤度関数と同じ値になります。
一方、上記のタイデータがある例では \(d_i = 2\) でしたから \(2! = 2\) 通りの事象を考えましたが、一般には各 \(t_i\) について \(d_i!\) 通りの事象を考える必要があります。これによって、 \(d_{i}\) に対応するタイデータについて \(d_i!\) 通りの場合を考慮した計算が必要になり、タイデータが多く存在する場合は推定が困難になります。
Breslow による近似法 (1974) と Efron による近似法 (1977)
Breslow による近似法[Breslow, 1974] では、部分尤度 \(L(\boldsymbol{\beta})\) を次のように定義します:
$$
L(\boldsymbol{\beta}) =\prod _{ i=1 }^{ p }{ \frac { \exp(\boldsymbol{\beta}^\top{ \boldsymbol{x}_{(i)+} }) }{ \left[ \sum _{ j\in R\left( t_{i} \right) }{ \exp({ \boldsymbol{\beta}^\top\boldsymbol{x}_j }) } \right] ^{ d_{ i } } } }
$$
一方、 Efron による近似法[Efron, 1977] では、部分尤度 \(L(\boldsymbol{\beta})\) を次のように定義します:
$$
L(\boldsymbol{\beta}) =\prod _{ i=1 }^{ p }{ \frac { \exp(\boldsymbol{\beta}^\top{ \boldsymbol{x}_{(i)+} }) }{ \prod _{ k=1 }^{ d_{ i } }{ \left[ \sum _{ j\in R\left( t_{ i } \right) }{ \exp({ \boldsymbol{\beta}^\top\boldsymbol{x}_j }) } -\frac { k-1 }{ { d }_{ i } } \sum _{ j\in D\left( t_{ i } \right) }{ \exp({ \boldsymbol{\beta}^\top\boldsymbol{x}_j }) } \right] } } }
$$
ここで、上記2つの式に現れた \(\boldsymbol{x}_{(i)+}\) は \(D(t_i)\) に含まれる \(d_{i}\) 人のすべての特徴量の和を表します:
$$
\boldsymbol{x}_{(i)+} = \sum _{j \in D\left( t_{i} \right)}{ \boldsymbol{x}_j }
$$
Breslow による近似法は、各時刻 \(t_i\) でのリスク集合に対するタイデータの割合 \(d_i \left/ R(t_i) \right.\) が小さければうまく機能することが知られていますが、この割合が大きくなると近似が悪くなる可能性があります。一方、 Efron による近似法ではこの割合が大きい場合に Breslow の近似法よりかなりいい近似が得られる [Kalbfleisch and Prentice, 2011] とされています。一方で、 Efron の近似法は Breslow の近似法に比べてやや計算が重いという弱点はあります。
2018年3月現在、生存時間解析のPythonライブラリ Lifelines のCox比例ハザードモデルの実装では Efron の近似法が実装されています。
部分尤度の最大化
さて、部分尤度が定まったのちには、この部分尤度を最大化することでパラメータ \( \boldsymbol{\beta} \) を推定します。部分尤度最大化のために用いる最適化手法は何でも構いませんが、 [Chen, 2006] や [Kalbfleisch and Prentice, 2011] ではニュートン・ラフソン法 (Newton-Raphson method) を利用する例が示されています。この場合には、部分尤度 \(L(\boldsymbol{\beta})\) の自然対数 \(\ln L(\boldsymbol{\beta})]\) の2階偏導関数まで計算した後に反復処理でパラメータを推定します
2. 基準ハザード関数の推定
Cox比例ハザードモデルを用いて、ある時間 \(t\) で、新しい特徴量をもつ被験者のイベント発生率を予測しようとする時、基準ハザード関数が必要になります。また、上では基準ハザード関数を考慮せずにパラメータ推定を行ったため、様々な関数を基準ハザード関数として採用することが可能になります。
指数分布 (exponential distribution), ワイブル分布 (Weibull distribution) などパラメトリックな分布を基準ハザード関数として指定し、最尤推定法でパラメータを推定するのも可能ですし、ノンパラメトリックな手法を使って分布を決めてしまうことも可能です。ただし、処理データに対する時間分布が分からないときは、ノンパラメトリックな手法で基準ハザード関数を推定することが一般的です1。
ノンパラメトリックな推定については、 [Kalbfleisch and Prentice, 2011], [Weng, 2007] などで示された手法が存在します。
今回は、生存時間分析パッケージ Lifelines で用いられている Nelson-Aalen 推定量 を変形した実装手法を紹介します。この手法では基準ハザード関数 \(h_0(t_i)\) は以下のように定義されます:
$$
\begin{equation}
h_0(t_i) = \frac{d_i}{\sum_{j \in R(t_i)}{\exp({\boldsymbol{\beta}^\top\boldsymbol{x}_j })}}
\end{equation}
$$
この手法では、トレーニングデータを基準に各時間 \(t\) でのイベント発生率を計算し、それを時間 \(t\) での基準ハザード値としてみなすことで基準ハザード関数を推定していると解釈することができます。
おわりに
本記事ではCox比例ハザードモデルとその学習の過程を紹介しました。
次の記事では、さらに、Cox比例ハザードモデルを拡張したモデルと、それを実データに適用した結果を紹介する予定です。
参考文献
[Breslow, 1974] N. E. Breslow "Covariance analysis of censored survival data." Biometrics, 30, 1974, 89-100.
[Chen, 2006] Zehua Chen, "Regression Models for Survival Data." ST3242: Introduction to Survival Analysis, 2006. 1-16.
[Cox, 1992] David R. Cox. "Regression Models and Life-Tables." Breakthroughs in statistics. Springer, New York, NY, 1992. 527-541.
[Efron, 1977] Bradley Efron. "The efficiency of Cox’s likelihood function for censored data". Journal of the American Statistical Association, 72. 1977, 557-565.
[Kalbfleisch and Prentice, 2011] John D. Kalbfleisch and Ross L. Prentice. "The Statistical Analysis of Failure Time Data.", John Wiley & Sons, 2011, 462.
[Klein and Moeschberger, 2006] John P. Klein and Melvin L. Moeschberger. "Survival Analysis: Techniques for Censored and Truncated Data." Springer Science & Business Media, 2006.
[Weng, 2007] Yi-Ping Weng. "Baseline Survival Function Estimators under Proportional Hazards Assumption." 2007.
[Zhang, 2005] Daowen Zhang. "Cox Proportional Hazards Regression Models (cont’d)." Analysis of Survival Data (ST745), 2005. 142-179.
脚注
- そもそも、基準ハザード関数をパラメトリックな分布に決めてしまうと、Cox比例ハザードモデルはパラメトリックな手法になります。 ↩