作者：赵天业；审核：龚志忠
已经知道要进行Cox回归，但是不知道如何使用R来操作的小伙伴，你们想要的教程来啦！本文将提供非常详细的教程帮助大家使用R来进行Cox回归。
PS. 如果对是否选择Cox回归分析还有疑虑，敬请阅读医咖会在2017年发布的文章《》
Cox模型也称Cox比例风险模型，它的假设是“比例风险”，即在研究的全程，不同患者发生结局事件的风险之比是恒定的。我们使用的协变量是研究开始时（基线）患者的特征，也就是说我们假设基线特征决定了患者发生结局事件的相对风险。在许多研究问题中，这种假设很可能是不成立的，这也是Cox回归一个重要的局限性。我们应该先确定协变量满足比例风险假设（Proportional hazards assumption，PH假定），再进行Cox回归。
但在实践中，我们可以先构建Cox模型，通过对Cox模型的评价来评估协变量是否满足PH假定。如果协变量满足PH假定，则模型是可用的，反之则应该对数据进行调整或更换统计分析方法。进行Cox回归前应观察生存曲线，如果不同分组的生存曲线发生交叉，则该变量很可能不满足PH假定。
Cox回归分析的输出是Cox模型的参数估计值和相关的统计检验。Cox模型的参数可以用于估计协变量对结局事件发生的影响，并估计HR（风险比）及其置信区间。HR>1表示协变量为“危险因素”，即协变量每增加一个单位（连续型变量）或者与参考水平相比（分类变量），结局事件发生的风险升高了（HR-1）×100%；反之，HR在0-1之间表示协变量为“保护因素”，即协变量每增加一个单位（连续型变量）或者与参考水平相比（分类变量），结局事件发生的风险降低了（1-HR）×100%。
一
准备工作
1 输入和输出
我们首先明确在R中进行Cox回归的输入和输出是什么。在R中使用Survival包进行Cox回归，最核心的两行代码如下：
cox_model <- coxph(Surv(time, status) ~ v1 + v2 + v3, data = data)
summary(cox_model)其中data是输入的数据，status（结局事件）、time（生存时间）和v（协变量）是data包含的变量（列），summary为输出Cox回归结果。
Cox回归的输入是一个数据文件，必须包含的变量有结局事件、生存时间和协变量。结局变量需要是数值型二分类变量，通常用数字1和0表示结局事件发生和未发生。生存时间为从随访开始至结局事件发生的时间间隔。对于未观察到结局事件的患者，生存时间计算至最后一个已知的、患者未发生结局事件的时间点，通常是末次随访，或者患者因其他原因死亡的时间。协变量可以是分类变量、连续变量或其他类型的变量。在多因素分析时，不同类型的协变量可以一起使用。有些协变量可能适合转化为其他类型的变量，例如将连续变量分为“高和低”，成为二分类变量，或者分为“高、中、低”，成为等级变量。
结局变量建议将文字编码为数字，例如用1和0表示死亡和生存；如果结局变量不是二分类，例如将未发生结局事件但因其他原因死亡的患者分为单独一类合并。
图. 用于Cox回归的数据示例
图中展示的是一个名为“mayo”的数据文件，censor为结局变量，time为生存时间，mayoscore5和mayoscore4分别是两个协变量。
二
单因素Cox回归
这里使用survivalROC包附带的数据集“mayo”进行演示，这种数据集属于软件包，只要加载软件包就可以调用，但为了模拟我们自己的数据，先通过以下代码导出为csv文件，再进行使用。如果没有安装survivalROC，可以先安装。
本文演示在R项目中工作（参见https://martinctc.github.io/blog/rstudio-projects-and-working-directories-a-beginner's-guide/）。此处文件路径中开头的点（.）的写法表示R项目中的路径，不是省略。
library(survivalROC)
data(mayo)
write.csv(mayo, "./Data/mayo.csv")进行Cox回归时，首先加载软件包，并读取数据。注：read.csv省略了一些默认参数，例如sep = ","，encoding = "UTF-8"，读取我们自己的数据时，可能需要调整这些参数。
library(survival)
data <- read.csv(file = "./Data/mayo.csv", header = TRUE)这里使用mayoscore5作为协变量，构建Cox模型。cox_model为生成的Cox模型， censor为结局变量，time为生存时间，Surv(time, censor) ~后连接协变量构成了“formula”。Summary函数的输入是Cox模型。
cox_model <- coxph(Surv(time, censor) ~ mayoscore5, data = data)
summary(cox_model)输出如下：
Call:
coxph(formula = Surv(time, censor) ~ mayoscore5, data = data)
n= 312, number of events= 125
coef exp(coef) se(coef) z Pr(>|z|)
mayoscore5 1.00086 2.72063 0.07134 14.03 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
mayoscore5 2.721 0.3676 2.366 3.129
Concordance= 0.843 (se = 0.02 )
Likelihood ratio test= 198.5 on 1 df, p=<2e-16
Wald test = 196.8 on 1 df, p=<2e-16
Score (logrank) test = 252.2 on 1 df, p=<2e-16得到输出之后，首先要检查结果中对数据的描述与数据的实际情况是否一致，例如患者是否为312人，发生结局事件的患者是否为125人。进行这种检查可以避免错误，并有助于寻找结果异常的原因。结果显示，coef为正数，表明随着协变量的增加，结局事件发生的风险增加（反之，如果水平为负数，则说明随着自变量水平的增加，结局事件发生的风险降低）；exp(coef)为2.721，即HR=2.72；lower .95和upper .95分别为2.366和3.129，即HR的95%置信区间为（2.37，3.13）；z为Wald检验的统计量，Pr(>|z|)<2e-16（2×10-16），即Wald检验的P值<0.01；HR的95%置信区间不含1，也提示协变量对结局事件的发生有显著影响。在解读结果时，可以认为mayoscore5每增加1，发生结局事件的风险升高172%。
在采用这些结果之前，还需要进行Schoenfeld残差检验，以评估协变量是否满足PH假定。为此使用survival包的cox.zph函数，输入是刚刚得到的cox模型。
cox.zph(cox_model)输出如下。
chisq df p
mayoscore5 0.541 1 0.46
GLOBAL 0.541 1 0.46结果显示，对协变量mayoscore5进行残差检验得到的χ2=0.541（chisq为χ2），P值为0.46，大于0.05，可以认为满足PH假定。
前述summary的结果还包含了其他可以根据研究目的报道的统计检验结果。
C（一致性）指数用来反应预测结果和实际情况的一致性。在Cox模型中，不同学者对什么是“预测结果和实际情况一致”以及如何评价有不同的观点，由此产生了多种估计C指数的方法，coxph函数提供其中一种C指数，即哈氏（Harrell's）C指数。Concordance= 0.843 (se = 0.02 )，即C指数为0.843，其SE（标准误）为0.02，可以用C指数±1.96×SE估计C指数的95%置信区间，即（0.804, 0.882）。
Likelihood ratio test（似然比检验）、Wald test和Score (logrank) test评价Cox模型总体的显著性，如果这三种检验不显著，提示模型中纳入的协变量均缺乏预测结局事件发生的能力。在这个例子中，三者均显著。
三
连续变量转为二分类变量
在进行Cox回归时，有时将连续变量转为分类变量可以得到临床意义更明确、更容易解读的结果。在确定截断值时，应参考相关文献和指南，如果文献中没有报道，可以用survminer包的surv_cutpoint函数通过Maximally Selected Test Statistics选择Log-rank检验的统计量最大（相应地，P值最小）的截断值，作为“最佳截断值”。这里继续使用“mayo”进行演示，将mayoscore5转为二分类变量。注：#开头的行为注释。
library(survminer)
sur.cut <- surv_cutpoint(data, time = "time",
event = "censor",
variables = "mayoscore5")
summary(sur.cut)
#对截断值两侧的数据分布可视化。
plot(sur.cut, "mayoscore5", palette = "npg")summary(sur.cut)的输出是：
cutpoint statistic
mayoscore5 6.627459 11.31263
图. plot输出区域的截图
即连续型变量“mayoscore5”的最佳截断值为6.63，标准化的Log-rank检验统计量为11.313。可以新建一个变量，利用survminer包的surv_categorize函数，使mayoscore5小于等于截断值时这个变量取值为0（较低得分），大于截断值取值为1（较高得分）。
#使用surv_categorize函数，根据截断值转为二分类变量。
sur.cat <- surv_categorize(sur.cut, labels = c(0, 1))注意此时数据框data已经不再是我们刚读取“mayo.csv”时的样子了。在处理数据的过程中，新增变量和调整变量的格式是很常见的做法，我们可以将处理后的数据保存，以便后续使用。再次使用时，可以直接读取，读取的数据框的名称仍为“data”。
#存储为.RData文件。
save(data, file = "./Data/mayo.RData")
#再次读取时，使用：
load(file = "./Data/mayo.RData")四
使用分类变量
这里继续使用刚刚创建的变量mayoscore5_cut演示在Cox回归中使用二分类变量。运行class(data)，输出 "data.frame"，即我们读取的数据被存储为数据框；运行class(data$mayoscore5_cut)，输出"numeric"，即这个变量的格式为数字（$表示前者的某一列），可以用as.factor函数转为factor（因子），并选定参考水平。
进行Cox回归时，其他取值将与参考水平比较，故输出将不包含参考水平的结果。对于分类变量，R选择参考水平的默认规则是最小的数字（如果以数字编码）或者最小的数字/顺序最靠前的字母（如果以文本编码），使用relevel函数可以指定一个格式为factor的变量的参考水平。
data$mayoscore5_cut <- as.factor(data$mayoscore5_cut)
#0为较低得分，1为较高得分，指定较高得分为参考水平。参考水平需要加引号。
data$mayoscore5_cut <- relevel(data$mayoscore5_cut, ref = "1")参考水平的选择会影响输出结果。
#按照刚刚的设置，以较高得分为参考水平。
cox_model_1 <- coxph(Surv(time, censor) ~ mayoscore5_cut, data = data)
summary(cox_model_1)
#将参考水平改为较低得分。
data$mayoscore5_cut <- relevel(data$mayoscore5_cut, ref = "0")
cox_model_2 <- coxph(Surv(time, censor) ~ mayoscore5_cut, data = data)
summary(cox_model_2)两次summary的输出为（摘录）：
#参考水平为较高得分。
exp(coef) exp(-coef) lower .95 upper .95
mayoscore5_cut0 0.09684 10.33 0.06603 0.142
#参考水平为较低得分。
exp(coef) exp(-coef) lower .95 upper .95
mayoscore5_cut1 10.33 0.09684 7.04 15.15可以看到，当参考水平为较低得分时，HR=10.33；当参考水平为较高得分时，HR=0.097（1/10.33）。可以解读为，mayoscore5较高得分的患者发生结局事件的风险是较低得分的患者的10.33倍，较低得分的患者发生结局事件的风险是较高得分的患者的0.097倍。
对于多分类变量，有时将其转为二分类变量，可能有助于得到更容易解释的结果。如果不适合转为二分类变量，可以转为哑变量，请参考《》。在进行Cox回归时，以factor格式存在的、有n个水平的多分类变量会被自动转为n-1个哑变量进入模型，这些哑变量的HR值的是与参考水平比较的结果。
五
多因素Cox回归与逐步回归
我们使用survival包附带的数据集“lung”演示多因素Cox回归。为了模拟我们自己的数据，仍然通过以下代码导出为csv文件，再进行使用。
library(survival)
data(lung)
write.csv(lung, "./Data/lung.csv")
#读取数据
data_1 <- read.csv(file = "./Data/lung.csv", header = TRUE)在这个数据集中，status为结局变量，time为生存时间。此时这个数据集的变量sex（性别）的格式为“integer”（整数），我们使用as.factor将其处理为factor。ph.ecog（一种评分）是等级变量，赋值为0-4五个等级。在处理等级变量时，如果我们基于数据的性质/基于专业知识/主观认为不同等级之间的差异是“等距”的，可以作为连续型变量使用，否则可以处理为多分类变量。在这个数据集中，ph.ecog的评分规则是“0=无症状，1=有症状但完全可以走动，2=每天<50%的时间在床上，3=在床上>50%的时间但没有卧床，4=卧床”，不同等级的差异看起来并不“等距”，我们将其作为多分类变量使用。
进行多因素Cox回归时，协变量之间用+相连。这里纳入age等五个变量进行多因素Cox回归。
data_1$sex <- as.factor(data_1$sex)
#以女性为参考水平
data_1$sex <- relevel(data_1$sex, ref = 2)
data_1$ph.ecog <- as.factor(data_1$ph.ecog)
#以ph.ecog=0为参考水平，因其为最小的数字，无需调整。
cox_model_3 <- coxph(Surv(time, status) ~ age + sex + ph.ecog + meal.cal + wt.loss, data = data_1)
summary(cox_model_3)
#检验PH假定
cox.zph(cox_model_3)输出如下。可以看到因为（任一变量）存在缺失值，有58个患者被从分析中排除。
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + meal.cal +
wt.loss, data = data_1)
n= 170, number of events= 123
(因为不存在，58个观察量被删除了)
coef exp(coef) se(coef) z Pr(>|z|)
age 7.511e-03 1.008e+00 1.136e-02 0.661 0.508550
sex1 5.346e-01 1.707e+00 1.998e-01 2.676 0.007452 **
ph.ecog1 3.358e-01 1.399e+00 2.333e-01 1.439 0.150148
ph.ecog2 1.002e+00 2.725e+00 2.926e-01 3.426 0.000613 ***
ph.ecog3 2.086e+00 8.050e+00 1.049e+00 1.989 0.046709 *
meal.cal -3.972e-05 1.000e+00 2.530e-04 -0.157 0.875255
wt.loss -1.194e-02 9.881e-01 7.720e-03 -1.547 0.121824
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0075 0.9925 0.9854 1.030
sex1 1.7068 0.5859 1.1538 2.525
ph.ecog1 1.3990 0.7148 0.8855 2.210
ph.ecog2 2.7246 0.3670 1.5355 4.835
ph.ecog3 8.0502 0.1242 1.0308 62.866
meal.cal 1.0000 1.0000 0.9995 1.000
wt.loss 0.9881 1.0120 0.9733 1.003
Concordance= 0.642 (se = 0.029 )
Likelihood ratio test= 23.46 on 7 df, p=0.001
Wald test = 23.86 on 7 df, p=0.001
Score (logrank) test = 25.74 on 7 df, p=6e-04
> cox.zph(cox_model_3)
chisq df p
age 6.70e-01 1 0.4130
sex 2.02e+00 1 0.1557
ph.ecog 5.49e+00 3 0.1391
meal.cal 7.20e+00 1 0.0073
wt.loss 7.18e-05 1 0.9932
GLOBAL 1.39e+01 7 0.0532结果显示，除了meal.cal，其他协变量均满足PH假定。结果中sex1表示以女性作为参考水平，HR为1.71，95%置信区间为（1.15，2.53），P值<0.01，提示在调整了多因素Cox模型中其他变量之后，男性发生结局事件的风险显著大于女性，男性发生结局事件的风险是女性的1.71倍。
作为多分类变量的ph.ecog，根据是否为1、是否为2、是否为3建立了三个哑变量（数据集中只有得分0-3分的患者，没有4分的患者），并分别计算了P值、估计了HR和置信区间。ph.ecog为1分的P值为0.15，95%置信区间（0.89，2.21）；2分和3分的P值均小于0.05，95%置信区间分别为（1.54，4.84）和（1.03，62.87），表明在调整了多因素Cox模型中其他变量之后，得1分的患者与得0分的患者相比，发生结局事件的风险没有显著差异，得2分和得3分的患者发生结局事件的风险显著高于得0分的患者。对于满足PH假定的多分类变量，可以使用multcomp包的glht函数进行多重比较，检验特定两组的回归系数之间的差异是否显著（由回归系数可计算HR）。
library(multcomp)
#比较方法使用Tukey法。
contr <- glht(cox_model_3, linfct = mcp(ph.ecog = "Tukey"))
summary(contr)输出显示得2分和得3分的回归系数没有显著差异（3 - 2 == 0，P=0.69），也即无法认为两组的HR有显著差异。
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + meal.cal +
wt.loss, data = data_1)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
1 - 0 == 0 0.3358 0.2333 1.439 0.43797
2 - 0 == 0 1.0023 0.2926 3.426 0.00256 **
3 - 0 == 0 2.0857 1.0487 1.989 0.16623
2 - 1 == 0 0.6666 0.2390 2.789 0.02202 *
3 - 1 == 0 1.7499 1.0341 1.692 0.29358
3 - 2 == 0 1.0834 1.0313 1.050 0.69081
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)对于我们自己的数据，可能包含很多协变量，需要进行预处理，对于缺失数据可以通过将含有缺失值的行删除，或者进行插补实现。
#调整status的编码
data_1$status1 <- ifelse(data_1$status==2,1,0)
#删除有缺失值的行
data_2 <- na.omit(data_1)这里介绍两种逐步回归的方法。第一种使用step函数，第一个参数为包含了所有计划纳入的协变量的Cox模型，第二个参数为逐步回归的方向，可以选择“both”（双向）、“backward”（向后）和“forward”（向前）。这里演示基于AIC（Akaike's Information Criterion）的向后回归法。
cox_model_4 <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno + pat.karno + meal.cal + wt.loss, data = data_2)
step_model <- step(cox_model_4, direction="backward")输出的底部给出了逐步回归后的模型，可以将其复制，粘贴到R程序中备用。
Step: AIC=1002.15
Surv(time, status) ~ sex + ph.ecog + ph.karno + wt.loss完成变量筛选后我们只将这些协变量纳入Cox回归，进行多因素Cox回归（输出略）。
cox_model_4_1 <- coxph(Surv(time, status) ~ sex + ph.ecog + ph.karno + wt.loss, data = data_2)
summary(cox_model_4_1)第二种方法是使用近年发布的StepReg包的stepwise函数。选择标准可以设置为AIC，也可以设置为AICc（Corrected Akaike's Information Criterion，校正后的AIC）。这里演示向前回归法。stepwise函数还可以使用更多选择标准，参阅StepReg包的vignette（https://cran.r-project.org/web/packages/StepReg/vignettes/StepReg.html）。
library(StepReg)
step <- stepwise(Surv(time, status) ~ age + sex + ph.ecog + ph.karno + pat.karno + meal.cal + wt.loss,
data = data_2,
type = "cox",
strategy = "forward",
metric = "AICc")
step输出为四个表格，第四个表格列出了筛选后纳入的协变量。
Table 4. Summary of coefficients for the selected model with Surv(time, status) under forward and AICc
‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗
Variable coef exp(coef) se(coef) z Pr(>|z|)
———————————————————————————————————————————————————————————————
ph.ecog1 0.650124 1.915779 0.280696 2.316116 0.020552
ph.ecog2 1.676944 5.349182 0.441824 3.795501 0.000147
ph.ecog3 2.88359 17.878348 1.121914 2.570242 0.010163
sex1 0.564688 1.758899 0.199667 2.828153 0.004682
wt.loss -0.012793 0.987289 0.007677 -1.666438 0.095626
ph.karno 0.018538 1.018711 0.011154 1.662025 0.096508
‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗‗相比单因素Cox回归，多因素Cox回归可以控制混杂偏倚，可以构建更加全面、准确的模型，还可能评估协变量之间的交互作用。然而，在进行多因素Cox回归时，仍有必要进行单因素Cox回归，因为需要对比单因素和多因素Cox回归的结果。如果一个协变量在单因素分析中显著，在多因素分析中不显著，提示它可能实际上不具有预测结局事件发生的能力。如果一个协变量的HR在单因素与多因素分析中有明显的差异，甚至效应方向相反，例如从危险因素变为保护因素，则很可能需要研究者对此进行讨论。
六
进行多因素Cox回归时
其他应考虑的事
进行多因素Cox回归时，还有一些事情需要考虑：
（1）在更高的预测能力与更简单的模型之间平衡。研究者可以将Cox回归的结果以列线图的形式展示，或者基于Cox回归的结果构建评分，来构建在临床工作中使用的工具。向Cox回归中引入更多的协变量，可能有助于提高模型的预测能力，起到“兼听则明”之效。然而，在模型中保留更多的协变量，特别是需要经不同方法测量才能得到的协变量，也意味着使用这个工具需要更多的前提与成本，降低其临床应用价值。
（2）协变量之间的多重共线性。如果两个协变量之间有较强的相关性，Cox模型会难以将它们的效应分离，导致回归的结果不准确。可以通过绘制协变量之间的散点图，或者分析变量之间的相关性，找到具有相关性的变量。还可以使用方差膨胀因子（variance inflation factor，VIF）评价多重共线性对Cox模型的影响。当存在相关性较强的协变量时，可以只将其中一个协变量纳入模型，或者考虑使用岭回归（Ridge Regression）等可以处理多重共线性的回归方法。
（3）数据质量。受限于研究方法与研究的实施质量，我们的数据可能存在错漏，同一个变量在不同患者的测量方式可能不尽相同。研究者在撰写论文/研究报告时，应该在研究方法提供必要的信息，以便读者解读与评价文章中Cox回归的结果是否可靠。
热门跟贴