R语言绘制Cox回归双CI森林图：单因素与多因素结果可视化

如果你在医学、公共卫生或生物统计领域工作一定会遇到这样的场景你完成了一项队列研究通过 Cox 回归模型分析了一组变量对生存结局的影响得到了单因素和多因素分析的结果。现在你需要向导师、审稿人或合作者展示这些结果。你面对的是一张密密麻麻的表格里面塞满了 HR、95% CI 和 P 值。如何让这些冰冷的数字变得直观、有力一眼就能看出哪些因素是保护性的哪些是危险因素并且同时展示单因素和多因素的结果对比答案就是Cox 回归森林图。它不仅是论文图表中的“颜值担当”更是高效传递信息的利器。一张好的森林图能让读者在几秒钟内抓住研究的核心发现。然而从原始数据到一张发表级的森林图中间隔着不少“坑”如何用 R 语言一次性提取单因素和多因素分析结果如何将两个模型的结果优雅地整合到同一张图上如何自定义坐标轴、字体、颜色以满足不同期刊的挑剔要求网上教程虽多但往往只讲单因素图或者代码复杂晦涩难以直接套用。本文将彻底解决这个问题。我将手把手带你使用 R 语言中强大的forestplot包绘制一张同时包含单因素和多因素 Cox 回归结果的森林图。我们不止步于“画出来”更要深入“为什么这么画”并分享一套可直接复用的、模块化的代码模板。无论你是 R 语言新手还是希望优化科研工作流的老手这篇文章都能让你在生存分析的可视化上迈出坚实的一步。1. 为什么你需要掌握“双CI”森林图在深入代码之前我们首先要理解为什么这种同时展示单因素和多因素结果的森林图如此重要。单因素分析 vs. 多因素分析单因素分析是“初筛”它单独考察每个变量与结局的关系但结果可能受到其他混杂因素的干扰。多因素分析则是“精炼”它把所有重要变量放在同一个模型中评估每个变量的“独立贡献”。审稿人最关心的往往是多因素分析的结果。传统做法的局限很多研究者会分别绘制单因素和多因素森林图或者只在表格中并列展示。前者浪费空间且不便于对比后者则不够直观无法快速识别效应值的变化趋势例如某个变量从单因素显著变为多因素不显著这本身就是一个有趣的故事。“双CI”森林图的优势将单因素和多因素的 HR 及其置信区间并排放在同一行用不同的符号如方块和菱形或颜色区分。这种设计让你可以一目了然地对比直接看到校正混杂因素前后每个变量效应值HR和精确度CI宽度的变化。高效利用空间在一张图上呈现所有信息符合顶级期刊对图表简洁高效的要求。讲述数据故事例如一个变量单因素分析时 HR2.5 (1.3-4.8)多因素分析后变为 HR1.2 (0.8-1.8)这强烈暗示其效应可能被其他变量所解释或混淆这本身就是分析讨论的一部分。因此掌握这项技能不是简单的“美化图表”而是提升你科研结果呈现专业度和洞察深度的关键一步。2. 核心工具与概念梳理在开始动手前我们需要明确几个核心概念和将要用到的工具。Cox 比例风险模型生存分析中最常用的回归模型用于分析一个或多个变量协变量对生存时间的影响。其结果以风险比Hazard Ratio, HR呈现。HR 1 表示该变量是危险因素增加事件发生风险HR 1 表示保护因素。森林图一种用于展示多个独立研究结果或同一研究中多个变量效应估计值如 HR, OR, RR及其置信区间的图表。每个变量占一行用一条水平线段置信区间和一个点点估计值如 HR表示。forestplot包R 语言中绘制发表级森林图的利器。相比基础的forestplot函数或survminer包forestplot包提供了无与伦比的灵活性可以轻松构建复杂的、多列的、高度自定义的森林图完美契合我们绘制“双CI”图的需求。本次任务的数据流数据准备一个包含生存时间、生存状态和若干协变量的数据集。模型拟合分别对每个变量进行单因素 Cox 回归再对所有变量进行多因素 Cox 回归。结果提取从模型结果中提取变量名、HR、置信区间和 P 值并整理成规整的数据框。绘图数据构建构建一个用于forestplot函数的列表结构包含文本标签、效应值估计和置信区间。图形绘制与美化使用forestplot函数绘图并调整所有视觉元素。接下来我们将按照这个流程一步步实现。3. 环境准备安装与加载必要的R包确保你的 R 环境已经就绪。我们将主要依赖以下包survival: 用于拟合 Cox 回归模型的核心包。forestplot: 用于绘制高度自定义的森林图。dplyr/tidyverse: 用于数据清洗和整理让代码更简洁。这里我们以dplyr为例。如果你还没有安装这些包请运行以下代码# 安装必要包 install.packages(c(survival, forestplot, dplyr))安装完成后在脚本开头加载它们# 加载包 library(survival) library(forestplot) library(dplyr)4. 从数据到模型拟合单因素与多因素Cox回归我们使用 R 内置的lung数据集肺癌患者数据进行演示。这个数据集包含生存时间、状态以及年龄、性别、体能评分等变量。# 加载示例数据 data(lung) # 查看数据结构 head(lung) # 处理数据将status转换为0/1格式生存包要求并处理缺失值 lung_data - lung %% mutate(status ifelse(status 2, 1, 0)) %% # 假设status2为死亡事件 select(time, status, age, sex, ph.ecog) %% # 选择要分析的变量 na.omit() # 删除含有缺失值的行实际分析中需谨慎处理缺失值 # 定义要分析的协变量列表 covariates - c(age, sex, ph.ecog)第一步批量进行单因素 Cox 回归分析我们不希望为每个变量单独写一遍coxph函数而是用循环或lapply函数高效完成。# 创建一个空列表来存储单因素模型结果 uni_models - list() # 循环拟合单因素Cox模型 for (covar in covariates) { formula - as.formula(paste(Surv(time, status) ~, covar)) uni_models[[covar]] - coxph(formula, data lung_data) }第二步进行多因素 Cox 回归分析将所有的协变量同时放入一个模型。# 拟合多因素Cox模型 multi_formula - as.formula(paste(Surv(time, status) ~, paste(covariates, collapse ))) multi_model - coxph(multi_formula, data lung_data)5. 结果提取与整理构建绘图数据框这是最关键的一步。我们需要从上面拟合的模型中提取出整齐的数据以便输入给forestplot函数。# 1. 提取单因素分析结果并整理 uni_results - lapply(names(uni_models), function(covar) { model - uni_models[[covar]] sum_model - summary(model) # 提取HR, 95% CI, P值 hr - round(sum_model$conf.int[, 1], 2) ci_low - round(sum_model$conf.int[, 3], 2) ci_high - round(sum_model$conf.int[, 4], 2) p_val - round(sum_model$coefficients[, 5], 3) # 将P值格式化为科学计数法通常用于小于0.001的值 p_val_format - ifelse(p_val 0.001, 0.001, as.character(p_val)) data.frame( variable covar, hr_uni hr, ci_low_uni ci_low, ci_high_uni ci_high, p_uni p_val_format, stringsAsFactors FALSE ) }) %% bind_rows() # 将列表合并为一个数据框 # 2. 提取多因素分析结果并整理 sum_multi - summary(multi_model) multi_results - data.frame( variable rownames(sum_multi$conf.int), hr_multi round(sum_multi$conf.int[, 1], 2), ci_low_multi round(sum_multi$conf.int[, 3], 2), ci_high_multi round(sum_multi$conf.int[, 4], 2), p_multi round(sum_multi$coefficients[, 5], 3) ) multi_results$p_multi - ifelse(multi_results$p_multi 0.001, 0.001, as.character(multi_results$p_multi)) # 3. 合并单因素和多因素结果 # 注意确保变量顺序一致这里按单因素结果的顺序合并 plot_data - uni_results %% left_join(multi_results, by variable) # 查看合并后的数据 print(plot_data)运行后plot_data数据框应该类似这样variable hr_uni ci_low_uni ci_high_uni p_uni hr_multi ci_low_multi ci_high_multi p_multi 1 age 1.02 1.00 1.04 0.028 1.02 1.00 1.04 0.041 2 sex 0.59 0.42 0.84 0.003 0.66 0.46 0.95 0.025 3 ph.ecog 1.75 1.36 2.26 0.001 1.71 1.32 2.22 0.0016. 构建forestplot所需的输入数据结构forestplot函数需要特定的输入格式一个列表或矩阵其中包含要在图上显示的文本标签、效应值估计和置信区间。# 1. 构建文本标签列表格的左侧部分 # 通常包括变量名、单因素分析的HR(CI)和P值、多因素分析的HR(CI)和P值 labeltext - list( # 第一列变量名 c(NA, plot_data$variable), # 第一行是表头“Variable”用NA占位 # 第二列单因素分析结果 c(Univariate\nHR (95% CI), paste0(plot_data$hr_uni, (, plot_data$ci_low_uni, -, plot_data$ci_high_uni, ))), # 第三列单因素P值 c(P Value, plot_data$p_uni), # 第四列多因素分析结果 c(Multivariate\nHR (95% CI), paste0(plot_data$hr_multi, (, plot_data$ci_low_multi, -, plot_data$ci_high_multi, ))), # 第五列多因素P值 c(P Value, plot_data$p_multi) ) # 2. 构建效应值估计和置信区间矩阵用于绘图部分 # forestplot要求一个矩阵每行对应一个变量每列对应一个模型这里有两列单因素和多因素 # 矩阵的每一行是一个列表包含mean点估计, lowerCI下限, upperCI上限 mean_values - matrix(c(plot_data$hr_uni, plot_data$hr_multi), ncol 2) lower_values - matrix(c(plot_data$ci_low_uni, plot_data$ci_low_multi), ncol 2) upper_values - matrix(c(plot_data$ci_high_uni, plot_data$ci_high_multi), ncol 2) # 将上述三个矩阵组合成一个列表的列表这是forestplot期望的格式 estimates_list - list() for (i in 1:nrow(plot_data)) { estimates_list[[i]] - list( mean c(mean_values[i, 1], mean_values[i, 2]), lower c(lower_values[i, 1], lower_values[i, 2]), upper c(upper_values[i, 1], upper_values[i, 2]) ) }7. 绘制并美化“双CI”森林图现在万事俱备只差绘图。我们将使用forestplot函数并添加大量参数来控制图形外观。# 绘制基础森林图 forestplot(labeltext labeltext, mean estimates_list, lower lapply(estimates_list, function(x) x$lower), upper lapply(estimates_list, function(x) x$upper), is.summary c(TRUE, rep(FALSE, nrow(plot_data))), # 第一行是汇总行表头 graph.pos 3, # 将森林图线条部分放在第3列即单因素P值和多因素HR列之间 xlog TRUE, # X轴使用对数刻度因为HR是对数尺度上的比值 xticks c(0.5, 1, 2), # 设置X轴刻度1是无效线 clip c(0.5, 3), # 限制置信区间的显示范围超出部分会被截断显示为箭头 col fpColors(box c(royalblue, darkred), # 单因素和多因素的点/方框颜色 lines c(royalblue, darkred), # 置信区间线条颜色 summary black), # 汇总行的颜色 boxsize 0.2, # 点/方框的大小 line.margin 0.1, # 行间距 colgap unit(4, mm), # 列间距 graphwidth unit(0.3, npc), # 森林图区域的宽度 txt_gp fpTxtGp(label gpar(cex0.8), # 文本大小 ticks gpar(cex0.8), xlab gpar(cex0.9)), legend c(Univariate, Multivariate), # 图例 legend_args fpLegend(pos list(x0.85, y0.98), # 图例位置 gp gpar(col#CCCCCC, fill#F9F9F9)), # 图例框样式 hrzl_lines list(2 gpar(lty1, lwd1)), # 在第2行后即表头后画一条横线 title Cox Regression Analysis of Factors Associated with Survival (Lung Cancer Data) )运行这段代码你将得到一张包含单因素和多因素结果的双列森林图。单因素结果用蓝色表示多因素结果用红色表示它们并排显示在每个变量行上。8. 进阶美化与自定义上面的图形已经可用但为了达到发表级别我们可能还需要进一步调整。调整X轴和对数刻度对于HR范围较大的情况需要调整xticks和clip参数。# 如果HR范围很大例如从0.1到10 xticks_custom - c(0.1, 0.2, 0.5, 1, 2, 5, 10) clip_custom - c(0.1, 10) # 在forestplot函数中替换对应的参数即可添加参考线并高亮显著结果我们可以通过后处理使用grid包函数或在构建labeltext时标记显著变量。# 方法在变量名前添加星号(*)来标记P0.05的变量 plot_data - plot_data %% mutate(variable_display ifelse(p_multi 0.05, paste0(variable, *), variable)) # 然后在构建labeltext时使用 variable_display 列分组合并变量如果你的变量有分组如“人口学特征”、“临床指标”可以通过在labeltext中插入空行和分组标题行并设置is.summary参数来实现。# 假设我们有三个组 group_labels - c(Demographics, Clinical Factors, Laboratory) # 在labeltext的变量名列表中相应位置插入组名 # 同时需要扩展 estimates_list在对应位置插入 NA # 并将 is.summary 中对应组名和空行的位置设为 TRUE # 这是一个更高级的操作需要仔细调整数据结构和索引保存高清图片使用png,pdf,tiff等函数将图形保存为文件。png(Cox_ForestPlot_Dual_CI.png, width 3200, height 1800, res 300) # 高分辨率PNG # 重新运行 forestplot 绘图代码 dev.off() # 或者保存为PDF矢量图无限缩放 pdf(Cox_ForestPlot_Dual_CI.pdf, width 12, height 7) # 重新运行 forestplot 绘图代码 dev.off()9. 常见问题与排查思路在实践过程中你可能会遇到以下问题问题现象可能原因排查方式解决方案图形空白或只显示部分estimates_list结构错误mean/lower/upper长度不一致。使用str(estimates_list)检查结构。确保每个子列表的mean,lower,upper都是长度为2的数值向量。仔细检查构建estimates_list的循环逻辑确保从plot_data中提取的数据维度正确。置信区间线条错位或缺失graph.pos参数设置错误或者labeltext的列数与图形位置不匹配。确认graph.pos的值是否在labeltext的列数范围内。例如有5列文本graph.pos可以是2到5。调整graph.pos参数。通常放在中间列如第3列视觉效果较好。X轴刻度标签重叠或显示不全xticks设置过于密集或者xlogTRUE时刻度值选择不当。尝试减少xticks中的刻度数量或使用pretty函数自动生成。手动设置一组更稀疏、更有代表性的刻度值如c(0.5, 1, 2, 4)。调整图形宽度 (graphwidth) 也可能有帮助。图例不显示或位置不对legend参数已设置但图例被画布边界截断。检查legend_args中的pos坐标是否超出了画布范围通常是0到1。调整pos list(x, y)的值例如x0.8, y0.95。确保forestplot函数调用中包含了legend参数。中文变量名或标签显示为乱码图形设备不支持中文字体。在保存为文件或RStudio中预览时检查。在绘图前设置中文字体。例如在保存为PNG时png(..., family“SimHei”)。在RStudio中可能需要修改图形设备的默认字体。P值“0.001”在图中显示为TRUE/FALSE在数据整理时p_val_format逻辑判断产生逻辑值而非字符。检查创建p_val_format的ifelse语句。确保yes和no参数都是字符型。使用as.character()明确转换p_val_format - ifelse(p_val 0.001, “0.001”, as.character(p_val))。10. 最佳实践与工程化建议将绘图代码脚本化、模块化能极大提升你的工作效率和结果的可重复性。1. 封装为函数将整个流程数据清洗、模型拟合、结果提取、绘图封装成一个函数。这样对于新的数据集你只需要更换数据源和变量名即可。draw_dual_ci_forest - function(data, time_var, status_var, covar_list, title Cox Regression Forest Plot) { # 函数体包含本文第4-7步的所有代码 # ... # 最终返回 forestplot 对象或直接绘图 } # 使用函数 draw_dual_ci_forest(lung_data, time, status, c(age, sex, ph.ecog), My Analysis)2. 参数化配置将颜色、图形尺寸、字体大小、刻度等视觉元素提取为函数的参数方便批量生成不同风格的图用于不同期刊。3. 结果输出自动化在函数末尾集成图形保存逻辑并自动生成结果摘要的 CSV 文件。# 在函数内部 write.csv(plot_data, cox_regression_results.csv, row.names FALSE) png(forestplot.png, width10, height6, unitsin, res300) print(forestplot(...)) dev.off()4. 处理分类变量上面的示例针对连续变量。对于分类变量如性别、种族Cox 回归会以某一类为参照生成多个HR。在整理plot_data时需要将因子变量展开并为每个水平创建单独的行。summary(coxph_model)中的coefficients行数会告诉你具体有多少个水平被估计。5. 考虑交互项如果你的多因素模型包含了交互项在提取结果和绘图时需要特别小心。交互项的结果通常不适合与主效应并排放在同一张简单的森林图中可能需要单独绘制或使用其他形式展示。掌握 Cox 回归双 CI 森林图的绘制远不止学会调用一个 R 包。它迫使你更清晰地理解从原始数据、统计模型到结果呈现的完整链条。当你能够游刃有余地生成这样一张信息密度高、视觉专业的图表时你向合作者和读者传递的不仅是数据更是严谨的分析逻辑和高效的沟通能力。本文提供的代码框架是一个强大的起点你可以根据自己的数据结构和审美偏好进行修改和扩展。建议将核心代码保存为脚本模板在未来的分析工作中反复调用和优化这将成为你生物统计或临床研究工具箱中一件不可或缺的利器。