生存時間分析の基礎４（Cox 比例ハザードモデル）
秋になり大分涼しくなってきましたが，皆さんお元気でしょうか？
前回から大分間が空いてしまいましたが，本シリーズも第 4 回目となりました．
本シリーズの内容の多くは，「エモリー大学クラインバウム教授の生存時間解析（サイエンティスト社）」にもとづいたものとなっていますので，より詳しく勉強されたい方は合わせてそちらを読んでみることをお勧めします．
さて，今回からモデリングのお話へとうつっていきます．
最終目標は「機械学習の生存時間分析モデルの紹介」ですので，まだまだ先は長そうです．こんなシリーズでもお付き合いいただける方がいましたら嬉しく思います．
本記事の項目は以下となります．
1. セミノンパラメトリックモデル
前回の記事では，2 群（介入群・対照群）における Kaplan - Meier 生存曲線（KM 生存曲線）に差があるかないかを検定するログランク検定をみてきました．
生存時間分析は大別すると，以下の 3 種類に分けることができます．
ノンパラメトリックモデル
確率分布を仮定せず，また共変量（説明変数）を使用しない．セミノンパラメトリックモデル
確率分布は仮定しないが，共変量は使用するモデル．確率分布を仮定しないことに伴う利点・欠点の両方がある．パラメトリックモデル
指数分布やワイブル分布といった確率分布を仮定し，共変量も使用するモデル．確率分布を仮定する関係上，分布の仮定ミスというリスクがあるため，あまり使用されることがない．
前回みてきたログランク検定は，確率分布は仮定していませんでしたし，共変量が生存時間にどう影響を与えるのかを数値的に表す分析ではありませんでした．
そのため，ログランク検定はここでいうところの「ノンパラメトリックモデル」に該当します．
今回見ていく Cox 比例ハザードモデル（Cox PH Model）はセミノンパラメトリックモデルに該当し，生存時間分析のモデルとしては王道の中の王道といえます．
セミノンパラメトリックモデル である Cox PH Model が他の 2 つのモデルと比較してよく使用される理由は，
確率分布の仮定ミスがなく，共変量で調整・回帰係数の推定を行うことができる
ハザード関数の式がその特性上好ましい（以降で説明します）
パラメトリックモデルである logistic 回帰と比較すると，生存時間の情報が利用可能（詳しくはエモリー大学クラインバウム教授の生存時間解析の p112 参照）
といったことが挙げられます．特に 1 つ目の特徴は，セミノンパラメトリックモデルが他の 2 つのモデルと比較して優れている点といえるでしょう．
では，Cox PH Model がどんな形をしているかをみていきましょう．
2. Cox 比例ハザードモデル
Cox PH Model はその名前に「ハザード（Hazard）」とあるように，生存関数ではなくハザード関数のモデリングを行います．
ハザード関数に関しては，本シリーズの第 1 回で紹介をしました．
忘れている方のために少し復習をしておきますと，ハザード関数 $${h(t)}$$ というのは
$${h(t)}$$: 「ある対象が時間 $${t}$$ まで生存している条件下」で，単位時間あたりにイベントが起きる瞬間的な可能性をあらわします．少し分かりにくく，あまり馴染みのない関数かもしれませんが，生存時間分析では必須となる関数です．確率を表す式ではありませんので，1 以上の値になる可能性もあります．
$$
\begin{array}{}
h(t) &=& \lim_{{\Delta}t \to 0} \frac{Pr(t < T < t + {\Delta}t | T \ge t)}{\Delta t} = \frac{f(t)}{S(t)} \; \; \; (continous) \\
\\
&=& Pr(T = t | T \ge t) = \frac{f(t)}{S(t-1)} \; \; \; (discrete)
\end{array}
$$
というものでした（notation など定義が思い出せない場合は，第 1 回記事を振り返ってみてください）
Cox PH Model の場合は，$${h(t)}$$ を少し拡張して $${h(t, X)}$$ と表します．つまり，ハザード関数が時間 $${t}$$ だけではなく，「時間と独立な共変量 $${X}$$ 」にも依存する形になります．これは Cox PH Model がセミノンパラメトリックモデルであることを思い出せば，自然な拡張になっています．
共変量 $${X}$$ が時間 $${t}$$ には依存しないと仮定している点は特に重要ですので頭に留めておきましょう．
ハザード関数 $${h(t, X)}$$ は以下の式で表されます．
$$
\begin{array}{}
h(t, X) &=& h_0(t) \cdot \exp({\sum^{p}_{i=1} \beta_i X_i}) \\
\\
X &=& (X_1, X_2, …, X_p)
\end{array}
$$
上式で表されるハザード関数に関して，いくつかポイントを述べます．
ポイント 1
上式の左辺にある $${h_0(t)}$$ は「基準ハザード関数」と呼ばれ，時間 $${t}$$ のみに依存し，共変量 $${X}$$ からは独立しています．共変量 $${X}$$ が 0 の値をとる時（$${p}$$ 個の共変量が全て 0 となる時）に $${h(t, X=0) = h_0(t)}$$ となりますので，$${h_0(t)}$$ が基準と呼ばれる所以が伺えます．
$${h_0(t)}$$ は（時間 $${t}$$ に依存するため厳密には異なりますが）線形回帰における切片と似たようなものと考えておくと理解がしやすいかもしれません．実際に，$${h(t, X)}$$ の定義式の両辺の対数をとってみると線形回帰の基本式と似たような形が得られます．
ポイント 2
このモデルが「比例」ハザードと呼ばれる理由ですが，異なる共変量同士のハザード関数の比が時間に依存しないことからきています．具体的にそのことをみてみましょう．
二つの共変量 $${X}$$ と $${X^*}$$ があった場合を考えます．それぞれのハザード関数とその比は以下の式で表すことができます．
$$
\begin{array}{}
h(t, X) &=& h_0(t) \cdot \exp({\sum^{p}_{i=1} \beta_i X_i}) \\
\\
h(t, X^*) &=& h_0(t) \cdot \exp({\sum^{p}_{i=1} \beta_i X^{*}_{i}}) \\
\\
\frac{h(t, X)}{h(t, X^*)} &=& \frac{\exp({\sum^{p}_{i=1} \beta_i X_i})}{\exp({\sum^{p}_{i=1} \beta_i X^{*}_{i}})} = constant \; (independent \; of \; time)
\end{array}
$$
上式のようにハザード関数の比をとると基準ハザード関数 $${h_0(t)}$$ は消えてしまい，時間 $${t}$$ には依らず共変量のみによって決まることが分かります．
この仮定は「比例ハザード仮定」と呼ばれ，Cox PH Model の前提となります（この仮定は実際のところかなり強い仮定となっており，その妥当性の検討に関しては次回以降で解説します）
この仮定に従うと，共変量 $${X}$$ と $${X^*}$$ をそれぞれもつ 2 つの対象のハザード関数 $${h(t, X)}$$ の比は時間に依らず一定なわけです．
例えば，2 つの異なる共変量をもつハザード関数を時間に対してプロットした場合，間違ってもどこか途中でハザード関数同士が交差するようなことがあってはいけません（この考え方は次記事の比例ハザード仮定の検討方法で重要となります）
ポイント 3
各共変量に対応する回帰係数 $${β_i}$$ と各共変量 $${X_i}$$ を含んだ線形結合の式が指数 $${e}$$ の乗数になっています．
そのため，線形結合の値が大きな負の値になっても，基準ハザード関数 $${h_0(t)}$$ が正の値であれば，ハザード関数 $${h(t, X)}$$ は負の値にはならないことが分かります（イベントが発生し，一度脱落した対象が復活したりしないことを示します）
ポイント 4
ハザード関数は確率分布を仮定していませんので，確率分布の設定ミスによる推定バイアスなどからは無縁です．このことはここまでみてきた数式から明らかですね．
Cox PH Model におけるハザード関数の式とその性質を確認できました．
では，次にこのモデルの推定方法をみていきましょう．
3. Cox 比例ハザードモデルでの最尤推定
Cox PH Model では最尤法を利用してパラメータ（回帰係数 $${β}$$）の推定を行うのが一般的です．
とはいえ Cox PH Model はセミノンパラメトリックなモデルでしたので，確率分布を仮定していません．これは確率分布をもとにした尤度では推定を行うことができないことを意味します．
ではどうすればよいでしょうか？
Cox PH Model の場合は，「Cox 尤度」と呼ばれる少し特殊な尤度を使いパラメータの推定を行うことになります．
下図の具体例をもとにみてみましょう．
本シリーズの第 1 回の時と同じように，くまがそれぞれの蜂蜜を食べ終えるまでを分析するケースを考えます．くま B は時間が 6 のタイミングで打ち切り（censoring）となっていますが，他のくま達に関しては蜂蜜を食べ終えるイベントが発生するまで観察を行うことができています．
（本記事では同時刻にイベントが複数の対象で発生するような場合は考えません）
まず，くま A ~ D に対応したハザード関数を表しておきます．
ここでは簡単のため，共変量は 1 つだけ $${X_{(A, B, C, D)}}$$ とします．
$$
\begin{array}{}
h_A(t, X_A) &=& h_0(t) \cdot \exp({\beta X_A}) \\
\\
h_B(t, X_B) &=& h_0(t) \cdot \exp({\beta X_B}) \\
\\
h_C(t, X_C) &=& h_0(t) \cdot \exp({\beta X_C}) \\
\\
h_D(t, X_D) &=& h_0(t) \cdot \exp({\beta X_D}) \\
\end{array}
$$
次に，このハザード関数をうまく利用して，くま達のイベント（蜂蜜を食べ終えるというイベント）が発生する「順番」を表現してあげることを考えます．
上図の場合，くま A が $${t=5}$$ で蜂蜜を食べ終えました．そして，他のくま達は蜂蜜を食べ終えていません．
その場合の確率（くま A に最初にイベントが発生する確率）は以下のようにあらわすことができます．
$$
L_{(t=5)} = \frac{h_{A} (t, \; X_A)}{h_{A} (t, \; X_A) + h_{B} (t, \; X_B) + h_{C} (t, \; X_C) + h_{D} (t, \; X_D)}
$$
ハザード関数は，
「ある対象が時間 $${t}$$ まで生存している条件下で単位時間あたりにイベントが起きる瞬間的な可能性」
であったことを思い出しましょう．
$${t=5}$$ において，各くまが蜂蜜を食べ終えるというイベントが発生する確率（ハザード関数）は，それぞれ $${h_A(t=5, X_A), …, h_D(t=5, X_D)}$$ でしたが，くま A だけ蜂蜜を食べ終えたというのが実際に起きたことです．
そう考えると，上式が $${t=5}$$ においてくま A だけが蜂蜜を食べ終えるという事象が発生した確率になっていることが分かると思います．
また，Cox 尤度では打ち切られた対象に関する確率は考慮しません．
例えば，時間が 5 以降におけるイベント発生確率を計算する場合を考えてみましょう．残っているのはくま B，C，D ですが，打ち切りのあるくま B を除き，くま C と D のみで以下のように計算をする形になります．
$$
\begin{array}{}
L_{(t=7)} &=& \frac{h_{C}(t, X_C)}{h_{C}(t, X_C) + h_{D}(t, X_D)} \\
\\
L_{(t=10)} &=& \frac{h_{D}(t, X_D)}{h_{D}(t, X_D)} &=& 1 \\
\end{array}
$$
但し，くま A のイベント発生確率の計算式の分母にはくま B のハザード関数が含まれていた点には注意が必要です（$${L_{(t=5)}}$$ の場合）
なぜならば，くま A のイベントが発生した時点では，くま B の打ち切りは発生していないからです．
打ち切りのある対象 $${i}$$でも，打ち切り時点以前における他の対象のイベント発生確率には，その打ち切り対象 $${i}$$ の生存情報（ハザード関数）が反映されているのです．
ここまでくればもうお分かりかと思いますが，念のため Cox 尤度 $${L}$$ を書き出してみると・・・
$$
\begin{array}{}
L &=& L_{(t=5)} \times L_{(t=7)} \times L_{(t=10)} \\
\\
&=& \frac{\exp({\beta X_A})}{\exp({\beta X_A}) + \exp({\beta X_B}) + \exp({\beta X_C}) + \exp({\beta X_D})} \times \frac{\exp({\beta X_C}) }{\exp({\beta X_C}) + \exp({\beta X_D})} \times 1 \\
\end{array}
$$
となりますが，ひとつ重要な点があります．
上の式の下段で，基準ハザード関数 $${h_{0}(t)}$$ が打ち消されて残っていない点です．
つまり，Cox 尤度においては時間に依存する項が消えるため，いつイベントが発生したかは重要ではないということです．
言い換えれば，重要なのはどの順番でイベントが発生したかということになります．
そのため，例えば下図のように 2 通りのイベントの発生の仕方があった場合，その分析結果（Cox 尤度による推定結果）は全く同じものとなります．
このように Cox 尤度はイベントの発生順だけにもとづくため，部分尤度とも呼ばれます．
では，実際に R で Cox PH Model の推定を行ってみましょう．
4. R による Cox 比例ハザードモデル
今回もこれまで同様に白血病の寛解期間（回復期間のようなもの）のデータである gehan データを分析することにします．
但し，今回は白血球数（White Blood Cells）の対数（logWBC）を共変量として新たにデータに加えることにします．
データは下記のリポジトリ内に用意しておきましたので，興味のある方は手元に pull して以降の分析手順を追ってみると良いかもしれません．
データを読み込み，カラム名をちょっといじると，下図のような形になっています．
# Step 1: Load Data ----
library(tidyverse)
gehan <- read_csv("./input/gehan-wbc.csv",
col_types = cols(treat = col_character())) %>%
mutate(treat = ifelse(treat == "6-MP", "6MP", "control"))
5 列目の logWBC カラムが先ほどお話をした白血球数の対数をとったものです．
gehan データでは白血病からの寛解期間を生存時間として分析しますが，白血病は症状として白血球数（WBC）が極端に大きくなります．そのため，寛解期間は白血球数の大きさに対して正の相関をもつことが予想されます．
また，対数をとる理由の一つとして，極端に大きな白血球数をもつ対象が存在する（右裾が厚い）ため，そういった異常値の影響を緩和するということが考えられます．
次に，前回までと同様に survival::Surv 関数で Surv オブジェクトを作成します．このオブジェクトは Cox 比例ハザードモデルを行う関数 survival::coxph の引数として使用することになります．
# Step 2: Create a survival object ----
library(survival)
surv.obj <- Surv(time = gehan$time, event = gehan$cens)
Cox 比例ハザードモデルの fitting を行ってみると，下のような結果を得ることができます．
# Step 3: Fits a Cox proportional hazards regression model ----
coxph(surv.obj ~ treat + logWBC,
data = gehan,
method = "efron")Call:
coxph(formula = surv.obj ~ treat + logWBC, data = gehan, method = "efron")
coef exp(coef) se(coef) z p
treatcontrol 1.3861 3.9991 0.4248 3.263 0.0011
logWBC 1.6909 5.4243 0.3359 5.034 4.8e-07
Likelihood ratio test=46.71 on 2 df, p=7.187e-11
n= 42, number of events= 30coxph 関数の option である method では，複数の対象が同一時刻にイベント発生（tie）となるような場合における Cox 尤度の近似計算方法を指定しています．Default では method=c("efron","breslow","exact") となっており，"efron" が選択されるようになっています．各 method の詳細（論文）を下記リンクで簡単に紹介していますので，興味のある方はご覧ください．
Mar. 2023 追記
どの version からかは把握していませんが，Mar. 2023 時点の最新 survival パッケージ 3.5.5 では，method option は引き続き指定可能なものの，ties option の方を推奨気味です．現時点では method = ties となっているため，どちらを使用しても問題はないようです．
上のリンク先で紹介している論文では，他の手法よりも「efron の近似手法」がバイアス・信頼区間（効率）の点で優れていると主張しています．
一方で，本記事の参考書籍「エモリー大学クラインバウム教授の生存時間解析」では「breslow の近似手法」を採用して結果を記述していますので，本記事の結果とは異なる値が記載されています．
続いて，coxph モデリングの結果の見方を簡単に解説しておきます．
coef: 各共変量の回帰係数 β の値（treat 共変量の場合: 1.3861）です．ハザード関数の定義式より分かる通り，ハザード関数値にするためには exp(coef) （= 3.9991）に変換する必要があります．
se(coef): 回帰係数 β の標準誤差（= 0.4248）です．
z: coef / se(coef) で得られる Wald 統計量（= 3.263）です．最尤推定量が漸近的に正規分布に分布収束することから，この z を各共変量の有意性（p = 0.0011 ***）の検定に使用します．
Likelihood ratio test: 尤度比検定の値（= 46.71 on 2 df）を示しています．null model（共変量が一切存在しないモデル，つまり切片項のみのモデル）に対する尤度比の値をもとに検定を行います．上の場合では，共変量が treat と logWBC の 2 つ存在するため，尤度比の値が漸近的に従う χ^2 分布の自由度（df）は 2 となっています．
p=7.187e-11 であることから，2 つの共変量を採用した方が null model よりはモデルのあてはまりが良いことが分かります．
さて，treat と logWBC の交互作用項を含む場合の Cox 比例ハザードモデルの fitting は以下のようになります．
# Step 4: Interaction term ----
coxph(surv.obj ~ treat + logWBC + treat:logWBC,
data = gehan,
method = "efron")Call:
coxph(formula = surv.obj ~ treat + logWBC + treat:logWBC, data = gehan,
method = "efron")
coef exp(coef) se(coef) z p
treatcontrol 2.3749 10.7500 1.7055 1.393 0.164
logWBC 1.8724 6.5040 0.4514 4.148 3.35e-05
treatcontrol:logWBC -0.3175 0.7280 0.5258 -0.604 0.546
Likelihood ratio test=47.07 on 3 df, p=3.356e-10
n= 42, number of events= 30交互作用項を含むことで，treat の Wald 統計量による p 値が大きくなってしまっています．また，treat と logWBC の交互作用項はどうやら有意とはいえないようです。上記の尤度比検定は null model に対するものですので，交互作用項を含む場合と含まない場合の尤度比検定を下記のように行ってみます．
# Step 5: Likelihood Ratio Test ----
m1 <- coxph(surv.obj ~ treat + logWBC,
data = gehan,
method = "efron")
m2 <- coxph(surv.obj ~ treat + logWBC + treat:logWBC,
data = gehan,
method = "efron")
anova(m1, m2)Analysis of Deviance Table
Cox model: response is surv.obj
Model 1: ~ treat + logWBC
Model 2: ~ treat + logWBC + treat:logWBC
loglik Chisq Df P(>|Chi|)
1 -69.828
2 -69.648 0.3594 1 0.5488尤度比検定の p 値は 0.5488 ですので，この結果からもどうやら treat と logWBC の交互作用項は有意とはいえないようです．
従って今回は，交互作用項を含まないモデルを採用することになります．
細かい結果をまとめてみたい場合は，総称的関数である summary を使用すれば良いでしょう．
# Step 6: Summary ----
summary(m1)Call:
coxph(formula = surv.obj ~ treat + logWBC, data = gehan, method = "efron")
n= 42, number of events= 30
coef exp(coef) se(coef) z Pr(>|z|)
treatcontrol 1.3861 3.9991 0.4248 3.263 0.0011 **
logWBC 1.6909 5.4243 0.3359 5.034 4.8e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
treatcontrol 3.999 0.2501 1.739 9.195
logWBC 5.424 0.1844 2.808 10.478
Concordance= 0.852 (se = 0.04 )
Likelihood ratio test= 46.71 on 2 df, p=7e-11
Wald test = 33.6 on 2 df, p=5e-08
Score (logrank) test = 46.07 on 2 df, p=1e-10最後に 調整生存曲線 を描いておきましょう．
「調整」とはここでは共変量の値によってハザード関数を条件付けていることを意味します．
生存関数はハザード関数の定義から分かるように下式の関係にあります．
$$
S(\bar{x}, t) = e^{{- \int^t_0 \; h(\bar{x}, u) \; du}}
$$
各共変量 $${X}$$ の回帰係数 $${β}$$ の推定値や標準誤差をもとに，累積ハザードを計算しその負値の指数をとることで，生存関数の推定値や信頼区間を計算することができます（信頼区間の求め方は何種類かありますが，ここでは詳しい説明は割愛します）．
共変量の値は一般的には全サンプルの平均値を使用するため，上式では $${\bar{x}}$$ としています．なお，次で示す R による計算でも平均値が採用されています．
R では以下のようにして調整生存曲線を描くことができます．
# Step 7: Draw Predicted Survival Function ----
library(survminer)
ggsurvplot(
fit = survfit(m1, conf.int = 0.95),
data = gehan,
conf.int = T,
ggtheme = ggplot2::theme_light(),
)survival::survfit 関数を使用して，既に fitting 済みの coxph オブジェクト m1 の生存曲線を計算しています．また，conf.int は信頼区間の幅を指定しています．
survminer::ggsurvplot 関数は第 2 回で紹介したように ggplot-like な生存曲線を描くためのものです．
推定値および信頼区間をもつ調整生存曲線を無事描くことができました．
ですが，投薬の有無による調整生存曲線を描いて見比べることができればもっと良さそうです．
そのためには，以下のようにしてあげます．
# Step 8: Draw Another Predicted Survival Function ----
new.data <- data.frame(
treat = c("control", "6MP"),
logWBC = rep(mean(gehan$logWBC), 2)
)
ggsurvplot(
fit = survfit(m1,newdata = new.data, conf.int = 0.95),
data = gehan,
conf.int = T,
legend.labs = c("control", "6MP"),
ggtheme = ggplot2::theme_light(),
)new.data は投薬の有無の差だけがある（logWBC は全サンプルの平均値を採用）2 行 2 列の擬似データです．
描画してみると上図のように，信頼区間を加味しても，投薬群（6MP）の生存率が非投薬群（control）より高そうな結果が得られました．
ここまでくると想像がつくかと思いますが，new.data として任意の共変量の値をもつ入力を複数設定することで，望みの複数の生存曲線を計算し描画することができます．
さて，少し長くなってしまいましたが，今回は Cox 比例ハザードモデルを確認しました．次回は「比例ハザード仮定の検証方法」の解説をしたいと思います．
Until next time!
5. Codes
今回使用したコードは以下となります．
6. References
この記事が気に入ったらサポートをしてみませんか？