在临床研究和生存分析中Cox比例风险模型是评估预后因素与生存时间关系的核心工具。然而仅仅计算出风险比HR和P值往往不够直观尤其是在需要同时展示多个变量或对比单因素与多因素分析结果时。森林图Forest Plot以其直观、信息量大的特点成为展示Cox回归结果特别是包含置信区间CI的HR值的标准可视化方式。本文将手把手教你如何从数据准备开始到最终绘制出专业的单因素、多因素Cox回归双CI森林图涵盖R语言完整代码、参数详解以及常见问题排查确保你能将这套方法直接应用到自己的研究项目中。1. 背景与核心概念为什么需要森林图在深入代码之前我们有必要厘清几个关键概念这能帮助你理解每一步操作背后的意义而不仅仅是复制代码。Cox比例风险回归模型这是一种半参数模型用于分析一个或多个“协变量”即我们关心的因素如年龄、治疗方案、基因表达量等对个体生存时间的影响。它的核心输出是风险比Hazard Ratio, HR。HR 1表示该因素是危险因素会增加事件如死亡、复发发生的风险HR 1则表示是保护因素HR 1则意味着该因素对风险无影响。置信区间Confidence Interval, CI我们得到的HR是一个点估计值存在抽样误差。95% CI给出了这个HR值可能的范围。如果CI包含了1则在统计学上通常认为该因素无显著意义P 0.05。在森林图中CI通常以一条水平线段表示。森林图它将多个研究或多个变量的统计分析结果通常是效应量及其CI整合在一张图上形似“森林”故得名。在Cox回归中森林图可以直观对比一眼看出哪些因素显著CI横线不跨越1以及效应强弱HR点距离1的远近。展示信息在一张图上集中展示变量名、HR值、95% CI、P值信息密度高。对比分析并排展示单因素和多因素分析结果可以观察在调整了其他混杂因素后某个因素的效应是否依然独立存在。单因素 vs. 多因素分析单因素分析将每个变量单独放入Cox模型。优点是简单能初步筛选变量。缺点是可能受到其他混杂因素的影响结果可能不准确。多因素分析将所有需要考察的变量同时放入一个Cox模型。优点是能评估每个因素的独立效应即在控制其他因素不变的情况下该因素的作用。结果是更可靠的结论依据。绘制“双CI森林图”即在一张图上为每个变量绘制两条CI横线一条代表单因素分析结果一条代表多因素分析结果是进行上述对比的绝佳方式。2. 环境准备与数据模拟我们将使用R语言完成全部操作因为它拥有最成熟和丰富的生存分析及绘图生态。2.1 软件与包准备首先确保你已安装R和RStudio推荐。我们将主要依赖以下R包请先安装它们# 如果尚未安装请运行以下命令 install.packages(c(survival, survminer, dplyr, ggplot2, forestplot))survival核心生存分析包提供coxph函数进行Cox回归。survminer基于ggplot2的生存分析可视化包其中的ggforest函数可以方便地绘制森林图但自定义“双CI”需要一些技巧。forestplot专门用于绘制高度定制化森林图的包非常适合我们制作复杂表格样式的森林图。dplyr数据整理神器。ggplot2绘图基础包。2.2 模拟一份临床数据集为了完整演示我们模拟一份包含生存时间、结局事件和若干协变量的数据集。# 设置随机种子保证结果可复现 set.seed(123) # 生成300个样本 n - 300 data - data.frame( # 生存时间月 time round(runif(n, min1, max120), 1), # 结局状态1事件发生如死亡0删失 status rbinom(n, size1, prob0.3), # 协变量1年龄连续变量 Age round(rnorm(n, mean60, sd10)), # 协变量2性别分类变量 (0女 1男) Gender factor(rbinom(n, size1, prob0.5), labelsc(Female, Male)), # 协变量3肿瘤分级有序分类变量 (I, II, III) Grade factor(sample(1:3, n, replaceTRUE, probc(0.2, 0.5, 0.3)), labelsc(I, II, III)), # 协变量4治疗方案分类变量 (A, B) Treatment factor(rbinom(n, size1, prob0.6), labelsc(A, B)), # 协变量5某个生物标志物表达量连续变量 Biomarker round(runif(n, min0.5, max3.5), 2) ) # 查看数据结构 head(data) str(data)运行后你会看到一个包含time,status,Age,Gender,Grade,Treatment,Biomarker共7列的数据框。这就是我们后续分析的基础。3. 核心分析执行单因素与多因素Cox回归我们的目标是得到每个变量在两个模型单因素、多因素下的HR、CI和P值。3.1 单因素Cox回归分析循环或使用lapply对每个变量单独进行Cox回归。library(survival) library(dplyr) # 定义要分析的变量名排除time和status covariates - c(Age, Gender, Grade, Treatment, Biomarker) # 初始化一个列表来存储单因素模型结果 uni_models - list() # 循环拟合单因素模型 for (cov in covariates) { formula - as.formula(paste(Surv(time, status) ~, cov)) uni_models[[cov]] - coxph(formula, data data) } # 提取单因素分析结果HR, 95% CI, P值 uni_results - lapply(uni_models, function(model) { # 获取模型摘要 sum_model - summary(model) # 提取关键信息 hr - sum_model$coefficients[, exp(coef)] # 风险比HR ci_lower - sum_model$conf.int[, lower .95] ci_upper - sum_model$conf.int[, upper .95] p_value - sum_model$coefficients[, Pr(|z|)] # 对于多水平的分类变量如Grade结果有多行需要处理变量名 var_name - rownames(sum_model$coefficients) # 返回一个数据框 data.frame( Variable var_name, HR_uni round(hr, 3), CI_lower_uni round(ci_lower, 3), CI_upper_uni round(ci_upper, 3), P_uni format.pval(p_value, digits3, eps0.001) # 格式化P值 ) }) # 将列表合并成一个数据框 uni_df - bind_rows(uni_results) print(uni_df)3.2 多因素Cox回归分析将所有变量同时放入一个模型。# 构建多因素模型公式 # 注意对于分类变量R会自动处理默认以第一水平为参照。 multi_formula - as.formula(Surv(time, status) ~ Age Gender Grade Treatment Biomarker) multi_model - coxph(multi_formula, data data) # 提取多因素分析结果 sum_multi - summary(multi_model) multi_df - data.frame( Variable rownames(sum_multi$coefficients), HR_multi round(sum_multi$coefficients[, exp(coef)], 3), CI_lower_multi round(sum_multi$conf.int[, lower .95], 3), CI_upper_multi round(sum_multi$conf.int[, upper .95], 3), P_multi format.pval(sum_multi$coefficients[, Pr(|z|)], digits3, eps0.001) ) print(multi_df)4. 数据整合为绘制双CI森林图做准备绘制森林图需要一份特定格式的数据。我们需要将单因素和多因素的结果按变量对齐合并。# 合并两个结果数据框 # 注意合并键是Variable但需要确保顺序和名称一致。 # 对于分类变量单因素结果中变量名可能是GenderMale多因素中也是可以直接合并。 # 但Grade这类多水平变量在结果中会展开为GradeII, GradeIII。 forest_data - merge(uni_df, multi_df, by Variable, all TRUE) # 为了绘图美观我们常常需要更清晰的变量标签。 # 创建一个用于显示的标签列 forest_data$Label - forest_data$Variable # 可以手动美化一下标签例如 forest_data$Label - gsub(GenderMale, Gender (Male vs Female), forest_data$Label) forest_data$Label - gsub(GradeII, Grade II vs I, forest_data$Label) forest_data$Label - gsub(GradeIII, Grade III vs I, forest_data$Label) forest_data$Label - gsub(TreatmentB, Treatment (B vs A), forest_data$Label) # Age和Biomarker保持不变或也可以加标签 forest_data$Label - ifelse(forest_data$Variable Age, Age (years), forest_data$Label) forest_data$Label - ifelse(forest_data$Variable Biomarker, Biomarker level, forest_data$Label) # 重新排序通常将分类变量的参照水平省略或按逻辑顺序排列。 # 这里我们简单按原始顺序你可以根据需要调整。 forest_data - forest_data[order(factor(forest_data$Variable, levels c(Age, GenderMale, GradeII, GradeIII, TreatmentB, Biomarker))), ] # 查看整理好的数据 print(forest_data)现在forest_data数据框中包含了每个变量对应的单因素HR/CI/P和多因素HR/CI/P。这是绘图的基础。5. 使用forestplot包绘制专业双CI森林图forestplot包提供了极大的灵活性可以绘制出类似医学期刊上发表的高质量森林图。5.1 构建绘图所需的表格文本forestplot函数的核心是接受一个矩阵或数据框作为“表格文本”并在其旁边绘制森林图。library(forestplot) # 1. 构建最左侧的文本标签列表格内容 tabletext - cbind( c(Variable, forest_data$Label), # 第一列变量标签 c(HR (95% CI)\nUnivariate, paste(forest_data$HR_uni, (, forest_data$CI_lower_uni, -, forest_data$CI_upper_uni, ), sep)), # 第二列单因素结果 c(P Value\nUni, forest_data$P_uni), # 第三列单因素P值 c(HR (95% CI)\nMultivariate, paste(forest_data$HR_multi, (, forest_data$CI_lower_multi, -, forest_data$CI_upper_multi, ), sep)), # 第四列多因素结果 c(P Value\nMulti, forest_data$P_multi) # 第五列多因素P值 ) # 2. 构建用于绘图的数值数据HR和CI # forestplot需要均值HR和下限、上限。 # 我们需要为单因素和多因素分别准备一组数据。 mean_uni - c(NA, forest_data$HR_uni) # 第一行是标题用NA填充 lower_uni - c(NA, forest_data$CI_lower_uni) upper_uni - c(NA, forest_data$CI_upper_uni) mean_multi - c(NA, forest_data$HR_multi) lower_multi - c(NA, forest_data$CI_lower_multi) upper_multi - c(NA, forest_data$CI_upper_multi) # 将两组数据组合成一个列表 mean_list - list(mean_uni, mean_multi) lower_list - list(lower_uni, lower_multi) upper_list - list(upper_uni, upper_multi)5.2 绘制基础双CI森林图# 绘制森林图 forestplot(labeltext tabletext, mean mean_list, # 传入列表绘制多条CI线 lower lower_list, upper upper_list, # 图形参数设置 graph.pos 3, # 森林图放在第3列后面即表格的第1、2列是文本第3列后开始画图 hrzl_lines list(2 gpar(lwd1, col#000000)), # 在第2行标题行后画一条黑色横线 txt_gp fpTxtGp(label gpar(cex0.9), # 标签字体大小 ticks gpar(cex0.8), # 刻度字体大小 xlab gpar(cex0.9)), # X轴标签字体大小 col fpColors(box c(blue, red), # 单因素和多因素HR点框的颜色 lines c(blue, red), # 对应CI线的颜色 summary c(darkblue, darkred)), # 此处未用用于汇总线 xlab Hazard Ratio (HR), # X轴标签 zero 1, # 无效线HR1的位置 lwd.zero 1.5, # 无效线的粗细 lwd.xaxis 1, # X轴线的粗细 grid TRUE, # 添加垂直网格线 # 图例 legend c(Univariate Analysis, Multivariate Analysis), legend_args fpLegend(pos list(topright), # 图例位置 gp gpar(col#CCCCCC, fill#F9F9F9)), # 图例框样式 # 箱线图HR点的样式 boxsize 0.2, # HR点的大小 line.margin 0.1, # 行间距 colgap unit(4, mm), # 列间距 graphwidth unit(0.3, npc) # 森林图部分的宽度占可用空间的30% )这段代码会生成一张森林图其中每个变量对应两条水平线段蓝色单因素红色多因素和一个方框代表HR点值。图例说明了颜色对应关系。X轴以1为中心如果CI横线不跨越1则说明该因素在该模型下显著。5.3 高级定制与美化你可能需要调整图形以满足特定出版或报告要求。# 更精细控制的示例 forestplot(labeltext tabletext, mean mean_list, lower lower_list, upper upper_list, graph.pos 3, is.summary c(TRUE, rep(FALSE, nrow(forest_data))), # 第一行是标题可以加粗等此处仅作标识 hrzl_lines list(2 gpar(lwd1.5, colblack), 13 gpar(lwd1.5, colblack)), # 在最后一行数据后再加一条线 txt_gp fpTxtGp(label list(gpar(fontfacebold), # 第一行标签加粗 gpar(cex0.85)), # 数据行标签 ticks gpar(cex0.8), xlab gpar(cex1, fontfacebold)), col fpColors(box c(#1C86EE, #EE3B3B), # 使用更柔和的颜色 lines c(#1C86EE, #EE3B3B)), xlab Hazard Ratio, zero 1, lwd.zero 2, lwd.xaxis 1.5, grid structure(c(0.5, 2, 5), # 在HR0.5, 2, 5处画网格线 gp gpar(lty2, col#CCCCCC)), xticks c(0.1, 0.5, 1, 2, 5, 10), # 设置X轴刻度 clip c(0.1, 10), # 限制绘图区域防止极端CI值导致图形变形 legend c(Univariate, Multivariate), legend_args fpLegend(pos list(x0.85, y0.98), titleAnalysis, r unit(0.1, snpc), gp gpar(colblack, lwd1)), boxsize 0.15, line.height unit(0.7, cm), colgap unit(5, mm), graphwidth unit(0.35, npc) )通过调整xticks,clip,grid,line.height等参数你可以获得更符合心意的图形。6. 使用survminer的ggforest进行快速绘制单模型survminer包的ggforest()函数能快速为单个Cox模型生成美观的森林图但它原生不支持并排绘制两个模型的CI。不过我们可以通过一些技巧来组合。library(survminer) library(ggplot2) library(patchwork) # 用于拼图 # 绘制单因素模型森林图以多因素模型为例 p_multi - ggforest(multi_model, data data, main Multivariate Cox Analysis) print(p_multi) # 如果你想并排展示单因素和多因素需要为每个变量手动创建一个数据框 # 这里展示一个更通用的方法用ggplot2从头构建 # 首先将long格式的数据准备好 plot_df - data.frame( Variable rep(forest_data$Label, 2), Model rep(c(Univariate, Multivariate), each nrow(forest_data)), HR c(forest_data$HR_uni, forest_data$HR_multi), Lower c(forest_data$CI_lower_uni, forest_data$CI_lower_multi), Upper c(forest_data$CI_upper_uni, forest_data$CI_upper_multi), P c(forest_data$P_uni, forest_data$P_multi) ) # 确保Variable的顺序 plot_df$Variable - factor(plot_df$Variable, levels rev(forest_data$Label)) # rev用于反转顺序使绘图从上到下 # 使用ggplot2绘制 ggplot(plot_df, aes(x HR, y Variable, color Model)) geom_vline(xintercept 1, linetype dashed, color grey50) geom_point(position position_dodge(width 0.6), size 3) # 错开点 geom_errorbarh(aes(xmin Lower, xmax Upper), height 0.2, position position_dodge(width 0.6)) scale_x_log10() # 由于HR通常呈对数尺度变化使用log10刻度更直观 labs(x Hazard Ratio (log scale), y , color Analysis Model) theme_minimal() theme(legend.position top, panel.grid.major.y element_blank(), panel.grid.minor.y element_blank())这种方法用ggplot2实现了双CI的绘制颜色区分模型点图并排错开更加灵活。7. 常见问题与排查思路在绘制过程中你可能会遇到以下问题问题现象常见原因解决思路错误:x和labels长度不一致用于绘图的数值向量mean,lower,upper与标签文本labeltext的行数不匹配。检查mean_uni等向量长度是否等于tabletext的行数。通常数值向量第一行是NA对应标题行。森林图中CI横线特别长图形变形某个变量的HR置信区间非常宽例如下限接近0上限几十导致X轴尺度被拉大。使用forestplot的clip参数如clipc(0.1, 10)限制绘图区域极端值会用箭头表示。分类变量的参照水平不见了在回归中分类变量的一个水平通常是第一个被设为参照其HR为1不会出现在结果中。这是正常现象。森林图中通常只显示与参照水平比较的其他水平如Grade II vs I。如果你想显示参照水平需要在结果数据框中手动添加一行HR1, CINA。P值显示为0.001或科学计数法默认的format.pval函数会将极小的P值格式化为0.001。如果你希望显示具体值可以使用formatC(p_value, format e, digits 2)显示科学计数法或调整format.pval的eps参数。图形中文字重叠或溢出变量标签太长或列宽太窄。调整forestplot的colgap列间距、graphwidth图宽度参数。缩短变量标签。在ggplot2中调整theme的plot.margin和axis.text.y。ggforest报错找不到数据ggforest函数需要原始数据data来重构模型矩阵。确保在ggforest(model, datayour_data)中传入了正确的原始数据框。多因素模型中某些变量P值很大如0.9可能存在多重共线性或者该变量在调整其他变量后确实没有独立预测作用。检查变量间的相关性cor()。考虑使用方差膨胀因子VIF检查共线性。这属于模型解释问题需结合专业知识判断。8. 最佳实践与工程建议将统计分析可视化不仅仅是跑通代码更要注意结果的准确性和图形的专业性。数据预处理是关键在运行Cox回归前务必检查数据。缺失值使用summary(data)或is.na()检查。根据情况选择删除、插补或作为单独类别。变量类型确保分类变量是factor类型连续变量是numeric类型。coxph会自动处理因子变量。比例风险假设Cox模型的核心假设。可以使用cox.zph()函数进行检验。如果假设被严重违反需要考虑时依协变量模型或其他模型。模型结果解读HR的解释对于连续变量如AgeHR表示该变量每增加一个单位如1岁风险变化的倍数。对于分类变量如Treatment B vs AHR表示B组相对于A组的风险比。CI的重要性永远同时报告HR和95% CI而不仅仅是P值。CI提供了效应大小的估计精度。单因素与多因素结果的差异如果某个变量单因素显著而多因素不显著说明它可能通过其他变量混杂因素产生影响其本身并非独立预后因素。反之如果多因素显著而单因素不显著则可能是抑制变量或交互作用。图形美化与输出保持一致性在同一篇文章或报告中森林图的样式颜色、字体、刻度应保持一致。高分辨率输出用于发表的图片需要高DPI通常300-600。在RStudio中使用ggsave()函数对于ggplot2图形或png()/pdf()设备。# 保存forestplot图形可能稍复杂推荐保存为PDF矢量图 pdf(Cox_ForestPlot_Dual_CI.pdf, width12, height8) forestplot(...) # 你的绘图命令 dev.off() # 保存ggplot2图形 ggsave(Cox_ForestPlot_ggplot2.png, plotlast_plot(), dpi300, width10, height6)添加必要注释在图形下方或图例中说明参照组、使用的模型、样本量等。代码可复现性在脚本开头使用set.seed()确保模拟数据可复现。使用相对路径或here包管理文件路径。将主要步骤封装成函数特别是如果你需要对多个数据集进行相同分析时。超越基础亚组分析森林图可以绘制不同亚组如不同癌症分期中同一治疗效果的森林图。网络Meta分析森林图比较多种干预措施。交互效应可视化如果存在显著的交互项可以绘制交互效应图而不仅仅是主效应森林图。掌握Cox回归双CI森林图的绘制不仅能提升你数据分析结果的呈现质量更能加深你对多变量生存分析模型的理解。从数据清理、模型拟合、结果提取到图形定制每一步都蕴含着对统计原理和编程技巧的运用。建议你用自己的实际数据替换本文的模拟数据从头到尾实践一遍遇到问题再回头查阅相关章节的解决方案。