【8.3】Cox模型
一、方法说明
Kaplan-Meier方法只能针对单一的变量进行分析，无法同时考察多个因素。当需要同时考察多个因素的影响时，这时我们可以使用Cox比例风险回归模型。
Cox比例风险回归模型(Cox’s proportional hazards regression model)，简称Cox模型，Cox来自提出者英国统计学家D.R.Cox的名字，主要用于肿瘤和其他疾病的预后分析。这个模型是一种半参数回归模型，因为它的公式中既包含参数模型又包含非参数模型。
其中
- t是生存时间，
- x1, x2到xp指的是具有预测效应的多个变量，
- b1,b2到bp则是每个变量对应的effect size，即效应量，可以理解为结果的影响程度。
- h(t)就是不同时间t的 hazard，即风险值，例如在观测死亡事件时，指的是研究对象从试验开始到某个特定时间t之前存活，但在t时间点发生死亡的概率。
- h0(t)是基准风险函数，也就是说在其他协变量x1, x2到xp都为0时，即不起作用时，衡量风险值的函数。
根据公式我们可以看到指数部分是参数模型，因为其参数个数有限，即b1,b2到bp，而基准风险函数h0(t)由于于其未确定性，可根据不同数据来使用不同的分布模型，因此是非参数模型。所以说, Cox模型是一种半参数模型。
从公式中我们可以看到，Cox模型能够把诸多可能影响生存率的因素都当作协变量引入到公式中去，在该公式中即x1, x2到xp，所以可以同时考察多个因素的影响。
我们的主要目标是通过一定方法来找到合适的h0(t)，以及所有协变量的系数b1,b2到bp。实际上cox模型是需要用到极大似然估计等计算方法，首先构建特定的似然函数，通过梯度下降等方法来求解模型的参数，使得函数求解值最大，这里不对细节进行解读。
假设我们已经通过计算得到了合适的h0(t)和协变量系数，如何去解读结果呢？我们可以比较某个协变量x1 在不同值时对应的不同风险比(hazard ratio)，这里 x1和x1+1，即若增加1个单位，增加前后的风险比实际上等于 exp(b1)。
假如x1指的是年龄，那么对于年龄 51岁 (x+1) 和年龄 50 岁 (x) 的人，可能死亡的风险比为 exp(b1）。如果b1>0，则 exp(b1)>1，意味着年龄+1，死亡风险增加；如果b1<0, 则 exp(b1)<1，意味着年龄+1，死亡风险降低；如果b1=0，exp(b1)=1，意味着年龄变化对死亡风险不起作用。从hazard ratio推导的结果看到，它是不包括时间t的。这是Cox模型可用的一个基本假设，即任意两人的风险比例是不随时间变化的。
二、计算与解读
研究者开发了方便进行生存分析的R包，survival和survminer。首先安装并加载这两个包：
install.packages(c("survival","survminer"))
library("survival")
library("survminer")
在survival包中提供了coxph()函数可以用来计算cox模型：
coxph(formula, data, method)
method默认为 “efron”，也可以是 “breslow”和“exact” 。以示例数据为例：
data("lung")
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
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
从结果中看到：sex对应的系数(coef)为-0.5310，小于0表示sex增加会降低风险，风险比(hazard ratio)为exp(coef) =0.588，该数值小于1，同样表明sex增加会导致风险增加，即女性比男性预后更好。
除了关注系数外，同时需要关注的是p value，即该参数估计是否具有统计学显著性，这里给出三种方法的结果，分别是Likelihood ratio test，Wald test和Score logrank test。
分析多个因素的影响：
res.cox <- coxph(Surv(time, status) ~ age+sex+ph.ecog, data = lung)
summary(res.cox)
最后是结果的可视化：
fit <- survfit(res.cox, data = lung)
ggsurvplot(fit, conf.int = TRUE, legend.labs=c("Sex=1", "Sex=2"),
ggtheme = theme_minimal())
以上是对生存分析中主要知识的一个整理，希望梳理清楚生存分析中的大多数概念，有助于大家在自己的工作中使用相关方法进行分析
参考资料
