1. 从临床问题到统计可视化：为什么亚组分析与森林图是黄金搭档

在临床研究、流行病学调查乃至任何涉及异质性人群的分析中，一个核心问题常常困扰着我们：某个干预措施或暴露因素的效果，在不同特征的人群中是否一致？比如，一种新药对年轻人和老年人的疗效一样好吗？某个风险因素对男性和女性的影响有差异吗？这就是亚组分析要回答的问题。而森林图，则是呈现这种分析结果最直观、最有力的工具，没有之一。它把每个亚组的效应估计值（比如风险比HR、比值比OR）及其置信区间，以点线图的形式并排展示，一目了然。在R语言生态里，实现从数据清洗、统计建模到精美森林图绘制的完整流程，已经成为研究者的一项必备技能。今天，我们就来彻底拆解这个过程，不仅告诉你每一步怎么操作，更会深入解释背后的统计逻辑和绘图美学，让你做出的图能直接用于顶级期刊的投稿。

2. 分析基石：数据准备与COX比例风险模型拟合

在进行任何花哨的可视化之前，扎实的统计建模是基础。对于生存数据（比如从诊断到死亡或复发的时间），COX比例风险回归模型是亚组分析最常用的工具。我们假设你手头有一个名为surv_data的数据框，至少包含以下变量：time（生存时间）、status（生存状态，1=事件发生，0=删失）、treatment（处理组，如“Drug” vs “Placebo”）、以及一系列用于定义亚组的协变量，如age_group（“<65”, “>=65”）、gender（“Male”, “Female”）、stage（“I”, “II”, “III”）等。

2.1 构建分层回归模型

核心思路是：我们不是为每个亚组单独跑一个模型，而是构建一个包含交互项的全局模型。这样做可以利用全部数据来估计基线风险，同时检验交互项的显著性，在统计上更为严谨。

# 加载生存分析包 library(survival) library(survminer) # 用于模型诊断和基础绘图 # 拟合包含交互项的COX模型 # 假设我们想研究治疗（treatment）在不同年龄组（age_group）和性别（gender）中的效果 cox_model <- coxph(Surv(time, status) ~ treatment * age_group + treatment * gender + stage, data = surv_data) summary(cox_model)

查看summary(cox_model)的输出，你需要重点关注treatmentDrug:age_group>=65和treatmentDrug:genderMale这样的交互项。如果交互项的p值小于0.05（或你设定的显著性水平），则提示治疗效应在该亚组变量上存在异质性，即亚组分析是有意义的。这里一个常见的误解是，只要主效应显著，就可以做亚组分析。实际上，交互项不显著意味着没有统计学证据表明效应大小在不同亚组间有差异，此时强行比较各亚组的HR并做出“对A组有效、对B组无效”的结论是危险的，容易导致假阳性发现。

2.2 提取各亚组的效应估计值

当交互项显著，或基于强有力的生物学假设必须展示亚组结果时，我们需要计算每个亚组内治疗 vs 对照的风险比（HR）。一种清晰的方法是使用emmeans或marginaleffects包进行边际效应估计，但更直接且符合传统文献报告习惯的是拟合分层模型。我们可以用循环或purrr包来优雅地实现。

library(dplyr) library(purrr) # 定义亚组变量列表 subgroup_vars <- c("age_group", "gender") # 使用map_dfr进行循环拟合，并提取结果 subgroup_results <- map_dfr(subgroup_vars, function(var) { # 按亚组变量拆分数据 surv_data %>% group_by(!!sym(var)) %>% # !!sym()用于将字符串转换为变量名 group_modify(~ { # 在每个亚组内拟合一个只包含treatment的简单COX模型 # 注意：这里省略了其他调整变量，实际分析中可能需要调整 fit <- coxph(Surv(time, status) ~ treatment, data = .x) # 提取HR和置信区间 hr_summary <- summary(fit) data.frame( Subgroup = unique(.x[[var]]), Variable = var, HR = hr_summary$conf.int["treatmentDrug", "exp(coef)"], CI_low = hr_summary$conf.int["treatmentDrug", "lower .95"], CI_high = hr_summary$conf.int["treatmentDrug", "upper .95"], P_value = hr_summary$coefficients["treatmentDrug", "Pr(>|z|)"] ) }) })

