Rで回帰分析:COX比例ハザードモデル

Rでデータサイエンス

COX比例ハザードモデル

生存関数\(S(t)\)

  • 被験者、機器、またはその他の対象物が特定の時間を超えて生存する確率を与える関数。
  • \(T\,\)を区間\(\,[0,\infty]\,\)上の累積分布関数\(\,F(t)\,\)を持つ連続確率変数としたとき、その生存関数(または信頼性関数)は式(1)の通り。
  • \(\,f(u)\,\)は時点\(\,u\,\)における死亡率(=時点\(\,t\,\)における死亡者数/サンプルサイズ)を表す関数であり、\(\,F(t)\,\)\(\,f(u)\,\)を積分した累積死亡率の時間的な変化を表す関数である。

\[ \begin{equation} S(t)=P\left(T>t\right)=\int_t^{\infty}f(u)\,du=1-F(t)\tag{1} \end{equation} \]

ハザード関数\(\lambda(t)\)

  • 被験者が時点\(\,t\,\)まで生存したという条件のもとで、その時間に死亡する確率。
  • 時点\(\,t\,\)における死亡者数(=\(\,f(t)\times\,\)サンプルサイズ)を同時点における生存者数(=\(\,S(t)\times\,\)サンプルサイズ)で除した瞬間死亡率。

\[ \begin{equation} \lambda(t)=\dfrac{f(t)\times N}{S(t)\times N}=\dfrac{f(t)}{S(t)}\tag{2} \end{equation} \]

生存関数とハザード関数

式(2)を変形すると、

\[ \begin{align} \lambda(t)&=\dfrac{f(t)}{S(t)}=\dfrac{1}{S(t)}\cdot\dfrac{d}{dt}F(t)\\&=\dfrac{1}{S(t)}\cdot\dfrac{d}{dt}(1-S(t))\\&=\dfrac{1}{S(t)}\cdot-\,\dfrac{d}{dt}S(t)\tag{3} \end{align} \]

式(3)の両辺を積分して、

\[ \begin{align} \int \lambda(t)dt & =\int\dfrac{1}{S(t)}\cdot\left(-\,\dfrac{d}{dt}S(t)\right)dt\\ & =-\int\dfrac{1}{S(t)}dS(t)\\ & =-\,\textrm{ln}(S(t))+C\tag{4} \end{align} \]

よって

\[ \begin{equation} S(t)=\textrm{exp}\left(C-\int\lambda(t)\,dt\right)\tag{5} \end{equation} \]

さらに

\[ \begin{equation} S(0)=\exp(C-0)=\textrm{exp}(C)=1\tag{6} \end{equation} \]

が成り立つには\(\,C=0\,\)である必要があるため、式(5)に\(\,C=0\,\)を代入して

\[ \begin{equation} S(t)=\textrm{exp}\left(-\,\int\lambda(t)dt\right)\tag{7} \end{equation} \]

となり

\[ \begin{equation} \textrm{ln}(S(t))=-\,\int\lambda(t)dt\tag{8} \end{equation} \]

となる。

つまり、生存関数\(\,S(t)\,\)はハザード関数\(\,\lambda(t)\,\)によって決定される。

そこで、ハザードを目的変数とし、生存率に影響する因子を説明変数とする重回帰分析により、説明変数が生存率に与える影響を求めることが出来る。

しかし、ハザード関数\(\,\lambda(t)\,\)を具体的に規定するのは困難であるため、ハザード関数と基準ハザード関数の比の自然対数を取った対数ハザード比を目的変数とした一般化線形モデル(比例ハザードモデル)を次節の通りに想定する。

なお、基準ハザード関数\(\,\lambda_0(t)\,\)とは全ての説明変数の値が\(\,0\,\)のときのハザード関数である。

比例ハザードモデル

対数ハザード比を目的変数とし、説明変数との間の線形関係を仮定して、リンク関数を式(9)の通りとした一般化線形モデルを考える。

\[ \begin{equation} \eta=\textrm{ln}\left(\dfrac{\lambda\left(t|x_1,\cdots,x_p\right)}{\lambda_0(t)}\right)=\beta_0+\beta_1x_1+\cdots+\beta_px_p+\epsilon\tag{9} \end{equation} \]

あるいは

\[ \begin{align} \eta_i&=\textrm{ln}\left(\dfrac{\lambda\left(t|x_i\right)}{\lambda_0(t)}\right)\\&=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}+\epsilon_i=\hat{\eta_i}+\epsilon_i=y_i+\epsilon_i\tag{10} \end{align} \]

行列式で表現すると

\[ \begin{equation} \pmb{\eta}=\textrm{ln}\left(\dfrac{\lambda(t)}{\lambda_0(t)}\right)=\textbf{x}\,\pmb{\beta}+\pmb{\epsilon}=\textbf{y}+\pmb{\epsilon}\tag{11} \end{equation} \]

となり、\(\,\textbf{y}\,\)

\[ \begin{equation} \begin{bmatrix}y_1\\\vdots\\y_n\end{bmatrix}=\begin{bmatrix}1&x_{11}&\cdots&x_{1p}\\\vdots&\vdots&&\vdots\\1&x_{n1}&\cdots&x_{np}\end{bmatrix}\begin{bmatrix}\beta_0\\\vdots\\\beta_p\end{bmatrix}\tag{12} \end{equation} \]

