重复测量置换方差分析的事后检验适配方法咨询

针对你用aovperm()构建的重复测量置换方差分析模型，无法用lsmeans()/glht()做事后检验的问题，以下是几种可行的替代方案：

1. 利用perm包内置的置换事后检验功能 #

aovperm()所属的perm包提供了permTest()函数，可以直接针对模型中的效应做置换-based的事后比较，适配嵌套/重复测量结构：

library(perm)

# 定义感兴趣的两两对比（以treatment的3组比较为例）
treatment_contrasts <- list(
  T1_vs_T2 = c(1, -1, 0),
  T1_vs_T3 = c(1, 0, -1),
  T2_vs_T3 = c(0, 1, -1)
)

# 针对species×treatment×month的交互效应做置换事后检验
# 指定blocks参数保留植株内的重复测量相关性
posthoc_perm <- permTest(area.aovp,
                         effect = "species:treatment:month",
                         contrasts = treatment_contrasts,
                         np = 2000,
                         blocks = growth$number)

# 查看检验结果
summary(posthoc_perm)

2. 切换到混合模型框架+置换事后检验 #

将重复测量数据用线性混合模型拟合，再结合emmeans和置换检验实现事后比较，能更好地处理个体嵌套结构：

library(lme4)
library(emmeans)
library(permute)

# 拟合线性混合模型：植株ID作为随机截距，固定效应包含所有交互项
lmer_model <- lmer(p_area ~ species*treatment*month + (1|number), data = growth)

# 获取各处理组合的边际均值
emm_obj <- emmeans(lmer_model, ~ species*treatment*month)

# 执行置换-based的两两比较，按植株ID分块保持相关性
set.seed(123)
pairwise_results <- pairwise.emmeans(emm_obj,
                                     adjust = "none",
                                     perm = TRUE,
                                     nperm = 2000,
                                     blocks = growth$number)

# 手动做Bonferroni校正控制多重比较误差
pairwise_results$contrasts$p.value <- p.adjust(pairwise_results$contrasts$p.value, method = "bonferroni")
print(pairwise_results$contrasts)

3. 非参数事后检验（适配小样本+方差不齐） #

如果置换检验的计算成本过高，可针对分组做非参数检验，结合多重比较校正：

library(dplyr)
library(rstatix)

# 按物种+处理分组，对同一植株的不同月份做Friedman检验（组内重复测量的非参数检验）
growth %>%
  group_by(species, treatment) %>%
  friedman_test(p_area ~ month | number) %>%
  # 事后配对Wilcoxon检验，用Bonferroni校正p值
  pairwise_wilcox_test(p_area ~ month, paired = TRUE, p.adjust.method = "bonferroni")

注意事项 #

• 置换次数np建议设置为10000以上，结果稳定性更好（需权衡计算时间）。
• 所有涉及置换的操作都要通过blocks = growth$number指定植株ID为块，避免破坏重复测量的个体内相关性。
• 多重比较校正方法可根据需求调整（如Holm法比Bonferroni更宽松）。

内容的提问来源于stack exchange，提问作者蔡譯禎