前言
在临床研究中，数据建模是最基本但却最容易被忽视的一环。很多研究者在分析数据时直接套用模型，却忽视了模型背后的理论基础。实际上，只有从理论到实践全面理解建模过程中的模型设置，才能更好地解读结果。今天，我们将以Cox比例风险模型（CoxPH）为例，详细介绍其实现和可视化方法。
一、理论简介
Cox比例风险模型是一种经典的生存分析方法，主要用于研究多个因素对时间相关事件（如疾病生存时间）的影响。它是一种半参数模型，既不完全依赖特定的分布假设，又能高效地分析高维数据。
模型的核心特点在于假设不同个体的风险比随时间保持恒定，这就是所谓的比例风险假设。它允许研究者引入多个协变量，并评估它们对风险的相对影响。在医疗、公共卫生等领域，这一模型被广泛应用于评估治疗效果、分析疾病进展以及预测患者生存率。
二、代码实现
我们将使用R语言中的survival和survminer包来实现CoxPH模型。以下是完整的代码实现步骤。
1. 数据加载与预处理
我们使用R自带的lung数据集进行演示：
# 安装或加载必要包
library(survival)
library(survminer)
# 查看数据
data(lung)
lung$sex_label <- ifelse(lung$sex == 1, "Male", "Female") # 添加性别标签
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
## sex_label
## 1 Male
## 2 Male
## 3 Male
## 4 Male
## 5 Male
## 6 Male
2. 构建CoxPH模型
我们构建一个包含多个协变量的CoxPH模型，并查看其结果：
# 构建CoxPH模型
cox_model <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
# 查看模型结果
summary(cox_model)
## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
##
## n= 227, number of events= 164
## (因为不存在，1个观察量被删除了)
##
## 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.025 )
## 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
输出结果包含协变量的回归系数、显著性检验结果以及模型整体的拟合情况。
3. 检验比例假设
比例假设是CoxPH模型的重要前提。我们使用cox.zph函数检验：
# 检验比例假设
cox.zph(cox_model)
## chisq df p
## age 0.188 1 0.66
## sex 2.305 1 0.13
## ph.ecog 2.054 1 0.15
## GLOBAL 4.464 3 0.22
若输出中的p > 0.05，则模型符合比例假设。
三、高质量可视化
为满足高质量发表需求，我们对图表的设计进行了优化，使其具备类似于Nature杂志的配色和字体格式。
1. 基准生存曲线
基准生存曲线展示了所有患者的总体生存概率：
# 基准生存曲线
surv_plot <- ggsurvplot(survfit(cox_model),
data = lung,
conf.int = TRUE,
ggtheme = theme_classic(base_size = 14, base_family = "serif") +
theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 14)),
palette = c("#4E79A7", "#F28E2B"),
title = "基准生存曲线",
xlab = "时间（天）",
ylab = "生存概率")
print(surv_plot)
2. 按性别分组的生存曲线
进一步分析性别对生存率的影响：
# 性别分组生存曲线
ggsurvplot(survfit(Surv(time, status) ~ sex_label, data = lung),
conf.int = TRUE,
pval = TRUE,
ggtheme = theme_classic(base_size = 14, base_family = "serif") +
theme(plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 14)),
palette = c("#4E79A7", "#F28E2B"),
title = "按性别分组的生存曲线",
xlab = "时间（天）",
ylab = "生存概率")
3. 协变量效应森林图
协变量效应森林图帮助我们直观比较不同变量对风险的影响：
# 协变量效应森林图
# 协变量效应森林图
ggforest(cox_model,
data = lung,
main = "协变量效应森林图",
fontsize = 1) # fontsize 控制整个图表的文字大小
四、总结
希望今天的更新让大家详细了解Cox比例风险模型从理论到实践的完整实现流程，特别是学会必要的可视化代码提高结果的美观。欢迎大家在自己的数据上练习本文的代码！
(代码定制，商务咨询等，可加Michael，注明来意)
感谢关注，你的支持是我不懈的动力！