ここで、ハザード関数\(\,\lambda(t)\,\)\(\,t\,\)とは無関係に一定であるとの仮定を置くと、

\[ \begin{equation} \lambda(t)=\lambda\tag{13} \end{equation} \]

式(7)は

\[ \begin{equation} S(t)=\textrm{exp}\left(-\,\int\lambda(t)dt\right)=\textrm{exp}(\lambda t)\tag{14} \end{equation} \]

となり、生存関数,\(S(t)\,\)は指数関数となって、理論的な生存時間は\(\,\infty\,\)となるが、\(\,S(t)\,\)が微小量\(\,\delta\,\)になったときに死亡すると仮定すると、生存時間とハザードモデルの間には以下の関係が成り立つ。

\[ \begin{equation} S(t)=\exp(-\lambda t)=\delta\tag{15} \end{equation} \]

\[ \begin{equation} \textrm{ln}(\delta)=-\lambda t\tag{16} \end{equation} \]

\[ \begin{equation} \lambda=-\dfrac{\textrm{ln}(\delta)}{t}\tag{17} \end{equation} \]

また、基準ハザード関数\(\,\lambda_0(t)\,\)

\[ \begin{equation} \lambda_0=-\dfrac{\textrm{ln}(\delta)}{t_0}\tag{18} \end{equation} \]

と表すこととし、式(17)および式(18)を式(11)に代入すると

\[ \begin{align} \eta&=\textrm{ln}\left(\dfrac{\lambda(t)}{\lambda_0(t)}\right) \\&=\textrm{ln}\left(\dfrac{-\dfrac{\textrm{ln}(\delta)}{t}}{-\dfrac{\textrm{ln}(\delta)}{t_0}}\right) \\&=\textrm{ln}\left(\dfrac{t_0}{t}\right) \\&=\textrm{ln}(t_0)-\textrm{ln}(t) \\&=\pmb{\beta}\textbf{x}+\pmb{\epsilon}\tag{19} \end{align} \]

よって、

\[ \begin{equation} \textrm{ln}(t)=\textrm{ln}(t_0)-\left(\pmb{\beta}\textbf{x}+\pmb{\epsilon}\right)\tag{20} \end{equation} \]

比例ハザードモデルは式(20)を解くことになる。

  • そのため推定された\(\,\pmb{\beta}\,\)の符号は反対になる。
  • 偏回帰係数のみを推定する。

サンプル

# イベントは全て発生とする(event = 1)。
# 誤差は正規分布に従うとする。
library(dplyr)
seed <- 20230707
set.seed(seed = seed)
n <- 200
b0 <- 2
b <- c(4, -2.5)
x1 <- runif(n = n)
x2 <- sample(x = c(1, 0), size = n, replace = T)
X <- c(x1, x2) %>%
    matrix(ncol = 2)
t <- exp(X %*% b + b0 + rnorm(n = n))
df <- data.frame(x1, x2, t, event = rep(1, n))
glimpse(df)
Rows: 200
Columns: 4
$ x1    <dbl> 0.81400512, 0.20429657, 0.02972173, 0.67730553, 0.01749885, 0.08…
$ x2    <dbl> 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1…
$ t     <dbl> 42.6773058, 10.7353750, 9.9645473, 3.0172411, 0.5855208, 0.26208…
$ event <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
library(survival)
result <- coxph(formula = Surv(time = t, event = event) ~ x1 + x2, data = df)
# 偏回帰係数
result %>%
    summary()
# 式(20)を解いているため、算出した偏回帰係数の符号は上記 b とは反転している。
Call:
coxph(formula = Surv(time = t, event = event) ~ x1 + x2, data = df)

  n= 200, number of events= 200 

       coef exp(coef) se(coef)      z            Pr(>|z|)    
x1 -3.76690   0.02312  0.30227 -12.46 <0.0000000000000002 ***
x2  2.43885  11.45991  0.20368  11.97 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

   exp(coef) exp(-coef) lower .95 upper .95
x1   0.02312   43.24586   0.01279   0.04182
x2  11.45991    0.08726   7.68787  17.08270

Concordance= 0.835  (se = 0.012 )
Likelihood ratio test= 235.8  on 2 df,   p=<0.0000000000000002
Wald test            = 202.2  on 2 df,   p=<0.0000000000000002
Score (logrank) test = 223.2  on 2 df,   p=<0.0000000000000002
# 定数項
mean(log(t)) - (mean(x1) * (-1 * coefficients(result)[1]) + mean(x2) * (-1 * coefficients(result)[2]))
      x1 
2.065826 
# 偏回帰係数の95%信頼区間
confint(object = result)
       2.5 %    97.5 %
x1 -4.359341 -3.174462
x2  2.039644  2.838066

最終更新

Sys.time()
[1] "2024-04-03 04:39:50 JST"

R、Quarto、Package

R.Version()$version.string
[1] "R version 4.3.3 (2024-02-29 ucrt)"
quarto::quarto_version()
[1] '1.4.542'
packageVersion(pkg = "tidyverse")
[1] '2.0.0'

著者