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)
![How to obtain the results of a Tukey HSD post-hoc test in a table showing grouped pairs?](https://cdn-ak-scissors.b.st-hatena.com/image/square/efe2a4b85715116c366abc9aba8c10b9f28dffb9/height=288;version=1;width=512/https%3A%2F%2Fcdn.sstatic.net%2FSites%2Fstats%2FImg%2Fapple-touch-icon%402.png%3Fv%3D344f57aa10cc)