(This article was first published on [R] tricks, and kindly contributed to R-bloggers)
Producing this kind of graphs (below) in R can be a pain in the a*s.

Here is a simple code that requires that data are presented in lists (see the example below).
multbar <- function(list.of.lists,...,condnames=0,pal=colorRampPalette(c('grey','cornsilk')),seriesnames=0,legendpos='topleft',legh=TRUE,do.pty='s')
{
par(pty=do.pty,mgp=c(1.9,0.8,0),oma=c(0,0,0,0),mar=c(4,3,2,1),bg='transparent',bty='o',tck=0.02,yaxs='i')
NofList <- length(list.of.lists)
NofSubList <- length(list.of.lists[[1]])
if(condnames==0){condnames=c(1:NofList)}
if(seriesnames==0){seriesnames=c(1:NofSubList)}
dim.mat.treat <- (NofList+1)*NofSubList
pos <- c(1:dim.mat.treat)
pos <- matrix(pos,nrow=(NofList+1))
nbreaks <- dim.mat.treat+1
par(xaxs='i',yaxs='i',bty='o')
plot(0,0,pch=NA,xlim=c(0,nbreaks),ylim=c(min(pretty(range(list.of.lists))),max(pretty(range(list.of.lists)))),xaxt='n',yaxt='n',...)
calc.ylab <- pretty(range(list.of.lists))
lab.ylab <- round(calc.ylab,2)
abline(h=pretty(range(list.of.lists)),lty=3,col='grey',lwd=2)
abline(h=0,col='darkgrey',lwd=1.4)
palette <- pal(NofList)
for(condition in 1:NofList)
{
for(treatment in 1:NofSubList)
{
CD <- list.of.lists[[condition]][[treatment]]
CDm <- mean(CD)
CDv <- sd(CD)
x.coord <- pos[condition,treatment]
rect(x.coord,0,x.coord+1,CDm,col=palette[condition])
arrows(x.coord+0.5,CDm+CDv,x.coord+0.5,CDm-CDv,angle=90,code=3,length=0.05)
}
}
axis(side=2,at=calc.ylab,labels=lab.ylab,tick=TRUE)
#abline(h=calc.ylab[1])
abline(h=pretty(range(list.of.lists))[1])
fills <- c(1:NofList)
legend(legendpos,legend=condnames,fill=palette[fills],bty='n',horiz=legh,cex=1)
cl.tab <- (c(1:NofSubList)-1)*NofSubList + NofList/2 +1
lab.cond <- NULL
int.lab.cond <- pos[c(1:(nrow(pos)-1)),]
for(COL in 1:ncol(int.lab.cond))
{
lab.cond[COL] <- mean(int.lab.cond[,COL])+0.5
}
axis(side=1,at=lab.cond,labels=seriesnames,tick=FALSE)
par(xaxs='r')
if(min(pretty(range(list.of.lists)))<0){abline(h=0,lwd=1.1)}
box()
}
list.1 <- list(rnorm(10,mean=3.4),rnorm(10,mean=6.8))
list.2 <- list(rnorm(10,mean=4),rnorm(10,mean=8.8))
list.3 <- list(rnorm(10,mean=3.1),rnorm(10,mean=7.4))
data <- list(list.1,list.2,list.3)
png()
par(family='serif')
multbar(data,main='Faked data',seriesnames=c('Control','Stress'),condnames=c('Ind. 1','Ind. 2','Ind. 3'),xlab='Treatment',ylab='Response')
dev.off()
To leave a comment for the author, please follow the link and comment on his blog: [R] tricks.
R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series,ecdf, trading) and more...

Zero Inflated Models and Generalized Linear Mixed Models with R.
Zuur, Saveliev, Ieno (2012).