如何在gtsummary中为ancova配置参数生成组间对比P值列

核心思路 #

gtsummary的ancova检验默认返回整体组间差异的P值，无法直接给出与参考组的两两对比P值。我们需要自定义函数，拟合调整协变量的线性模型（即ANCOVA），提取组2、组3与参考组（组1）的对比P值，再通过add_stat将结果整合到汇总表中。

步骤1：加载依赖包 #

library(gtsummary)
library(dplyr)
library(broom) # 提取模型结果
# 如果需要事后比较校正，可加载emmeans包
# library(emmeans)

步骤2：自定义P值计算函数 #

这个函数会拟合包含协变量的线性模型，提取组2、组3相对于组1的系数P值：

compare_ref_group <- function(data, variable, adj_vars) {
  # 构建模型公式：因变量 ~ 分组 + 协变量
  formula <- as.formula(paste(variable, "~ group +", paste(adj_vars, collapse = " + ")))
  # 拟合线性模型（ANCOVA）
  model <- lm(formula, data = data)
  # 提取group2和group3对应的P值
  p_vals <- tidy(model) %>%
    filter(stringr::str_detect(term, "^group")) %>%
    pull(p.value)
  
  # 返回两个P值，对应组2 vs组1、组3 vs组1
  tibble(
    p_group2_vs_1 = p_vals[1],
    p_group3_vs_1 = p_vals[2]
  )
}

若需要对P值做多重比较校正（如Bonferroni），可改用emmeans实现：

compare_ref_group_emmeans <- function(data, variable, adj_vars) {
  formula <- as.formula(paste(variable, "~ group +", paste(adj_vars, collapse = " + ")))
  model <- lm(formula, data = data)
  # 计算边际均值并指定对比方式
  emm <- emmeans(model, ~ group)
  contrasts <- contrast(emm, 
                        list(group2_vs_1 = c(-1,1,0), group3_vs_1 = c(-1,0,1)),
                        adjust = "bonferroni") # 多重比较校正
  # 提取校正后的P值
  p_vals <- tidy(contrasts) %>% pull(p.value)
  
  tibble(
    p_group2_vs_1 = p_vals[1],
    p_group3_vs_1 = p_vals[2]
  )
}

步骤3：生成汇总表并添加自定义P值列 #

final_table <- df %>%
  select(group, y1, y2, x1, x2) %>%
  # 生成基础分组汇总表
  tbl_summary(by = group, missing = "no") %>%
  add_n() %>%
  # 添加自定义的参考组对比P值列
  add_stat(
    fns = list(
      # 为y1、y2分别指定计算函数，传入协变量
      y1 ~ ~compare_ref_group(data = ., variable = "y1", adj_vars = c("x1", "x2")),
      y2 ~ ~compare_ref_group(data = ., variable = "y2", adj_vars = c("x1", "x2"))
      # 若用emmeans版本，替换为compare_ref_group_emmeans
    ),
    # 定义新列的标签
    labels = list(p_group2_vs_1 = "Group 2 vs 1 (P)", p_group3_vs_1 = "Group 3 vs 1 (P)")
  ) %>%
  # 调整列顺序，将P值列放在末尾
  modify_column_order(c(stat_1, stat_2, stat_3, n, p_group2_vs_1, p_group3_vs_1)) %>%
  # 格式化P值为手稿常用样式（如<0.001）
  modify_fmt_fun(c(p_group2_vs_1, p_group3_vs_1) ~ style_pvalue)

# 查看最终表格
final_table

关键知识点解释 #

1. add_stat的作用：gtsummary中用于添加自定义统计量的核心函数，支持返回多列结果，完美解决替换/新增P值列的需求。
2. ANCOVA的本质：代码中拟合的lm(y ~ group + x1 + x2)就是ANCOVA模型，通过控制协变量x1、x2的影响，得到分组变量对因变量的净效应P值。
3. 参考组的设定：因为group因子的水平是c("1","2","3")，模型默认以group1为参考组，所以提取的系数P值直接对应组2、组3与组1的对比。

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