R语言之生信⑧Cox比例风险模型(多因素)
目录
R语言之生信①差异基因分析1
R语言之生信②差异基因分析2
R语言之生信③差异基因分析3
R语言之生信④TCGA生存分析1
R语言之生信⑤TCGA生存分析2
R语言之生信⑥TCGA生存分析3
R语言之生信⑦Cox比例风险模型(单因素)
R语言之生信⑧Cox比例风险模型(多因素)
======================================
正文
现在，我们想描述这些因素如何共同影响生存。 为了回答这个问题，我们将进行多变量Cox回归分析。 由于变量ph.karno在单变量Cox分析中不显着，我们将在多变量分析中跳过它。 我们将3个因素（性别，年龄和ph.ecog）纳入多变量模型。
多变量Cox回归分析
时间常数协变量的时间到死亡的Cox回归规定如下：
> res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
> summary(res.cox)
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=1e-06
Wald test = 29.93 on 3 df, p=1e-06
Score (logrank) test = 30.5 on 3 df, p=1e-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％，这不是一个重要的贡献。
可视化估计的生存时间分布
将Cox模型拟合到数据后，可以在特定风险组的任何给定时间点可视化预测的存活比例。 函数survfit（）估计生存比例，默认为协变量的平均值。
# Create the new data
sex_df <- with(lung,
data.frame(sex = c(1, 2),
age = rep(mean(age, na.rm = TRUE), 2),
ph.ecog = c(1, 1)
)
)
sex_df
sex age ph.ecog
1 1 62.44737 1
2 2 62.44737 1
绘制图片
fit <- survfit(Surv(time, status) ~sex, data = lung)
ggsurvplot(fit, conf.int = TRUE, legend.labs=c("Sex=1", "Sex=2"),
ggtheme = theme_minimal())