这样我们就得到了一个数据框subgroup_results，它包含了每个亚组的名称、所属变量、HR值、95%置信区间上下限和P值。这是绘制森林图的直接输入数据。实操心得：在提取结果时，务必确认你提取的系数名称（如treatmentDrug）与模型中的因子水平完全一致。R默认以因子第一个水平为参照，如果你不确定，最好先用levels(your_data$treatment)检查一下。

3. 森林图绘制进阶：从forestplot到forestploter的跃迁

有了结果数据，就可以绘图了。基础的survminer::ggforest()或forestplot包能快速出图，但定制化程度有限，尤其在处理多层级分组、复杂注释时显得力不从心。而forestploter包，正如其名，是一个“绘图器”，它把森林图拆解成一个个可以灵活编辑的表格单元格，实现了数据与视觉元素的完美分离，功能强大到令人惊叹。

3.1 准备forestploter所需的表格数据

forestploter的核心思想是先创建一个“空白”表格（一个数据框），然后在指定位置填入文本、点估计值和置信区间。表格的每一行对应森林图的一行。

library(forestploter) # 假设我们已将subgroup_results整理为如下格式的`df_plot`： # Subgroup | Variable | HR | CI_low | CI_high | P_value # <65 | age_group| 0.65| 0.45 | 0.94 | 0.023 # >=65 | age_group| 1.10| 0.78 | 1.55 | 0.580 # Male | gender | 0.80| 0.60 | 1.07 | 0.134 # Female | gender | 0.50| 0.32 | 0.78 | 0.002 # 1. 创建基础表格 # 我们需要合并“亚组变量”和“亚组水平”两列，并创建用于绘图的占位列 df_plot$` ` <- paste(rep(" ", 20), collapse = "") # 创建一个空列，用于放置森林图 df_plot$`HR (95% CI)` <- sprintf("%.2f (%.2f to %.2f)", df_plot$HR, df_plot$CI_low, df_plot$CI_high) df_plot$`P Value` <- ifelse(df_plot$P_value < 0.001, "<0.001", sprintf("%.3f", df_plot$P_value)) # 2. 确定表格列的顺序 dt <- df_plot[, c("Variable", "Subgroup", " ", "HR (95% CI)", "P Value")] # 3. 将占位列转换为forestploter可识别的格式 # 这一列将存储每个亚组的效应量(HR)和置信区间信息 dt$` ` <- sprintf("%.2f (%.2f to %.2f)", dt$HR, dt$CI_low, dt$CI_high) # 但注意，上面我们为了显示创建了文本列，绘图需要的是数值列。 # 更标准的做法是保留数值列，在绘图函数中指定。 # 让我们重构一下： dt_for_plot <- df_plot[, c("Variable", "Subgroup", "HR", "CI_low", "CI_high", "P_value")] # 添加一个空的森林图列 dt_for_plot$` ` <- NA # 添加显示用的文本列 dt_for_plot$`HR (95% CI)` <- sprintf("%.2f (%.2f to %.2f)", dt_for_plot$HR, dt_for_plot$CI_low, dt_for_plot$CI_high) dt_for_plot$`P Value` <- ifelse(dt_for_plot$P_value < 0.001, "<0.001", sprintf("%.3f", dt_for_plot$P_value))

3.2 使用forestploter绘制与高度定制

现在，我们可以使用forestploter::forest()函数进行绘制。其强大之处在于est、lower、upper参数直接指定点估计和区间，而is_summary参数可以轻松标记汇总行（如交互作用P值行）。

