这次我们来看一个在医学统计和生存分析中非常实用的技能如何绘制 Cox 回归的单因素和多因素分析森林图并且同时展示双置信区间。对于从事临床研究、流行病学或数据分析的朋友来说这几乎是论文和报告中的“标配”图表。它不仅能直观展示每个变量的风险比还能清晰呈现其统计显著性是评估预后因素的核心工具。这个项目的重点不是教你复杂的统计学理论而是直接解决“怎么画出来”的问题。我们将聚焦于使用 R 语言中的survival和forestplot等核心包从数据整理、模型构建到图形美化一步步实现一个专业级的森林图。整个过程不依赖商业软件完全开源可复现适合需要批量分析或集成到自动化报告流程中的场景。本文会带你完成从零到一的完整流程首先准备一个模拟的生存数据集然后分别进行单因素和多因素 Cox 回归分析最后使用forestplot包绘制包含双 CI例如 95% 和 99% 置信区间的森林图。我们会重点关注代码的实用性、图形的可定制性以及如何避免常见的绘图错误确保你读完就能在自己的数据上跑通。1. 核心能力速览能力项说明核心功能绘制 Cox 比例风险回归模型的单因素及多因素分析森林图支持展示多重置信区间。主要工具R 语言核心依赖包survival,forestplot(或ggplot2结合ggforest)。输入要求包含生存时间 (time)、生存状态 (status) 及多个协变量的数据框。输出成果高清可发表的森林图 (PDF/PNG)图表元素字体、颜色、区间高度可定制。适合场景临床研究论文、生存分析报告、批量生成多组数据的分析图表、教学演示。硬件门槛极低。普通电脑即可运行分析过程主要消耗 CPU 和内存与数据集大小相关。自动化潜力高。可通过 R 脚本批量处理多个数据集或亚组并输出统一格式的图表。2. 适用场景与使用边界这个工具适合谁临床研究人员与医学生需要分析患者生存数据寻找独立的预后因素并在论文中展示规范的结果。流行病学与公共卫生分析师从事队列研究评估多种暴露因素与结局事件的风险关联。数据科学家与统计顾问为客户提供生存分析服务需要生成直观、专业的可视化报告。高校教师与学生用于统计学、医学统计学等课程的教学案例与实操练习。能解决什么问题可视化风险比较将 Cox 回归得到的风险比及其置信区间以图形方式呈现一目了然。区分单/多因素结果清晰对比仅考虑单一变量和校正了混杂因素后变量效应的变化。展示统计不确定性通过置信区间的长短和是否跨越无效线HR1判断结果的精确性与显著性。提升报告效率通过脚本化绘图避免手动在图形软件中拼接便于重复分析和更新。不适合什么场景数据不符合 Cox 比例风险假设时需先进行检验。样本量过小或事件数过少导致模型不稳定或无法拟合。需要绘制非常复杂的、包含多层亚组分析的森林图时可能需要更专业的图形包或手动调整。合规与伦理边界所使用的生存数据必须符合伦理审查要求获得必要的知情同意并经过脱敏处理。在发表或报告结果时需对统计方法进行准确描述避免误导性解读。森林图是统计结果的展示其解释应结合专业背景知识不能仅凭图形下结论。3. 环境准备与前置条件为了顺利运行后续代码你需要准备好以下环境操作系统Windows, macOS 或 Linux 均可。R 环境安装最新版的 R 语言环境。建议版本 4.0.0。R 集成开发环境推荐使用 RStudio它提供了友好的代码编辑、运行和图形查看界面。必要的 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, mean60, sd10)) sex - factor(sample(c(Male, Female), n, replaceTRUE, probc(0.6, 0.4))) stage - factor(sample(c(I, II, III, IV), n, replaceTRUE, probc(0.2, 0.3, 0.3, 0.2)), levelsc(I, II, III, IV)) treatment - factor(sample(c(Drug_A, Drug_B), n, replaceTRUE, probc(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(boxroyalblue, linedarkblue, summaryroyalblue), txt_gp fpTxtGp(label gpar(cex0.8), # 标签字体大小 xlab gpar(cex0.9), # X轴标签字体大小 title gpar(cex1.1)) # 标题字体大小 )运行这段代码你将得到一个基础的单因素森林图。图中每个变量对应一条水平线和一个方块方块的位置代表点估计值HR水平线的长度代表95%置信区间。X轴为对数刻度竖线为 HR1 的无效线。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), ncol3), matrix(c(mean_multi, lower_99_multi, upper_99_multi), ncol3) ) # 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(boxdarkred, lines c(darkred, darkgreen), # 为两条CI线指定不同颜色 summarydarkred), fn.ci_norm c(fpDrawNormalCI, fpDrawCircleCI), # 使用不同的图形绘制CI线和点 lwd.ci c(2, 1), # 设置两条CI线的粗细 ci.vertices FALSE, # 不显示CI线两端的顶点 txt_gp fpTxtGp(label gpar(cex0.8), xlab gpar(cex0.9), title gpar(cex1.1)), legend list(95% CI list(coldarkred, lwd2), 99% CI list(coldarkgreen, lwd1)), legend_args fpLegend(pos list(x0.85, y0.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), ncol3), matrix(c(mean_multiv, lower_multiv, upper_multiv), ncol3) ) # 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(boxc(blue, red), # 单因素蓝色多因素红色 lines c(blue, red), summaryc(blue, red)), fn.ci_norm c(fpDrawNormalCI, fpDrawNormalCI), lwd.ci c(1.5, 1.5), ci.vertices FALSE, txt_gp fpTxtGp(label gpar(cex0.75), xlab gpar(cex0.9), title gpar(cex1.1)), legend list(Univariable list(colblue, lwd1.5), Multivariable list(colred, lwd1.5)), legend_args fpLegend(pos list(x0.85, y0.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(boxsteelblue, linesteelblue, summaryroyalblue, hrz_lines #444444), # 水平分隔线颜色 lwd.ci 2, ci.vertices TRUE, # 显示CI线两端的短横线 ci.vertices.height 0.1, # 短横线高度 txt_gp fpTxtGp(label gpar(fontfamily sans, cex0.9), ticks gpar(cex0.8), xlab gpar(cex0.9, fontfacebold), title gpar(cex1.1, fontfacebold)), hrzl_lines list(2 gpar(lwd1, lty2, colgray)), # 在第2行后加虚线 grid TRUE # 添加垂直网格线 )通过调整graph.pos,colgap,lineheight,col,hrzl_lines,grid等参数你可以控制图表的几乎所有视觉元素使其完全符合目标期刊或报告的格式要求。10. 常见问题与排查方法在绘制森林图的过程中你可能会遇到以下问题问题现象可能原因排查方式解决方案错误object ‘Surv’ not found未加载survival包。检查是否运行了library(survival)。在分析前加载survival包。错误could not find function “forestplot”未安装或加载forestplot包。检查是否运行了install.packages(“forestplot”)和library(forestplot)。安装并加载forestplot包。图形中所有置信区间都特别宽或异常模型可能未收敛或某个组的样本量/事件数太少。查看coxph模型的summary()输出检查系数是否异常大标准误是否巨大。检查数据考虑合并某些分类水平或增加样本量。确保模型满足比例风险假设。因子变量的参照组也出现在图中在整理结果时错误地将参照组HR1也纳入了数据框。检查用于绘图的结果数据框univ_results或multiv_results。在 Cox 回归中因子变量的第一水平是参照组其 HR 为 1。在整理结果时通常只纳入非参照组的水平。确保你的结果数据框中不包含 HR 为 1 且 CI 为 NA 的行。森林图的线条或方块位置错乱提供给forestplot的mean,lower,upper向量长度与labeltext行数不匹配。使用length()函数检查各个向量的长度。确保labeltext的行数等于mean等向量的长度。仔细核对数据准备步骤。通常labeltext的第一行是表头对应的效应值向量第一个元素是NA。确保这种对应关系一致。图形显示不全或超出边界图形设备绘图窗口太小或者变量太多行高不够。尝试调整 RStudio 的绘图窗口大小或使用pdf(),png()等函数设置更大的输出尺寸。使用pdf(“forestplot.pdf”, width10, height12)在绘图前打开一个 PDF 设备绘图完成后用dev.off()关闭。调整width和height参数。想改变颜色、线型不生效fpColors或col参数设置不正确。查看?forestplot和?fpColors帮助文档确认参数格式。col参数应是一个由fpColors()函数创建的对象或一个颜色向量。确保你传递给lines参数的颜色数量与 CI 的组数一致。如何添加亚组分析数据结构和绘图矩阵组织更复杂。规划好标签文本矩阵将亚组标题行对应的is.summary设为TRUE。在准备labeltext和效应值矩阵时插入亚组标题行。将该行对应的效应值设为NA并在is.summary向量中将该行标记为TRUE它就会以加粗或不同格式显示。11. 最佳实践与使用建议数据清洗先行在运行 Cox 回归前务必处理缺失值、检查变量类型数值型、因子型。对于因子变量使用factor()并设置正确的参照水平。模型假设检验使用cox.zph()函数检验比例风险假设。如果假设被严重违反考虑使用参数模型或分层 Cox 模型。结果备份将coxph()模型对象和summary()结果保存为.RData文件便于后续复查或绘制其他图形。save(multiv_model, multiv_summary, file “cox_multiv_model_results.RData”)脚本模块化将数据准备、模型拟合、结果提取和图形绘制写成独立的函数或 R 脚本文件。这样便于对不同数据集进行批量分析。图形输出设置对于发表或报告使用高分辨率输出。png(“my_forestplot.png”, width3200, height2400, res300) # 高分辨率PNG # 在这里运行你的 forestplot() 代码 dev.off()版本控制使用 Git 等工具管理你的分析脚本特别是当分析需要多次迭代或合作完成时。解释与报告在报告中不仅要提供森林图还应附上相应的表格列出具体的 HR、CI 和 P 值。对图中任何重要的发现如 CI 跨越 1或校正前后 HR 方向改变进行文字解释。掌握 Cox 回归森林图的绘制能极大提升你处理生存分析数据的效率和成果展示的专业性。从单因素到多因素从基础图形到包含双 CI 的复杂图表整个流程都可以通过 R 脚本自动化完成。建议你从本文的模拟数据代码开始逐步替换成自己的真实数据并尝试调整各种图形参数直到生成完全符合你需求的图表。 
长河君 · 2026/9/3 6:10:37 · 阅读 1.2k R语言绘制Cox回归森林图：单因素/多因素分析与双置信区间可视化
这次我们来看一个在医学统计和生存分析中非常实用的技能如何绘制 Cox 回归的单因素和多因素分析森林图并且同时展示双置信区间。对于从事临床研究、流行病学或数据分析的朋友来说这几乎是论文和报告中的“标配”图表。它不仅能直观展示每个变量的风险比还能清晰呈现其统计显著性是评估预后因素的核心工具。这个项目的重点不是教你复杂的统计学理论而是直接解决“怎么画出来”的问题。我们将聚焦于使用 R 语言中的survival和forestplot等核心包从数据整理、模型构建到图形美化一步步实现一个专业级的森林图。整个过程不依赖商业软件完全开源可复现适合需要批量分析或集成到自动化报告流程中的场景。本文会带你完成从零到一的完整流程首先准备一个模拟的生存数据集然后分别进行单因素和多因素 Cox 回归分析最后使用forestplot包绘制包含双 CI例如 95% 和 99% 置信区间的森林图。我们会重点关注代码的实用性、图形的可定制性以及如何避免常见的绘图错误确保你读完就能在自己的数据上跑通。1. 核心能力速览能力项说明核心功能绘制 Cox 比例风险回归模型的单因素及多因素分析森林图支持展示多重置信区间。主要工具R 语言核心依赖包survival,forestplot(或ggplot2结合ggforest)。输入要求包含生存时间 (time)、生存状态 (status) 及多个协变量的数据框。输出成果高清可发表的森林图 (PDF/PNG)图表元素字体、颜色、区间高度可定制。适合场景临床研究论文、生存分析报告、批量生成多组数据的分析图表、教学演示。硬件门槛极低。普通电脑即可运行分析过程主要消耗 CPU 和内存与数据集大小相关。自动化潜力高。可通过 R 脚本批量处理多个数据集或亚组并输出统一格式的图表。2. 适用场景与使用边界这个工具适合谁临床研究人员与医学生需要分析患者生存数据寻找独立的预后因素并在论文中展示规范的结果。流行病学与公共卫生分析师从事队列研究评估多种暴露因素与结局事件的风险关联。数据科学家与统计顾问为客户提供生存分析服务需要生成直观、专业的可视化报告。高校教师与学生用于统计学、医学统计学等课程的教学案例与实操练习。能解决什么问题可视化风险比较将 Cox 回归得到的风险比及其置信区间以图形方式呈现一目了然。区分单/多因素结果清晰对比仅考虑单一变量和校正了混杂因素后变量效应的变化。展示统计不确定性通过置信区间的长短和是否跨越无效线HR1判断结果的精确性与显著性。提升报告效率通过脚本化绘图避免手动在图形软件中拼接便于重复分析和更新。不适合什么场景数据不符合 Cox 比例风险假设时需先进行检验。样本量过小或事件数过少导致模型不稳定或无法拟合。需要绘制非常复杂的、包含多层亚组分析的森林图时可能需要更专业的图形包或手动调整。合规与伦理边界所使用的生存数据必须符合伦理审查要求获得必要的知情同意并经过脱敏处理。在发表或报告结果时需对统计方法进行准确描述避免误导性解读。森林图是统计结果的展示其解释应结合专业背景知识不能仅凭图形下结论。3. 环境准备与前置条件为了顺利运行后续代码你需要准备好以下环境操作系统Windows, macOS 或 Linux 均可。R 环境安装最新版的 R 语言环境。建议版本 4.0.0。R 集成开发环境推荐使用 RStudio它提供了友好的代码编辑、运行和图形查看界面。必要的 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, mean60, sd10)) sex - factor(sample(c(Male, Female), n, replaceTRUE, probc(0.6, 0.4))) stage - factor(sample(c(I, II, III, IV), n, replaceTRUE, probc(0.2, 0.3, 0.3, 0.2)), levelsc(I, II, III, IV)) treatment - factor(sample(c(Drug_A, Drug_B), n, replaceTRUE, probc(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(boxroyalblue, linedarkblue, summaryroyalblue), txt_gp fpTxtGp(label gpar(cex0.8), # 标签字体大小 xlab gpar(cex0.9), # X轴标签字体大小 title gpar(cex1.1)) # 标题字体大小 )运行这段代码你将得到一个基础的单因素森林图。图中每个变量对应一条水平线和一个方块方块的位置代表点估计值HR水平线的长度代表95%置信区间。X轴为对数刻度竖线为 HR1 的无效线。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), ncol3), matrix(c(mean_multi, lower_99_multi, upper_99_multi), ncol3) ) # 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(boxdarkred, lines c(darkred, darkgreen), # 为两条CI线指定不同颜色 summarydarkred), fn.ci_norm c(fpDrawNormalCI, fpDrawCircleCI), # 使用不同的图形绘制CI线和点 lwd.ci c(2, 1), # 设置两条CI线的粗细 ci.vertices FALSE, # 不显示CI线两端的顶点 txt_gp fpTxtGp(label gpar(cex0.8), xlab gpar(cex0.9), title gpar(cex1.1)), legend list(95% CI list(coldarkred, lwd2), 99% CI list(coldarkgreen, lwd1)), legend_args fpLegend(pos list(x0.85, y0.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), ncol3), matrix(c(mean_multiv, lower_multiv, upper_multiv), ncol3) ) # 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(boxc(blue, red), # 单因素蓝色多因素红色 lines c(blue, red), summaryc(blue, red)), fn.ci_norm c(fpDrawNormalCI, fpDrawNormalCI), lwd.ci c(1.5, 1.5), ci.vertices FALSE, txt_gp fpTxtGp(label gpar(cex0.75), xlab gpar(cex0.9), title gpar(cex1.1)), legend list(Univariable list(colblue, lwd1.5), Multivariable list(colred, lwd1.5)), legend_args fpLegend(pos list(x0.85, y0.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(boxsteelblue, linesteelblue, summaryroyalblue, hrz_lines #444444), # 水平分隔线颜色 lwd.ci 2, ci.vertices TRUE, # 显示CI线两端的短横线 ci.vertices.height 0.1, # 短横线高度 txt_gp fpTxtGp(label gpar(fontfamily sans, cex0.9), ticks gpar(cex0.8), xlab gpar(cex0.9, fontfacebold), title gpar(cex1.1, fontfacebold)), hrzl_lines list(2 gpar(lwd1, lty2, colgray)), # 在第2行后加虚线 grid TRUE # 添加垂直网格线 )通过调整graph.pos,colgap,lineheight,col,hrzl_lines,grid等参数你可以控制图表的几乎所有视觉元素使其完全符合目标期刊或报告的格式要求。10. 常见问题与排查方法在绘制森林图的过程中你可能会遇到以下问题问题现象可能原因排查方式解决方案错误object ‘Surv’ not found未加载survival包。检查是否运行了library(survival)。在分析前加载survival包。错误could not find function “forestplot”未安装或加载forestplot包。检查是否运行了install.packages(“forestplot”)和library(forestplot)。安装并加载forestplot包。图形中所有置信区间都特别宽或异常模型可能未收敛或某个组的样本量/事件数太少。查看coxph模型的summary()输出检查系数是否异常大标准误是否巨大。检查数据考虑合并某些分类水平或增加样本量。确保模型满足比例风险假设。因子变量的参照组也出现在图中在整理结果时错误地将参照组HR1也纳入了数据框。检查用于绘图的结果数据框univ_results或multiv_results。在 Cox 回归中因子变量的第一水平是参照组其 HR 为 1。在整理结果时通常只纳入非参照组的水平。确保你的结果数据框中不包含 HR 为 1 且 CI 为 NA 的行。森林图的线条或方块位置错乱提供给forestplot的mean,lower,upper向量长度与labeltext行数不匹配。使用length()函数检查各个向量的长度。确保labeltext的行数等于mean等向量的长度。仔细核对数据准备步骤。通常labeltext的第一行是表头对应的效应值向量第一个元素是NA。确保这种对应关系一致。图形显示不全或超出边界图形设备绘图窗口太小或者变量太多行高不够。尝试调整 RStudio 的绘图窗口大小或使用pdf(),png()等函数设置更大的输出尺寸。使用pdf(“forestplot.pdf”, width10, height12)在绘图前打开一个 PDF 设备绘图完成后用dev.off()关闭。调整width和height参数。想改变颜色、线型不生效fpColors或col参数设置不正确。查看?forestplot和?fpColors帮助文档确认参数格式。col参数应是一个由fpColors()函数创建的对象或一个颜色向量。确保你传递给lines参数的颜色数量与 CI 的组数一致。如何添加亚组分析数据结构和绘图矩阵组织更复杂。规划好标签文本矩阵将亚组标题行对应的is.summary设为TRUE。在准备labeltext和效应值矩阵时插入亚组标题行。将该行对应的效应值设为NA并在is.summary向量中将该行标记为TRUE它就会以加粗或不同格式显示。11. 最佳实践与使用建议数据清洗先行在运行 Cox 回归前务必处理缺失值、检查变量类型数值型、因子型。对于因子变量使用factor()并设置正确的参照水平。模型假设检验使用cox.zph()函数检验比例风险假设。如果假设被严重违反考虑使用参数模型或分层 Cox 模型。结果备份将coxph()模型对象和summary()结果保存为.RData文件便于后续复查或绘制其他图形。save(multiv_model, multiv_summary, file “cox_multiv_model_results.RData”)脚本模块化将数据准备、模型拟合、结果提取和图形绘制写成独立的函数或 R 脚本文件。这样便于对不同数据集进行批量分析。图形输出设置对于发表或报告使用高分辨率输出。png(“my_forestplot.png”, width3200, height2400, res300) # 高分辨率PNG # 在这里运行你的 forestplot() 代码 dev.off()版本控制使用 Git 等工具管理你的分析脚本特别是当分析需要多次迭代或合作完成时。解释与报告在报告中不仅要提供森林图还应附上相应的表格列出具体的 HR、CI 和 P 值。对图中任何重要的发现如 CI 跨越 1或校正前后 HR 方向改变进行文字解释。掌握 Cox 回归森林图的绘制能极大提升你处理生存分析数据的效率和成果展示的专业性。从单因素到多因素从基础图形到包含双 CI 的复杂图表整个流程都可以通过 R 脚本自动化完成。建议你从本文的模拟数据代码开始逐步替换成自己的真实数据并尝试调整各种图形参数直到生成完全符合你需求的图表。