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

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

需要進(jìn)行多元統(tǒng)計(jì)建模

在臨床研究中,，有許多情況下，幾個(gè)已知量（稱為協(xié)變量）可能影響患者,。

統(tǒng)計(jì)模型是一個(gè)經(jīng)常使用的工具,，可以同時(shí)分析多個(gè)因素的生存情況。另外,，統(tǒng)計(jì)模型提供了每個(gè)因素的效應(yīng)大小,。

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

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

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

h（t）=h0（t）×e×p（b1X1+b2X2+,。。,。+bpXp）H（?）=H0（?）×?Xp（b1X1+b2X2+,。。,。+bpXp）

換言之,，高于1的風(fēng)險(xiǎn)比指示與事件概率正相關(guān)的協(xié)變量，并因此與生存期的長(zhǎng)短負(fù)相關(guān),。

綜上所述,，

HR = 1：沒(méi)有效果

HR <1：減少危險(xiǎn)

HR> 1：危害增加

請(qǐng)注意，在癌癥研究中：

危險(xiǎn)比> 1（即：b> 0）的協(xié)變量被稱為不良預(yù)后因素

危險(xiǎn)比<1（即：b <0）的協(xié)變量被稱為良好的預(yù)后因子

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

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

換句話說(shuō),，如果一個(gè)人在某個(gè)初始時(shí)間點(diǎn)有死亡風(fēng)險(xiǎn)，是另一個(gè)人的兩倍,，那么在所有的晚些時(shí)候,，死亡風(fēng)險(xiǎn)仍然是兩倍。

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

安裝并加載所需的R包

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

生存計(jì)算生存分析

幸存者可視化生存分析結(jié)果

安裝軟件包

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

加載包

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

R函數(shù)來(lái)計(jì)算Cox模型：coxph（）

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

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

coxph(formula,data,method)

公式：以生存對(duì)象為響應(yīng)變量的線性模型,。Survival對(duì)象是使用Surv（）函數(shù)創(chuàng)建的，如下所示：Surv（time,，event）,。

數(shù)據(jù)：包含變量的數(shù)據(jù)框

方法：用于指定如何處理關(guān)系。默認(rèn)是'efron',。其他選項(xiàng)是“breslow”和“確切”,。默認(rèn)的“efron”通常比一度流行的“breslow”方法更受歡迎?！按_切”的方法計(jì)算密集得多,。

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

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

data("lung")head(lung)

inst：機(jī)構(gòu)代碼

時(shí)間：以天為單位的生存時(shí)間

狀態(tài)：審查狀態(tài)1 =審查,，2 =死亡

年齡：年齡

性別：男= 1女= 2

ph.ecog：ECOG表現(xiàn)評(píng)分（0 =好5 =死）

ph.karno：Karnofsky表現(xiàn)評(píng)分（bad = 0-好= 100）由醫(yī)師評(píng)定

pat.karno：Karnofsky表現(xiàn)評(píng)分由患者評(píng)估

膳食：餐時(shí)消耗的卡路里

wt.loss：過(guò)去六個(gè)月的體重下降

計(jì)算Cox模型

我們將使用以下協(xié)變量進(jìn)行Cox回歸：年齡，性別,，ph.ecog和wt.loss,。

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

單變量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      psex -0.531    0.588    0.167 -3.18 0.0015Likelihood ratio test=10.6  on 1 df, p=0.00111n= 228, number of events= 165

Cox模型的函數(shù)summary（）產(chǎn)生更完整的報(bào)告：

summary(res.cox)

Call:coxph(formula = Surv(time, status) ~ sex, data = lung)n= 228, number of events= 165coef 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 ' ' 1exp(coef) exp(-coef) lower .95 upper .95sex    0.588      1.701    0.4237    0.816Concordance= 0.579  (se = 0.022 )Rsquare= 0.046  (max possible= 0.999 )Likelihood ratio test= 10.63  on 1 df,  p=0.001111Wald test            = 10.09  on 1 df,  p=0.001491Score (logrank) test = 10.33  on 1 df,  p=0.001312

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

統(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ì)評(píng)估是否beta（ββ）系數(shù)在統(tǒng)計(jì)上顯著不同于0,。從上面的輸出,，我們可以得出結(jié)論，變量性別具有高度統(tǒng)計(jì)學(xué)意義的系數(shù),。

