Cox 迴歸分析 (Cox Regression)

在進行存活分析時，可以透過 Log-rank test 來比較組別之間的存活曲線是否有顯著差異。但當自變項是連續變項或自變項超過2個以上時的分析時，就需要透過 Cox Regression (Cox Proportional Hazard Model) 來達成。Cox Regression Model 是使用有母數分析的預測方式，用以分析影響死亡率的變數(或影響存活率的重要因子)。在討論 cox regression 時有二個部分需要先釐清，即 存活函式 與 風險函式:

存活函式 (Survival Function):

目的是呈現在某特定時間點下，個案可以活過此特定時間點的機率為何。

\[ S(t) = P(T > t) = \int_t^{\infty}f(u)du \]

其中 \(t\) 為時間點，\(S(0)=1\), \(S(\infty)=0\)。

風險函式 (Hazard Function):

\[ h(t) = \frac{\frac{-d[S(t)]}{dt}}{S(t)} = \frac{f(t)}{S(t)} \\ S(t) = exp(-\int_0^th(u)du) \]

其中 \(t\) 為時間點，風險函式 \(h(t)\) 值與存活時間有關，此代表在某一時間點下，事件機率(死亡)除以生存函式的值。若時間越久 (\(t\) 越大)，存活函式 \(S(t)\) 越小 (存活函式為遞減函式)，假設死亡機率密度 \(f(t)\) 對任何時間點都一樣，則風險函式 \(h(t)\) 值越大，預期存活時間短。

在經典的 Cox Regression 的等比例風險假設(Proportional hazard assumption, PH assumption)假設下，針對某一危險因子而言，其風險比，不能隨著時間而有所改變，必須要固定。。

由上可之，推導 Cox Regression 可得：

\[ h(t) = \lambda \\ S(t) = exp(-\int_0^th(u)du) = exp(-\lambda t) = e^{-\lambda t} \\ f(t) = h(t) * S(t) = \lambda * e^{-\lambda t} \]

而因 \(\lambda\) 為一定值，故可以推出 Cox Regression 為:

\[ log(HR(x))=log(\frac{h(t|x)}{h_0(t)}) = \beta_1{x_1} + \beta_2{x_2} + ... + \beta_k{x_k} + \varepsilon \]

此即為迴歸方程式。Hazard Ratio (\(HR\)) 表示某個時間下會發生事件的風險比，而 \(HR(x)\) 表示在給定 x 下會發生事件的風險比，x 即是自變項的數值，如年齡90歲便是一個數值 x。\(h(t|x)\) 表示在第 \(t\) 個時間點時，給定 x 值時的風險。\(h_0(t)\) 表示在第 \(t\) 個時間點時的基礎風險。而 Hazard Ratio (\(HR\)) 計算範例為當僅有一變項 \(x_1\)，\(x_1=1\) 表示為男性，\(x_1=0\) 表示為女性，則 \(exp(\beta_1)\) 為男性相對於女性的危險比率，且 \(HR=\frac{exp(\beta_0 + \beta_1)}{exp(\beta_0)}\)。

Cox Regression 不僅可以用來分析病因與死亡率之間的關係，更可以延伸至企業倒閉風險分析等議題上。

套件安裝

於 R 實作中，可以透過套件 survival 來實現存活分析。

packageName <- c("survival")
for(i in 1:length(packageName)) {
  if(!(packageName[i] %in% rownames(installed.packages()))) {
    install.packages(packageName[i])
  }
}
lapply(packageName, require, character.only = TRUE)

## Loading required package: survival

## [[1]]
## [1] TRUE

使用函式介紹

函式名稱	函式用途
Surv	創造一個存活分析的物件，通常在迴歸公式中代表應變數(response variable)
coxph	Fit Proportional Hazards Regression Model
coxph.control	控制 coxph 擬合的輔助參數
cox.zph	測試迴歸模型中的比例風險假設
cox.detail	Cox 迴歸模型的細節，此函式會回傳在每個獨立時間下的 1 階, 2 階的導數矩陣
strata	用來分層分析變數

資料準備

使用套件 survival 中內建的資料 ovarian，此為子宮癌生存時間的資料集。

data("ovarian")
head(ovarian, 10)

##    futime fustat     age resid.ds rx ecog.ps
## 1      59      1 72.3315        2  1       1
## 2     115      1 74.4932        2  1       1
## 3     156      1 66.4658        2  1       2
## 4     421      0 53.3644        2  2       1
## 5     431      1 50.3397        2  1       1
## 6     448      0 56.4301        1  1       2
## 7     464      1 56.9370        2  2       2
## 8     475      1 59.8548        2  2       2
## 9     477      0 64.1753        2  1       1
## 10    563      1 55.1781        1  2       2

