翻譯自：http://www.sthda.com/english/wiki/cox-proportional-hazards-model
Cox比例風(fēng)險模型（考克斯鸟廓，1972年）是常用的統(tǒng)計在醫(yī)學(xué)研究調(diào)查的患者和一個或多個預(yù)測變量的存活時間之間的關(guān)聯(lián)回歸模型技竟。
在上一章 生存分析基礎(chǔ) 中避归，我們描述了生存分析的基本概念以及生存數(shù)據(jù)的分析和匯總方法鼓鲁，包括：
- 風(fēng)險和生存函數(shù)的定義，
- 不同患者群體的Kaplan-Meier生存曲線的構(gòu)建
- 對數(shù)秩檢驗（Log-Rank test）创肥，用于比較兩個或多個生存曲線
上述方法-Kaplan-Meier曲線和logrank檢驗-是單變量分析的示例矿酵。他們根據(jù)調(diào)查中的一個因素描述了生存情況浪箭，但忽略了其他因素的影響婴梧。
此外下梢，僅當(dāng)預(yù)測變量為分類變量時（例如：治療A與治療B；男性與女性）塞蹭，Kaplan-Meier曲線和對數(shù)秩檢驗才有用孽江。對于定量預(yù)測指標(biāo)（例如基因表達，體重或年齡）番电，它們并不容易工作岗屏。
一種替代方法是Cox比例風(fēng)險回歸分析，它既適用于定量預(yù)測變量也適用于類別變量。此外这刷，Cox回歸模型擴展了生存分析方法涎跨，可以同時評估幾種風(fēng)險因素對生存時間的影響。
在本文中崭歧，我們將描述Cox回歸模型并提供使用R軟件的實際示例。
內(nèi)容
多元統(tǒng)計建模的需求
在臨床研究中撞牢，有許多情況率碾，其中幾個已知量（稱為 協(xié)變量covariates）可能會影響患者的預(yù)后。
例如屋彪，假設(shè)比較了兩組患者：有和沒有特定基因型的患者所宰。如果其中一組還包含較年長的個體，則生存率的任何差異都可能歸因于基因型或年齡畜挥，或兩者都有仔粥。因此，在調(diào)查與任何一個因素相關(guān)的生存率時蟹但，通常需要針對其他因素的影響進行調(diào)整躯泰。
統(tǒng)計模型是一種常用工具，可以同時分析多個因素的生存率华糖。此外麦向，統(tǒng)計模型還提供了每個因素的影響大小。
考克斯比例風(fēng)險模型是用于對生存分析數(shù)據(jù)進行建模的最重要方法之一客叉。下一節(jié)介紹Cox回歸模型的基礎(chǔ)诵竭。
Cox比例風(fēng)險模型的基礎(chǔ)
該模型的目的是同時評估幾個因素對生存的影響。換句話說兼搏，它允許我們檢查特定因素如何影響特定時間點特定事件（例如卵慰，感染，死亡）的發(fā)生率佛呻。該比率通常稱為風(fēng)險比率裳朋。預(yù)測變量（或因子）在生存分析文獻中通常稱為協(xié)變量 covariates。
Cox模型由 h(t) 表示的風(fēng)險函數(shù)表示吓著。簡而言之再扭，危險函數(shù)可以解釋為在時間t死亡的風(fēng)險∫勾＃可以估計如下：
其中：
- t 表示生存時間
- h(t) 是由一組p協(xié)變量（x1泛范，x2，…紊撕，xp）確定的風(fēng)險函數(shù)
- 系數(shù)（b1罢荡，b2，…，bp）衡量協(xié)變量的影響（即效應(yīng)大星浴）（可以理解為風(fēng)險權(quán)重）
- 術(shù)語 h0 稱為基線危險惭缰。如果所有的系數(shù)等于零（既exp(0)等于1），則對應(yīng)于風(fēng)險的值笼才。h(t)中的“t”提醒我們漱受，危險可能隨時間而變化。
Cox模型可以被寫為變量x(i)的危險對數(shù)的多元線性回歸骡送，而基線危險是隨時間變化的“截距”項昂羡。
系數(shù) bi 稱為危險比率（HR，hazard ratio）摔踱。bi 值大于零虐先，或相當(dāng)于風(fēng)險比率大于1，表明隨著第 i 個協(xié)變量值的增加派敷，事件風(fēng)險增加蛹批，因此生存時間縮短。
換句話說篮愉，風(fēng)險比大于1表示協(xié)變量與事件概率正相關(guān)腐芍，因此與存活時間負相關(guān)。
總之试躏，
HR=1：無影響
HR<1：危害降低
HR>1：危險增加
在癌癥研究中：
- 風(fēng)險比 > 1（即：b > 0）的協(xié)變量稱為不良預(yù)后因素
- 風(fēng)險比 < 1（即：b < 0）的協(xié)變量被稱為良好的預(yù)后因素
Cox模型的關(guān)鍵假設(shè)是觀察組（或患者）的危險曲線應(yīng)成比例甸赃，并且不能交叉。
假設(shè)兩個x值不同的患者k和k'冗酿。相應(yīng)的風(fēng)險函數(shù)可以簡單地寫成如下：
-
患者k的風(fēng)險函數(shù)：
-
患者k'的風(fēng)險函數(shù)：
-
這兩名患者的危險比 與時間t無關(guān)埠对。
因此，Cox 模型是一個比例風(fēng)險模型：任何一組事件的風(fēng)險都是其他任何一組事件風(fēng)險的常數(shù)倍裁替。這一假設(shè)意味著项玛，如上所述，各組的危險曲線應(yīng)成比例弱判，不能交叉襟沮。
換言之，如果一個人在某個初始時間點的死亡風(fēng)險是另一個人的兩倍昌腰，那么在以后的任何時候开伏，死亡風(fēng)險仍然是另一個人的兩倍。
這種比例風(fēng)險的假設(shè)應(yīng)該得到檢驗遭商。我們將在本系列的下一篇文章中討論評估比例性的方法：Cox模型假設(shè)固灵。
在R中計算Cox模型
安裝并加載所需的R包
我們將使用兩個R包：
survival 用于計算存活分析
survminer 用于可視化生存分析結(jié)果
安裝軟件包
install.packages(c("survival", "survminer"))
- 加載包
library("survival")
library("survminer")
用于計算Cox模型的R函數(shù)：coxph（）
函數(shù)coxph（）[在survival包中]可用于計算R中的Cox比例風(fēng)險回歸模型。
簡化格式如下：
coxph(formula, data, method)
- formula：是以生存對象為響應(yīng)變量的線性模型劫流。使用函數(shù)Surv（）創(chuàng)建生存對象巫玻，如下所示：Surv（time丛忆，event）。
- data：包含變量的數(shù)據(jù)框
- method：用于指定如何處理
tie仍秤。默認值為“efron”熄诡。其他選項有“breslow”和“exact”。默認的“efron”通常優(yōu)先于曾經(jīng)流行的“breslow”方法诗力』烁。“exact”方法的計算量要大得多。
示例數(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
- 機構(gòu)：機構(gòu)代碼
- 時間：以天為單位的生存時間
- 狀態(tài)：審查狀態(tài)1 =審查袜茧，2 =失效
- 年齡：歲
- 性別：男= 1女= 2
- ph.ecog：ECOG成績得分（0 =好5 =死）
- ph.karno：醫(yī)師對Karnofsky績效評分（不良= 0-良好= 100）
- pat.karno：Karnofsky表現(xiàn)評分，按患者評分
- meal.cal：用餐時消耗的卡路里
- wt.loss：最近六個月的體重減輕
計算考克斯模型
我們將使用以下協(xié)變量來擬合Cox回歸：年齡圈澈，性別，ph.ecog和wt.loss尘惧。
我們首先為所有這些變量計算單變量Cox分析康栈。然后我們將使用兩個變量來擬合多元Cox分析，以描述這些因素如何共同影響生存喷橙。
單變量Cox回歸
單變量Cox分析的計算公式如下：
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)生更完整的報告：
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é)果可以解釋為：
統(tǒng)計意義啥么。標(biāo)記為“ z”的列給出了Wald統(tǒng)計值。它對應(yīng)于每個回歸系數(shù)與其標(biāo)準(zhǔn)誤差的比率（z = coef / se（coef））贰逾。Wald統(tǒng)計量會評估beta給定變量的系數(shù)在統(tǒng)計上顯著不同于0悬荣。從上面的輸出中，我們可以得出結(jié)論疙剑，變量性別具有很高的統(tǒng)計學(xué)上顯著的系數(shù)氯迂。
回歸系數(shù)。Cox模型結(jié)果中要注意的第二個特征是回歸系數(shù)（coef）的符號言缤。正號表示該變量值較高的受試者的危險（死亡風(fēng)險）較高嚼蚀，因此預(yù)后較差。變量sex被編碼為數(shù)字向量管挟。1：男性轿曙，2：女性。Cox模型的 R summary 提供了第二組相對于第一組（即女性與男性）的危險比（HR）僻孝。在這些數(shù)據(jù)中导帝，性別的beta系數(shù)= -0.53表示女性比男性具有較低的死亡風(fēng)險（較低的生存率）。
危險率穿铆。指數(shù)系數(shù)（exp（coef）= exp（-0.53）= 0.59）您单，也稱為危險比，給出了協(xié)變量的影響大小荞雏。例如睹限，女性（性別= 2）可以將危險降低0.59或41％譬猫。女性是有良好預(yù)后的。
危險比的置信區(qū)間羡疗。摘要輸出還給出了危險比（exp（coef））的上限和下限95％置信區(qū)間染服，下限95％= 0.4237，上限95％= 0.816叨恨。
該模型的全局統(tǒng)計意義柳刮。最后，輸出為模型的總體重要性提供了三個備選檢驗的p值：似然比檢驗痒钝，Wald檢驗和得分對數(shù)秩統(tǒng)計秉颗。這三種方法漸近等效。對于足夠大的N送矩，它們將給出相似的結(jié)果蚕甥。對于較小的N，它們可能會有所不同栋荸。對于小樣本量菇怀，似然比測試具有更好的行為，因此通常是首選晌块。
要將單變量coxph函數(shù)一次應(yīng)用于多個協(xié)變量爱沟，請輸入以下命令：
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
上面的輸出顯示了每個變量相對于總生存率的回歸beta系數(shù)，效應(yīng)大写冶场（以危險比給出）和統(tǒng)計顯著性呼伸。通過單獨的單變量Cox回歸評估每個因素。
從上面的輸出中钝尸，
性別括享，年齡和ph.ecog變量具有極高的統(tǒng)計顯著性系數(shù)，而wt.loss的系數(shù)則不顯著珍促。
年齡和ph.ecog的β系數(shù)為正奶浦，而性別的系數(shù)為負。因此踢星，年齡更大和ph.ecog較高與較差的生存率有關(guān)澳叉，而女性（性別= 2）則與較好的生存率有關(guān)。
現(xiàn)在沐悦，我們要描述這些因素如何共同影響生存成洗。為了回答這個問題，我們將執(zhí)行多元Cox回歸分析藏否。由于變量ph.karno在單變量Cox分析中不重要瓶殃，因此在多變量分析中將其跳過。我們將3個因素（性別副签，年齡和ph.ecog）納入多元模型遥椿。
多元Cox回歸分析
時間常數(shù)協(xié)變量的死亡時間的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
所有三個總體測試（似然性基矮，Wald和得分）的p值均顯著，表明該模型具有顯著性冠场。這些測試評估了所有beta的綜合零假設(shè)為0家浇。在上面的示例中，檢驗統(tǒng)計量非常一致碴裙，并且完全拒絕了綜合零假設(shè)磁椒。
在多變量Cox分析中务傲，協(xié)變量性別和ph.ecog保持顯著性（p <0.05）坷随。但是饱岸，協(xié)變量年齡不顯著（p = 0.23，大于0.05）载慈。
性別的p值為0.000986惭等，危險比HR = exp（coef）= 0.58，表明患者的性別與死亡風(fēng)險降低之間有很強的關(guān)系办铡。協(xié)變量的危險比可解釋為對危險的倍增效應(yīng)辞做。例如，保持其他協(xié)變量不變（女性（性別= 2））可將危險降低0.58或42％料扰。我們得出結(jié)論凭豪，成為女性與良好的預(yù)后相關(guān)焙蹭。
同樣晒杈，ph.ecog的p值為4.45e-05，危險比HR = 1.59孔厉，表明ph.ecog值與死亡風(fēng)險增加之間有很強的關(guān)系拯钻。保持其他協(xié)變量不變，ph.ecog的值越高撰豺，生存率越低粪般。
相比之下，年齡的p值現(xiàn)在為p = 0.23污桦。危險比HR = exp（coef）= 1.01亩歹，95％置信區(qū)間為0.99至1.03。由于HR的置信區(qū)間為1凡橱，因此這些結(jié)果表明小作，在調(diào)整phog值和患者的性別之后，年齡對HR差異的貢獻較小稼钩，并且僅趨于顯著顾稀。例如，在其他協(xié)變量保持不變的情況下坝撑，再增加一歲會引起每日死亡危險静秆，其系數(shù)為expβ= 1.01或1％粮揉，這并不是一個重要的貢獻。
可視化生存時間的估計分布
將Cox模型擬合到數(shù)據(jù)后抚笔，就可以可視化特定風(fēng)險組在任何給定時間點的預(yù)測生存率扶认。函數(shù)survfit（）估計生存比例，默認情況下為協(xié)變量的平均值塔沃。
# Plot the baseline survival function
ggsurvplot(survfit(res.cox), color = "#2E9FDF",
ggtheme = theme_minimal())
我們不妨展示估計的生存率如何取決于目標(biāo)協(xié)變量的值蝠引。
考慮到這一點，我們想評估性別對估計生存率的影響蛀柴。在這種情況下螃概，我們用兩行構(gòu)造一個新的數(shù)據(jù)幀，每一行代表性別鸽疾。其他協(xié)變量固定為其平均值（如果是連續(xù)變量）或最低水平（如果它們是離散變量）吊洼。對于偽協(xié)變量，平均值為數(shù)據(jù)集中編碼為1的比例制肮。該數(shù)據(jù)幀通過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回歸模型，用于同時評估多種風(fēng)險因素與患者生存時間之間的關(guān)系豺鼻。我們演示了如何使用生存包計算Cox模型综液。此外，我們描述了如何使用survminer軟件包來可視化分析結(jié)果儒飒。
參考文獻
- 考克斯博士（1972）谬莹。回歸模型和壽命表（有討論）桩了。JR Statist Soc B 34：187–220
- MJ Bradburn附帽，TG Clark，SB Love和DG Altman井誉。生存分析第二部分：多元數(shù)據(jù)分析–概念和方法介紹蕉扮。英國癌癥雜志（2003）89，431 – 436
