Introduction
一次要看很多个效应量的时候：十个亚组各自的 OR、三个模型对同一个暴露的估计、二十项研究的 meta 分析。每行一个估计，一条线一个区间，竖着排开跟同一条无效线比。
只有一两个效应量就别用了——两行的森林图不如把数字直接写进正文。
读的时候按这个顺序：
- 先找无效线。 比值型指标（HR / OR / RR）在 1，差值型（beta / MD）在 0。
- 看区间碰没碰到它。 碰到了是"这份数据没能排除无效"，不是"证明了无效"。
- 看区间的宽度。 宽是样本少或事件少，跟效应大小无关。一条又宽又偏的线说的是"不知道"，不是"很强"。
- 别横向比亚组之间的显著性。 A 组显著、B 组不显著，不等于两组效应不同。
Example Data
森林图要的输入是：每行一个效应量——点估计 + 置信区间上下限，外加要一起显示的文本列。
这是这张图唯一麻烦的地方，画的部分反而是机械的。forestploter::forest() 要的不是一张表，而是拆开的四样东西：
| 要什么 | 说明 |
|---|---|
data |
只放文本列的 data.frame，逐列变成表格里的文字 |
| 占位列 | data 里必须有一列是全空格，CI 图就画在它的位置；空格数量决定 CI 面板的宽度 |
ci_column |
占位列是第几列 |
est / lower / upper |
三个独立的数值向量，按行和 data 对齐 |
还有两样是 forest() 不管、得自己造的：HR (95% CI) 那一列文本（sprintf 拼），以及 p 值的两个版本——数值的用来判断要不要加粗，字符串的用来显示（带 <0.001 截断）。
plot_forest() 把后面这些全接管了：只给文本列和三个向量，它自动插入占位列和 CI 文本列。所以它的 data 比 forest() 的 data 少两列，ci_column 指的也是插入之前的位置。这是两者最容易混的地方。
forest_subgroups 和 forest_models 分开保存两种分析结果。前者是 9 行亚组估计，后者是 BMI 三个水平 × 三种校正模型的 9 行 Cox 结果。CSV 只放分析数据；亚组标题与缩进在绘图时根据 group 生成。
library(forestploter)
library(grid)
library(biopalette)
# 公开地址，和数据集页上「下载 CSV」给的是同一个 —— 不写仓库相对路径：
# 那个目录不进仓库，读者 clone 下来也没有这个文件，这段代码就跑不了
subgroups <- read.csv(
"https://assets.evanzhou.org/tessera/csv/forest_subgroups.csv"
)
models <- read.csv(
"https://assets.evanzhou.org/tessera/csv/forest_models.csv"
)
# CSV 没有展示专用的标题行。Overall 直接保留；其余每个 group 在这里插入
# 一行没有效应量的标题，标题行展示该组的交互作用 p 值。
rows <- do.call(rbind, lapply(unique(subgroups$group), function(g) {
members <- subgroups[subgroups$group == g, ]
leaves <- data.frame(
label = members$level,
depth = if (g == "Overall") 0L else 1L,
n = members$n,
p_value = members$p_value,
p_interaction = NA_real_,
est = members$HR,
lower = members$CI_lower,
upper = members$CI_upper
)
if (g == "Overall") return(leaves)
header <- data.frame(
label = g, depth = 0L, n = NA_integer_, p_value = NA_real_,
p_interaction = unique(members$p_interaction)[1L],
est = NA_real_, lower = NA_real_, upper = NA_real_
)
rbind(header, leaves)
}))
# 文本列和三个数值向量分开 —— 这就是上面那张表说的形状
sg <- data.frame(
Subgroup = rows$label,
`No. of patients` = ifelse(is.na(rows$n), "",
format(rows$n, big.mark = ",", trim = TRUE)),
p_value = rows$p_value,
p_interaction = rows$p_interaction,
check.names = FALSE # 列名里有空格和点，不让 R 改写
)
est <- rows$est
lower <- rows$lower
upper <- rows$upper
# 缩进从刚才插入的标题行推，CSV 本身不承担展示逻辑
ind <- rows$depth
nrow(sg)
#> 12Palettes
亚组图用 heat_light 的暖红标显著、冷蓝标不显著；模型比较图用 three_body 的三个颜色依次表示未校正、年龄与性别校正及完全校正。颜色全部来自已有色板，不在 recipe 里另写一套 HEX。
按显著性给 CI 上色是有争议的：它把一个连续量（p 值）压成了两档，而 p = 0.049 和 p = 0.051 之间没有实质区别。用它的理由只有一个——图上行数多到一眼扫不完时，颜色帮读者先定位。行数少就别上色，黑色一种更干净。
斑马底纹的作用是横向对齐：行多的时候，眼睛从最左边的标签扫到最右边的 p 值容易串行。它刻意是无彩色的浅灰，不取自任何色板——底纹一旦带上色相，就会被当成第三种编码。
heat_light <- get_palette("heat_light", type = "qualitative")
SIG <- heat_light[1] # 暖红，p < 0.05
NONSIG <- heat_light[2] # 冷蓝
three_body <- get_palette("three_body", type = "qualitative")
MODEL_COL <- setNames(
three_body,
c("Unadjusted", "Age and sex adjusted", "Fully adjusted")
)
STRIPE <- "#F0F0F0" # 斑马底纹，浅到几乎看不见就对了Recipe
| No. | Method | Input Data | Palettes |
|---|---|---|---|
| 1 | forestploter |
d（凑好形状的表） |
— |
| 2 | ukbflow |
sg |
— |
| 3 | ukbflow |
sg |
SIG / NONSIG |
| 4 | ukbflow |
models（9 行原表） |
three_body |
1 · 手写 forest()
这一段的全部内容就是把 sg 改造成 forest() 认的形状：造 CI 文本列、造空格占位列、把 p 值格式化成字符串、再按顺序排好列。
#| fig: manual
#| fig-units: cm
#| fig-width: 24
#| fig-height: 10
d <- sg
# ① CI 文本列 —— forest() 不会替你生成
d$`HR (95% CI)` <- ifelse(is.na(est), "",
sprintf("%.2f (%.2f, %.2f)", est, lower, upper))
# ② 空格占位列 —— 空格数量就是 CI 面板的宽度，不是随手写的 22。
# 行数一变最容易忘掉调它，结果 CI 被挤成一小段
d$` ` <- strrep(" ", 22)
# ③ p 值转成显示用的字符串（数值那份留着判断加粗）
d$p_value <- ifelse(is.na(sg$p_value), "",
ifelse(sg$p_value < 0.001, "<0.001",
sprintf("%.3f", sg$p_value)))
d$p_interaction <- ifelse(is.na(sg$p_interaction), "",
sprintf("%.3f", sg$p_interaction))
# ④ 排好列序，占位列在第 3 位
d <- d[, c("Subgroup", "No. of patients", " ", "HR (95% CI)",
"p_value", "p_interaction")]
names(d) <- c("Subgroup", "No. of patients", "", "HR (95% CI)",
"P value", "P interaction")
p_manual <- forest(
d,
est = est,
lower = lower,
upper = upper,
ci_column = 3, # 指的是上面那个 " " 列
ref_line = 1,
xlim = c(0.7, 2.6),
# 不给 ticks_at 会在 xlim 上五等分：c(0.7, 2.6) 得到 0.7 / 1.175 / 1.65 / 2.125 / 2.6,
# 无效线 1 甚至不在刻度上
ticks_at = c(0.7, 1, 1.5, 2, 2.5)
)
# forest() 按文字自动算出的列宽偏松；这里按这六列的实际内容收紧。
p_manual$widths <- grid::unit(c(3, 35, 35, 52, 40, 22, 30, 3), "mm")
p_manual <- add_border(
p_manual, part = "header", row = 1, where = "top",
gp = gpar(lwd = 3)
)
p_manual <- add_border(
p_manual, part = "header", row = 1, where = "bottom",
gp = gpar(lwd = 3)
)
p_manual <- add_border(
p_manual, part = "body", row = nrow(d), where = "bottom",
gp = gpar(lwd = 3)
)
plot(p_manual)
