ここまでで紹介してきた生存関数の推定やハザード関数の推定、ログランク検定といったノンパラメトリックな手法は、一標本の生存時間データ解析、二群間の比較に有用な方法になります。
一方で、実際の医学研究では、個々の被験者について、生存・死亡といったアウトカムの情報だけでなく、健診項目がどう変動したか、性別や年齢など生存時間に影響を与える補助的な情報も得られます。これらの値は説明変数と呼ばれ、生存時間データをモデル化することで、観測時点の起点における各患者の説明変数の値がその患者グループの生存にどのような影響を与えているか特定することができます。
生存時間解析では、ハザード関数をモデル化の対象（目的変数）とします。
このモデル化には大きく2つの目的があります。
- 潜在的な説明変数のどの変数がハザード関数の形状に影響を与えているか特定する
- ある被験者のハザード関数・生存関数そのものの推定値を得る
得られた説明変数の影響、生存関数の推定値などは治療計画の立案や患者の予後に関するカウンセリングに有用な情報となり得ます。
セミパラメトリック回帰モデル
上記で前触れもなく「ノンパラメトリック」というワードを使用しました。
パラメータ依存の分布を仮定したモデルを「パラメトリックモデル」、パラメータ依存の分布を仮定しないモデルを「ノンパラメトリックモデル」といいます。
パラメトリックなモデルは、大きく下記の用途で使用されます。
- 目的変数が従う分布を表現する
- 説明変数の関数として、分布がどのように変化するか特徴づける
パラメトリックなモデルは多くの分析で利用されてきましたが、ある状況においては後者のみが重要であり、前者に必要な制約が非現実的となる場合があります。
そのような場合に採用されるのがセミパラメトリックモデルです 。
セミパラメトリック回帰モデルでは、一部にパラメトリックなモデルを採用することで完全な回帰構造を担保する一方、その他の部分を不特定のままモデル化することで説明変数と目的変数の関係を表現したまま精度の高いモデリングを可能としています。
ハザード関数について、時間に依存する部分を不特定（$h_0(t)$）としたセミパラメトリック回帰モデルを定義すると、下記のように表せます。
$$ h(t, x, \beta) = h_0(t)r(x, \beta) $$
ここで、$h_0(t)$は、$r(x=0, \beta)$のときのハザード関数とみなせることから、ベースラインハザードと呼ばれます。
ハザード関数をセミパラメトリック回帰モデルでモデリングすると、ハザード比をパラメトリックなモデルで表すことができます。 例えば、共変量 $x_0, x_1$のときのハザード比を考えると、ベースラインが分母分子で打ち消されるため、下記となります。
$$ \begin{align} Hazard\ Ratio(t, x_0, x_1) & = \frac{h(t, x_1, \beta)}{h(t, x_0, \beta)} \\ & = \frac{r(x_1, \beta)}{r(x_0, \beta)} \end{align} $$
Cox (1972) は、$r(x, \beta) = \exp(x\beta)$ とした下記のモデルを提案しました。
$$ h(t, x, \beta) = h_0(t)\exp(x\beta) $$
このモデルは、Cox比例回帰モデル、Cox比例ハザードモデルなどと呼ばれます。
Cox比例ハザードモデル
一般化したCox比例ハザードモデルは、下記のように表せます。
$$ \begin{align} h(t|\beta_1, \beta_2, \cdots, \beta_n) & = h_0(t) \exp{(\beta_1x_1+\beta_2x_2+\cdots\beta_nx_n)} \ \end{align} $$
上記で触れたように、Cox比例ハザードモデルのハザード比はパラメトリックなモデルで表すことができ、
$$ Hazard\ Ratio(t, x_1, x_0) = e^{\beta(x_1 - x_0)} $$
となります。
例として、説明変数が性別のように二値変数の場合（男性が1、女性が0等）、
$$ Hazard\ Ratio(t, x_1, x_0) = e^{\beta} $$
となります。
ハザード比が乗法的に関与するような前提を置くことから、比例ハザードモデルと名づけらたようです。
生存関数
以前の記事で触れたように、累積ハザード関数と生存関数は下記のような関係にあります。
$$ S(t, x, \beta) = e^{-H(t, x, \beta)} $$
累積ハザード関数は、ハザード関数の積分であるため、
$$ \begin{align} H(t, x, \beta) & = \int_0^t{h(u, x, \beta)}du \\ & = r(x, \beta)\int_0^t{h_0(u)}du \\ & = r(x, \beta)H_0(t) \end{align} $$
よって、パラメトリックハザード関数の生存関数は、
$$ \begin{align} S(t, x, \beta) & = e^{-r(x, \beta)}H_0(t) \\ & = [e^{-H_0(t)}]^{r(x, \beta)} \\ & = [S_0(t)]^{r(x, \beta)} \\ \end{align} $$
となります。
ここで、$S_0(t)$はベースラインハザードに対応するベースライン生存関数となります。
以上より、Cox比例ハザードモデルにおける生存関数は
$$ S(t, x, \beta) = [S_0(t)]^{exp(x, \beta)} $$
と表せます。
Cox比例ハザードモデルのあてはめ
モデルの線形成分の回帰係数（$\beta_1, \beta_2, \cdots \beta_n$）及びベースラインハザード$h_0$を推定しなければなりません。
この2つは別々に推定することができ、回帰係数（$\beta_1, \beta_2, \cdots \beta_n$）を推定した後にこれらの値を使ってベースラインハザード$h_0$推定する手順となります。
最尤推定
最尤推定については各所で丁寧な解説記事があるのでそちらを参考にしていただければと思いますので軽めに。
標本データの尤度は観測データの同時確率であり、仮定するモデルにおける未知のパラメータの関数とみなされます。この関数を尤度関数と呼びます。最尤推定法では、この尤度関数を最大化することによって尤もらしいパラメータの値を推定することになります。
この最大化では、計算上尤度関数の対数を最大化する方が都合が良いため、一般に対数尤度関数を最大化することでパラメータの推定値を求めます。
また、最尤推定法に依って求まったパラメータ推定値の分散の近似値は尤度関数の二回導関数から求めることができます（クラメール・ラオの限界）。
Cox比例ハザードモデルにおける尤度関数
Coxは比例ハザードモデルの尤度関数が以下のように与えられることを示しました（尤度関数が下記であらわせることの妥当性は後日）。
$$ L(\boldsymbol{\beta}) = \prod_{j=1}^r{\frac{\exp{(\boldsymbol{\beta}'\boldsymbol{x}_{(j)})}}{\sum_{l\in R(t(j))}\exp{(\boldsymbol{\beta}'\boldsymbol{x}_l)}}} $$
$x_{(j)}$：$j$番目の死亡時刻$t_{(j)}$で死亡した被験者の共変量ベクトル
$R(t(j))$：$j$番目の死亡時刻$t_{(j)}$におけるリスク集合
$\delta_j$を打ち切りが生じた場合は1、生じていない場合は0とする指示変数とすると、 打ち切りを考慮した尤度関数は下記のように表せます。
$$ L(\boldsymbol{\beta}) = \prod_{j=1}^r\Bigg\{{\frac{\exp{(\boldsymbol{\beta}'\boldsymbol{x}_{(j)})}}{\sum_{l\in R(t(j))}\exp{(\boldsymbol{\beta}'\boldsymbol{x}_l)}}}\Bigg\}^{\delta_j} $$
依って、対数尤度関数は下記のように表せます。
$$ \log{L(\boldsymbol{\beta})} = \sum_{j=1}^r\delta_j\Bigg\{{\boldsymbol{\beta}'\boldsymbol{x}_{(j)} - \log{\sum_{l\in R(t(j))}\exp{(\boldsymbol{\beta}'\boldsymbol{x}_l)}}}\Bigg\} $$
$\boldsymbol{\beta}$の最尤推定値は数値計算によってこの関数を最大化することによって求めることができます。 通常この尤度関数の最大化にはニュートン・ラフソン法を用います（後日説明）。
lifelinesでCox比例ハザードモデルによる推定
lifelinesパッケージで提供されているテストデータ（rossi）を用いてCox比例ハザードモデルを用いた推定を行います。
rossiには、432個の観測値が含まれています。week の列は期間であり、arrest の列はイベント（再逮捕）が発生したかどうかを示しています。他の列を説明変数に使用して、各説明変数とイベント発生の関係性をモデル化します。
>>> 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
Cox比例ハザードモデルへの当てはめを行う際は、CoxPHFitterを用います。
難しいこと考える必要なしに.fitメソッドのduration_colに期間に該当する列名、event_colにイベントのあり/なしに該当する列名を指定することで対数尤度関数を最大化し、係数の推定を行ってくれます。
from lifelines import CoxPHFitter cph = CoxPHFitter() cph.fit(rossi, duration_col='week', event_col='arrest')
.plot()メソッドにて、推定された係数（対数ハザード比）の値とその95%信頼区間をプロットできます。
cph.plot()
print_summary関数でモデルについて推定の結果求まった係数と関連する統計値を表形式で出力することができます。
>>> cph.print_summary() <lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations> duration col = 'week' event col = 'arrest' baseline estimation = breslow number of observations = 432 number of events observed = 114 partial log-likelihood = -658.75 time fit was run = 2020-11-01 07:37:03 UTC --- coef exp(coef) se(coef) coef lower 95% coef upper 95% exp(coef) lower 95% exp(coef) upper 95% covariate fin -0.38 0.68 0.19 -0.75 -0.00 0.47 1.00 age -0.06 0.94 0.02 -0.10 -0.01 0.90 0.99 race 0.31 1.37 0.31 -0.29 0.92 0.75 2.50 wexp -0.15 0.86 0.21 -0.57 0.27 0.57 1.30 mar -0.43 0.65 0.38 -1.18 0.31 0.31 1.37 paro -0.08 0.92 0.20 -0.47 0.30 0.63 1.35 prio 0.09 1.10 0.03 0.04 0.15 1.04 1.16 z p -log2(p) covariate fin -1.98 0.05 4.40 age -2.61 0.01 6.79 race 1.02 0.31 1.70 wexp -0.71 0.48 1.06 mar -1.14 0.26 1.97 paro -0.43 0.66 0.59 prio 3.19 <0.005 9.48 --- Concordance = 0.64 Partial AIC = 1331.50 log-likelihood ratio test = 33.27 on 7 df -log2(p) of ll-ratio test = 15.37
まずはprio（過去の犯罪件数）に着目してみます。この説明変数の係数は0.09と推測されていることから、犯罪件数に関するハザード比の値は$\exp(0.09) = 1.094$ であることがが分かります。これは、過去の犯罪件数が1上がるごとに再犯率が約10%上昇することを意味しています。
また、mar（既婚の場合 1、未婚の場合0の） に着目すると、この説明変数の係数は-0.43と推測されていることから、結婚に関連するハザード比の値は、$\exp(-0.43) = 0.650$であることが分かります。
よって、既婚の人の再犯率は未婚の人の65%であることが分かります。
他のすべてを等しく保ちながら単一の共変量を変化させたときに、生存曲線がどのように変化するかをプロットできます。 過去の逮捕数が2回、4回、6回、8回、10回の場合の生存関数をプロットすると、
rossi = load_rossi() cph = CoxPHFitter() 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')
50週間後の再犯率は、過去の犯罪数が2回の場合は約80%であるのに対し、過去の犯罪数が10回の場合は60%を下回ることが分かります。
このように、セミパラメトリックモデルを構築すると各説明変数と目的変数の関係性を把握することが可能となります。
おまけ
収束後の対数尤度関数の値は.log_likelihood_から分散行列の値は.variance_matrix_からそれぞれ取得することができます。
# 収束後の対数尤度関数 >>> cph.log_likelihood_ -658.7476594460861 # 分散行列 >>> cph.variance_matrix_ covariate fin age race wexp mar paro prio covariate fin 0.036626 -0.000259 -0.003744 0.000066 0.001877 0.002661 -0.000202 age -0.000259 0.000484 -0.000303 -0.001437 -0.001018 0.000474 -0.000019 race -0.003744 -0.000303 0.094860 0.001942 0.010422 -0.001894 0.000505 wexp 0.000066 -0.001437 0.001942 0.045039 -0.014642 -0.001623 0.001700 mar 0.001877 -0.001018 0.010422 -0.014642 0.145823 -0.004449 -0.000484 paro 0.002661 0.000474 -0.001894 -0.001623 -0.004449 0.038321 0.000843 prio -0.000202 -0.000019 0.000505 0.001700 -0.000484 0.000843 0.000821
