Post-hoc test on mixed model - function for groups display

Question

I want to run post-hoc tests for a linear mixed model. I’ve initally used lsmeans∷lsmeans, which works fine. Problem is: with lsmeans, the result is not “sorted” in groups (say: a ab b c etc.), at least not that I know of. So I’ve sorted the groups manually. As I have many treatments in my experimental design, and many post-hoc tests to run, this is somehow painful (and error-prone).

I’ve looked at agricolae∷HSD.test, which does give groups for post-hoc tests. Downside is that HSD.test does not run with mixed models (neither do multcomp::cld, multcomp::hsd, nor multcomp::multcompLetters for that matter).

EDIT: ANSWER

# here, the groups are sorted so that the highest mean is "a"
library(lsmeans)
comp_2_by_2<- lsmeans(model, pairwise ~ FactorA:FactorB)
groups_HSD<- cld(comp_2_by_2, alpha=0.05, Letters=c("a","b","c","d"), reverse=TRUE)