之前文章介绍了Kaplan-Meier生存曲线分析，Kaplan-Meier模型除了展示预后状况，也可以用log-rank法检测是否分组预后有显著差异。cox比例风险模型则适合衡量具体某一因素对生存的影响程度，用HR(hazard ratio)值体现，HR是某一因素影响生存的比率。cox模型公式如下。
-
t- 时间 -
h(t)- 该时间点死亡风险 -
x- 因素 -
- 某因素HR
HR值对应含义如下
- HR = 1: 无影响
- HR < 1: 降低风险
- HR > 1: 增加风险
不过我们不只看HR值，还要看95%CI即95%置信度区间，如果95%CI跨越了1，一般就不认为该因素对生存有显著影响。另外要提示的是HR值受输入值的规模(身高用CM为单位还是M为单位)影响，所以如果得到非常巨大/小的HR结果，要思考自己数据缩放问题。
下面用 lung 数据集分别展示单因素、多因素cox分析。
> library(survival, quietly = TRUE)
> library(survminer, quietly = TRUE)
> head(lung, n = 3)
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1 3 306 2 74 1 1 90 100 1175 NA
2 3 455 2 68 1 0 90 90 1225 15
3 3 1010 1 56 1 0 90 90 NA 15
使用 coxph 函数进行cox回归分析。我们看性别对肺癌预后有多大影响。
> reg <- coxph(Surv(time, status) ~ sex, data = lung)
> summary(reg)
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.021 )
Likelihood ratio test= 10.63 on 1 df, p=0.001
Wald test = 10.09 on 1 df, p=0.001
Score (logrank) test = 10.33 on 1 df, p=0.001
用 summary 函数查看结果，其中 coef 是系数beta， exp(coef) 就是HR值，在这里是 0.5880，95%CI是 0.4237 ~ 0.816. P值给3个，如果样本少适合"Likelihood ratio test"，如果样本量大，3个方法P值不会差异太多。
如果有多个因素的数据，进行多因素cox回归是适合的。也同样使用 coxph 和 summary 函数。
> reg <- coxph(Surv(time, status) ~ sex + age + `ph.ecog` + `ph.karno`, data = lung)
> summary(reg)
Call:
coxph(formula = Surv(time, status) ~ sex + age + ph.ecog + ph.karno,
data = lung)
n= 226, number of events= 163
(2 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
sex -0.572802 0.563943 0.169222 -3.385 0.000712 ***
age 0.012868 1.012951 0.009404 1.368 0.171226
ph.ecog 0.633077 1.883397 0.176034 3.596 0.000323 ***
ph.karno 0.012558 1.012637 0.009514 1.320 0.186842
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
sex 0.5639 1.7732 0.4048 0.7857
age 1.0130 0.9872 0.9945 1.0318
ph.ecog 1.8834 0.5310 1.3338 2.6594
ph.karno 1.0126 0.9875 0.9939 1.0317
Concordance= 0.632 (se = 0.025 )
Likelihood ratio test= 31.27 on 4 df, p=3e-06
Wald test = 30.61 on 4 df, p=4e-06
Score (logrank) test = 31.06 on 4 df, p=3e-06
可以看到性别因素sex是显著影响生存预后的，但是年龄因素age(HR:1.013, P:0.171, 95%CI:0.9945~1.0318)不是。然后使用 ggforest 函数直接画出森林图。
> survminer::ggforest(reg, data = lung)
图片如下，效果有点复古。也可以自己手动画森林图，95%CI用线段(geom_segment)表示，HR值用点(geom_point)然后颜色表示P值即可。
参考
Cox Proportional-Hazards Model - Easy Guides - Wiki - STHDA
Drawing Survival Curves using 'ggplot2' • survminer
