Cox比例ハザード解析

（説明）
異なる説明変数ベクトル$\vec x_i, \vec x_j$に対して、そのハザード比$\dfrac{h(t;\vec x_j)}{h(t;\vec x_i)}$は、
$$\begin{aligned} \dfrac{h(t;\vec x_j)}{h(t;\vec x_i)} &= \dfrac {h_0 (t) \cdot \exp \{ \beta_1 x_{j,1} + \cdots + \beta_d x_{j,d} \}} {h_0 (t) \cdot \exp \{ \beta_1 x_{i,1} + \cdots + \beta_d x_{i,d} \}} \\[10px] &= \exp \{ \beta_1 (x_{j,1}-x_{i,1}) + \cdots + \beta_d (x_{j,d}-x_{i,d}) \} \\[10px] &= \small{(Const)} \end{aligned}$$となります。




つまり、$h(t;\vec x_j)$は$h(t;\vec x_i)$の時間$t$に依らない定数倍になっているということで、これを『比例ハザード性』と言います。

（補足）
比例ハザード性が成立するのは、ひとえにベースラインハザード$h_0 (t)$が共通しているためです。




またここで、
$$\begin{aligned} \vec x_i &=(x_{1},\cdots, x_{k-1}, x_{k}, x_{k+1},\cdots, x_{d})^{\top} \\ \vec x_j &=(x_{1},\cdots, x_{k-1}, x_{k} \textcolor{red}{+ 1}, x_{k+1},\cdots, x_{d})^{\top} \end{aligned}$$と説明変数ベクトルの$k$番目が1だけずれた場合について、そのハザード比は、
$$\begin{aligned} \dfrac{h (t; \vec x_j)}{h (t; \vec x_i)} &= \exp \{ \beta_1 (x_1-x_1) + \cdots + \beta_{k-1} (x_{k-1}-x_{k-1}) + \beta_k (x_k \textcolor{red}{+ 1}-x_k) + \\ &\beta_{k+1} (x_{k+1}-x_{k+1}) + \cdots + \beta_d (x_d-x_d) \} \\ &= \exp (\beta_k) \end{aligned}$$となります。




つまり$k$番目の説明変数が1だけ大きくなる時、ハザード関数は$\exp (\beta_k)$倍となります。この事実は回帰係数$\beta_k$の解釈として重要です。

（部分尤度の導出）
まず尤度関数$L(\vec \beta, h_0)$は、
$$\begin{aligned} L(\vec \beta, h_0) &\propto \prod_{i=1}^{n} \{ h(t_i) \}^{E_i} \cdot \exp \{ – H(t_i) \} ~~~~~ \mathrm{(*)} \\[10px] &= \prod_{i=1}^{n} \{ h_0 (t_i) \cdot \exp (\vec \beta^{\top} \vec x_i) \}^{E_i} \cdot \exp \{ -H_0 (t_i) \cdot \exp (\vec \beta^{\top} \vec x_i) \} \\[10px] &= \prod_{j=1}^{D} \{ h_0 (t_j) \cdot \exp (\vec \beta^{\top} \vec x_j) \} \cdot \exp \{ – \sum_{i=1}^{n} H_0 (t_i) \cdot \exp (\vec \beta^{\top} \vec x_i) \} \\ &{\scriptsize (t_jはイベント発生時点を表す)} \end{aligned}$$となります。

$(*)$：<生存時間解析の基本概念> > 3.1では、$C_i = \vec 1 \{ V_i \leqq X_i \}$ （打ち切りがなければ$0$、あれば$1$）、を用いて尤度$L$は、
$$\begin{aligned} L &\propto \prod_{i=1}^{n} f (t_i)^{1-C_i} S(t_i)^{C_i} \end{aligned}$$と紹介しました。




今回は$C_i$ではなく$E_i$を用いてますが、
$$\begin{aligned} C_i = 1-E_i \end{aligned}$$の関係にあります。

すたどく

これは、最後まで観察された症例も形上打ち切りが生じたとみなすと、すべての症例が打ち切りあり（イベント発生なし）、または、打ち切りなし（イベント発生あり）、に割り振ることができるからですね。

これと、<生存時間解析の基本概念> > 2 > おまけ1/2、で紹介したとおり、

・$H(t) = – \log S(t)$
・$f(t) = h(t) S(t)$

であることを用いると以下の様になります。
$$\begin{aligned} L &\propto \prod_{i=1}^{n} f (t_i)^{1-C_i} S(t_i)^{C_i} \\[10px] &= \prod_{i=1}^{n} f (t_i)^{E_i} S(t_i)^{1-E_i} \\[10px] &= \prod_{i=1}^{n} \{ h(t_i) S(t_i) \}^{E_i} S(t_i)^{1-E_i} \\[10px] &= \prod_{i=1}^{n} h(t_i)^{E_i} S(t_i) \\[10px] &= \prod_{i=1}^{n} h(t_i)^{E_i} \cdot \exp\{ -H(t_i) \} \end{aligned}$$

