# One-way ANOVA, clustering levels using Tukey Kramer HSD

Suppose I am doing a one-way fixed effect ANOVA on an unbalanced data set. Let the number of levels of the one way factor be "a".

Suppose the omnibus null hypothesis is rejected.

So I do a post hoc test, say Tukey Kramer HSD.

This will tell me which pairwise means are different. My query is the following: Can I form groups of levels that are NOT significantly different in pairwise tests, and conclude that they form ONE group?

The precise algo for this would be:

For i in 2 : a
For j in 1 : i-1
if mean of level i is not significantly different to the
mean of level j,then put i and j in the same cluster.
After the first time mean of level i is not different
to the mean of level j , just goto the next i , no
need to compare with remaining j's.


Can I use the above to cluster the levels in a data-driven way? If not, can someone tell me what test or package in R will be appropriate for clustering?

Alternately we run the following algo:

For i in 2 : a
For j in 1: n-1
compute the contrast that mean of level i is the same
as the mean of level j.If they are the same then put them in
the same cluster, skip all remaining comparisons in the
inner loop and goto the next iteration in the outer loop.


To control the family wise error rate, since we are doing AT MOST 1 + 2 + ... + (n-1 ) = (n-1)(n)/2 we can do a bonferroni or scheffe correction. What correction will be the best given we have AT MOST (n-1)(n)/2 comparisons?

• This strikes me as a sub-optimal idea since your final conclusions are going to be data-driven. It would be better to use theory or past studies to guide any re-categorisation. Commented Dec 31, 2017 at 15:11

Are you just trying to display the results as a compact letter display? Where, means sharing a letter are not significantly different.

If so, this can be accomplished with the multcompLetters function in the multcompView package.

But you may want to use the emmeans package for multiple comparisons.

Note that Tukey's test controls for the familywise error rate, so further corrections to the p-values are not needed.

if(!require(multcompView)){install.packages("multcompView")}
if(!require(emmeans)){install.packages("emmeans")}

fm1 = aov(breaks ~ wool + tension, data = warpbreaks)

TH = TukeyHSD(fm1, "tension", ordered = TRUE)

TH

P.value = TH\$tension[,4]

library(multcompView)

multcompLetters(P.value)

###   M   L   H
### "a" "b" "a"

library(emmeans)

marginal = emmeans(fm1, ~ tension)

cld(marginal, Letters=letters)

###  tension   emmean       SE df lower.CL upper.CL .group
###  H       21.66667 2.738184 50 16.16686 27.16647  a
###  M       26.38889 2.738184 50 20.88908 31.88869  a
###  L       36.38889 2.738184 50 30.88908 41.88869   b
###
### Results are averaged over the levels of: wool
### Confidence level used: 0.95
### P value adjustment: tukey method for comparing a family of 3 estimates
### significance level used: alpha = 0.05

• Dear Sal, Many thanks for your reply. I have one more query. Suppose in your example,one level is assigned to 2 groups,say H was in group a, M was in group ab and L was in group b. How would I then cluster the levels into 2 groups ? Where would I classify M, in group a or in group b? Commented Jan 1, 2018 at 5:10
• The meaning of the letters is "Means sharing a letter are not significantly different." So in this new case, H is not different from M, and M is not different L, but H is different from L. This can be difficult to summarize in words. Sometimes I will say something like "Treatment I [a] and II [a] and III [ab] were among the highest", but in any case including the mean separation letters are probably the best way to show the results. Commented Jan 1, 2018 at 14:49
• Dear Sal, the cld is necessary but not sufficient for my needs. It will tell me which 2 pairs of means are significantly different or not. It will not say anything about which sets of means belong together. How can I partition the a levels into groups where the means in each group are similar vs the means in different groups are significantly different ? In the special case where each mean belongs to only one group, the cld will suffice since we are able to see the grouping of means. If one mean belongs to 2 or more groups we wont be able to cleanly partition the means into groups. Commented Jan 2, 2018 at 6:13
• I think this approach will not do what you want for your data, then. It sounds like you are trying to force your results into boxes they don't fit in to. Of course, there are other approaches to cluster data that result in distinct categories. Commented Jan 2, 2018 at 11:58
• Okay thank you. Can you please name the other approaches to cluster data that result in distinct categories? Commented Jan 2, 2018 at 16:43