Cox模型是干什么的？
Cox模型的目的是同时评估几个因素对生存的影响。换句话说，它使我们能够检查特定因素在特定时间点如何影响特定事件（例如，感染，死亡）的发生率。这个速度通常被称为危险率。预测变量（或因子）在生存分析文献中通常被称为协变量。
Cox模型由h（t）表示的危险函数表示。简而言之，危险函数可以解释为在时间t死亡的风险。可以估计如下：
h(t)=h0(t)×exp(b1 x1+b2x2+...+bpxp)
- t代表生存时间
- h(t) is the hazard function determined by a set of p covariates (x1,x2,...,xp)
换言之，高于1的风险比指示与事件概率正相关的协变量，并因此与生存期的长短负相关。
综上所述，
HR = 1：没有效果
HR<1：减少危险
HR> 1：危害增加
请注意，在癌症研究中：
危险比(hazard ratio) > 1（即：b> 0）的协变量被称为不良预后因素
危险比(hazard ratio) <1（即：b <0）的协变量被称为良好的预后因子
Cox模型的一个关键假设是，观察组（或患者）的危险曲线应该是成比例的并且不能交叉。
因此，考克斯模型是一个比例 - 危险模型：任何一组中事件的危险性是其他危险的常数倍数。这一假设意味着，如上所述，各组的危险曲线应该是成比例的，不能交叉。
换句话说，如果一个人在某个初始时间点有死亡风险，是另一个人的两倍，那么在所有的晚些时候，死亡风险仍然是两倍。
怎么用R来实现Cox回归模型？
没有包先装包：
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/survminer", build_vignettes = TRUE)
装完包，载入包
library("survival")
library("survminer")
cox模型用到的function就是：coxph()[in survival package]，主要语法结构如下：
coxph(formula, data, method)
载入数据
这两个包里自带了示例数据，载入示例数据演示一下：
data("lung")
head(lung)
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
4 5 210 2 57 1 1 90 60 1150 11
5 1 883 2 60 1 0 100 90 NA 0
6 12 1022 1 74 1 1 50 80 513 0
- inst: Institution code
- time: Survival time in days
- status: censoring status 1=censored, 2=dead
- age: Age in years
- sex: Male=1 Female=2
- ph.ecog: ECOG performance score (0=good 5=dead)
- ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
- pat.karno: Karnofsky performance score as rated by patient
- meal.cal: Calories consumed at meals
- wt.loss: Weight loss in last six months
Cox模型可以分为单变量和多变量，单变量就是单独看每一个变量的关联。
单变量Cox回归 Univariate Cox regression
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox
Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)
coef exp(coef) se(coef) z p
sex -0.531 0.588 0.167 -3.18 0.0015
Likelihood ratio test=10.6 on 1 df, p=0.00111
n= 228, number of events= 165
The function summary() for Cox models produces a more complete report:
summary(res.cox)
# summary() 结果：
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
Cox回归结果可以解释如下：
-
统计显着性。标记为“z”的列给出了Wald统计值。它对应于每个回归系数与其标准误差的比率（z = coef / se（coef））。wald统计评估是否beta（β）系数在统计上显着不同于0。从上面的输出，我们可以得出结论，变量性别具有高度统计学意义的系数。
-
回归系数（The regression coefficients）。Cox模型结果中要注意的第二个特征是回归系数（coef）的符号。正值（+）意味着危险（死亡风险）较高，因此预后更差。1：男，2：女。Cox模型的R总结给出了第二组相对于第一组，即女性与男性的风险比（HR）。性别的β系数= -0.53表明在这些数据中，女性的死亡风险（低存活率）低于男性。
-
危害比例（Hazard ratios）。指数系数（exp（coef）= exp（-0.53）= 0.59）也称为风险比，给出协变量的效应大小。例如，女性（性别= 2）将危害降低了0.59倍，即41％。女性与预后良好相关。
-
风险比的置信区间。总结结果还给出了风险比（exp（coef））的95％置信区间的上限和下限，下限95％界限= 0.4237，上限95％界限= 0.816。
-
模型的全局显著性。最后，输出为模型的总体显着性提供了三个替代测试的p值：likelihood-ratio，Wald测试和得分logrank统计。这三种方法是渐近等价的。对于足够大的N，他们会得到相似的结果。对于小N来说，它们可能有所不同。似然比检验对于小样本量具有更好的表现，所以通常是优选的。
有多个单因素，想一次性跑完？没问题，看这里：
# 把你想跑的单因素都放到这个vector里
covariates <- c("age", "sex", "ph.karno", "ph.ecog", "wt.loss")
# 用sapply创建一个循环，把每个因素的公式写好
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(time, status)~', x)))
# 循环运行上一步写好的所有公式
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
# Extract data
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
p.value<-signif(x$wald["pvalue"], digits=2)
wald.test<-signif(x$wald["test"], digits=2)
beta<-signif(x$coef[1], digits=2);#coeficient beta
HR <-signif(x$coef[2], digits=2);#exp(beta)
HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
HR <- paste0(HR, " (",
HR.confint.lower, "-", HR.confint.upper, ")")
res<-c(beta, HR, wald.test, p.value)
names(res)<-c("beta", "HR (95% CI for HR)", "wald.test",
"p.value")
return(res)
#return(exp(cbind(coef(x),confint(x))))
})
# 这个例子还很贴心的把所有结果给你搞成一个dataframe：
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
看看最后的结果吧，很工整有木有：
beta HR (95% CI for HR) wald.test p.value
age 0.019 1 (1-1) 4.1 0.042
sex -0.53 0.59 (0.42-0.82) 10 0.0015
ph.karno -0.016 0.98 (0.97-1) 7.9 0.005
ph.ecog 0.48 1.6 (1.3-2) 18 2.7e-05
wt.loss 0.0013 1 (0.99-1) 0.05 0.83
多变量Cox回归分析 Multivariate Cox regression analysis
现在，我们想描述这些因素如何共同影响生存。 为了回答这个问题，我们将进行多变量Cox回归分析。 由于变量ph.karno在单变量Cox分析中不显着，我们将在多变量分析中跳过它。 我们将3个因素（性别，年龄和ph.ecog）纳入多变量模型。
res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
summary(res.cox)
#summary()结果：
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
n= 227, number of events= 164
(1 observation deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.011067 1.011128 0.009267 1.194 0.232416
sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***
ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0111 0.9890 0.9929 1.0297
sex 0.5754 1.7378 0.4142 0.7994
ph.ecog 1.5900 0.6289 1.2727 1.9864
Concordance= 0.637 (se = 0.026 )
Rsquare= 0.126 (max possible= 0.999 )
Likelihood ratio test= 30.5 on 3 df, p=1.083e-06
Wald test = 29.93 on 3 df, p=1.428e-06
Score (logrank) test = 30.5 on 3 df, p=1.083e-06
- 所有三个总体测试（可能性，Wald和得分）的p值都很显着，表明该模型是显着的。这些测试评估了所有beta（ββ）为0的综合零假设。在上面的例子中，测试统计数据非常接近，并且完全无效的假设被完全拒绝。
- 在多变量Cox分析中，协变量性别和ph.ecog仍然显着（p <0.05）。然而，协变量年龄不显着（p = 0.23，其大于0.05）。
性别的p值为0.000986，风险比HR = exp（coef）= 0.58，表明患者性别与死亡风险降低之间存在密切关系。协变量的风险比可解释为对危害的乘法效应。例如，保持其他协变量不变，为女性（性别= 2）可将危险降低0.58倍，即42％。我们得出结论，女性与预后良好有关。 - 类似地，ph.ecog的p值为4.45e-05，风险比HR = 1.59，表明ph.ecog值与死亡风险增加之间存在密切关系。保持其他协变量不变，较高的ph.ecog值与较差的存活率相关。
- 相比之下，年龄的p值现在为p = 0.23。风险比HR = exp（coef）= 1.01,95％置信区间为0.99至1.03。因为HR的置信区间包括1，所以这些结果表明，在调整ph.ecog值和患者的性别之后，年龄对HR的差异做出较小的贡献，并且仅趋向于显着性。例如，保持其他协变量不变，另外一年的年龄会导致每日死亡危险因子为exp（beta）= 1.01或1％，这不是一个重要的贡献。
森林图作图
分析完了，来张图可视化一下吧
ggforest(res.cox,main="hazard ratio",cpositions=c(0.02,0.22,0.4),fontsize=0.8,refLabel="reference",noDigits=2)
forest plot
survival analysis 可视化
# 建立一个表达式
fit <- survfit(Surv(time, status) ~ ph.ecog, data = lung)
# 画图
ggsurvplot(fit)
生存曲线
想看更多请关注公众号：bioinfo-c
想看更多请关注
网友评论