R语言绘制Cox回归森林图：单因素/多因素分析与双置信区间可视化

这次我们来看一个在医学统计和生存分析中非常实用的技能：如何绘制 Cox 回归的单因素和多因素分析森林图，并且同时展示双置信区间。对于从事临床研究、流行病学或数据分析的朋友来说，这几乎是论文和报告中的“标配”图表。它不仅能直观展示每个变量的风险比，还能清晰呈现其统计显著性，是评估预后因素的核心工具。

这个项目的重点不是教你复杂的统计学理论，而是直接解决“怎么画出来”的问题。我们将聚焦于使用 R 语言中的survival和forestplot等核心包，从数据整理、模型构建到图形美化，一步步实现一个专业级的森林图。整个过程不依赖商业软件，完全开源可复现，适合需要批量分析或集成到自动化报告流程中的场景。

本文会带你完成从零到一的完整流程：首先准备一个模拟的生存数据集，然后分别进行单因素和多因素 Cox 回归分析，最后使用forestplot包绘制包含双 CI（例如 95% 和 99% 置信区间）的森林图。我们会重点关注代码的实用性、图形的可定制性，以及如何避免常见的绘图错误，确保你读完就能在自己的数据上跑通。

1. 核心能力速览

2. 适用场景与使用边界

这个工具适合谁？

• 临床研究人员与医学生：需要分析患者生存数据，寻找独立的预后因素，并在论文中展示规范的结果。
• 流行病学与公共卫生分析师：从事队列研究，评估多种暴露因素与结局事件的风险关联。
• 数据科学家与统计顾问：为客户提供生存分析服务，需要生成直观、专业的可视化报告。
• 高校教师与学生：用于统计学、医学统计学等课程的教学案例与实操练习。

能解决什么问题？

1. 可视化风险比较：将 Cox 回归得到的风险比及其置信区间以图形方式呈现，一目了然。
2. 区分单/多因素结果：清晰对比仅考虑单一变量和校正了混杂因素后，变量效应的变化。
3. 展示统计不确定性：通过置信区间的长短和是否跨越无效线（HR=1），判断结果的精确性与显著性。
4. 提升报告效率：通过脚本化绘图，避免手动在图形软件中拼接，便于重复分析和更新。

不适合什么场景？

• 数据不符合 Cox 比例风险假设时（需先进行检验）。
• 样本量过小或事件数过少，导致模型不稳定或无法拟合。
• 需要绘制非常复杂的、包含多层亚组分析的森林图时，可能需要更专业的图形包或手动调整。

合规与伦理边界：

• 所使用的生存数据必须符合伦理审查要求，获得必要的知情同意，并经过脱敏处理。
• 在发表或报告结果时，需对统计方法进行准确描述，避免误导性解读。
• 森林图是统计结果的展示，其解释应结合专业背景知识，不能仅凭图形下结论。

3. 环境准备与前置条件

为了顺利运行后续代码，你需要准备好以下环境：

1. 操作系统：Windows, macOS 或 Linux 均可。
2. R 环境：安装最新版的 R 语言环境。建议版本 >= 4.0.0。
3. R 集成开发环境：推荐使用 RStudio，它提供了友好的代码编辑、运行和图形查看界面。
4. 必要的 R 包：我们将主要使用survival和forestplot包。survival是生存分析的标准包，forestplot在绘制森林图方面功能强大且灵活。

安装与检查步骤：打开 R 或 RStudio，在控制台执行以下命令来安装和加载必要的包。

# 安装必需的包（如果尚未安装） install.packages(c("survival", "forestplot", "dplyr", "tidyr")) # 加载包到当前会话 library(survival) library(forestplot) library(dplyr) # 用于数据整理 library(tidyr) # 用于数据整理

确保安装过程没有报错。forestplot包是绘制我们目标图形的关键。

4. 数据准备与模拟数据集

在实战前，我们首先创建一个模拟的生存数据集。这个数据集将包含生存时间、生存状态以及几个常见的临床协变量，例如年龄、性别、肿瘤分期等。