回歸系數(shù),。Cox模型結(jié)果中要注意的第二個(gè)特征是回歸系數(shù)（coef）的符號(hào)。一個(gè)積極的信號(hào)意味著危險(xiǎn)（死亡風(fēng)險(xiǎn)）較高,，因此對(duì)于那些變量值較高的受試者,，預(yù)后更差。變量性被編碼為數(shù)字向量,。1：男,，2：女。Cox模型的R總結(jié)給出了第二組相對(duì)于第一組,，即女性與男性的風(fēng)險(xiǎn)比（HR）,。性別的β系數(shù)= -0.53表明在這些數(shù)據(jù)中,，女性的死亡風(fēng)險(xiǎn)（低存活率）低于男性。

危害比例,。指數(shù)系數(shù)（exp（coef）= exp（-0.53）= 0.59）也稱為風(fēng)險(xiǎn)比,，給出協(xié)變量的效應(yīng)大小。例如,，女性（性別= 2）將危害降低了0.59倍,，即41％。女性與預(yù)后良好相關(guān),。

風(fēng)險(xiǎn)比的置信區(qū)間,。總結(jié)結(jié)果還給出了風(fēng)險(xiǎn)比（exp（coef））的95％置信區(qū)間的上限和下限,，下限95％界限= 0.4237,，上限95％界限= 0.816。

全球統(tǒng)計(jì)學(xué)意義的模型,。最后,，輸出為模型的總體顯著性提供了三個(gè)替代測(cè)試的p值：可能性比率測(cè)試，Wald測(cè)試和得分logrank統(tǒng)計(jì),。這三種方法是漸近等價(jià)的,。對(duì)于足夠大的N，他們會(huì)得到相似的結(jié)果,。對(duì)于小N來(lái)說(shuō),，它們可能有所不同。似然比檢驗(yàn)對(duì)于小樣本量具有更好的表現(xiàn),，所以通常是優(yōu)選的,。

要一次性將單變量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 datauniv_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 betaHR<-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.valueage      0.019            1 (1-1)       4.1  0.042sex      -0.53  0.59 (0.42-0.82)        10  0.0015ph.karno -0.016      0.98 (0.97-1)      7.9  0.005ph.ecog    0.48        1.6 (1.3-2)        18 2.7e-05wt.loss  0.0013        1 (0.99-1)      0.05    0.83

上面的輸出顯示了與總體生存相關(guān)的每個(gè)變量的回歸β系數(shù),，效應(yīng)大?。ㄗ鳛槲：Ρ冉o出）和統(tǒng)計(jì)顯著性。每個(gè)因素通過(guò)單獨(dú)的單變量Cox回歸評(píng)估,。

從上面的輸出中,，

變量性別，年齡和ph.ecog具有高度統(tǒng)計(jì)學(xué)意義的系數(shù),，而ph.karno的系數(shù)不顯著,。

年齡和ph.ecog有正的beta系數(shù)，而性別有負(fù)系數(shù),。因此,，年齡越大和ph.ecog越高，生存率越差，而女性（性別= 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.232416sex    -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 ' ' 1exp(coef) exp(-coef) lower .95 upper .95age        1.0111    0.9890    0.9929    1.0297sex        0.5754    1.7378    0.4142    0.7994ph.ecog    1.5900    0.6289    1.2727    1.9864Concordance= 0.637  (se = 0.026 )Rsquare= 0.126  (max possible= 0.999 )Likelihood ratio test= 30.5  on 3 df,  p=1.083e-06Wald test            = 29.93  on 3 df,  p=1.428e-06Score (logrank) test = 30.5  on 3 df,  p=1.083e-06

所有三個(gè)整體測(cè)試（可能性，Wald和得分）的p值都是顯著的,，表明該模型是顯著的,。這些測(cè)試評(píng)估了所有的beta（ββ）為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é)變量的風(fēng)險(xiǎn)比可以解釋為對(duì)風(fēng)險(xiǎn)的倍增效應(yīng),。例如，保持其他協(xié)變量不變,，女性（性別= 2）將危害降低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。風(fēng)險(xiǎn)比HR = exp（coef）= 1.01，95％置信區(qū)間為0.99至1.03,。由于HR的置信區(qū)間為1,，這些結(jié)果表明，年齡在調(diào)整了ph值和患者性別后對(duì)HR的差異的貢獻(xiàn)較小,，并且僅趨向顯著性,。例如，保持其他協(xié)變量不變,，額外的年齡會(huì)導(dǎo)致每日死亡危險(xiǎn)因素為exp（beta）= 1.01或1％,，這不是一個(gè)重要的貢獻(xiàn)。

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

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

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

我們可能希望展示估計(jì)的生存如何取決于感興趣的協(xié)變量的值,。

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

# Create the new datasex_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.ecog1  1 62.44737      12  2 62.44737      1

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