forest

数据

亚组分析的典型形状：总体一行，然后每个亚组一个标题行加若干水平行。标题行的 est 是 NA。

library(forestploter)
library(grid)

sg <- data.frame(
  Subgroup = c("Overall",
               "Age", "<60", ">=60",
               "Sex", "Male", "Female",
               "BMI", "<25", "25-30", ">=30"),
  `No. of patients` = c("4,521", "", "2,140", "2,381", "", "2,260", "2,261",
                        "", "1,502", "1,808", "1,211"),
  p_value = c(0.003, NA, 0.012, 0.041, NA, 0.008, 0.19,
              NA, 0.44, 0.021, 0.0004),
  check.names = FALSE
)

est   <- c(1.42, NA, 1.31, 1.55, NA, 1.61, 1.18, NA, 1.09, 1.38, 1.87)
lower <- c(1.13, NA, 1.06, 1.02, NA, 1.13, 0.92, NA, 0.88, 1.05, 1.32)
upper <- c(1.78, NA, 1.62, 2.36, NA, 2.30, 1.51, NA, 1.35, 1.81, 2.65)

# 标题行缩进 0，水平行缩进 1
ind <- c(0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1)

一、手写：把四样东西凑齐

这一段的全部内容就是把上面那张表改造成 forest() 认的形状：造 CI 文本列、造空格占位列、把 p 值格式化成字符串、然后按顺序排好列。

#| fig: manual
#| fig-units: cm
#| fig-width: 18
#| fig-height: 12
d <- sg

# ① CI 文本列 —— forest() 不会替你生成
d$`OR (95% CI)` <- ifelse(is.na(est), "",
                          sprintf("%.2f (%.2f, %.2f)", est, lower, upper))

# ② 空格占位列 —— 空格数量就是 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)))

# ④ 排好列序，占位列在第 3 位
d <- d[, c("Subgroup", "No. of patients", " ", "OR (95% CI)", "p_value")]

plot(forest(
  d,
  est       = est,
  lower     = lower,
  upper     = upper,
  ci_column = 3,          # 指的是上面那个 " " 列
  ref_line  = 1,
  xlim      = c(0.5, 3),
  ticks_at  = c(0.5, 1, 2, 3)
))

出来的是能用的图，但表头还顶着 p_value 这种原始列名，缩进、加粗、配色、三线边框全都没有——那些是 edit_plot() 和 add_border() 的活，再写三四十行。

二、一行版

同一份数据交给 plot_forest()。注意 data 传的是 sg 本身——没有占位列、没有 CI 文本列，那两列由它插入。

#| fig: wrapped
#| fig-units: cm
#| fig-width: 22
#| fig-height: 20
library(ukbflow)

plot(plot_forest(
  data      = sg,
  est       = est,
  lower     = lower,
  upper     = upper,
  ci_column = 3,                     # 插入之前的位置
  p_cols    = "p_value",             # 这一列按 p 值格式化并加粗
  header    = c("Subgroup", "No. of patients", "",
                "OR (95% CI)", "P value"),
  xlim      = c(0.5, 3),
  ticks_at  = c(0.5, 1, 2, 3),
  save      = FALSE
))

三线边框、斑马底纹、p < 0.05 加粗、列对齐都是默认就有的。header 的长度是 ncol(sg) + 2，因为那两列已经插进去了；占位列那一位传空字符串。

三、亚组图

加上缩进和按显著性配色。ci_col 接受一个和行数等长的向量，所以显著性配色就是一句 ifelse()。

#| fig: subgroup
#| fig-units: cm
#| fig-width: 22
#| fig-height: 20
plot(plot_forest(
  data      = sg,
  est       = est,
  lower     = lower,
  upper     = upper,
  ci_column = 3,
  p_cols    = "p_value",
  indent    = ind,                   # 标题行不缩进，水平行缩进一级
  ci_col    = ifelse(!is.na(sg$p_value) & sg$p_value < 0.05,
                     "#B2182B", "#4575B4"),
  header    = c("Subgroup", "No. of patients", "",
                "OR (95% CI)", "P value"),
  xlim      = c(0.5, 3),
  ticks_at  = c(0.5, 1, 2, 3),
  arrow_lab = c("Lower risk", "Higher risk"),
  save      = FALSE
))

bold_label 不用传：它默认从 indent 推——缩进为 0 的是父行，自动加粗。

四、结果表直通

上面那三张都是先有效应量、再画图。但效应量本来就是回归跑出来的，那张结果表里已经有点估计和上下限了。

plot_forest() 认这种表：est / lower / upper 都不给的时候，它去看列名，认出来就自己派生。识别规则很简单——有 CI_lower 和 CI_upper，再加上 HR / OR / beta / SHR 里的任意一个。连参考线放哪都跟着效应类型走：比值型放 1，beta 放 0。

所以不必真的跑一遍回归，手搭一张同样形状的表就能验证：

#| fig: assoc
#| fig-units: cm
#| fig-width: 28
#| fig-height: 18
res <- data.frame(
  model    = rep(c("Unadjusted", "Age- and sex-adjusted", "Fully adjusted"),
                 each = 3),
  term     = rep(c("Overweight", "Obese", "Severely obese"), times = 3),
  n_events = rep(c(412, 268, 97), times = 3),
  HR       = c(1.18, 1.52, 1.94, 1.14, 1.44, 1.81, 1.09, 1.36, 1.72),
  CI_lower = c(1.02, 1.27, 1.48, 0.99, 1.20, 1.37, 0.94, 1.12, 1.28),
  CI_upper = c(1.37, 1.82, 2.54, 1.32, 1.73, 2.39, 1.27, 1.65, 2.31),
  p_value  = c(0.028, 0.0002, 0.00004, 0.067, 0.0004, 0.0002,
               0.25, 0.0021, 0.0004)
)

plot(plot_forest(
  res,
  xlim     = c(0.9, 3),      # 自动推的范围留白过多，比值型指标上会让刻度没法读
  ticks_at = c(1, 1.5, 2, 3),
  save     = FALSE
))

这条路省掉的正是本页开头那一整节：输入的形状根本不用你去凑，est / lower / upper 和表头都是派生的。

一点要注意：自动推的 xlim 会为了留白把下界压得很低，比值型指标上刻度会挤成一团（这张表实测推出的下界是 0.0035，而 CI 最小才 0.17），所以 xlim 和 ticks_at 通常还是自己给。

版式要再改，就回到前面几节——把 label_col / show_cols 或者干脆 est / lower / upper 显式传进去。任何一个显式给了的参数都会盖掉派生值。

存成文件

plot_forest(save = TRUE, dest = "forest") 会一次写出 png / pdf / jpg / tiff 四份，300 dpi。宽度默认 20 cm，高度不给的话按 行数 × 0.9 + 3 估。

手写路线就是常规的 ggplot2::ggsave() 或者开图形设备再 plot()——森林图是个 gtable，不是 ggplot 对象，ggsave() 要显式把它当 plot 传。

两种做法怎么选

1. 常规亚组图、多模型对比 → plot_forest()。 占位列、CI 文本、p 值格式、加粗、斑马纹、三线边框全在里面，剩下的只是填数据。
2. 要对数轴、多列 CI、或者版式很特殊 → 手写 forest()。 x_trans = "log" 这类参数封装没暴露，而比值型指标本来就该用对数轴。