# 定义需要绘图的数值列 estimate <- dt_for_plot$HR low <- dt_for_plot$CI_low high <- dt_for_plot$CI_high # 插入一行空白行，用于分隔不同的亚组变量，并添加“交互作用P值”行 # 首先，我们需要在数据中标记哪些是分组标题行或汇总行 dt_for_plot$is_summary <- c(FALSE, FALSE, TRUE, FALSE, FALSE) # 假设我们的数据顺序是：年龄<65，年龄>=65，（空白/交互P值行），男性，女性 # 第三行我们打算放“Age Group Interaction p = 0.032” # 在第三行插入汇总信息 dt_for_plot$Variable[3] <- "Interaction P value" dt_for_plot$Subgroup[3] <- "0.032" dt_for_plot$`HR (95% CI)`[3] <- "" dt_for_plot$`P Value`[3] <- "" estimate[3] <- NA low[3] <- NA high[3] <- NA # 绘制森林图 p <- forest(dt_for_plot[, c("Variable", "Subgroup", " ", "HR (95% CI)", "P Value")], est = estimate, lower = low, upper = high, ci_column = 3, # 森林图绘制在第3列（即我们预留的空列） ref_line = 1, # 在HR=1处画垂直线 xlim = c(0, 2), # X轴范围 ticks_at = c(0.5, 1, 1.5, 2), # X轴刻度 title = "亚组分析：治疗X在不同人群中的风险比", theme = tm) # tm是一个通过tm_*()函数定义的主题 # 使用tm_*函数进行精细美化 tm <- forest_theme(base_size = 10, core = list(bg_params=list(fill = c("white", "gray95"))), # 隔行换色 summary_fill = "lightblue", # 汇总行背景色 summary_col = "black") # 汇总行文字色 print(p)

关键技巧与避坑指南：

1. 列对齐问题：forestploter对表格的列宽非常敏感。如果出现文字溢出或错位，请检查数据框各列的内容，确保没有异常长的字符串。可以使用colwidths参数手动调整每列宽度。
2. 缺失值处理：对于汇总行或空白行，其est、lower、upper必须设置为NA，否则绘图会出错。
3. 多层级分组：对于嵌套分组（如先按“人口学特征”，下面再分“年龄”、“性别”），可以通过在数据框中插入带有缩进空格的行标题来实现，例如Variable列填写" Age Group"。更高级的做法是构建一个包含分组信息的data.frame，并利用group_by参数。
4. 保存高清图：使用ggsave()保存forestploter输出的对象时，需要指定宽度和高度，特别是当行数很多时，要增加height参数，否则文字会重叠。例如：ggsave("forestplot.png", plot = p, width = 12, height = 8, dpi = 300)。

4. 结果解读与报告：超越图形本身

绘制出漂亮的森林图只是第一步，正确解读和报告结果更为关键。森林图上每一个“方块”和“横线”都讲述着一个故事。

4.1 如何解读森林图中的信息

1. 点估计值（方块）：代表该亚组的HR。方块通常大小与样本量或权重成比例（forestploter可通过size参数设置）。方块在垂直参考线（HR=1）左侧，表示治疗降低风险（HR<1）；在右侧则表示增加风险（HR>1）。
2. 置信区间（横线）：横线越长，表示不确定性越大（标准误大）。这是解读的重中之重：
   • 如果某个亚组的置信区间横线完全跨过了HR=1的垂直线，那么即使点估计值看起来很有希望（比如HR=0.7），我们也不能认为在该亚组中效应是统计学显著的。因为区间包含了1（无效值）。
   • 相反，如果横线没有跨过HR=1，则说明在该亚组内，效应是显著的。
   • 比较不同亚组时，不能只看点估计值的位置。如果两个亚组的置信区间有大量重叠，那么即使它们的点估计值一左一右，我们也不能武断地说效应存在差异。正式的检验应依赖于模型中的交互项P值，它通常会被标注在森林图下方或对应分组行。

4.2 在论文中报告亚组分析结果的规范

仅仅贴一张图是不够的，在论文的方法和结果部分需要清晰描述：

