生存時間解析における多変量解析である「Cox比例ハザードモデル」について解説し、生存時間解析のPythonパッケージである lifelines を用いて、Cox比例ハザードモデルによる生存回帰を行います。本記事は以下の記事の内容に沿っていますが、「比例ハザード性の検証」や「層化」についてはまた別の機会に扱いたいと思います。
関連記事として、「Kaplan-Meier推定量」や「Nelson-Aalen推定量」に関する記事のリンクを載せておきます。
Cox比例ハザードモデルの概要
■ 生存関数とハザード関数
確率変数Xの累積分布関数を$F(x)$とするとき、Xが寿命を表すような場合には、$S(x)=1-F(x)$は「時刻xにまだ生きている確率」を表すという意味で、生存関数 (survival function) と呼ばれます。また寿命Xが連続確率変数として、
$$h(x)=\frac{f(x)}{1-F(x)}=(-logS(x))'$$
をハザード関数 (hazard function) と呼びます。ハザード関数は時刻xにおいて生存している者のうち、その後短時間に死亡する者の率に対応します。
生存時間解析の主要な目的は「生存関数モデルの構築」にありますが、そのためにハザード関数からアプローチしていくのが一般的です。
■ 比例ハザードモデル
比例ハザードモデルでは、対数ハザードが要因（共変量）の線形結合になっていると仮定し、$\alpha, \beta$を推定します。
$$h(t;x)=exp(\alpha+\beta_1x_1+...\beta_px_p)=h_0exp(\beta_1x_1+...\beta_px_p)$$
■ Cox比例ハザードモデル
比例ハザードモデルにおける$h_0$を時間依存性のある$h_0(t)$に置き換えたものがCox比例ハザードモデルです。$h_0(t)=h(t;x=0)$なので、$h_0(t)$は基準ハザード (baseline hazard) と呼ばれます。要因の影響の大きさは、ハザード比やその信頼区間によって評価できます。
$$h(t;x)=h_0(t)exp(\beta_0+\beta_1x_1+...\beta_px_p)$$
生存時間を解析するための要因として1変数しか利用できないカプラン・マイヤー (Kaplan-Meier) 法に対して、Cox比例ハザードモデルでは複数の要因を解析することができます。
■ 比例ハザード性
比例ハザードモデルのもとでは、以下で表されるようにハザード比は時間に依存しません。これを比例ハザード性といい、比例ハザードモデルを適用する際にはこの比例ハザード性を確認する必要があります。
$$\frac{h(t;x)}{h(t;x*)}=exp((x-x*)^{T}\beta)$$
扱うデータセットについて
1970年に米国メリーランド州の刑務所から釈放され、釈放後に一年間追跡調査された受刑者に関するデータです。釈放された受刑者のうち半数は実験的措置としての財政援助を受け、残りの半数は受けませんでした。
| カラム名 | 説明 | |
|---|---|---|
| week | 釈放後の最初の逮捕までの週数、または打ち切り | 生存時間 |
| arrest | 逮捕の有無 (1:有) | 状態変数 |
| fin | 財政援助の有無 (1:有) | 共変量 |
| age | 釈放時の年齢 | 共変量 |
| race | 黒人かどうか (1:黒人) | 共変量 |
| wexp | 労働経験の有無 (1:有) | 共変量 |
| mar | 釈放時の婚姻状況 (1:結婚) | 共変量 |
| paro | 仮釈放かどうか (1:仮釈放) | 共変量 |
| prio | 過去の受刑回数 | 共変量 |
week を生存期間、arrest を状態変数とし、その他のカラムは回帰のための説明変数 (共変量) とします。fin (財政援助の有無) や age (釈放時の年齢) といった共変量が状態変数 arrest にどのように影響しているか解析します。
from lifelines.datasets import load_rossi
rossi = load_rossi()
rossi.head()
| week | arrest | fin | age | race | wexp | mar | paro | prio | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 20 | 1 | 0 | 27 | 1 | 0 | 0 | 1 | 3 |
| 1 | 17 | 1 | 0 | 18 | 1 | 0 | 0 | 1 | 8 |
| 2 | 25 | 1 | 0 | 19 | 0 | 1 | 0 | 1 | 13 |
| 3 | 52 | 0 | 1 | 23 | 1 | 1 | 1 | 1 | 1 |
| 4 | 52 | 0 | 0 | 19 | 0 | 1 | 0 | 1 | 3 |
lifelinesパッケージを用いたCox回帰分析
■ モデルの適合
lifelinesでは CoxPHFitter にCoxモデルが実装されています。fit() でモデルをデータセットに適合し、print_summary() 関数を呼び出すことでテーブル形式で係数やその他の統計情報を確認できます。
from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi
rossi = load_rossi()
cph = CoxPHFitter()
cph.fit(rossi, duration_col='week', event_col='arrest')
cph.print_summary()
カラム coef が推定されたパラメータ$\beta$に相当します。例えば、prio (過去の受刑回数) は係数約0.09と推定されているため、prio が1増加すると基準ハザードは$exp(0.09)=1.10$倍増加することになります。この$exp(0.09)$はハザード比と呼ばれ、比例ハザードモデルの下では時間に依らず一定となります。
労働経験 wexp や財政援助 fin は負の係数をとっていることから、「有」の場合の方がハザードは小さく、すなわち再逮捕されるリスクは低くなり、これは感覚と矛盾しません。ただし、wexp, fin とも95%信頼区間が比較的広くなっているため、解釈には注意が必要です。
plot() メソッドを用いることで、各共変量の係数に対する推定値を信頼区間とともにプロットできます。
cph.plot()
■ 生存関数の予測
生存時間解析におけるゴールは生存関数モデルを構築することです。生存関数の予測には predict_survival_function() メソッドを用います。共変量をデータフレーム形式で引数として渡します。
X = rossi
result = cph.predict_survival_function(X)
rossiデータ先頭2行の受刑者の生存関数を比較してみましょう。共変量の値は以下のようになっています。ageとprioに違いが見られますが、1行目の方が年齢が高く、受刑回数が少ないことから、生存関数は高い値を示すと予想されます。（実際、再逮捕までの期間は1行目の受刑者の方が長くなっています）
rossi.iloc[0:2]
| week | arrest | fin | age | race | wexp | mar | paro | prio | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 20 | 1 | 0 | 27 | 1 | 0 | 0 | 1 | 3 |
| 1 | 17 | 1 | 0 | 18 | 1 | 0 | 0 | 1 | 8 |
sample0 = result.iloc[:, 0]
sample1 = result.iloc[:, 1]
import matplotlib.pyplot as plt
plt.plot(sample0.index, sample0)
plt.plot(sample1.index, sample1)
予想通りの結果となりました。受刑回数が8回ある受刑者は、受刑回数が0回で他の共変量が等しい受刑者と比べて再逮捕のリスクは$exp(8×0.09)=2.05$倍になります。予測された生存関数から、この受刑者が1年間再逮捕されない確率は50%に満たないと考えられます。
■ Formula
v0.25.0 からは、formulaを指定できるようになっています。例えば、次のようなモデルを仮定する場合、formulaを "fin + wexp + age * prio" のように指定します。
$$\beta_1×fin+\beta_2×wexp+\beta_3×age+\beta_4×prio+\beta_5×(age×prio)$$
cph.fit(rossi, duration_col='week', event_col='arrest', formula="fin + wexp + age * prio")
cph.print_summary()
■ 共変量の効果を可視化する
モデルの適合後、他の共変量を一定にしたままある共変量を変化させたときに、生存関数がどのように変化するかを可視化することができます。これには plot_partial_effects_on_outcome() メソッドを使用し、「関心のある共変量」と「値のリスト」を引数として渡します。
cph.fit(rossi, duration_col='week', event_col='arrest')
cph.plot_partial_effects_on_outcome(covariates='prio', values=[0, 2, 4, 6, 8, 10], cmap='coolwarm')
prio の値が大きいほど、生存関数が下に位置していることが視覚的に確認できます。
