你肯定见过那种密密麻麻、信息量巨大的森林图Forest Plot尤其是在医学、流行病学或生存分析的论文里。它把一堆风险比Hazard Ratio, HR和置信区间Confidence Interval, CI直观地铺开一眼就能看出哪些因素是保护性的哪些是危险的。但当你自己动手特别是用R语言处理Cox比例风险模型时从跑出coxph()结果到画出一张既专业又清晰、同时包含单因素和多因素分析结果的森林图中间往往隔着好几个让人头疼的“坑”。比如单因素分析结果是一堆变量多因素分析结果是筛选后的另一堆变量怎么把它们优雅地整合到一张图里置信区间的横线怎么调整才不显得拥挤变量名的标签太长怎么办那些代表统计显著性的星号*和P值是该放在图里还是图外的表格里更实际的是当你把代码从自己电脑搬到合作者的环境或者准备把图嵌入报告、论文时字体、尺寸、导出格式又是一堆问题。这远不是调用一个plot()函数那么简单它涉及数据整理、图形美学和科研表达的交叉。今天我们不谈复杂的统计学原理假设你已经了解Cox回归的基本思想而是聚焦于一个非常实际的工程问题如何将Cox单因素和多因素分析的结果系统性地、可复现地整合成一张用于发表级报告的森林图。我会基于常见的survival包和forestplot包的工作流拆解从数据准备、图形绘制到细节打磨的全过程。你会发现真正的难点不在于画图本身而在于对分析流程的掌控和对结果呈现的深度定制。1. 理解任务核心不是“画图”而是“结果整合与呈现”在开始写任何代码之前我们必须明确一点绘制森林图本身只是一个可视化步骤它之前有更重要的数据准备步骤。对于Cox回归的森林图尤其是结合单因素和多因素分析时我们通常要完成以下数据流单因素分析 (Univariate Analysis)将每个感兴趣的变量单独放入Cox模型得到一组HR和CI。目的是初步筛选可能与生存结局相关的变量。多因素分析 (Multivariate Analysis)将单因素分析中显著或基于学科知识的变量同时放入一个Cox模型得到调整了其他因素后的HR和CI。目的是评估每个变量的独立效应。结果整合将两步分析的结果变量名、HR值、CI上下限、P值整理成一个结构化的数据框Data Frame。可视化利用这个数据框绘制森林图并清晰地标注出哪些结果来自单因素分析哪些来自多因素分析。很多教程只讲第4步但前3步的数据处理才是决定森林图是否准确、清晰的关键。你的主判断应该是一张好的Cox森林图其质量70%取决于前期数据整理的严谨与清晰30%才是图形参数的调整。2. 构建分析流程从数据到结果数据框让我们用一个模拟的生存数据集来演示。假设我们有一个名为my_surv_data的数据框包含生存时间time、生存状态status以及若干个协变量如年龄age连续、性别sex分类、肿瘤分级grade分类、治疗方案treatment分类等。2.1 单因素Cox回归循环我们不想对每个变量手动运行coxph那样效率低且容易出错。更通用的方法是编写一个循环或使用lapply。library(survival) library(broom) # 用于整洁地提取模型结果 # 假设这是你的数据 # my_surv_data - read.csv(your_data.csv) # 定义要分析的变量名列表 univar_vars - c(age, sex, grade, treatment, another_var) # 初始化一个列表来存储每个单因素模型的结果 univar_results - list() for (var in univar_vars) { # 构建公式Surv(time, status) ~ variable formula - as.formula(paste(Surv(time, status) ~, var)) # 拟合Cox模型 cox_model - coxph(formula, data my_surv_data) # 使用broom::tidy提取关键结果并加上变量名 result_df - broom::tidy(cox_model, exponentiate TRUE, conf.int TRUE) result_df$variable - var result_df$analysis - Univariate # 存储 univar_results[[var]] - result_df } # 将列表合并成一个数据框 univar_df - do.call(rbind, univar_results) # 查看整理后的单因素结果 head(univar_df)broom::tidy()函数非常有用它把模型输出变成一个整洁的数据框默认包含term模型项、estimateHR因为exponentiateTRUE、std.error、statistic、p.value、conf.low、conf.high等列。2.2 多因素Cox回归多因素分析需要你事先决定放入模型的变量。这里假设我们基于单因素结果或临床意义选择age,grade,treatment进入多因素模型。# 拟合多因素Cox模型 multivar_formula - Surv(time, status) ~ age grade treatment multivar_model - coxph(multivar_formula, data my_surv_data) # 提取多因素结果 multivar_df - broom::tidy(multivar_model, exponentiate TRUE, conf.int TRUE) multivar_df$analysis - Multivariate # 注意多因素模型的term列已经是变量名我们不需要再添加variable列但为了与单因素数据框结构一致可以重命名或新增 # 这里我们简单处理将term复制到variable列 multivar_df$variable - multivar_df$term # 查看多因素结果 multivar_df2.3 整合单因素与多因素结果这是关键步骤。我们需要一个最终的数据框能够清晰地对应每个变量在单因素和多因素分析中的结果。通常我们会把单因素和多因素的结果行并排放在一起或者用子图区分。# 选择我们需要展示的列 cols_to_keep - c(variable, analysis, estimate, conf.low, conf.high, p.value) univar_for_plot - univar_df[univar_df$variable %in% c(age, grade, treatment), cols_to_keep] multivar_for_plot - multivar_df[, cols_to_keep] # 合并 combined_df - rbind(univar_for_plot, multivar_for_plot) # 为了绘图时顺序正确我们可以设定因子水平 combined_df$variable - factor(combined_df$variable, levels c(age, grade, treatment)) combined_df$analysis - factor(combined_df$analysis, levels c(Univariate, Multivariate)) # 按变量和分析类型排序 combined_df - combined_df[order(combined_df$variable, combined_df$analysis), ] # 生成森林图所需的标签文本 # 通常包括变量名、HR(95% CI)、P值 combined_df$label - paste0( combined_df$variable, (, combined_df$analysis, )\n, sprintf(%.2f, combined_df$estimate), (, sprintf(%.2f, combined_df$conf.low), -, sprintf(%.2f, combined_df$conf.high), )\n, P, sprintf(%.3f, combined_df$p.value) ) print(combined_df)现在combined_df数据框包含了绘制森林图所需的所有核心数据每个估计值HR及其置信区间以及我们自定义的标签。3. 使用forestplot包进行高级绘图虽然R基础绘图或survminer包的ggforest()也能画森林图但forestplot包在定制化方面更灵活尤其适合处理像我们这样整合了多种分析的数据结构。3.1 基础森林图绘制首先安装并加载包install.packages(forestplot)。library(forestplot) library(dplyr) # 为forestplot准备数据矩阵 # forestplot需要几个部分 # 1. 文本标签tabletext # 2. 均值估计mean # 3. 置信区间下限lower # 4. 置信区间上限upper # 提取数据 mean - combined_df$estimate lower - combined_df$conf.low upper - combined_df$conf.high # 创建文本标签矩阵。通常第一列是变量/分析标签后面是HR(95%CI)和P值。 tabletext - cbind( c(Variable (Analysis), combined_df$label), # 第一列标签 c(HR (95% CI), paste0(sprintf(%.2f, combined_df$estimate), (, sprintf(%.2f, combined_df$conf.low), -, sprintf(%.2f, combined_df$conf.high), ))), c(P Value, sprintf(%.3f, combined_df$p.value)) ) # 绘制基础森林图 forestplot(labeltext tabletext, mean c(NA, mean), # 第一行是标题所以用NA lower c(NA, lower), upper c(NA, upper), is.summary c(TRUE, rep(FALSE, nrow(combined_df))), # 第一行是汇总行标题 xlog TRUE, # Cox模型的HR通常取对数刻度显示 boxsize 0.2, col fpColors(box royalblue, line darkblue), xticks c(0.5, 1, 2, 4), # 根据你的HR范围设置 graph.pos 2) # 森林图放在第几列文本后面这段代码会生成一张包含所有结果的森林图。xlog TRUE非常重要因为HR的尺度是对称的HR1表示无效应1表示保护因素1表示危险因素对数刻度能让图形更直观。3.2 区分单因素与多因素分析上面的图把所有结果混在一起了。为了更清晰我们通常希望用视觉元素区分单因素和多因素结果。# 方法使用不同的颜色或形状 # 首先为不同分析类型定义颜色 analysis_colors - c(Univariate #E69F00, Multivariate #0072B2) # 橙色和蓝色 # 创建颜色向量对应每一行数据不包括标题行 box_colors - analysis_colors[combined_df$analysis] # 扩展is.summary除了标题行我们还可以将变量名所在行设为“汇总行”以加粗显示 # 这里我们创建一个逻辑向量标记每个变量第一次出现的行为TRUE作为分组标题 var_first_occurrence - !duplicated(combined_df$variable) is_summary_vec - c(TRUE, var_first_occurrence) # 加上第一行标题 # 绘制 forestplot(labeltext tabletext, mean c(NA, mean), lower c(NA, lower), upper c(NA, upper), is.summary is_summary_vec, xlog TRUE, boxsize 0.2, col fpColors(box box_colors, line box_colors), # 按分析类型着色 fn.ci_norm fpDrawCircleCI, # 用圆圈代替方块可能更清晰 vertices TRUE, xticks c(0.25, 0.5, 1, 2, 4), graph.pos 2, txt_gp fpTxtGp(label gpar(cex0.8), # 调整标签字体大小 ticks gpar(cex0.7), xlab gpar(cex0.9)), hrzl_lines list(2 gpar(lty2)), # 在第二行第一个标题行后加虚线 mar unit(c(4,1,4,1), mm)) # 调整图形边距通过col参数和自定义的box_colors向量单因素和多因素的结果点现在用不同颜色显示。is.summary参数将每个变量的第一行加粗起到了视觉分组的作用。4. 进阶定制与避坑指南一张能直接用于论文或报告的森林图还需要考虑许多细节。4.1 处理分类变量和参照组在Cox模型中分类变量如grade II,grade III会以参照组如grade I为基础生成多个哑变量。在整合结果时你需要明确展示出每个级别与参照组的比较。在数据整理阶段broom::tidy()提取的结果中term列会显示为gradeII,gradeIII。你需要更清晰的标签例如“Grade II vs I”。在标签文本tabletext中第一列应该清晰地标明比较对象。你可能需要手动构建一个更易读的标签向量而不是直接使用变量名。参照线务必确保森林图的垂直参照线在xlogTRUE时位于x1的位置。4.2 控制图形尺寸与导出在RStudio的预览窗口里看起来不错的图导出为PDF或TIFF用于投稿时可能会走样。# 保存为高分辨率PDF pdf(Cox_ForestPlot.pdf, width 10, height 6) # 宽度通常需要大一些以容纳文本 forestplot(...你的绘图参数...) # 重新运行绘图命令 dev.off() # 保存为高分辨率TIFF适合投稿 tiff(Cox_ForestPlot.tiff, width 10, height 6, units in, res 300) forestplot(...你的绘图参数...) dev.off()注意在脚本中pdf()和dev.off()之间的所有绘图命令都会输出到文件。确保你的图形在屏幕显示时就已经布局合理否则保存后问题会更明显。如果标签被截断尝试增加width参数或减小txt_gp中的字体大小cex。4.3 常见错误排查置信区间异常宽或窄检查数据是否存在共线性、模型是否收敛、或生存数据是否存在极端值。使用coxph()后运行summary(model)查看输出确认没有警告信息如“Loglik converged before variable X”可能预示问题。图形中HR或CI值显示为NA检查combined_df中的estimate,conf.low,conf.high列是否存在NA或Inf值。这通常是由于某个变量在某个亚组中事件数为0导致的完全分离。需要考虑合并类别或使用其他统计方法。标签错位或重叠forestplot的labeltext参数接受矩阵。确保你构建的矩阵行数与mean、lower、upper参数的长度一致考虑标题行。使用cex参数调整字体大小或考虑将过长的变量名缩写。P值格式对于非常小的P值如0.001在表格中通常表示为“P0.001”而不是具体的科学计数法值。你可以在构建tabletext时用ifelse语句处理ifelse(p.value 0.001, 0.001, sprintf(%.3f, p.value))。4.4 创建可复现的分析脚本最好的实践是将整个流程封装在一个R脚本或R Markdown文档中。结构如下# 1. 加载库与数据 library(survival) library(broom) library(forestplot) library(dplyr) my_data - read.csv(data.csv) # 2. 定义分析变量 univar_vars - c(age, sex, grade, treatment) multivar_formula - Surv(time, status) ~ age grade treatment # 3. 单因素分析函数 run_univariate - function(vars, data) { ... } # 4. 多因素分析函数 run_multivariate - function(formula, data) { ... } # 5. 结果整合函数 combine_results - function(univar_df, multivar_df) { ... } # 6. 绘图函数 draw_forest_plot - function(combined_df) { ... } # 7. 执行主流程 univar_res - run_univariate(univar_vars, my_data) multivar_res - run_multivariate(multivar_formula, my_data) final_df - combine_results(univar_res, multivar_res) draw_forest_plot(final_df) # 8. 导出图形 ggsave(final_forestplot.png, width10, height7, dpi300) # 如果使用ggplot2系 # 或使用pdf()/tiff()这样的脚本确保了从原始数据到最终图形的全过程可复现也便于你后续更新数据或调整变量。绘制Cox回归的森林图尤其是整合单因素与多因素结果是一个典型的“数据分析管道”任务。它考验的不仅仅是你对某个绘图包函数的熟悉程度更是你对整个统计分析流程的理解和数据操作能力。核心在于构建一个清晰、准确、包含所有必要信息的中间结果数据框。一旦这个数据框准备妥当无论你是用forestplot、ggplot2还是其他工具绘图都变成了相对简单的参数调整问题。因此下次当你需要绘制这样的图时请把至少一半的时间和精力分配给数据整理和验证。先用View()或print()仔细检查combined_df里的每一个HR、CI和P值是否合理确认分类变量的处理方式然后再进入绘图阶段。这张图最终会成为你研究结论的视觉基石值得你投入时间把它打磨精确。