如何为多组混合设计ANOVA调整alpha值或执行Bonferroni校正？

我正在分析一项3×2×2混合因子设计的实验数据：

• 3水平被试内变量：Viewing number
• 另外两个2水平变量（一个被试内、一个被试间）
• 分析需求：
  1. 先运行全模型3×2×2 ANOVA
  2. 拆分Viewing number的不同水平组合，运行一系列2×2×2 ANOVA
  3. 需对多组ANOVA的结果进行多重检验校正，但不清楚如何在现有工具中设置调整后的alpha值或校正方法

jmv::anovaRM() #

我用该函数分析过分差分数，但改用原始数据后无法拆分Viewing number的水平，且无法手动指定校正后的alpha：

model.diff2.1 <- anovaRM(data = dat1, 
                 rm = list(list(label = 'Armed',
                                levels = c('armeddiff2.1', 'unarmdiff2.1'))), 
                 rmCells = list(list(measure = 'armeddiff2.1', cell = 'armeddiff2.1'),
                                list(measure = 'unarmdiff2.1', cell = 'unarmdiff2.1')),
                 rmTerms = list('Armed'),
                 bs = 'condition',
                 bsTerms = list('condition'),
                 effectSize = c('partEta'),
                 leveneTest = TRUE,
                 spherTests = TRUE,
                 spherCorr = c('none','GG'), #no correction and gg correction
                 postHoc = list('Armed', 'condition'),
                 postHocCorr = list('holm','tukey'),
                 emMeans = ~ Armed + condition + Armed:condition,
                 emmTables = T)

ez::ezANOVA() #

该函数指定多被试内变量更简便，但同样无法清晰设置多重检验校正，也不知道如何修改alpha值：

anova_results <- ezANOVA(
  data = dat1long,
  dv = .(danger),
  wid = .(index),
  within = .(armed, view),
  between = .(condition),
  detailed = TRUE
)

1. 手动对多组ANOVA的p值进行校正 #

因为jmv和ezANOVA都没有内置的多ANOVA校正参数，你可以先批量运行所有需要的ANOVA，收集关键p值后用R内置的p.adjust()函数做校正：

步骤示例：

• 拆分Viewing number的水平组合（比如组合1：View1&View2；组合2：View1&View3；组合3：View2&View3）
• 对每个子集运行ANOVA，提取主效应和交互效应的p值
• 用p.adjust()进行校正（支持holm、bonferroni、fdr等方法）

# 假设dat1long是原始数据长格式，view是Viewing number变量
# 定义需要的Viewing number组合
view_combinations <- list(
  c("View1", "View2"),
  c("View1", "View3"),
  c("View2", "View3")
)

# 存储所有ANOVA的p值
all_p_values <- c()

# 批量运行ANOVA
for (comb in view_combinations) {
  subset_data <- dat1long[dat1long$view %in% comb, ]
  anova_res <- ezANOVA(
    data = subset_data,
    dv = .(danger),
    wid = .(index),
    within = .(armed, view),
    between = .(condition),
    detailed = TRUE
  )
  # 提取所有效应的p值（比如主效应和交互效应）
  p_vals <- anova_res$ANOVA$p
  all_p_values <- c(all_p_values, p_vals)
}

# 对所有p值进行校正（这里用holm方法）
adjusted_p <- p.adjust(all_p_values, method = "holm")

# 可以将原始p和校正后p合并查看
cbind(original_p = all_p_values, adjusted_p = adjusted_p)

2. 用事后检验/简单效应替代多组ANOVA #

其实你的需求可以通过全模型ANOVA + 针对Viewing number的简单效应检验实现，更符合统计规范，也能直接用内置校正：

用emmeans包实现（配合ezANOVA的全模型）

library(emmeans)

# 先跑全模型3×2×2 ANOVA
full_anova <- ezANOVA(
  data = dat1long,
  dv = .(danger),
  wid = .(index),
  within = .(armed, view),
  between = .(condition),
  detailed = TRUE
)

# 拟合线性混合模型（用于emmeans分析）
library(lme4)
lmer_model <- lmer(danger ~ armed * view * condition + (1 + armed + view | index), data = dat1long)

# 针对Viewing number的不同水平组合，做简单效应检验并校正
# 比如比较每个Viewing number水平下的armed×condition交互
emm <- emmeans(lmer_model, ~ armed * condition | view)
pairwise_emmeans <- pairs(emm, adjust = "holm") # 用holm校正
print(pairwise_emmeans)

3. 修改alpha阈值的替代方式 #

如果一定要手动指定调整后的alpha（比如Bonferroni校正后alpha=0.05/3），你可以在提取ANOVA结果后，自行对比p值和调整后的alpha：

# 示例：Bonferroni校正，假设跑3组ANOVA，调整后alpha=0.05/3≈0.0167
adjusted_alpha <- 0.05 / length(view_combinations)

# 提取某组ANOVA的p值后，判断是否显著
if (anova_res$ANOVA$p[1] < adjusted_alpha) {
  cat("主效应显著")
}

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