翻譯自：http://www.sthda.com/english/wiki/cox-proportional-hazards-model

Cox比例風(fēng)險(xiǎn)模型（考克斯，1972年）是常用的統(tǒng)計(jì)在醫(yī)學(xué)研究調(diào)查的患者和一個(gè)或多個(gè)預(yù)測(cè)變量的存活時(shí)間之間的關(guān)聯(lián)回歸模型。

在上一章 生存分析基礎(chǔ) 中，我們描述了生存分析的基本概念以及生存數(shù)據(jù)的分析和匯總方法，包括：

• 風(fēng)險(xiǎn)和生存函數(shù)的定義，
• 不同患者群體的Kaplan-Meier生存曲線的構(gòu)建
• 對(duì)數(shù)秩檢驗(yàn)（Log-Rank test），用于比較兩個(gè)或多個(gè)生存曲線

上述方法-Kaplan-Meier曲線和logrank檢驗(yàn)-是單變量分析的示例。他們根據(jù)調(diào)查中的一個(gè)因素描述了生存情況，但忽略了其他因素的影響。

此外，僅當(dāng)預(yù)測(cè)變量為分類變量時(shí)（例如：治療A與治療B；男性與女性），Kaplan-Meier曲線和對(duì)數(shù)秩檢驗(yàn)才有用。對(duì)于定量預(yù)測(cè)指標(biāo)（例如基因表達(dá)，體重或年齡），它們并不容易工作。

一種替代方法是Cox比例風(fēng)險(xiǎn)回歸分析，它既適用于定量預(yù)測(cè)變量也適用于類別變量。此外，Cox回歸模型擴(kuò)展了生存分析方法，可以同時(shí)評(píng)估幾種風(fēng)險(xiǎn)因素對(duì)生存時(shí)間的影響。

在本文中，我們將描述Cox回歸模型并提供使用R軟件的實(shí)際示例。

內(nèi)容

多元統(tǒng)計(jì)建模的需求

在臨床研究中，有許多情況，其中幾個(gè)已知量（稱為 協(xié)變量covariates）可能會(huì)影響患者的預(yù)后。

例如，假設(shè)比較了兩組患者：有和沒(méi)有特定基因型的患者。如果其中一組還包含較年長(zhǎng)的個(gè)體，則生存率的任何差異都可能歸因于基因型或年齡，或兩者都有。因此，在調(diào)查與任何一個(gè)因素相關(guān)的生存率時(shí)，通常需要針對(duì)其他因素的影響進(jìn)行調(diào)整。

統(tǒng)計(jì)模型是一種常用工具，可以同時(shí)分析多個(gè)因素的生存率。此外，統(tǒng)計(jì)模型還提供了每個(gè)因素的影響大小。

考克斯比例風(fēng)險(xiǎn)模型是用于對(duì)生存分析數(shù)據(jù)進(jìn)行建模的最重要方法之一。下一節(jié)介紹Cox回歸模型的基礎(chǔ)。

Cox比例風(fēng)險(xiǎn)模型的基礎(chǔ)

該模型的目的是同時(shí)評(píng)估幾個(gè)因素對(duì)生存的影響。換句話說(shuō)，它允許我們檢查特定因素如何影響特定時(shí)間點(diǎn)特定事件（例如，感染，死亡）的發(fā)生率。該比率通常稱為風(fēng)險(xiǎn)比率。預(yù)測(cè)變量（或因子）在生存分析文獻(xiàn)中通常稱為協(xié)變量 covariates。

Cox模型由 h(t) 表示的風(fēng)險(xiǎn)函數(shù)表示。簡(jiǎn)而言之，危險(xiǎn)函數(shù)可以解釋為在時(shí)間t死亡的風(fēng)險(xiǎn)?？梢怨烙?jì)如下：

其中：

• t 表示生存時(shí)間
• h(t) 是由一組p協(xié)變量（x1，x2，…，xp）確定的風(fēng)險(xiǎn)函數(shù)
• 系數(shù)（b1，b2，…，bp）衡量協(xié)變量的影響（即效應(yīng)大?。梢岳斫鉃轱L(fēng)險(xiǎn)權(quán)重）
• 術(shù)語(yǔ) h0 稱為基線危險(xiǎn)。如果所有的系數(shù)等于零（既exp(0)等于1），則對(duì)應(yīng)于風(fēng)險(xiǎn)的值。h(t)中的“t”提醒我們，危險(xiǎn)可能隨時(shí)間而變化。