# 设置随机种子以保证结果可重复 set.seed(123) # 生成模拟数据 n <- 200 # 样本量 patient_id <- 1:n age <- round(rnorm(n, mean=60, sd=10)) sex <- factor(sample(c("Male", "Female"), n, replace=TRUE, prob=c(0.6, 0.4))) stage <- factor(sample(c("I", "II", "III", "IV"), n, replace=TRUE, prob=c(0.2, 0.3, 0.3, 0.2)), levels=c("I", "II", "III", "IV")) treatment <- factor(sample(c("Drug_A", "Drug_B"), n, replace=TRUE, prob=c(0.5, 0.5))) # 生成生存时间（假设基线风险与年龄、分期相关） base_hazard <- 0.01 * exp(0.03*(age-60) + 0.8*(as.numeric(stage)-1)) survival_time <- rexp(n, rate = base_hazard) # 生成删失时间，假设研究随访期为5年 censor_time <- runif(n, 1, 5) # 确定观察到的生存时间和状态 time <- pmin(survival_time, censor_time) status <- as.numeric(survival_time <= censor_time) # 组合成数据框 sim_data <- data.frame( patient_id, age, sex, stage, treatment, time, status ) # 查看数据结构 head(sim_data) str(sim_data)

运行后，你将看到一个名为sim_data的数据框，包含 200 行观测和 7 个变量。这就是我们后续分析的基础。

5. Cox 单因素与多因素回归分析

5.1 单因素 Cox 回归分析

单因素分析是指将每个协变量单独放入 Cox 模型，评估其与生存结局的粗关联。

# 定义我们感兴趣的协变量列表 covariates <- c("age", "sex", "stage", "treatment") # 初始化一个列表来存储每个单因素模型的结果 univ_models <- list() univ_results <- data.frame() # 循环对每个变量进行单因素 Cox 回归 for (covar in covariates) { # 构建公式 formula <- as.formula(paste("Surv(time, status) ~", covar)) # 拟合 Cox 模型 cox_model <- coxph(formula, data = sim_data) # 将模型存入列表 univ_models[[covar]] <- cox_model # 提取模型摘要中的关键信息 model_summary <- summary(cox_model) # 获取风险比、置信区间和P值 hr <- exp(coef(cox_model)) ci_low <- exp(confint(cox_model))[,1] ci_high <- exp(confint(cox_model))[,2] p_value <- coef(model_summary)[,5] # 如果是因子变量，会有多行结果，需要处理变量名 var_levels <- names(hr) if (length(var_levels) > 1) { # 对于因子变量，我们通常以第一级为参照，这里我们保留所有级别 var_name <- rep(covar, length(var_levels)) result_df <- data.frame( variable = var_name, level = var_levels, hr = round(hr, 2), ci_low_95 = round(ci_low, 2), ci_high_95 = round(ci_high, 2), p_value = sprintf("%.3f", p_value) ) } else { # 对于连续变量 result_df <- data.frame( variable = covar, level = NA, hr = round(hr, 2), ci_low_95 = round(ci_low, 2), ci_high_95 = round(ci_high, 2), p_value = sprintf("%.3f", p_value) ) } univ_results <- rbind(univ_results, result_df) } # 查看单因素分析结果 print(univ_results)

这段代码会输出一个表格，包含每个变量（或因子变量的每个水平）的风险比、95%置信区间和P值。

5.2 多因素 Cox 回归分析

多因素分析是将所有感兴趣的协变量同时放入一个模型，以校正它们之间的相互影响，评估每个变量的独立效应。

# 构建多因素模型的公式 # 注意：对于因子变量，R会自动处理哑变量，以第一个水平为参照 multiv_formula <- as.formula("Surv(time, status) ~ age + sex + stage + treatment") # 拟合多因素 Cox 模型 multiv_model <- coxph(multiv_formula, data = sim_data) # 查看模型摘要 multiv_summary <- summary(multiv_model) print(multiv_summary) # 提取多因素分析结果用于绘图 multiv_coef <- coef(multiv_summary) multiv_hr <- exp(multiv_coef[, 1]) multiv_ci_low_95 <- exp(multiv_coef[, 3]) # 默认是95% CI的下限 multiv_ci_high_95 <- exp(multiv_coef[, 4]) # 默认是95% CI的上限 multiv_p_value <- multiv_coef[, 5] # 获取变量名（对于因子变量，会展开为多个哑变量） var_names <- rownames(multiv_coef) # 整理成数据框 multiv_results <- data.frame( variable = var_names, hr = round(multiv_hr, 2), ci_low_95 = round(multiv_ci_low_95, 2), ci_high_95 = round(multiv_ci_high_95, 2), p_value = sprintf("%.3f", multiv_p_value) ) print(multiv_results)

现在，我们得到了单因素和多因素分析的结果数据框，这是绘制森林图的原料。