（部分尤度の導出〜第一段階（開始））
$L(\vec \beta, h_0)$を最大化することを考えるにあたり、まずは$\vec \beta$をfixした上で$h_0 (,H_0)$について最大化することを考えます。




$H_0$が小さいほど$L(\vec \beta, h_0)$が大きくなるので、イベントが発生する時点を除いて$h_0 (t) = 0$が成立する時に$L(\vec \beta, h_0)$は最大化されます。




そこで、
$$\begin{aligned} h_{0j} = h_0 (t_j) ~~~~~ \small{(j=1, \ldots, D)} \end{aligned}$$とすると、
$$\begin{aligned} H_0 (t_i) &= \sum_{t \leqq t_i} h_{0} (t) \\ &= \sum_{j:t_{j} \leqq t_i} h_{0j} \end{aligned}$$と書くことができます。




よって、
$$\begin{aligned} L_{\vec \beta} (h_{01}, \ldots, h_{0D}) &\propto \prod_{j=1}^{D} \{ h_{0j} \cdot \exp (\vec \beta^{\top} \vec x_j) \} \cdot \exp \{ – \sum_{i=1}^{n} \sum_{j:t_{j} \leqq t_i} h_{0j} \cdot \exp (\vec \beta^{\top} \vec x_i) \} \\[10px] &= \prod_{j=1}^{D} \{ h_{0j} \cdot \exp (\vec \beta^{\top} \vec x_j) \} \cdot \exp \{ – \sum_{j=1}^{D} \sum_{i \in R_j} h_{0j} \cdot \exp (\vec \beta^{\top} \vec x_i) \} \\ &{\scriptsize(ただしR_jはt_j直前においてイベントが発生しうる個体の集合を表す)} \\[10px] &= \prod_{j=1}^{D} \{ h_{0j} \cdot \exp (\vec \beta^{\top} \vec x_j) \} \cdot \exp \{ – \sum_{j=1}^{D} h_{0j} \cdot \sum_{i \in R_j} \exp (\vec \beta^{\top} \vec x_i) \} ~~~~~ \mathrm{(A)} \end{aligned}$$


$\log$をとると、
$$\begin{aligned} \log L_{\vec \beta} (h_{01}, \ldots, h_{0D}) = \sum_{j=1}^{D} \{ \log h_{0j} + \vec \beta^{\top} \vec x_j \} – \sum_{j=1}^{D} h_{0j} \cdot \sum_{i \in R_j} \exp (\vec \beta^{\top} \vec x_i) \end{aligned}$$


これを$h_{0j} ~~ \small{(j=1,\ldots,D)}$について最大化すると（上式右辺を$h_{0j}$で偏微分して$=0$としたものを解くと）、
$$\begin{aligned} \hat{h}_{0j} &= \dfrac {1} {\sum_{i \in R_j} \exp (\vec \beta^{\top} \vec x_i)} \\ \Rightarrow \hat{H}_0 (t) &= \sum_{j:t_{j} \leqq t_i} \dfrac {1} {\sum_{i \in R_j} \exp (\vec \beta^{\top} \vec x_i)} \end{aligned}$$
（部分尤度の導出〜第一段階（終了））


（部分尤度の導出〜第二段階（開始））
つづいて上記のとおり求めた$\hat{h}_0$を代入して$\vec \beta$について最大化することを考えます。


$$\begin{aligned} L (\vec \beta, \hat{h}_0) &\propto \prod_{j=1}^{D} \{ \hat{h}_{0j} \cdot \exp (\vec \beta^{\top} \vec x_j) \} \cdot \underbrace{ \prod_{j=1}^{D} \exp \{ – \textcolor{red}{\hat{h}_{0j} \cdot \sum_{i \in R_j} \exp (\vec \beta^{\top} \vec x_i)} \}}_{\exp内の赤文字部分は上記より1になるため、ここは定数} \\ &{\scriptsize(\mathrm{(A)}を少々変形した)} \\[10px] &\propto \prod_{j=1}^{D} \dfrac {\exp (\vec \beta^{\top} \vec x_j)} {\sum_{i \in R_j} \exp (\vec \beta^{\top} \vec x_i)} ~~~~~ \mathrm{(B)} \\ &{\scriptsize(上記\hat{h}_0を代入した)} \end{aligned}$$


この時$\mathrm{(B)}$は冒頭の形と一致しています。
（部分尤度の導出〜第二段階（終了））

すたどく