Cox模型可以被寫為變量x(i)的危險(xiǎn)對(duì)數(shù)的多元線性回歸，而基線危險(xiǎn)是隨時(shí)間變化的“截距”項(xiàng)。

系數(shù) bi 稱為危險(xiǎn)比率（HR，hazard ratio）。bi 值大于零，或相當(dāng)于風(fēng)險(xiǎn)比率大于1，表明隨著第 i 個(gè)協(xié)變量值的增加，事件風(fēng)險(xiǎn)增加，因此生存時(shí)間縮短。

換句話說(shuō)，風(fēng)險(xiǎn)比大于1表示協(xié)變量與事件概率正相關(guān)，因此與存活時(shí)間負(fù)相關(guān)。
總之，
HR=1：無(wú)影響
HR<1：危害降低
HR>1：危險(xiǎn)增加

在癌癥研究中：

• 風(fēng)險(xiǎn)比 > 1（即：b > 0）的協(xié)變量稱為不良預(yù)后因素
• 風(fēng)險(xiǎn)比 < 1（即：b < 0）的協(xié)變量被稱為良好的預(yù)后因素

Cox模型的關(guān)鍵假設(shè)是觀察組（或患者）的危險(xiǎn)曲線應(yīng)成比例，并且不能交叉。

假設(shè)兩個(gè)x值不同的患者k和k'。相應(yīng)的風(fēng)險(xiǎn)函數(shù)可以簡(jiǎn)單地寫成如下：

• 患者k的風(fēng)險(xiǎn)函數(shù)：

  
  

• 患者k'的風(fēng)險(xiǎn)函數(shù)：

  
  

• 這兩名患者的危險(xiǎn)比 與時(shí)間t無(wú)關(guān)。

  
  

因此，Cox 模型是一個(gè)比例風(fēng)險(xiǎn)模型：任何一組事件的風(fēng)險(xiǎn)都是其他任何一組事件風(fēng)險(xiǎn)的常數(shù)倍。這一假設(shè)意味著，如上所述，各組的危險(xiǎn)曲線應(yīng)成比例，不能交叉。

換言之，如果一個(gè)人在某個(gè)初始時(shí)間點(diǎn)的死亡風(fēng)險(xiǎn)是另一個(gè)人的兩倍，那么在以后的任何時(shí)候，死亡風(fēng)險(xiǎn)仍然是另一個(gè)人的兩倍。

這種比例風(fēng)險(xiǎn)的假設(shè)應(yīng)該得到檢驗(yàn)。我們將在本系列的下一篇文章中討論評(píng)估比例性的方法：Cox模型假設(shè)。

在R中計(jì)算Cox模型

安裝并加載所需的R包

我們將使用兩個(gè)R包：

• survival 用于計(jì)算存活分析

• survminer 用于可視化生存分析結(jié)果

• 安裝軟件包

install.packages(c("survival", "survminer"))

• 加載包

library("survival")
library("survminer")

用于計(jì)算Cox模型的R函數(shù)：coxph（）

函數(shù)coxph（）[在survival包中]可用于計(jì)算R中的Cox比例風(fēng)險(xiǎn)回歸模型。

簡(jiǎn)化格式如下：

coxph(formula, data, method)

• formula：是以生存對(duì)象為響應(yīng)變量的線性模型。使用函數(shù)Surv（）創(chuàng)建生存對(duì)象，如下所示：Surv（time，event）。
• data：包含變量的數(shù)據(jù)框
• method：用于指定如何處理tie。默認(rèn)值為“efron”。其他選項(xiàng)有“breslow”和“exact”。默認(rèn)的“efron”通常優(yōu)先于曾經(jīng)流行的“breslow”方法。“exact”方法的計(jì)算量要大得多。

示例數(shù)據(jù)集

我們將在生存R數(shù)據(jù)包中使用肺癌數(shù)據(jù)。

data("lung")
head(lung)

  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