6. 绘制单因素分析森林图（基础版）

我们先从基础的单因素森林图开始，使用forestplot包。

# 准备森林图所需的数据 # forestplot 需要提供一个矩阵，其中包含要显示的文本标签和效应值 # 1. 准备标签文本 # 第一列通常是变量名和水平 label_text <- cbind( c("Variable", as.character(univ_results$variable)), # 变量名 c("HR (95% CI)", paste0(univ_results$hr, " (", univ_results$ci_low_95, "-", univ_results$ci_high_95, ")")) # HR和CI ) # 2. 准备效应值（HR）和置信区间矩阵 # 矩阵的行数 = 变量数 + 1（表头），列数 = 3（均值，下限，上限） mean_values <- c(NA, univ_results$hr) # 第一行是表头，用NA lower_values <- c(NA, univ_results$ci_low_95) upper_values <- c(NA, univ_results$ci_high_95) # 组合成矩阵 forest_data <- matrix(c(mean_values, lower_values, upper_values), ncol = 3) colnames(forest_data) <- c("mean", "lower", "upper") # 3. 绘制森林图 forestplot(labeltext = label_text, mean = forest_data[, "mean"], lower = forest_data[, "lower"], upper = forest_data[, "upper"], is.summary = c(TRUE, rep(FALSE, nrow(univ_results))), # 第一行是汇总行（表头） xlog = TRUE, # X轴使用对数刻度，因为HR是对称的 title = "Univariable Cox Regression Analysis", xlab = "Hazard Ratio (HR)", boxsize = 0.2, # 中间方块的大小 col = fpColors(box="royalblue", line="darkblue", summary="royalblue"), txt_gp = fpTxtGp(label = gpar(cex=0.8), # 标签字体大小 xlab = gpar(cex=0.9), # X轴标签字体大小 title = gpar(cex=1.1)) # 标题字体大小 )

运行这段代码，你将得到一个基础的单因素森林图。图中每个变量对应一条水平线和一个方块，方块的位置代表点估计值（HR），水平线的长度代表95%置信区间。X轴为对数刻度，竖线为 HR=1 的无效线。

7. 绘制多因素分析森林图并添加双 CI

现在，我们来绘制更高级的多因素森林图，并展示双置信区间（例如，同时显示 95% CI 和 99% CI）。这需要我们先计算 99% 的置信区间。

# 计算多因素模型中各系数的 99% 置信区间 conf_int_99 <- exp(confint(multiv_model, level = 0.99)) multiv_ci_low_99 <- conf_int_99[,1] multiv_ci_high_99 <- conf_int_99[,2] # 更新多因素结果数据框 multiv_results$ci_low_99 <- round(multiv_ci_low_99, 2) multiv_results$ci_high_99 <- round(multiv_ci_high_99, 2) # 准备森林图数据（包含双CI） # 1. 标签文本 label_text_multi <- cbind( c("Variable", multiv_results$variable), c("HR (95% CI)", paste0(multiv_results$hr, " (", multiv_results$ci_low_95, "-", multiv_results$ci_high_95, ")")), c("HR (99% CI)", paste0(multiv_results$hr, " (", multiv_results$ci_low_99, "-", multiv_results$ci_high_99, ")")) ) # 2. 效应值与CI矩阵 # 我们需要为每个CI范围准备一组数据。forestplot可以接受多组CI。 # 第一组：95% CI mean_multi <- c(NA, multiv_results$hr) lower_95_multi <- c(NA, multiv_results$ci_low_95) upper_95_multi <- c(NA, multiv_results$ci_high_95) # 第二组：99% CI lower_99_multi <- c(NA, multiv_results$ci_low_99) upper_99_multi <- c(NA, multiv_results$ci_high_99) # 组合成一个列表，每个元素是一个矩阵（均值，下限，上限） forest_data_multi <- list( matrix(c(mean_multi, lower_95_multi, upper_95_multi), ncol=3), matrix(c(mean_multi, lower_99_multi, upper_99_multi), ncol=3) ) # 3. 绘制带双CI的森林图 forestplot(labeltext = label_text_multi, mean = forest_data_multi[[1]][, 1], # 使用第一组数据的均值 lower = list(forest_data_multi[[1]][, 2], forest_data_multi[[2]][, 2]), # 两组下限 upper = list(forest_data_multi[[1]][, 3], forest_data_multi[[2]][, 3]), # 两组上限 is.summary = c(TRUE, rep(FALSE, nrow(multiv_results))), xlog = TRUE, title = "Multivariable Cox Regression Analysis with 95% & 99% CI", xlab = "Hazard Ratio (HR)", boxsize = 0.25, col = fpColors(box="darkred", lines = c("darkred", "darkgreen"), # 为两条CI线指定不同颜色 summary="darkred"), fn.ci_norm = c(fpDrawNormalCI, fpDrawCircleCI), # 使用不同的图形绘制CI（线和点） lwd.ci = c(2, 1), # 设置两条CI线的粗细 ci.vertices = FALSE, # 不显示CI线两端的顶点 txt_gp = fpTxtGp(label = gpar(cex=0.8), xlab = gpar(cex=0.9), title = gpar(cex=1.1)), legend = list("95% CI" = list(col="darkred", lwd=2), "99% CI" = list(col="darkgreen", lwd=1)), legend_args = fpLegend(pos = list(x=0.85, y=0.95)) # 图例位置 )

