R中如何对MICE插补数据执行Tukey、Gammel-Howell事后检验

核心逻辑：事后检验需要在每个插补数据集上单独计算统计量，再通过Rubin规则合并结果，默认的pool()方法仅能处理线性模型原始系数的合并，无法直接输出校正后的两两比较结果。

首先修正你原代码的语法问题：两处函数调用均缺失闭合右括号，修正后的基础流程框架如下：

# 加载基础依赖包
library(mice)
library(multcomp)

方案1：Tukey HSD检验（方差齐性假设满足时使用） #

Tukey检验可直接通过multcomp包的glht函数适配mice的with()语法，无需手动遍历数据集，输出结果自动完成多重比较校正，直接覆盖所有组别两两对比：

# 基础插补流程（补全原代码缺失的括号，增加seed保证结果可复现）
IMP <- mice(data, m=5, maxit=10, seed = 42)
IMP_long <- complete(IMP, include=TRUE, action = 'long')
# 在此处添加总分计算逻辑，例：IMP_long$total_score <- rowSums(IMP_long[, paste0("item_", 1:5)])
IMP_mids <- as.mids(IMP_long, .imp='.imp', .id = '.id')

# 在每个插补数据集上拟合带Tukey对比的线性模型
fit_tukey <- with(IMP_mids, {
  lm_fit <- lm(total_score ~ GroupingVariable)
  glht(lm_fit, linfct = mcp(GroupingVariable = "Tukey"))
})

# 合并插补结果，直接输出三组两两比较的统计量、校正p值
tukey_result <- summary(pool(fit_tukey))
print(tukey_result)

输出结果中会直接列出GroupingVariable: 2 - 1、GroupingVariable: 3 - 1、GroupingVariable: 3 - 2三组比较的结果，对应你需要的1vs2、1vs3、2vs3对比。

方案2：Games-Howell检验（方差齐性假设不满足时使用） #

Games-Howell检验不要求组间方差齐性，无法直接通过glht适配mice的合并逻辑，需要手动遍历插补数据集提取统计量，再按Rubin规则合并：

# 加载额外依赖
library(dplyr)
library(userfriendlyscience)

# 提取所有完成插补的数据集
imp_datasets <- complete(IMP_mids, action = "all")

# 遍历每个插补数据集计算Games-Howell检验统计量
gh_raw <- lapply(imp_datasets, function(dat) {
  gh_test <- oneway(y = dat$total_score, x = dat$GroupingVariable, posthoc = "games-howell")
  posthoc_tab <- gh_test$intermediate$posthoc
  data.frame(
    comp = rownames(posthoc_tab),
    mean_diff = posthoc_tab$diff,
    se = posthoc_tab$se,
    df = posthoc_tab$df
  )
})

# 按Rubin规则合并所有插补集的结果
gh_pooled <- bind_rows(gh_raw, .id = ".imp_id") %>%
  group_by(comp) %>%
  summarise(
    pooled_diff = mean(mean_diff),
    # 计算合并标准误
    within_var = mean(se^2),
    between_var = var(mean_diff),
    m = n_distinct(.imp_id),
    pooled_se = sqrt(within_var + between_var * (1 + 1/m)),
    # 计算检验统计量与校正自由度
    t_val = pooled_diff / pooled_se,
    lambda = ((1 + 1/m) * between_var) / (pooled_se^2),
    df_old = mean(df),
    pooled_df = (df_old + 1)/(df_old + 3) * df_old * (1/lambda),
    # 计算双侧p值
    p_val = 2 * pt(abs(t_val), df = pooled_df, lower.tail = FALSE)
  )

print(gh_pooled)

注意事项 #

• 不要通过提前设置因子对比的方式用默认pool()输出结果，该方法不会对多重比较做p值校正，结果假阳性率会偏高。
• 总分计算必须在长格式的插补数据集上完成，再转回mids对象，不要直接在原始mids对象上计算总分，会导致插补数据逻辑错误。
• 若分组变量不是因子类型，需提前用as.factor()转换，否则事后检验会将其按连续变量处理，无法输出两两比较结果。

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