线性混合模型的非参数ANOVA及功效计算相关技术咨询

1. 支持bootstrap ANOVA及事后比较的R包与方法 #

可以通过以下工具组合实现需求：

• lmeresampler + car + emmeans
  利用lmeresampler生成模型的bootstrap样本，结合car包的Anova()计算bootstrap版本的ANOVA统计量，再通过emmeans的内置bootstrap功能完成事后两两比较：

  # 加载依赖包
  library(lmeresampler)
  library(car)
  library(emmeans)
  
  # 对原lmer模型执行残差bootstrap
  boot_mod <- bootstrap(mylmm, type = "residual", B = 1000)
  
  # 计算bootstrap ANOVA结果
  boot_anova_list <- lapply(boot_mod$bootstraps, function(x) Anova(x, type = 3))
  boot_pvals <- sapply(boot_anova_list, function(x) x$`Pr(>Chisq)`)
  orig_anova <- Anova(mylmm, type = 3)
  # 计算bootstrap p值：统计量大于原模型的比例
  boot_p <- rowMeans(t(boot_pvals) > orig_anova$`Pr(>Chisq)`)
  
  # 事后两两比较的bootstrap推断
  emm_obj <- emmeans(mylmm, ~ group * exam)
  boot_emm <- bootEM(emm_obj, B = 1000)
  # 两两比较并做多重比较校正
  pairs(boot_emm, adjust = "tukey")
  

• boot包自定义实现：如果需要更灵活的控制，可使用boot包编写自定义函数，手动提取ANOVA统计量或事后比较的差异值，再计算置信区间或p值。

2. bootstrap/非参数ANOVA的功效计算 #

需要计算功效，小样本场景下，功效能直观反映检验检测真实效应的概率。推荐采用模拟法实现：

library(boot)
library(lmerTest)

# 定义数据模拟函数：基于原数据分布生成带指定效应的新数据
simulate_data <- function(orig_data, effect_magnitude) {
  new_data <- orig_data
  # 为指定group-exam组合添加效应（可按需修改）
  new_data$score <- new_data$score + effect_magnitude * 
    (new_data$group == "B" & new_data$exam == "Exam1")
  return(new_data)
}

# 模拟循环计算功效
n_simulations <- 1000
power_results <- numeric(n_simulations)

for (i in 1:n_simulations) {
  sim_dat <- simulate_data(mydata, effect_magnitude = 0.5)
  sim_lmm <- lmer(score ~1+group+exam+group*exam+(1|participant), data = sim_dat)
  # 执行bootstrap ANOVA
  boot_mod <- bootstrap(sim_lmm, type = "residual", B = 500)
  boot_anova_list <- lapply(boot_mod$bootstraps, function(x) Anova(x, type = 3))
  orig_anova <- Anova(sim_lmm, type = 3)
  boot_p <- rowMeans(t(sapply(boot_anova_list, function(x) x$`Pr(>Chisq)`)) > orig_anova$`Pr(>Chisq)`)
  # 记录检验是否显著
  power_results[i] <- as.integer(boot_p["group:exam"] < 0.05)
}

# 计算最终功效
final_power <- mean(power_results)

也可使用nparcomp等非参数工具，但自定义模拟更贴合你的数据结构。

3. 事后两两比较的功效计算（simr + emmeans） #

simr无法直接计算事后比较的功效，必须结合emmeans定义检验函数来实现：

library(simr)
library(emmeans)

# 假设已拟合好lme模型mod
# 扩展模型的模拟能力（若需调整样本量可修改n参数）
mod <- extend(mod, n = 10)

# 定义事后比较的检验函数：返回两两比较的p值
emm_test_func <- function(model) {
  emm_obj <- emmeans(model, ~ group * exam)
  comps <- pairs(emm_obj, adjust = "tukey")
  comps$p.value
}

# 模拟计算功效
power_out <- powerSim(mod, test = emm_test_func, nsim = 1000)
# 查看功效结果
summary(power_out)

simr负责生成带效应的模拟数据集，emmeans处理事后比较的统计推断，两者结合即可得到两两比较的功效。

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