A boxplot replacement that I call the "Jelly Bean Plot".

jelly.beans.R

jelly.beans <- function(data, brewerset = "Set1", cols = brewer.pal(length(data),brewerset),
    diffs = rep(NA,length(data)), notches = FALSE,
    xlab = NA, ylab = NA, las = rep(0,2), bean.width = c(3,8,3,1,2), ...) {
  require(RColorBrewer)
  # data = list(grp1=vector(),grp2=vector(),grp3=vector())
  # bean.width = c(0.95-lines,0.5-lines,median-line,grid-line)
  # diffs = c("A", "A", "B") eg to present TukeyHSD values
  do.call(rbind,lapply(1:length(data),function(i) c(i,data[[i]])))
  num.grps <- length(data)
  grp.names <- names(data)
  mu <- lapply(data,mean,na.rm=TRUE)
  hpd <- lapply(data,quantile,probs=c(0.025,0.25,0.5,0.75,0.975),na.rm=TRUE)
  yrange <- range(hpd)
  plot(1,type="n",xlim=c(0,num.grps+1),ylim=yrange,axes=FALSE,xlab=NA,ylab=NA,...)
  if (notches) {
    for (i in 1:num.grps) {
      iqr <- hpd[[i]][3]+c(-1,1)*(hpd[[i]][4]-hpd[[i]][2])*1.58/sqrt(length(data[[i]]))
      rect(0.5,iqr[1],num.grps+0.5,iqr[2],col=paste(cols[i],"33",sep=""),border=NA)
      lines(c(0.5,num.grps+0.5),rep(iqr[1],2),col=cols[i],lty=2)
      lines(c(0.5,num.grps+0.5),rep(iqr[2],2),col=cols[i],lty=2)
     }
  }
  for (i in 1:num.grps) {
    lines(c(i,i),yrange,col=rgb(.5,.5,.5),lwd=bean.width[4])
    lines(c(i,i),hpd[[i]][c(1,5)],col=cols[i],lwd=bean.width[1])
    lines(c(i,i),hpd[[i]][c(2,4)],col=cols[i],lwd=bean.width[2])
    if (notches) {
      lines(i+c(-.2,.2),rep(hpd[[i]][3],2),col=cols[i],lwd=bean.width[3]*1.1)
    }
    lines(i+c(-.2,.2),rep(hpd[[i]][3],2),col="white",lwd=bean.width[3])
  }
  axis.at <- 1:num.grps
  axis(1,at=axis.at,labels=grp.names,las=1,lwd=0,pos=yrange[1],cex.axis=.8)
  axis(3,at=axis.at,labels=diffs,las=1,lwd=0,pos=yrange[2],cex.axis=.8)
  x.tics <- axTicks(2)
  lines(c(0,0),c(min(yrange,x.tics),max(yrange,x.tics)))
  axis(2,at=x.tics,pos=0,lwd=0,lwd.ticks=1,las=1,tcl=-0.3)
  if (! is.na(xlab)) mtext(xlab,side=1,line=2,las=las[1])
  if (! is.na(ylab)) mtext(ylab,side=2,line=2,las=las[2])
}

library(multcompView)
data <- lapply(1:4,rnorm,n=20,sd=1)
names(data) <- c("Monkey","Possum","Armadillo","Tardigrade")
data.df <- data.frame(do.call(rbind,lapply(1:4,function(i)
                      cbind(Animal=i,Awesomeness=data[[i]]))))
data.df$Animal <- factor(data.df$Animal)
mo <- aov(Awesomeness~Animal,data=data.df)
mo.hsd <- TukeyHSD(mo)
diffs <- rev(multcompLetters2(Awesomeness~Animal,mo.hsd$Animal[,"p adj"],
             data.df,reversed=TRUE)$Letters)

pdf(file="jellyBeans.pdf",width=12,height=6)
par(mfrow=c(1,2))
jelly.beans(data,ylab="Awesomeness",xlab="Animal",notches=FALSE,
    las=c(1,0),bean.width=c(10,20,10,1,5),diffs=diffs)
jelly.beans(data,ylab="Awesomeness",xlab="Animal",notches=TRUE,
    las=c(1,0),bean.width=c(10,20,10,1,5),diffs=diffs)
dev.off()

which version of R did you develop this in? I am having trouble with the multcompView package.

> install.packages(multcompView)
Error in install.packages : object 'multcompView' not found