其中 futime 為生存或是資料不再更新(如刪除，不再追蹤)的時間，fustat 為是否有持續追蹤的指標(censoring status)，age 為年齡，resid.ds 是否目前仍有疾病 (1為無,2為有)，rx 為治療或控制組(1 為控制組 placebo, 2 為治療組 thiopeta)，ecog.ps 為美國東岸癌症臨床研究合作組織定義之日常體能狀態(1 為較佳，分數越高越差)。

迴歸分析

# the prototype of Surv
# time 為與時間並進的變數
# event 為指標，需考慮 type 為何；預設為 0 為存活，1 為死亡
# type 為 censoring 的型態，此會與 event 之值有關
Surv(time, time2, event,
    type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'))

# the prototype of coxph
# formula: 需注意先將應變數透過函式(Surv)轉為物件
# init: 反覆計算中使用的初始值向量，預設皆為 0
# control: 其他用於反覆計算之參數選項，可接受 coxph.control 定義之物件
# ties: 預設為 Efron，處理死亡時間更精確
# singular.ok: 如何處理模型矩陣中的共線性
# robust: 若為真，將可得到一穩定的變異(Variance)估計
coxph(formula, data, weights, subset, 
      na.action, init, control, 
      ties=c("efron","breslow","exact"), 
      singular.ok=TRUE, robust=FALSE, 
      model=FALSE, x=FALSE, y=TRUE, tt, ...)

# the prototype of coxph.control
coxph.control(eps = 1e-09, toler.chol = .Machine$double.eps^0.75,
      iter.max = 20, toler.inf = sqrt(eps), outer.max = 10, timefix=TRUE)

ova.fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps, data = ovarian)
ova.fit

## Call:
## coxph(formula = Surv(futime, fustat) ~ age + ecog.ps, data = ovarian)
## 
##           coef exp(coef) se(coef)    z      p
## age     0.1615    1.1753   0.0499 3.24 0.0012
## ecog.ps 0.0187    1.0188   0.5991 0.03 0.9751
## 
## Likelihood ratio test=14.3  on 2 df, p=0.000787
## n= 26, number of events= 12

Cox regression 結果中顯示，coef 為變數 age 或 ecog.ps 的係數，exp(coef) 是一個生存比率，se(coef) 為係數的標準差，而 Z 表示 Z-score，P 為 P-value，可以透過 Z-score 來推求出 P-value。

驗證風險假設

# the prototype of cox.zph
# transform: 指定在驗證假設前將存活時間進行轉換的方法
# global: 除了單變數檢驗外，是否需要進行全域卡方檢驗
cox.zph(fit.model, transform="km", global=TRUE)

ova.test <- cox.zph(ova.fit, transform="km", global=TRUE)
print(ova.test)

##            rho chisq     p
## age     -0.243 0.856 0.355
## ecog.ps  0.520 2.545 0.111
## GLOBAL      NA 3.195 0.202

由上可以看出，以 age 與 ecog.ps 二變數的 Cox 模型，其模型擬合結果並不顯著。

繪迴歸曲線結果

plot(ova.test)

迴歸模型細節

透過 coxph.detail 回傳許多計算細節，將有利於導入新的方法。

coxph.detail(ova.fit)