• 機(jī)構(gòu)：機(jī)構(gòu)代碼
• 時(shí)間：以天為單位的生存時(shí)間
• 狀態(tài)：審查狀態(tài)1 =審查，2 =失效
• 年齡：歲
• 性別：男= 1女= 2
• ph.ecog：ECOG成績(jī)得分（0 =好5 =死）
• ph.karno：醫(yī)師對(duì)Karnofsky績(jī)效評(píng)分（不良= 0-良好= 100）
• pat.karno：Karnofsky表現(xiàn)評(píng)分，按患者評(píng)分
• meal.cal：用餐時(shí)消耗的卡路里
• wt.loss：最近六個(gè)月的體重減輕

計(jì)算考克斯模型

我們將使用以下協(xié)變量來(lái)擬合Cox回歸：年齡，性別，ph.ecog和wt.loss。

我們首先為所有這些變量計(jì)算單變量Cox分析。然后我們將使用兩個(gè)變量來(lái)擬合多元Cox分析，以描述這些因素如何共同影響生存。

單變量Cox回歸

單變量Cox分析的計(jì)算公式如下：

res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox

Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)
      coef exp(coef) se(coef)     z      p
sex -0.531     0.588    0.167 -3.18 0.0015
Likelihood ratio test=10.6  on 1 df, p=0.00111
n= 228, number of events= 165 

Cox模型的功能摘要（）產(chǎn)生更完整的報(bào)告：

summary(res.cox)

Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)
  n= 228, number of events= 165 
       coef exp(coef) se(coef)      z Pr(>|z|)   
sex -0.5310    0.5880   0.1672 -3.176  0.00149 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    exp(coef) exp(-coef) lower .95 upper .95
sex     0.588      1.701    0.4237     0.816
Concordance= 0.579  (se = 0.022 )
Rsquare= 0.046   (max possible= 0.999 )
Likelihood ratio test= 10.63  on 1 df,   p=0.001111
Wald test            = 10.09  on 1 df,   p=0.001491
Score (logrank) test = 10.33  on 1 df,   p=0.001312

Cox回歸結(jié)果可以解釋為：

1. 統(tǒng)計(jì)意義。標(biāo)記為“ z”的列給出了Wald統(tǒng)計(jì)值。它對(duì)應(yīng)于每個(gè)回歸系數(shù)與其標(biāo)準(zhǔn)誤差的比率（z = coef / se（coef））。Wald統(tǒng)計(jì)量會(huì)評(píng)估beta給定變量的系數(shù)在統(tǒng)計(jì)上顯著不同于0。從上面的輸出中，我們可以得出結(jié)論，變量性別具有很高的統(tǒng)計(jì)學(xué)上顯著的系數(shù)。

2. 回歸系數(shù)。Cox模型結(jié)果中要注意的第二個(gè)特征是回歸系數(shù)（coef）的符號(hào)。正號(hào)表示該變量值較高的受試者的危險(xiǎn)（死亡風(fēng)險(xiǎn)）較高，因此預(yù)后較差。變量sex被編碼為數(shù)字向量。1：男性，2：女性。Cox模型的 R summary 提供了第二組相對(duì)于第一組（即女性與男性）的危險(xiǎn)比（HR）。在這些數(shù)據(jù)中，性別的beta系數(shù)= -0.53表示女性比男性具有較低的死亡風(fēng)險(xiǎn)（較低的生存率）。

3. 危險(xiǎn)率。指數(shù)系數(shù)（exp（coef）= exp（-0.53）= 0.59），也稱為危險(xiǎn)比，給出了協(xié)變量的影響大小。例如，女性（性別= 2）可以將危險(xiǎn)降低0.59或41％。女性是有良好預(yù)后的。

4. 危險(xiǎn)比的置信區(qū)間。摘要輸出還給出了危險(xiǎn)比（exp（coef））的上限和下限95％置信區(qū)間，下限95％= 0.4237，上限95％= 0.816。

5. 該模型的全局統(tǒng)計(jì)意義。最后，輸出為模型的總體重要性提供了三個(gè)備選檢驗(yàn)的p值：似然比檢驗(yàn)，Wald檢驗(yàn)和得分對(duì)數(shù)秩統(tǒng)計(jì)。這三種方法漸近等效。對(duì)于足夠大的N，它們將給出相似的結(jié)果。對(duì)于較小的N，它們可能會(huì)有所不同。對(duì)于小樣本量，似然比測(cè)試具有更好的行為，因此通常是首選。

