用R语言进行Cox回归生存分析

在生存分析中，探究生存时间的影响因素是一个重要的研究内容，通过KM和log-rank test检验的方法，只能够处理单个二分类因素的生存数据。当想探究多个因素或者离散型变量对生存时间的影响时，我们就需要借助于cox回归方法。cox回归的全称如下

cox proportional hazards regression model

称之为cox等比例风险回归模型， 对应的公式如下

将上述公式进行log转换，可以变换成以下格式

这个公式和逻辑回归的公式就非常的接近了, cox回归其实是在线性回归和逻辑回归的基础上延伸而出的一种方法，将影响生存的多个因素当做回归方程中的自变量，将风险函数h(t)和h0(t）的比值当做因变量。

每个自变量对应的系数，如b1，b2这类的 ，称之为偏回归系数。当偏回归系数大于0时， 随着该自变量值的增加，风险增加，生存时间减少，当系数小于0时，则相反；等于0时，没有影响。

将exp(b)称之为hazard ratio, 简称HR。将偏回归系数转换成HR, 对应的关系如下

1. HR = 1， 没有影响
2. HR > 1, 风险增加
3. HR < 1, 风险降低

在临床上，将HR>1的自变量称之为坏的预后因子，将HR<1的自变量称之为好的预后因子。通过survivalR包，可以轻松的实现cox回归分析，步骤如下

1. 准备生存数据

对于每个个体而言，其生存数据会出现两种情况，第一种是观测到生存时间，通常用1表示，第二种则是删失。通常用0表示。survival自带了一个测试数据lung, 内容如下所示

每一行代表一个样本，time表示生存时间，status表示删失情况，这里只有1和2两种取值，默认排序后的第一个level对应的值为删失，这里则为1表示删失。其他列为样本对应的性别，年龄等基本信息。

2. cox回归分析

代码如下

可以看到，cox回归的适用范围更广，以最后一个回归分析为例，结果如下所示

首先查看likehood ration test , wald test, logrank test三种检验方法的p值，p值小于0.05, 这个回归方程是统计学显著的。说明在这么多自变量中包含了对生存时间具有影响的因素。

然后查看每个自变量的p值，可以看到sex和ph.ecog这两个变量的p值小于0.05，而age的p值大于0.05, 说明sex和ph,ecog这两个变量对生存时间的影响更加显著。

最后查看自变量的coef等指标，coef就是偏回归系数，exp(coef)就是HR。sex的HR值小于1，该数据集中1=male, 2= female, HR表示的是数值大的风险/数值小的风险，在这里就是female/ male， 说明female死亡的相对较低。HR的值约为0.58, 说明female的死亡风险只占了male的58%， 相比male, female的死亡风险降低了42%。ph.ecog的HR值大于1， 说明随着ph.ecog数值的增加，死亡风险会增加。

3. 结果可视化

以sex为例， 可视化的代码如下

效果图如下