Graphics with error bars

Next we will demonstrate how to use R’s powerful graphics functions to add error bars to a plot. The example uses the Elashoff example discussed earlier. In this example we will briefly show how visual representations compliment the statistical tests. We use R’s jpg() graphics driver to generate a graph that can be viewed by a web browser. The command syntax may appear intimidating for beginners, but it is worth the increased efficiency in the long run.
You can also use the postscript() graphics driver to generate presentation-quality plots. PostScript files can be transformed into PDF format so that nowadays the graphs generated by R can be viewed and printed by virtually anyone.¹¹
Typically the graphs are first generated interactively with drivers like X11(), then the commands are saved and edited into a script file. A command syntax script eliminates the need to save bulky graphic files.
First we start the graphics driver jpg() and name the file where the graph(s) will be saved.

attach(Ela.uni)
jpeg(file = "ElasBar.jpg")

Then we find the means, the standard errors, and the 95% confidence bounds of the means.

tmean <- tapply(effect, IND = list(gp, drug, dose), mean)
tse   <- tapply(effect, IND = list(gp, drug, dose), se)
tbarHeight <- matrix(tmean, ncol=3)
dimnames(tbarHeight) <- list(c("gp1dr1","gp2dr1","gp1dr2","gp2dr2"), 
          c("dose1","dose2" ,"dose3")) 
tn <- tapply(effect, IND = list(gp, drug, dose), length)
tu <- tmean + qt(.975, df=tn-1) * tse    # upper bound of 95% confidence int.
tl <- tmean + qt(.025, df=tn-1) * tse    # lower bound
tcol <- c("blue", "darkblue", "yellow", "orange")  # color of the bars

After all the numbers are computed, we start building the barplot. First we plot the bars without the confidence intervals, axes, labels, and tick marks. Note that the barplot() function returns the x-axis values at where the center of the bars are plotted. Later we will use the values in tbars to add additional pieces.

tbars <- barplot(height=tbarHeight, beside=T, space=c(0, 0.5), 
                 ylab="", xlab="", axes=F, names.arg=NULL, ylim=c(-15, 40), 
                 col=tcol)

Then we add the 95% confidence intervals of the means to the bars.

segments(x0=tbars, x1=tbars, y0=tl, y1=tu)
segments(x0=tbars-.1, x1=tbars+0.1, y0=tl, y1=tl)    # lower whiskers
segments(x0=tbars-.1, x1=tbars+0.1, y0=tu, y1=tu)    # upper whiskers

The axes labels are added.

axis(2, at = seq(0, 40, by=10), labels = rep("", 5), las=1)
tx <- apply(tbars, 2, mean)   # center positions for 3 clusters of bars

We plot the horizontal axis manually so that we can ask R to put things at exactly where we want them.

segments(x0=0, x1=max(tbars)+1.0, y0=0, y1=0, lty=1, lwd = 2)
text(c("Dose 1", "Dose 2", "Dose 3"), x = tx, y = -1.5, cex =1.5)
mtext(text=seq(0,40,by=10), side = 2, at = seq(0,40,by=10), 
      line = 1.5, cex =1.5, las=1)
mtext(text="Drug Effectiveness", side = 2, line = 2.5, at = 20, cex =1.8)

Finally we want to plot the legend of the graph manually. R also has a legend() function, although less flexible.

tx1 <- c(0, 1, 1, 0)
ty1 <- c(-15, -15, -13, -13)
polygon(x=tx1, y=ty1, col=tcol[1])
polygon(x=tx1, y=ty1 + 2.5, col=tcol[2])  # 2nd, moved 2.5 points up
polygon(x=tx1, y=ty1 + 5, col=tcol[3])    # 3rd
polygon(x=tx1, y=ty1 + 7.5, col=tcol[4])  # last

Finally, we add the legend labels for the filled rectangles.

text(x = 2.0, y = -14, labels="group 1, drug 1",  cex = 1.5, adj = 0)
text(x = 2.0, y = -11.5, labels="group 2, drug 1",  cex = 1.5, adj = 0)
text(x = 2.0, y = -9, labels="group 1, drug 2",  cex = 1.5, adj = 0)
text(x = 2.0, y = -6.5, labels="group 2, drug 2",  cex = 1.5, adj = 0)

The greatest effectiveness is attained by subjects in group 2 when drug 2 is given. This suggests a group by drug interaction, which is confirmed by the aov() results outlined earlier. It also indicates an increasing effectiveness from dose 1 to 3, which is also confirmed by the statistics.

7.8.8  Another way to do error bars using plotCI()

The gplots library has a function plotCI(), which does confidence intervals for plots. Here is an example, which is Figure 1 in http://journal.sjdm.org/06137/jdm06137.htm. Note the use of a small horizontal offset (−.01) so that the error bars do not overlap. The font “NimbusSan” is supposed to fit well with Times Roman.

library("gplots")
m.pg <- c(-2.6400, 3.6000, 6.0000, 3.6800, 5.4400)
se.pg <- c(1.71938, 1.86548, 1.74738, 1.94484, 1.83492)
m.pl <- c(-4.9600, -3.7600, -2.3200, -.1600, 6.5600)
se.pl <- c(1.47024, 1.72170, 1.79139, 1.36587, 1.56852)

postscript(file="fig1.eps",family="NimbusSan",
           width=8,height=8,horiz=F,pointsize=18,paper="special")
plotCI(y=c(m.pg,m.pl),x=c(c(1:5)-.01,c(1:5)+.01),uiw=c(se.pg,se.pl),
       ylim=c(-6,8),ylab="Net IGT score",xlab="Block",lty=rep(c(1,2),c(5,5)))
lines(y=m.pg,x=c(1:5)-.01,lty=1)
lines(y=m.pl,x=c(1:5)+.01,lty=2)
legend(3.6,-3.7,legend=c("Prior gain","Prior loss"),lty=c(1,2))
dev.off()