• 方法部分：说明亚组分析是预先设定的还是探索性的。报告用于定义亚组的变量，以及检验交互作用的统计方法（如COX模型中的似然比检验或Wald检验）。
• 结果部分：
  • 首先报告整体人群的效应估计值（主效应）。
  • 然后报告交互作用的检验P值。例如：“治疗与年龄分组之间的交互作用具有统计学意义（交互作用P=0.032）。”
  • 再展示森林图，并附上表格，列出每个亚组的样本量、事件数、HR及其95% CI。
  • 文字描述应聚焦于有显著交互作用的亚组，避免对每一个无显著差异的亚组进行过度解读。例如：“亚组分析显示，治疗在年龄<65岁的患者中显著降低死亡风险（HR 0.65， 95% CI 0.45-0.94），而在年龄≥65岁的患者中未观察到显著获益（HR 1.10， 95% CI 0.78-1.55；交互作用P=0.032）。”
• 讨论部分：对发现的任何异质性进行生物学或临床上的合理解释，同时必须指出探索性亚组分析的局限性，强调其结论需要未来研究验证。

5. 高级应用与扩展：当数据变得更复杂

现实世界的数据分析需求往往更复杂。以下是一些进阶场景及其在R中的处理思路。

5.1 连续变量的亚组分析

当亚组变量是连续的（如年龄、血压），将其武断地二分法会损失信息并可能引入偏倚。更好的方法是：

1. 使用交互项：在COX模型中直接纳入连续变量与治疗变量的乘积项。例如：coxph(Surv(time, status) ~ treatment * age + gender + stage, data)。如果交互项显著，说明治疗效果随年龄变化。
2. 可视化：使用ggplot2绘制限制性立方样条图，展示治疗效应（如HR）随连续变量变化的平滑曲线及其置信带。这比简单的森林图更能揭示趋势。
3. 分层展示：如果仍需要类似森林图的表格，可以按临床常用的切点（如每10岁一个区间）将连续变量离散化，然后按上述流程进行。但务必在报告中说明切点的选择依据。

5.2 使用metafor包进行Meta分析的亚组分析与森林图

如果你的数据来自多项研究（Meta分析），那么亚组分析是探索异质性的核心工具。metafor包是这方面的权威。

library(metafor) # 假设有数据框`meta_data`，包含：study, subgroup, yi(效应量如logHR), vi(效应量方差) # 首先拟合随机效应模型 res <- rma(yi, vi, data=meta_data, method="REML") # 然后按亚组进行元回归（即检验亚组变量是否能解释异质性） res_subgroup <- rma(yi, vi, mods = ~ subgroup, data=meta_data, method="REML") # 查看模型结果，其中`subgroup`水平的系数检验即亚组差异检验 summary(res_subgroup) # 绘制森林图 forest(res, slab = meta_data$study, # 研究标签 xlab = "Hazard Ratio (log scale)", header = c("Study", "HR [95% CI]"), atransf = exp) # 将对数尺度转换回HR尺度 # 在图中添加亚组汇总 addpoly(res_subgroup, row=-1, mlab="Overall Subgroup Difference")

metafor的forest()函数功能也非常强大，可以自定义字体、颜色、布局等，并能轻松添加汇总菱形。

5.3 自动化报告与可重复性

将整个分析流程（数据清洗、模型拟合、结果提取、绘图）封装在一个R Markdown文档或Shiny应用中，是保证分析可重复、结果可追溯的最佳实践。你可以创建参数化报告，只需更改数据源或亚组变量定义，即可一键生成所有分析和图表。这不仅提高了效率，也最大限度地减少了人为操作错误。

最后，我想分享一个自己踩过的坑：早期做亚组分析时，我曾热衷于在森林图上用星号(*)或不同颜色高亮“显著”的亚组，并据此大做文章。后来才深刻理解，亚组分析中，单个亚组内效应的“显著性”与亚组间差异的“显著性”（即交互作用的显著性）是完全不同的两回事。前者可能因样本量小、多重比较而出现假阳性或假阴性；后者才是判断治疗效应是否真正存在异质性的依据。因此，现在我的森林图一定会把交互作用的P值放在最醒目的位置，并在图注中强调“亚组间差异的检验基于模型中的交互项”。这一个小小的改变，让你的分析在审稿人眼中立刻显得更加专业和严谨。