如何在gtsummary中实现重复测量方差分析与Friedman检验？

问题背景 #

在临床试验中收集了患者Day 0、3、7的空腹血糖重复测量数据，希望通过gtsummary包执行Friedman检验（替代默认的Kruskal-Wallis检验，因为是重复测量数据），同时为另一组数据执行重复测量方差分析，但自定义Friedman检验函数时出现报错。

错误原因分析 #

原自定义函数存在3个关键问题：

• friedman.test的公式写法错误，未正确引用分组变量和随机效应变量
• 数据过滤时遗漏了random变量，导致检验无法按受试者分组
• add_p中传递额外参数的方式不符合gtsummary的要求，需通过test.args传递

修正后的Friedman检验实现 #

步骤1：修正自定义检验函数 #

library(gtsummary)
library(dplyr)

# 修正后的Friedman检验函数
my_friedman <- function(data, variable, by, random, ...) {
  # 保留需要的变量并过滤缺失值
  data <- data %>% 
    dplyr::select(all_of(c(variable, by, random))) %>% 
    dplyr::filter(complete.cases(.))
  
  # 正确构建friedman.test的公式
  friedman.test(
    formula = as.formula(paste0(variable, " ~ ", by, " | ", random)),
    data = data
  )
}

步骤2：调用gtsummary生成表格并添加检验结果 #

# 使用样本数据执行分析
smpl %>% 
  dplyr::select(Participant.ID, value, Day) %>% 
  tbl_summary(by = Day, include = -Participant.ID) %>% 
  add_p(
    test = list(value ~ "my_friedman"),
    test.args = list(value = list(random = "Participant.ID"))
  )

执行后将得到和直接调用friedman.test(value ~ Day | Participant.ID, data = smpl)一致的结果：Friedman chi-squared = 1.2667, df = 2, p-value = 0.5308。

重复测量方差分析的实现 #

如果需要对符合正态性的重复测量数据执行方差分析，可以自定义基于混合效应模型的检验函数：

步骤1：加载所需包并定义函数 #

library(lmerTest)

my_repeated_anova <- function(data, variable, by, random, ...) {
  data <- data %>% 
    dplyr::select(all_of(c(variable, by, random))) %>% 
    dplyr::filter(complete.cases(.))
  
  # 构建混合效应模型（随机截距）
  model <- lmerTest::lmer(
    formula = as.formula(paste0(variable, " ~ ", by, " + (1 | ", random, ")")),
    data = data
  )
  
  # 提取方差分析的p值
  anova_result <- anova(model)
  p_val <- anova_result[by, "Pr(>F)"]
  
  # 返回符合gtsummary要求的结果列表
  list(p.value = p_val, statistic = anova_result[by, "F value"], df = anova_result[by, "NumDF"])
}

步骤2：调用函数生成结果 #

smpl %>% 
  dplyr::select(Participant.ID, value, Day) %>% 
  tbl_summary(by = Day, include = -Participant.ID) %>% 
  add_p(
    test = list(value ~ "my_repeated_anova"),
    test.args = list(value = list(random = "Participant.ID"))
  )

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