要將單變量coxph函數(shù)一次應(yīng)用于多個(gè)協(xié)變量，請(qǐng)輸入以下命令：

covariates <- c("age", "sex",  "ph.karno", "ph.ecog", "wt.loss")
univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(time, status)~', x)))

univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
# Extract data 
univ_results <- lapply(univ_models,
                       function(x){ 
                          x <- summary(x)
                          p.value<-signif(x$wald["pvalue"], digits=2)
                          wald.test<-signif(x$wald["test"], digits=2)
                          beta<-signif(x$coef[1], digits=2);#coeficient beta
                          HR <-signif(x$coef[2], digits=2);#exp(beta)
                          HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
                          HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
                          HR <- paste0(HR, " (", 
                                       HR.confint.lower, "-", HR.confint.upper, ")")
                          res<-c(beta, HR, wald.test, p.value)
                          names(res)<-c("beta", "HR (95% CI for HR)", "wald.test", 
                                        "p.value")
                          return(res)
                          #return(exp(cbind(coef(x),confint(x))))
                         })
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)

           beta HR (95% CI for HR) wald.test p.value
age       0.019            1 (1-1)       4.1   0.042
sex       -0.53   0.59 (0.42-0.82)        10  0.0015
ph.karno -0.016      0.98 (0.97-1)       7.9   0.005
ph.ecog    0.48        1.6 (1.3-2)        18 2.7e-05
wt.loss  0.0013         1 (0.99-1)      0.05    0.83

上面的輸出顯示了每個(gè)變量相對(duì)于總生存率的回歸beta系數(shù)，效應(yīng)大小（以危險(xiǎn)比給出）和統(tǒng)計(jì)顯著性。通過(guò)單獨(dú)的單變量Cox回歸評(píng)估每個(gè)因素。

從上面的輸出中，

• 性別，年齡和ph.ecog變量具有極高的統(tǒng)計(jì)顯著性系數(shù)，而wt.loss的系數(shù)則不顯著。

• 年齡和ph.ecog的β系數(shù)為正，而性別的系數(shù)為負(fù)。因此，年齡更大和ph.ecog較高與較差的生存率有關(guān)，而女性（性別= 2）則與較好的生存率有關(guān)。

現(xiàn)在，我們要描述這些因素如何共同影響生存。為了回答這個(gè)問(wèn)題，我們將執(zhí)行多元Cox回歸分析。由于變量ph.karno在單變量Cox分析中不重要，因此在多變量分析中將其跳過(guò)。我們將3個(gè)因素（性別，年齡和ph.ecog）納入多元模型。

多元Cox回歸分析

時(shí)間常數(shù)協(xié)變量的死亡時(shí)間的Cox回歸指定如下：

res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)
summary(res.cox)

Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
  n= 227, number of events= 164 
   (1 observation deleted due to missingness)
             coef exp(coef)  se(coef)      z Pr(>|z|)    
age      0.011067  1.011128  0.009267  1.194 0.232416    
sex     -0.552612  0.575445  0.167739 -3.294 0.000986 ***
ph.ecog  0.463728  1.589991  0.113577  4.083 4.45e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
        exp(coef) exp(-coef) lower .95 upper .95
age        1.0111     0.9890    0.9929    1.0297
sex        0.5754     1.7378    0.4142    0.7994
ph.ecog    1.5900     0.6289    1.2727    1.9864
Concordance= 0.637  (se = 0.026 )
Rsquare= 0.126   (max possible= 0.999 )
Likelihood ratio test= 30.5  on 3 df,   p=1.083e-06
Wald test            = 29.93  on 3 df,   p=1.428e-06
Score (logrank) test = 30.5  on 3 df,   p=1.083e-06

所有三個(gè)總體測(cè)試（似然性，Wald和得分）的p值均顯著，表明該模型具有顯著性。這些測(cè)試評(píng)估了所有beta的綜合零假設(shè)為0。在上面的示例中，檢驗(yàn)統(tǒng)計(jì)量非常一致，并且完全拒絕了綜合零假設(shè)。