実は、今回ここで導出したものは『profile尤度』と呼ばれるものです。あるパラメータについて尤度を最大化して、また別のパラメータについて尤度を最大化して、、、としていく導出を『ヒューリスティックな導出』と言いますが、それにより導出されるのがこの『profile尤度』なのです。
『部分尤度』自体はCoxさんが定義したものであり（興味ある方は原著を参照ください）、profile尤度という別の視点からも同じ形がでてくることを確認したわけです。

(説明)
$h(t; \vec x_i) = h_0 (t) \cdot \exp (\vec \beta^{\top} \vec x_i)$を前提とします。




すると、生存関数$S(t; \vec x_i)$について、
$$\begin{aligned} S(t; \vec x_i) &= \exp [-\int_0^t h(s; \vec x_i) ds] \\ &{\scriptsize(h(t) = -\frac{d}{dt} \log S(t)より変形)} \\[10px] &= \exp [-\int_0^t h_0 (s) \cdot \exp (\vec \beta^{\top} \vec x_i) ds] \\ &{\scriptsize(h(t; \vec x_i) = h_0 (t) \cdot \exp (\vec \beta^{\top} \vec x_i))} \\[10px] &= \{ S_0 (t) \}^{\exp (\vec \beta^{\top} \vec x_i)} \\ &{\scriptsize(ただしS_0 (t) := \exp[-\int_0^t h_0 (s) ds])} \end{aligned}$$となり、これより、
$$\begin{aligned} \log ( -\log S (t; \vec x_i) ) = \log ( -\log S_0 (t) ) + \vec \beta^{\top} \vec x_i \end{aligned}$$となります。




よって、異なる$\vec x_i, \vec x_j$について、
$$\begin{aligned} \log ( -\log S (t; \vec x_i) ) &= \textcolor{red}{\log ( -\log S_0 (t) )} + \vec \beta^{\top} \vec x_i \\[10px] \log ( -\log S (t; \vec x_j) ) &= \textcolor{red}{\log ( -\log S_0 (t) )} + \vec \beta^{\top} \vec x_j \end{aligned}$$となります。（赤文字部分は共通）

つまり、比例ハザード性が成立していれば、$\log ( -\log S (t; \vec x_i) ), \log ( -\log S (t; \vec x_j) )$は、$\vec \beta^{\top} \vec x_i-\vec \beta^{\top} \vec x_j$という定数分だけずれることになります。


これはプロット上においては、$\log ( -\log S (t; \vec x_i) )$と$\log ( -\log S (t; \vec x_j) )$が平行であることと同じです。

学習者

生存関数$S(t)$は<カプランマイヤー推定量・曲線とログランク検定>であつかった通りプロットすることができましたが、ここでは$\log ( -\log S (t) )$を層別にプロットして比例ハザード性を確認するわけですね。ちなみに、比例ハザード性を確認する検定はあるのでしょうか？

すたどく

はい。比例ハザード性を確認する検定として『Schoenfeld残差検定』などがありますが、この様な検定では帰無仮説が『比例ハザード性が成立する』と置かれてるため、帰無仮説が棄却された場合には積極的に『比例ハザード性が不成立』と主張できます。しかし、見た目には比例ハザード性が成立していそうでも、サンプルサイズやイベント発生数が大きくなると小さな差をも検出してしまい、帰無仮説が棄却される可能性が高くなります。


比例ハザード性の確認にあたってはこの様な検定に言及した医学論文は数少ないため、基本的にはプロットをして比例ハザード性が成立してそうなことを確認すればよいのかと思います（私見です）。

学習者

わかりました。ところでプロット的に明らかに比例ハザード性が成立してない時にはどのように解析すればよいのでしょうか？

すたどく

比例ハザード性が成立しない場合にはCox比例ハザードモデルを仮定した解析は成立しません。その場合には以下の様な対応が可能ですが、いずれが望ましいかについて一定の見解はない様です。


①時間依存性の説明変数をモデルに組み込む
②層別に（ある説明変数についてそのカテゴリーごとに）、異なるベースラインハザードを用いる
③$h(t; \vec x_i) = h_0 (t) \cdot \exp { \beta_1 x_{i,1} + \cdots + \beta_d x_{i,d} }$の、”$\beta_1 x_{i,1} + \cdots + \beta_d x_{i,d}$”に対して非線形変換を施す

まとめ.

• Cox比例ハザードモデルに基づく解析においては、部分尤度を最大化するパラメータを求める。
  
  
• $k$番目の説明変数が1だけ大きくなる時、ハザード関数は$\exp (\beta_k)$倍となる。
  
  
• Cox比例ハザードモデルは『比例ハザード性』の仮定の上に成立しており、このモデルを使用する際には$\log ( -\log S (t) )$を層別にプロットして平行であることを確認することで、その比例ハザード性が成立することを確認できる。