如何用gtsummary生成含混合模型p值与交互项的交叉试验分析表？

需求背景 #

开展交叉试验，为患者随机分配treatment 1和treatment 2，需用tbl_summary生成包含以下内容的表格：

• 处理效应列（treatment1 - treatment2的估计值）
• 混合模型组间比较的p值
• 序列交互作用的p值

已通过lme4构建混合模型，emmeans完成组间比较，现需将结果整合到目标表格中。

完整实现代码 #

# 加载所需包
library(emmeans)
library(lme4)
library(lmerTest)
library(gtsummary)
library(gt)
library(glue)

# 构建数据集
df <- data.frame (
  record_id  = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12),
  treatment = c(1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2),
  treatment_sequence = c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
  treatment_response = c(-43.5, 135.0, 8.4, -7.2, 99.0, 159.0, 12.0, -27.0, 3.0, 12.0, -15.0, 91.5, 6.0, -9.0, 177.0, 27.0, 52.8, -54.0, -50.7, 63.0, -9.0, 186.0, -72.0, 15.0)
)

# 构建混合模型
df_mm <- lmer(treatment_response ~ as.factor(treatment)*treatment_sequence + (1|record_id), data=df)

# 提取所需统计量
## 提取组间比较的处理效应和p值
emm_result <- emmeans(df_mm, list(pairwise ~ treatment), adjust = "bonferroni")
treatment_effect <- round(emm_result$`pairwise differences of treatment`$estimate, 1)
treatment_p <- round(emm_result$`pairwise differences of treatment`$p.value, 4)

## 提取序列交互作用的p值
anova_result <- anova(df_mm)
sequence_interaction_p <- round(anova_result$`Pr(>F)`[3], 2)

# 构建tbl_summary表格
tbl <- df %>%
  select(treatment_response) %>%
  tbl_summary(
    statistic = list(all_continuous() ~ "{mean} ({sd})"),
    label = list(treatment_response ~ "Treatment Response")
  ) %>%
  # 添加处理效应列
  add_stat(
    fns = list(all_continuous() ~ function(x) glue("{treatment_effect}")),
    label = "Treatment Effect (Trt1 - Trt2)"
  ) %>%
  # 添加组间比较p值列
  add_stat(
    fns = list(all_continuous() ~ function(x) glue("{treatment_p}")),
    label = "Group Comparison p-value"
  ) %>%
  # 添加序列交互p值列
  add_stat(
    fns = list(all_continuous() ~ function(x) glue("{sequence_interaction_p}")),
    label = "Sequence Interaction p-value"
  ) %>%
  # 优化表头显示
  modify_header(stat_0 ~ "Mean (SD)") %>%
  # 确保统计量仅对应目标行
  modify_table_body(
    ~ .x %>%
      mutate(
        across(c(stat_1, stat_2, stat_3), ~ ifelse(variable == "treatment_response", ., NA))
      )
  ) %>%
  # 转换为gt对象支持更多样式自定义
  as_gt()

# 输出表格
tbl

代码说明 #

1. 统计量提取：从emmeans结果中提取treatment1与treatment2的差值（处理效应）和组间比较p值；从anova结果中提取处理与序列交互项的p值。
2. 表格构建：
   • 先用tbl_summary生成基础描述性统计表格
   • 通过add_stat依次添加处理效应、组间比较p值、序列交互p值三列
   • 用modify_header和modify_table_body优化表格显示逻辑，保证统计量对应正确行
3. 样式扩展：转换为gt对象后，可根据需求进一步调整表格字体、边框、对齐方式等样式。

内容的提问来源于stack exchange，提问作者Kristoffer Berg Hansen