在这张图中，你将看到每个变量有两条置信区间线：一条较粗的深红色线代表 95% CI，一条较细的深绿色线代表 99% CI。99% CI 总是比 95% CI 更宽。图例说明了两种颜色的含义。这种展示方式能更细致地呈现估计的不确定性。

8. 合并单因素与多因素结果的森林图

在论文中，经常将单因素和多因素分析的结果并列展示。我们可以通过精心组织数据来实现。

# 目标：创建一个表格，左边是变量名，中间是单因素HR(CI)，右边是多因素HR(CI) # 1. 整理单因素结果（以多因素结果的变量顺序为基准） # 假设我们只关心出现在多因素模型中的这些变量 # 注意：多因素模型中的变量是“sexMale”，“stageII”等，我们需要匹配 univ_for_merge <- univ_results # 简化处理：这里我们创建一个匹配键。实际应用中可能需要更复杂的合并逻辑。 # 为了演示，我们假设单因素结果已按变量名整理好。 # 2. 创建合并的标签和效应值矩阵 # 变量名 all_vars <- multiv_results$variable # 单因素 HR(CI) 字符串 univ_hr_ci <- paste0(univ_for_merge$hr, " (", univ_for_merge$ci_low_95, "-", univ_for_merge$ci_high_95, ")") # 多因素 HR(CI) 字符串 multiv_hr_ci <- paste0(multiv_results$hr, " (", multiv_results$ci_low_95, "-", multiv_results$ci_high_95, ")") label_combined <- cbind( c("Variable", all_vars), c("Univariable HR (95% CI)", c(univ_hr_ci)), c("Multivariable HR (95% CI)", c(multiv_hr_ci)) ) # 3. 创建合并的效应值矩阵（用于绘图） # 我们需要两组数据：单因素的CI和多因素的CI mean_univ <- c(NA, univ_for_merge$hr) lower_univ <- c(NA, univ_for_merge$ci_low_95) upper_univ <- c(NA, univ_for_merge$ci_high_95) mean_multiv <- c(NA, multiv_results$hr) lower_multiv <- c(NA, multiv_results$ci_low_95) upper_multiv <- c(NA, multiv_results$ci_high_95) # 将两组数据放入一个列表 forest_data_combined <- list( matrix(c(mean_univ, lower_univ, upper_univ), ncol=3), matrix(c(mean_multiv, lower_multiv, upper_multiv), ncol=3) ) # 4. 绘制合并森林图 forestplot(labeltext = label_combined, mean = cbind(forest_data_combined[[1]][,1], forest_data_combined[[2]][,1]), # 两列均值 lower = cbind(forest_data_combined[[1]][,2], forest_data_combined[[2]][,2]), # 两列下限 upper = cbind(forest_data_combined[[1]][,3], forest_data_combined[[2]][,3]), # 两列上限 is.summary = c(TRUE, rep(FALSE, length(all_vars))), xlog = TRUE, title = "Comparison: Univariable vs. Multivariable Cox Regression", xlab = "Hazard Ratio (HR)", boxsize = 0.2, col = fpColors(box=c("blue", "red"), # 单因素蓝色，多因素红色 lines = c("blue", "red"), summary=c("blue", "red")), fn.ci_norm = c(fpDrawNormalCI, fpDrawNormalCI), lwd.ci = c(1.5, 1.5), ci.vertices = FALSE, txt_gp = fpTxtGp(label = gpar(cex=0.75), xlab = gpar(cex=0.9), title = gpar(cex=1.1)), legend = list("Univariable" = list(col="blue", lwd=1.5), "Multivariable" = list(col="red", lwd=1.5)), legend_args = fpLegend(pos = list(x=0.85, y=0.95)) )

