Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I would love to perform a TukeyHSD post-hoc test after my two-way Anova with R, obtaining a table containing the sorted pairs grouped by significant difference. (Sorry about the wording, I'm still new with statistics.)

I would like to have something like this:

enter image description here

So, grouped with stars or letters.

Any idea? I tested the function HSD.test() from the agricolae package, but it seems it doesn't handle two-way tables.

share|improve this question
up vote 11 down vote accepted

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:

yield <- read.dta("")
tx <- with(yield, interaction(fertilizer, irrigation))
amod <- aov(yield ~ tx, data=yield)
HSD.test(amod, "tx", group=TRUE)

This gives the results shown below:

Groups, Treatments and means
a        2.1     51.17547 
ab       4.1     50.7529 
abc      3.1     47.36229 
 bcd     1.1     45.81229 
  cd     5.1     44.55313 
   de    4.0     41.81757 
    ef   2.0     38.79482 
    ef   1.0     36.91257 
     f   3.0     36.34383 
     f   5.0     35.69507 

They match what we would obtain with the following commands:

. webuse yield
. regress yield fertilizer##irrigation
. pwcompare fertilizer#irrigation, group mcompare(tukey)

                      |                           Tukey
                      |     Margin   Std. Err.   Groups
fertilizer#irrigation |
                 1 0  |   36.91257   1.116571    AB    
                 1 1  |   45.81229   1.116571      CDE 
                 2 0  |   38.79482   1.116571    AB    
                 2 1  |   51.17547   1.116571         F
                 3 0  |   36.34383   1.116571    A     
                 3 1  |   47.36229   1.116571       DEF
                 4 0  |   41.81757   1.116571     BC   
                 4 1  |    50.7529   1.116571        EF
                 5 0  |   35.69507   1.116571    A     
                 5 1  |   44.55313   1.116571      CD  
Note: Margins sharing a letter in the group label are
      not significantly different at the 5% level.

The multcomp package also offers symbolic visualization ('compact letter displays', see Algorithms for Compact Letter Displays: Comparison and Evaluation for more details) of significant pairwise comparisons, although it does not present them in a tabular format. However, it has a plotting method which allows to conveniently display results using boxplots. Presentation order can be altered as well (option decreasing=), and it has lot more options for multiple comparisons. There is also the multcompView package which extends those functionalities.

Here is the same example analyzed with glht:

tuk <- glht(amod, linfct = mcp(tx = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))

Treatment sharing the same letter are not significantly different, at the chosen level (default, 5%).

enter image description here

Incidentally, there is a new project, currently hosted on R-Forge, which looks promising: factorplot. It includes line and letter-based displays, as well as a matrix overview (via a level plot) of all pairwise comparisons. A working paper can be found here: factorplot: Improving Presentation of Simple Contrasts in GLMs

share|improve this answer
Thank you so much for this exhaustive answer! I will be trying those different methods as soon as I get a few minutes. Cheers! – chtfn Jul 7 '12 at 13:51
I tried the multcomp package function, put when I use the 'cld()' function I get the error 'Error: sapply(split_names, length) == 2 is not all TRUE' Any idea why? – chtfn Jul 11 '12 at 10:03
@chtfn There seems to be a problem with variable labels. A quick look at the source code indicates that this error message comes from insert_absorb() which tries to extract pair of treatments. You can perhaps try to change the separator you used for coding levels of your interaction term? Without a working example, it's hard to tell what happened. – chl Jul 11 '12 at 10:14
I figured it out: I had points in my genotypes and treatments' names, and as qlht() uses a point to divide the pair names, it freaked out. Thank you so much for all your help, chl! :) – chtfn Jul 11 '12 at 15:02
I noticed today that I now have to add console=TRUE in HSD.test() in order to get the tables, in case someone tries this and sees no result. Probably an update of agricolae. – chtfn Nov 18 '13 at 6:07

There's a function called TukeyHSD that, according to the help file, calculates a set of confidence intervals on the differences between the means of the levels of a factor with the specified family-wise probability of coverage. The intervals are based on the Studentized range statistic, Tukey's "Honest Significant Difference" method. Does this do what you want?

share|improve this answer
Thank you for your response. Yes, I tried this function, but it gives me raw lists of comparisons. What I would like is to see them grouped like in the image in my question, to have a clear view of which group differs to which group, and eventually add the group names on my graphs (for example: a, ab, abc, bc, c) – chtfn Jul 5 '12 at 7:18

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.