## $time
##  [1]  59 115 156 268 329 353 365 431 464 475 563 638
## 
## $means
##          age  ecog.ps
## 59  69.09819 1.455425
## 115 68.46900 1.544049
## 156 66.17359 1.751348
## 268 66.13888 1.721811
## 329 59.69610 1.507553
## 353 59.77545 1.509985
## 365 59.28316 1.439937
## 431 58.40780 1.596546
## 464 58.75868 1.583439
## 475 58.91087 1.548638
## 563 56.39392 1.688778
## 638 56.53567 1.652492
## 
## $nevent
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $nrisk
##  [1] 26 25 24 23 22 21 20 17 15 14 12 11
## 
## $hazard
##  [1] 0.01207211 0.01442131 0.01991626 0.02228211 0.03944350 0.03963251
##  [7] 0.04529802 0.06142322 0.06738846 0.07301832 0.12124602 0.13538256
## 
## $score
##             age    ecog.ps
## 59    3.2333105 -0.4554246
## 115   6.0242037 -0.5440489
## 156   0.2922075  0.2486520
## 268   8.3652187  0.2781892
## 329 -16.5590971 -0.5075530
## 353   3.4437534  0.4900149
## 365   5.1415416 -0.4399369
## 431  -8.0681010 -0.5965458
## 464  -1.8216798  0.4165610
## 475   0.9439309  0.4513619
## 563  -1.2158159  0.3112218
## 638   0.2205275  0.3475083
## 
## $imat
## , , 59
## 
##                age    ecog.ps
## age     49.5808416 -0.4278646
## ecog.ps -0.4278646  0.2480130
## 
## , , 115
## 
##               age    ecog.ps
## age     56.798874 -0.1688140
## ecog.ps -0.168814  0.2480597
## 
## , , 156
## 
##              age   ecog.ps
## age     59.34410 1.4915098
## ecog.ps  1.49151 0.1868242
## 
## , , 268
## 
##               age  ecog.ps
## age     66.382212 1.659029
## ecog.ps  1.659029 0.200800
## 
## , , 329
## 
##                age    ecog.ps
## age     22.1041739 -0.2359377
## ecog.ps -0.2359377  0.2499430
## 
## , , 353
## 
##                age    ecog.ps
## age     20.8898416 -0.2775354
## ecog.ps -0.2775354  0.2499003
## 
## , , 365
## 
##                age    ecog.ps
## age     21.9383999 -0.5929218
## ecog.ps -0.5929218  0.2463924
## 
## , , 431
## 
##                age    ecog.ps
## age     21.4326674 -0.2817995
## ecog.ps -0.2817995  0.2406789
## 
## , , 464
## 
##                age    ecog.ps
## age     21.4157990 -0.3733844
## ecog.ps -0.3733844  0.2430379
## 
## , , 475
## 
##                age    ecog.ps
## age     22.9045464 -0.3358857
## ecog.ps -0.3358857  0.2476343
## 
## , , 563
## 
##                age   ecog.ps
## age     19.3536405 0.9661201
## ecog.ps  0.9661201 0.2143628
## 
## , , 638
## 
##               age   ecog.ps
## age     21.417710 1.1280254
## ecog.ps  1.128025 0.2267463
## 
## 
## $varhaz
##  [1] 0.0001457358 0.0002079742 0.0003966576 0.0004964923 0.0015557897
##  [6] 0.0015707358 0.0020519106 0.0037728117 0.0045412045 0.0053316755
## [11] 0.0147005982 0.0183284376
## 
## $y
##       time status
## 1  -1   59      1
## 2  -1  115      1
## 3  -1  156      1
## 22 -1  268      1
## 23 -1  329      1
## 24 -1  353      1
## 25 -1  365      1
## 26 -1  377      0
## 4  -1  421      0
## 5  -1  431      1
## 6  -1  448      0
## 7  -1  464      1
## 8  -1  475      1
## 9  -1  477      0
## 10 -1  563      1
## 11 -1  638      1
## 12 -1  744      0
## 13 -1  769      0
## 14 -1  770      0
## 15 -1  803      0
## 16 -1  855      0
## 17 -1 1040      0
## 18 -1 1106      0
## 19 -1 1129      0
## 20 -1 1206      0
## 21 -1 1227      0
## 
## $x
##        age ecog.ps
## 1  72.3315       1
## 2  74.4932       1
## 3  66.4658       2
## 22 74.5041       2
## 23 43.1370       1
## 24 63.2192       2
## 25 64.4247       1
## 26 58.3096       1
## 4  53.3644       1
## 5  50.3397       1
## 6  56.4301       2
## 7  56.9370       2
## 8  59.8548       2
## 9  64.1753       1
## 10 55.1781       2
## 11 56.7562       2
## 12 50.1096       1
## 13 59.6301       2
## 14 57.0521       1
## 15 39.2712       1
## 16 43.1233       2
## 17 38.8932       2
## 18 44.6000       1
## 19 53.9068       1
## 20 44.2055       1
## 21 59.5890       2

分層分析

經典的存活分析建立在等比例風險假設(Proportional hazard assumption, PH assumption)。但當風險因子會隨著時間改變時，原本的 Cox 迴歸模型就需要修正，其中最常使用的便是分層 Cox 迴歸模型。分層 Cox 迴歸模型將這些不符合等比例風險假設的因子當作分層變數，表示不同層的基礎危險函式(baseline hazard function)是不同的，但需注意分層數若太多可能會降低統計檢定力(statistical power)。在 R 實作中可透過函式 strata 針對某一變數進行分層分析。

strata.data <- data.frame(
  time = c(4,3,1,1,2,2,3),
  status = c(1,1,1,0,1,1,0),
  x = c(0,2,1,1,1,0,0),
  sex = c(0,0,0,0,1,1,1)
)
summary(coxph(Surv(time,status) ~ x + strata(sex), data=strata.data))

## Call:
## coxph(formula = Surv(time, status) ~ x + strata(sex), data = strata.data)
## 
##   n= 7, number of events= 5 
## 
##     coef exp(coef) se(coef)     z Pr(>|z|)
## x 0.8023    2.2307   0.8224 0.976    0.329
## 
##   exp(coef) exp(-coef) lower .95 upper .95
## x     2.231     0.4483    0.4451     11.18
## 
## Concordance= 0.667  (se = 0.362 )
## Rsquare= 0.144   (max possible= 0.669 )
## Likelihood ratio test= 1.09  on 1 df,   p=0.2971
## Wald test            = 0.95  on 1 df,   p=0.3293
## Score (logrank) test = 1.05  on 1 df,   p=0.3053