如何用gtsummary针对第三分类变量的不同水平做配对t检验

问题背景 #

需要针对desg变量的两个水平（x01/x02），分别在group的两个分组（I/F）内执行配对t检验，但当前代码输出结果不符合预期，还出现boundary (singular) fit警告。

错误原因 #

• 配对ID生成错误：当前代码按group分组后用row_number()生成id，破坏了原始的配对关系（真实配对依据是原始ID列，每个ID对应desg的两个水平）。
• 分析逻辑混乱：tbl_summary的by = group是用来对比I/F组的差异，但实际需求是在每个group内部，对比desg两个水平的配对差异，参数设置完全偏离目标。

修正后的代码 #

library(dplyr)
library(gtsummary)

df <- tibble(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,
                     13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19,
                     20, 20, 21, 21, 22, 22, 23, 23, 24, 24),
              desg = c(rep(c("x01", "x02"), 24)),
              group = c(rep(c("I", "I", "F", "F"), 12)),
              values = c(8, 9, 10,  9, 10,  9,  3, 10,  9,  8,  8,  8, 10,  9, 10,
                         7,  9,  6, 10,  8, 10, 10, 10, 10, 10,  8, 10, 10, 10,
                         10, 8, 9, 10, 9, 10, 9, 7, 8, 9, 6, 7, 8, 6, 6, 5, 7,
                         9, 8)
)

# 确保变量类型正确
df <- df %>%
  mutate(
    desg = as.factor(desg),
    group = as.factor(group),
    ID = as.factor(ID)
  )

# 按group分层，每层内针对desg的两个水平做配对t检验
df %>%
  # 确保每个ID都有desg的两个水平（原始数据已满足，此步可验证）
  group_by(ID) %>%
  filter(n() == 2) %>%
  ungroup() %>%
  # 按group分层分析
  tbl_strata(
    strata = group,
    .tbl_fun = function(data) {
      data %>%
        tbl_summary(
          by = desg,  # 对比desg的两个水平
          include = values,
          statistic = values ~ "{mean} ({sd})"  # 显示均值标准差，适合t检验
        ) %>%
        add_p(
          test = values ~ "paired.t.test",  # 指定配对t检验
          group = ID  # 指定配对依据是原始ID
        )
    }
  )

结果说明 #

修正后的代码会生成两个分层表格（分别对应group=I和group=F），每个表格展示desg两个水平下values的描述性统计，以及对应的配对t检验p值，完全符合“针对desg两个水平分别在group分组内做配对检验”的需求，同时消除了之前的警告。

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