在多變量Cox分析中，協(xié)變量性別和ph.ecog保持顯著性（p <0.05）。但是，協(xié)變量年齡不顯著（p = 0.23，大于0.05）。

性別的p值為0.000986，危險(xiǎn)比HR = exp（coef）= 0.58，表明患者的性別與死亡風(fēng)險(xiǎn)降低之間有很強(qiáng)的關(guān)系。協(xié)變量的危險(xiǎn)比可解釋為對(duì)危險(xiǎn)的倍增效應(yīng)。例如，保持其他協(xié)變量不變（女性（性別= 2））可將危險(xiǎn)降低0.58或42％。我們得出結(jié)論，成為女性與良好的預(yù)后相關(guān)。

同樣，ph.ecog的p值為4.45e-05，危險(xiǎn)比HR = 1.59，表明ph.ecog值與死亡風(fēng)險(xiǎn)增加之間有很強(qiáng)的關(guān)系。保持其他協(xié)變量不變，ph.ecog的值越高，生存率越低。

相比之下，年齡的p值現(xiàn)在為p = 0.23。危險(xiǎn)比HR = exp（coef）= 1.01，95％置信區(qū)間為0.99至1.03。由于HR的置信區(qū)間為1，因此這些結(jié)果表明，在調(diào)整phog值和患者的性別之后，年齡對(duì)HR差異的貢獻(xiàn)較小，并且僅趨于顯著。例如，在其他協(xié)變量保持不變的情況下，再增加一歲會(huì)引起每日死亡危險(xiǎn)，其系數(shù)為expβ= 1.01或1％，這并不是一個(gè)重要的貢獻(xiàn)。

可視化生存時(shí)間的估計(jì)分布

將Cox模型擬合到數(shù)據(jù)后，就可以可視化特定風(fēng)險(xiǎn)組在任何給定時(shí)間點(diǎn)的預(yù)測(cè)生存率。函數(shù)survfit（）估計(jì)生存比例，默認(rèn)情況下為協(xié)變量的平均值。

# Plot the baseline survival function
ggsurvplot(survfit(res.cox), color = "#2E9FDF",
           ggtheme = theme_minimal())

考克斯比例風(fēng)險(xiǎn)模型

我們不妨展示估計(jì)的生存率如何取決于目標(biāo)協(xié)變量的值。

考慮到這一點(diǎn)，我們想評(píng)估性別對(duì)估計(jì)生存率的影響。在這種情況下，我們用兩行構(gòu)造一個(gè)新的數(shù)據(jù)幀，每一行代表性別。其他協(xié)變量固定為其平均值（如果是連續(xù)變量）或最低水平（如果它們是離散變量）。對(duì)于偽協(xié)變量，平均值為數(shù)據(jù)集中編碼為1的比例。該數(shù)據(jù)幀通過(guò)newdata參數(shù)傳遞給survfit（）：

# Create the new data  
sex_df <- with(lung,
               data.frame(sex = c(1, 2), 
                          age = rep(mean(age, na.rm = TRUE), 2),
                          ph.ecog = c(1, 1)
                          )
               )
sex_df

  sex      age ph.ecog
1   1 62.44737       1
2   2 62.44737       1

# Survival curves
fit <- survfit(res.cox, newdata = sex_df)
ggsurvplot(fit, conf.int = TRUE, legend.labs=c("Sex=1", "Sex=2"),
           ggtheme = theme_minimal())

概要

在本文中，我們描述了Cox回歸模型，用于同時(shí)評(píng)估多種風(fēng)險(xiǎn)因素與患者生存時(shí)間之間的關(guān)系。我們演示了如何使用生存包計(jì)算Cox模型。此外，我們描述了如何使用survminer軟件包來(lái)可視化分析結(jié)果。

參考文獻(xiàn)

• 考克斯博士（1972）?；貧w模型和壽命表（有討論）。JR Statist Soc B 34：187–220
• MJ Bradburn，TG Clark，SB Love和DG Altman。生存分析第二部分：多元數(shù)據(jù)分析–概念和方法介紹。英國(guó)癌癥雜志（2003）89，431 – 436