作者：赵天业；审核：龚志忠
已经知道要进行Cox回归，但是不知道如何使用R来操作的小伙伴，你们想要的教程来啦！本文将提供非常详细的教程帮助大家使用R来进行Cox回归。
PS. 如果对是否选择Cox回归分析还有疑虑，敬请阅读医咖会在2017年发布的文章《如何使用R来做Cox回归？来看详细教程！》
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值的是与参考水平比较的结果。
确认删除