The agricolae::HSD.test function does exactly that, but you will need to let it know that you are interested in an interaction term. Here is an example with a Stata dataset: library(foreign) yield <- read.dta("http://www.stata-press.com/data/r12/yield.dta") tx <- with(yield, interaction(fertilizer, irrigation)) amod <- aov(yield ~ tx, data=yield) library(agricolae) HSD.test(amod, "tx", group=TRUE)