If you want to stick with the
aov() function you can use the
emmeans package which can handle
aovlist (and many other) objects.
# set orthogonal contrasts
options(contrasts = c("contr.sum", "contr.poly"))
aov_velocity <- aov(Velocity ~ Material + Error(Subject / Material), data = scrd)
After creating an
emmGrid object as follows
emm <- emmeans(aov_velocity, ~ Material)
it is very easy to get all (post hoc) pairwise comparisons using the
pairs() function or any desired contrast using the
contrast() function of the
emmeans package. Multiple-testing adjustments can be achieved via the
adjust argument of these functions:
pairs(emm) # adjust argument not specified -> default p-value adjustment in this case is "tukey"
For more information on this I found the detailed emmeans vignettes and the documentation to be very helpful.
Also, you can find a complete (reproducible) example including a description on how to get the correct contrast weights in my answer here.
Note, however, that using a univariate model for the post hoc tests can result in anti-conservative p-values if sphericity is violated.