см. Например, ?errbar
в упаковке Hmisc
group <- factor(sample(1:10,100,T))
y <- (1:10)[group] + rnorm(100)
grmean <- tapply(y,group,mean)
lims <- tapply(y,list(group),FUN = sd)*2
errbar( 1:10, grmean, grmean + lims, grmean - lims , xaxt="n" )
axis(3)
дает: