使用dplyr和gtsummary计算多结局变量的相对风险及置信区间

一、使用dplyr批量计算RR、置信区间与P值 #

基于你已有的数据处理流程，通过分组嵌套+批量映射可一次性处理所有结局变量，无需逐个筛选子集：

library(tidyverse)
library(epitools)

# 沿用示例数据生成代码
reprex.1 <- iris %>% 
  mutate(
    sepal.length.group = case_when(
      Sepal.Length < 5 ~ "< 5",
      between(Sepal.Length, 5, 6) ~ "5 - 6",
      between(Sepal.Length, 6, 7) ~ "6 - 7",
      Sepal.Length >= 7 ~ "≥ 7"
    ) %>% factor(levels = c("< 5", "5 - 6", "6 - 7", "≥ 7")),
    outcome.var.1 = Sepal.Width < 3,
    outcome.var.2 = Petal.Length < 3,
    outcome.var.3 = Petal.Length < 1
  )

# 按分组统计各结局的true/false数量
reprex.2 <- reprex.1 %>% 
  group_by(sepal.length.group) %>% 
  summarise(
    across(starts_with("outcome.var"), 
           list(n.true = ~sum(.x, na.rm = TRUE), n.false = ~sum(!.x, na.rm = TRUE)))
  ) %>% 
  pivot_longer(!sepal.length.group, 
               names_to = c("outcome.variable", ".value"),
               names_pattern = "(.*)_(.*)") %>% 
  arrange(outcome.variable, sepal.length.group)

# 批量计算所有结局的RR、CI和P值
rr_batch <- reprex.2 %>%
  group_by(outcome.variable) %>%
  nest() %>% # 按结局变量嵌套分组数据
  mutate(
    rr_results = map(data, function(df) {
      # 转换为epitab要求的矩阵格式
      tab_matrix <- df %>% select(n.false, n.true) %>% data.matrix()
      # 计算RR及相关统计量
      epitab(tab_matrix, method = "riskratio")$tab %>%
        as.data.frame() %>%
        select(riskratio, rr.lower = lower, rr.upper = upper, p.value)
    })
  ) %>%
  unnest(c(data, rr_results)) %>% # 展开嵌套结果合并为完整表格
  ungroup()

# 查看批量计算结果
rr_batch

关键逻辑说明 #

1. group_by(outcome.variable) %>% nest()：将每个结局变量对应的分组统计数据单独嵌套为子数据框，实现按结局批量处理的基础。
2. map(data, function(df) {...})：遍历每个嵌套的子数据框，调用epitools::epitab()计算RR、置信区间和P值。
3. unnest()：将嵌套的计算结果展开，合并为包含所有结局统计量的完整数据框。

二、编写适配gtsummary::add_stat()的工具函数 #

gtsummary的add_stat()要求工具函数接收data、variable、by三个核心参数，返回包含格式化统计量的数据框。以下是适配实现：

library(gtsummary)

# 定义适配gtsummary的RR计算函数
add_riskratio <- function(data, variable, by, conf.level = 0.95) {
  # 按分组统计结局的true/false数量
  df_summary <- data %>%
    group_by({{by}}) %>%
    summarise(
      n.true = sum({{variable}}, na.rm = TRUE),
      n.false = sum(!{{variable}}, na.rm = TRUE)
    ) %>%
    ungroup()
  
  # 计算RR、置信区间和P值
  tab_matrix <- df_summary %>% select(n.false, n.true) %>% data.matrix()
  rr_results <- epitab(tab_matrix, method = "riskratio", conf.level = conf.level)$tab %>%
    as.data.frame() %>%
    select(riskratio, rr.lower = lower, rr.upper = upper, p.value)
  
  # 格式化为gtsummary要求的输出格式
  df_summary %>%
    bind_cols(rr_results) %>%
    mutate(
      # 格式化RR与置信区间，参照组标记为(ref)
      stat = case_when(
        is.na(riskratio) ~ "1.00 (ref)",
        TRUE ~ sprintf("%.2f (%.2f-%.2f)", riskratio, rr.lower, rr.upper)
      ),
      # 格式化P值，参照组留空
      p.value = ifelse(is.na(p.value), "", sprintf("%.3f", p.value))
    ) %>%
    select({{by}}, stat, p.value)
}

# 使用示例：生成包含RR的gtsummary表格
tbl_rr <- tbl_summary(
  reprex.1,
  by = sepal.length.group, # 指定分组变量
  include = starts_with("outcome.var"), # 指定需要分析的结局变量
  type = all_logical() ~ "categorical" # 将逻辑变量按分类变量处理
) %>%
  add_stat(
    fns = everything() ~ add_riskratio, # 对所有结局应用RR计算函数
    by = sepal.length.group
  ) %>%
  modify_header(
    stat_1 ~ "**< 5**", 
    stat_2 ~ "**5 - 6**", 
    stat_3 ~ "**6 - 7**", 
    stat_4 ~ "**≥ 7**"
  ) %>%
  modify_footnote(everything() ~ "RR (95% CI); P值")

# 查看表格
tbl_rr

关键逻辑说明 #

1. 使用{{variable}}和{{by}}的tidy eval语法，支持gtsummary动态传递变量。
2. 自定义stat列的格式化规则，确保输出符合医学统计报告的常见格式，参照组明确标记。
3. 直接在tbl_summary()之后调用add_stat()，即可将RR统计量整合到汇总表格中。

内容的提问来源于stack exchange，提问作者obgynfetalmedDK