这张图将单因素和多因素分析的结果并排显示，方便直接比较校正前后风险比的变化。蓝色代表单因素分析，红色代表多因素分析。

9. 图形美化与高级定制

forestplot包提供了极高的定制自由度。以下是一些常见的美化技巧：

# 示例：创建一个更美观、更接近发表要求的单因素森林图 # 假设我们使用最初的单因素结果 `univ_results` # 1. 准备更精美的标签（包含P值） label_pretty <- cbind( c("Variable", paste(univ_results$variable, ifelse(is.na(univ_results$level), "", paste0(" (", univ_results$level, ")")))), c("HR", univ_results$hr), c("95% CI", paste0("[", univ_results$ci_low_95, ", ", univ_results$ci_high_95, "]")), c("P Value", univ_results$p_value) ) # 2. 绘制图形 forestplot(labeltext = label_pretty, mean = c(NA, univ_results$hr), lower = c(NA, univ_results$ci_low_95), upper = c(NA, univ_results$ci_high_95), is.summary = c(TRUE, rep(FALSE, nrow(univ_results))), graph.pos = 3, # 将森林图放在第3列（标签后） xlog = TRUE, title = "Univariable Cox Regression Analysis (Formatted)", xlab = "Hazard Ratio", colgap = unit(5, "mm"), # 列间距 lineheight = unit(0.7, "cm"), # 行高 boxsize = 0.25, line.margin = unit(0.1, "cm"), col = fpColors(box="steelblue", line="steelblue", summary="royalblue", hrz_lines = "#444444"), # 水平分隔线颜色 lwd.ci = 2, ci.vertices = TRUE, # 显示CI线两端的短横线 ci.vertices.height = 0.1, # 短横线高度 txt_gp = fpTxtGp(label = gpar(fontfamily = "sans", cex=0.9), ticks = gpar(cex=0.8), xlab = gpar(cex=0.9, fontface="bold"), title = gpar(cex=1.1, fontface="bold")), hrzl_lines = list("2" = gpar(lwd=1, lty=2, col="gray")), # 在第2行后加虚线 grid = TRUE # 添加垂直网格线 )

通过调整graph.pos,colgap,lineheight,col,hrzl_lines,grid等参数，你可以控制图表的几乎所有视觉元素，使其完全符合目标期刊或报告的格式要求。

10. 常见问题与排查方法

在绘制森林图的过程中，你可能会遇到以下问题：

11. 最佳实践与使用建议

1. 数据清洗先行：在运行 Cox 回归前，务必处理缺失值、检查变量类型（数值型、因子型）。对于因子变量，使用factor()并设置正确的参照水平。
2. 模型假设检验：使用cox.zph()函数检验比例风险假设。如果假设被严重违反，考虑使用参数模型或分层 Cox 模型。
3. 结果备份：将coxph()模型对象和summary()结果保存为.RData文件，便于后续复查或绘制其他图形。

   save(multiv_model, multiv_summary, file = “cox_multiv_model_results.RData”)

4. 脚本模块化：将数据准备、模型拟合、结果提取和图形绘制写成独立的函数或 R 脚本文件。这样便于对不同数据集进行批量分析。
5. 图形输出设置：对于发表或报告，使用高分辨率输出。

   png(“my_forestplot.png”, width=3200, height=2400, res=300) # 高分辨率PNG # 在这里运行你的 forestplot() 代码 dev.off()

6. 版本控制：使用 Git 等工具管理你的分析脚本，特别是当分析需要多次迭代或合作完成时。
7. 解释与报告：在报告中，不仅要提供森林图，还应附上相应的表格，列出具体的 HR、CI 和 P 值。对图中任何重要的发现（如 CI 跨越 1，或校正前后 HR 方向改变）进行文字解释。

掌握 Cox 回归森林图的绘制，能极大提升你处理生存分析数据的效率和成果展示的专业性。从单因素到多因素，从基础图形到包含双 CI 的复杂图表，整个流程都可以通过 R 脚本自动化完成。建议你从本文的模拟数据代码开始，逐步替换成自己的真实数据，并尝试调整各种图形参数，直到生成完全符合你需求的图表。