[R] Tukey HSD plot with lines indicating (non-)significance

[Karl Ove Hufthammer]

Dear list members,

I'm running some tests looking at differences between means for various 
levels of a factor, using Tukey's HSD method.

I would like to plot the data as boxplots or dotplots, with horizontal 
significance lines indicating which groups are statistically 
significantly different, according to Tukey HSD. Here's a nice image 
showing an example of such a graphical display:

   http://www.biomedcentral.com/1471-2148/11/99/figure/F2

Is there a R function for this? I have tried searching, and found 
several nice plotting fuctions for Tukey's HSD, but they all either hide 
the actual data, or display the results only with letters indicating 
(non-)signifiance.

Examples:

# Fit an one-way ANOVA.
library(DAAG)
l=aov(ShootDryMass ~ trt, data=rice)
summary(l)


# Calculate the Tukey HSD tests and CIs.
# The plot shows all 6*5/2 = 15 confidence intervals,
# but not the data or the actual group means.
l.hsd=TukeyHSD(l, "trt", ordered = TRUE)
l.hsd
plot(l.hsd)


# I rather like this one. It shows the group means in a nice way,
# but hides the data. And checking which groups are different
# isn't always easy (e.g., for NH4Cl - NH4NO3), as one
# has to compare differences to a non-aligned scale (bar).
onewayPlot(l)


# One function that gives a results closer to what I want,
# is the compact letter display (CLD) of the multcomp package.
# This shows boxplots (but it is easy to overlay the actual
# observations, using the points() function) and uses letters
# to indicate non-significance.
#
# (Note that (if I have understood everything correctly) this
# doesn't actually calculcate Tukey HSD values (the "Tukey"
# in the glht() call only selects Tukey *contrasts*), but
# the results should be *very* similar.)
library(multcomp)
l.glht=glht(l, linfct = mcp(trt = "Tukey"))
summary(l.glht)
l.cld=cld(l.glht)
old.par <- par( mai=c(1,1,2,1))
plot(l.cld)
par(old.par)


# Basically, I want a similar display as the one above (or
# preferably a ggplot2- or lattice-based one), but with lines
# instead of letters. Of course, for this, the levels need to be
# reordered by group means, and it is only guaranteed to
# work for balanced data. Example, with the lines missing:
library(ggplot2)
d=rice
d$trt=reorder(d$trt, -d$ShootDryMass)
ggplot(d, aes(x=trt, y=ShootDryMass)) + geom_boxplot(outlier.colour=NA) 
+
   geom_jitter(col="red", size=3, position=position_jitter(width=.1))

Is there such a function available in R?

-- 
Karl Ove Hufthammer