企业官网建设流程全解析

如果你在医学、公共卫生或生物统计领域工作，一定会遇到这样的场景：你完成了一项队列研究，通过 Cox 回归模型分析了一组变量对生存结局的影响，得到了单因素和多因素分析的结果。现在，你需要向导师、审稿人或合作者展示这些结果。你面对的是一张密密麻麻的表格，里面塞满了 HR、95% CI 和 P 值。如何让这些冰冷的数字变得直观、有力，一眼就能看出哪些因素是保护性的，哪些是危险因素，并且同时展示单因素和多因素的结果对比？

答案就是Cox 回归森林图。它不仅是论文图表中的“颜值担当”，更是高效传递信息的利器。一张好的森林图，能让读者在几秒钟内抓住研究的核心发现。

然而，从原始数据到一张发表级的森林图，中间隔着不少“坑”：如何用 R 语言一次性提取单因素和多因素分析结果？如何将两个模型的结果优雅地整合到同一张图上？如何自定义坐标轴、字体、颜色，以满足不同期刊的挑剔要求？网上教程虽多，但往往只讲单因素图，或者代码复杂晦涩，难以直接套用。

本文将彻底解决这个问题。我将手把手带你，使用 R 语言中强大的forestplot包，绘制一张同时包含单因素和多因素 Cox 回归结果的森林图。我们不止步于“画出来”，更要深入“为什么这么画”，并分享一套可直接复用的、模块化的代码模板。无论你是 R 语言新手，还是希望优化科研工作流的老手，这篇文章都能让你在生存分析的可视化上，迈出坚实的一步。

1. 为什么你需要掌握“双CI”森林图？

在深入代码之前，我们首先要理解，为什么这种同时展示单因素和多因素结果的森林图如此重要。

单因素分析 vs. 多因素分析：单因素分析是“初筛”，它单独考察每个变量与结局的关系，但结果可能受到其他混杂因素的干扰。多因素分析则是“精炼”，它把所有重要变量放在同一个模型中，评估每个变量的“独立贡献”。审稿人最关心的，往往是多因素分析的结果。

传统做法的局限：很多研究者会分别绘制单因素和多因素森林图，或者只在表格中并列展示。前者浪费空间且不便于对比；后者则不够直观，无法快速识别效应值的变化趋势（例如，某个变量从单因素显著变为多因素不显著，这本身就是一个有趣的故事）。

“双CI”森林图的优势：将单因素和多因素的 HR 及其置信区间并排放在同一行，用不同的符号（如方块和菱形）或颜色区分。这种设计让你可以：

1. 一目了然地对比：直接看到校正混杂因素前后，每个变量效应值（HR）和精确度（CI宽度）的变化。
2. 高效利用空间：在一张图上呈现所有信息，符合顶级期刊对图表简洁高效的要求。
3. 讲述数据故事：例如，一个变量单因素分析时 HR=2.5 (1.3-4.8)，多因素分析后变为 HR=1.2 (0.8-1.8)，这强烈暗示其效应可能被其他变量所解释或混淆，这本身就是分析讨论的一部分。

因此，掌握这项技能，不是简单的“美化图表”，而是提升你科研结果呈现专业度和洞察深度的关键一步。

2. 核心工具与概念梳理

在开始动手前，我们需要明确几个核心概念和将要用到的工具。

Cox 比例风险模型：生存分析中最常用的回归模型，用于分析一个或多个变量（协变量）对生存时间的影响。其结果以风险比（Hazard Ratio, HR）呈现。HR > 1 表示该变量是危险因素（增加事件发生风险），HR < 1 表示保护因素。

森林图：一种用于展示多个独立研究结果或同一研究中多个变量效应估计值（如 HR, OR, RR）及其置信区间的图表。每个变量占一行，用一条水平线段（置信区间）和一个点（点估计值，如 HR）表示。

forestplot包：R 语言中绘制发表级森林图的利器。相比基础的forestplot函数或survminer包，forestplot包提供了无与伦比的灵活性，可以轻松构建复杂的、多列的、高度自定义的森林图，完美契合我们绘制“双CI”图的需求。

本次任务的数据流：

1. 数据准备：一个包含生存时间、生存状态和若干协变量的数据集。
2. 模型拟合：分别对每个变量进行单因素 Cox 回归，再对所有变量进行多因素 Cox 回归。
3. 结果提取：从模型结果中提取变量名、HR、置信区间和 P 值，并整理成规整的数据框。
4. 绘图数据构建：构建一个用于forestplot函数的列表结构，包含文本标签、效应值估计和置信区间。
5. 图形绘制与美化：使用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)) %>% # 假设status=2为死亡事件 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.001

6. 构建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（点估计）, lower（CI下限）, upper（CI上限） 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(cex=0.8), # 文本大小 ticks = gpar(cex=0.8), xlab = gpar(cex=0.9)), legend = c("Univariate", "Multivariate"), # 图例 legend_args = fpLegend(pos = list(x=0.85, y=0.98), # 图例位置 gp = gpar(col="#CCCCCC", fill="#F9F9F9")), # 图例框样式 hrzl_lines = list("2" = gpar(lty=1, lwd=1)), # 在第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时标记显著变量。

# 方法：在变量名前添加星号(*)来标记P<0.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. 常见问题与排查思路

在实践过程中，你可能会遇到以下问题：

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", width=10, height=6, units="in", res=300) print(forestplot(...)) dev.off()

4. 处理分类变量：上面的示例针对连续变量。对于分类变量（如性别、种族），Cox 回归会以某一类为参照生成多个HR。在整理plot_data时，需要将因子变量展开，并为每个水平创建单独的行。summary(coxph_model)中的coefficients行数会告诉你具体有多少个水平被估计。

5. 考虑交互项：如果你的多因素模型包含了交互项，在提取结果和绘图时需要特别小心。交互项的结果通常不适合与主效应并排放在同一张简单的森林图中，可能需要单独绘制或使用其他形式展示。

掌握 Cox 回归双 CI 森林图的绘制，远不止学会调用一个 R 包。它迫使你更清晰地理解从原始数据、统计模型到结果呈现的完整链条。当你能够游刃有余地生成这样一张信息密度高、视觉专业的图表时，你向合作者和读者传递的，不仅是数据，更是严谨的分析逻辑和高效的沟通能力。本文提供的代码框架是一个强大的起点，你可以根据自己的数据结构和审美偏好进行修改和扩展。建议将核心代码保存为脚本模板，在未来的分析工作中反复调用和优化，这将成为你生物统计或临床研究工具箱中一件不可或缺的利器。