你可以做的是指定模型，lme然后glht从multcomp包中使用来做你想做的事。但是，lme 给出的 F 值与标准 ANOVA 略有不同（另请参阅我最近的问题here）。

lme_velocity = lme(Velocity ~ Material, data=scrd, random = ~1|Subject)
anova(lme_velocity)

require(multcomp)
summary(glht(lme_velocity, linfct=mcp(Material = "Tukey")), test = adjusted(type = "bonferroni"))

对于其他对比，然后是 bonferroni，请参见包作者的书multcomp。

您可能还想在 R-mailing list 上看到这篇文章，以及这篇关于在 R 中指定重复测量 ANOVA 的博客文章。

但是，正如我的这个问题所示，我不确定这种方法是否与 ANOVA 相同。此外， glht 只报告z 值而不是通常的t或F值。这似乎也不常见。

到目前为止，我还没有遇到过另一种方法。

如果您想坚持使用该aov()功能，您可以使用emmeans可以处理aovlist（和许多其他）对象的包。

library("emmeans")
# set orthogonal contrasts
options(contrasts = c("contr.sum", "contr.poly"))

aov_velocity <- aov(Velocity ~ Material + Error(Subject / Material), data = scrd)

创建emmGrid对象后如下

emm <- emmeans(aov_velocity, ~ Material)

使用该函数很容易获得所有（事后）成对比较，或者使用包的函数获得pairs()任何所需的对比。通过这些函数的参数可以实现多次测试调整：contrast()emmeansadjust

pairs(emm)  # adjust argument not specified -> default p-value adjustment in this case is "tukey"  

有关这方面的更多信息，我发现详细的emmeans vignettes和文档非常有帮助。

此外，您可以在我的答案中找到一个完整的（可重现的）示例，其中包括有关如何获得正确对比度权重的说明。

但是请注意，如果违反球形度，则使用单变量模型进行事后检验可能会导致反保守的p值。

如果满足球形度，则可以运行双向 ANOVA：

aov_velocity = aov(Velocity~Material+Subject, data=scrd)
posthoc      = TukeyHSD(aov_velocity, 'Material', conf.level=0.95).