R Barplot with multiple grouping levels

Question

I'm stuck with trying to create a somewhat complex grouped barplot.

Here is a snippet of my data file

"location"|"region"|"treatment"|"GB"
"Georgia"|"Keys"|"pre"|354
"Georgia"|"Keys"|"pre"|183
"Georgia"|"Keys"|"pre"|182
"Georgia"|"Keys"|"pre"|133
"Georgia"|"North East"|"pre"|44
"Georgia"|"North East"|"pre"|19
"Georgia"|"North East"|"pre"|70
"Georgia"|"North East"|"pre"|66
"Georgia"|"North West"|"pre"|102
"Georgia"|"North West"|"pre"|33
"Georgia"|"North West"|"pre"|106
"Georgia"|"North West"|"pre"|61
"Georgia"|"North West"|"pre"|101
"Georgia"|"Texas"|"pre"|150
"Georgia"|"Texas"|"pre"|187
"Georgia"|"Texas"|"pre"|152
"Georgia"|"Texas"|"pre"|148
"Georgia"|"Texas"|"pre"|100
"Maryland"|"Keys"|"pre"|637
"Maryland"|"Keys"|"pre"|52
"Maryland"|"Keys"|"pre"|43
"Maryland"|"Keys"|"pre"|156
"Maryland"|"Keys"|"pre"|38
"Maryland"|"North East"|"pre"|166
"Maryland"|"North East"|"pre"|91
"Maryland"|"North East"|"pre"|167
"Maryland"|"North East"|"pre"|104
"Maryland"|"North East"|"pre"|113
"Maryland"|"North West"|"pre"|370
"Maryland"|"North West"|"pre"|895
"Maryland"|"North West"|"pre"|198
"Maryland"|"North West"|"pre"|137
"Maryland"|"North West"|"pre"|168
"Maryland"|"Texas"|"pre"|95
"Maryland"|"Texas"|"pre"|331
"Maryland"|"Texas"|"pre"|163
"Maryland"|"Texas"|"pre"|90
"North Carolina"|"Keys"|"pre"|217
"North Carolina"|"Keys"|"pre"|91
"North Carolina"|"Keys"|"pre"|148
"North Carolina"|"Keys"|"pre"|208
"North Carolina"|"Keys"|"pre"|18
"North Carolina"|"North East"|"pre"|49
"North Carolina"|"North East"|"pre"|60
"North Carolina"|"North East"|"pre"|167
"North Carolina"|"North East"|"pre"|82
"North Carolina"|"North East"|"pre"|31
"North Carolina"|"North West"|"pre"|47
"North Carolina"|"North West"|"pre"|10
"North Carolina"|"North West"|"pre"|207
"North Carolina"|"North West"|"pre"|70
"North Carolina"|"North West"|"pre"|214
"North Carolina"|"Texas"|"pre"|183
"North Carolina"|"Texas"|"pre"|162
"North Carolina"|"Texas"|"pre"|94
"North Carolina"|"Texas"|"pre"|102
"South Carolina"|"Keys"|"pre"|101
"South Carolina"|"Keys"|"pre"|155
"South Carolina"|"Keys"|"pre"|85
"South Carolina"|"Keys"|"pre"|67
"South Carolina"|"Keys"|"pre"|60
"South Carolina"|"North East"|"pre"|173
"South Carolina"|"North East"|"pre"|148
"South Carolina"|"North East"|"pre"|575
"South Carolina"|"North East"|"pre"|96
"South Carolina"|"North West"|"pre"|51
"South Carolina"|"North West"|"pre"|86
"South Carolina"|"North West"|"pre"|34
"South Carolina"|"North West"|"pre"|67
"South Carolina"|"Texas"|"pre"|124
"South Carolina"|"Texas"|"pre"|155
"South Carolina"|"Texas"|"pre"|183
"South Carolina"|"Texas"|"pre"|101
"Georgia"|"Keys"|"post"|344
"Georgia"|"Keys"|"post"|241
"Georgia"|"Keys"|"post"|486
"Georgia"|"Keys"|"post"|191
"Georgia"|"North East"|"post"|128
"Georgia"|"North East"|"post"|14
"Georgia"|"North East"|"post"|192
"Georgia"|"North East"|"post"|298
"Georgia"|"North West"|"post"|540
"Georgia"|"North West"|"post"|236
"Georgia"|"North West"|"post"|172
"Georgia"|"North West"|"post"|87
"Georgia"|"Texas"|"post"|357
"Georgia"|"Texas"|"post"|221
"Georgia"|"Texas"|"post"|131
"Georgia"|"Texas"|"post"|55
"Maryland"|"Keys"|"post"|196
"Maryland"|"Keys"|"post"|85
"Maryland"|"Keys"|"post"|90
"Maryland"|"Keys"|"post"|530
"Maryland"|"North East"|"post"|477
"Maryland"|"North East"|"post"|447.253
"Maryland"|"North East"|"post"|509
"Maryland"|"North East"|"post"|64
"Maryland"|"North West"|"post"|1204
"Maryland"|"North West"|"post"|756
"Maryland"|"North West"|"post"|635
"Maryland"|"North West"|"post"|948
"Maryland"|"Texas"|"post"|740
"Maryland"|"Texas"|"post"|567
"Maryland"|"Texas"|"post"|549
"Maryland"|"Texas"|"post"|271
"North Carolina"|"Keys"|"post"|173
"North Carolina"|"Keys"|"post"|114
"North Carolina"|"Keys"|"post"|1159
"North Carolina"|"Keys"|"post"|113
"North Carolina"|"North East"|"post"|176
"North Carolina"|"North East"|"post"|187
"North Carolina"|"North East"|"post"|279
"North Carolina"|"North East"|"post"|182
"North Carolina"|"North West"|"post"|103
"North Carolina"|"North West"|"post"|230
"North Carolina"|"North West"|"post"|117
"North Carolina"|"North West"|"post"|143
"North Carolina"|"Texas"|"post"|358
"North Carolina"|"Texas"|"post"|458
"North Carolina"|"Texas"|"post"|102

The levels I have are 'treatment' (pre and post), 'region' [where samples were collected](Keys, northwest, north east and texas) and 'location' [where the experiments were performed](Maryland, Georgia, North Carolina and South Carolina). I was measuring 'GB' for each sample.

I want to get a plot that shows pre and post GB at each sample region within each experimental location

I can make a plot showing mean GB for pre and post treatment at region with no issue.

When I add location to the script as follows

data <- tapply(dat$GB, list(treat,regi,loci), mean)

I can get R to calculate the mean GB for pre and post for each region at each location. Happy days!!

But then when I try to plot this in to a barplot I get an error in R saying 'height' must be a vector or a matrix'

Here is the script I have written. You will see I ordered the data for each level and then used this new 'ordered' data to create means and SE's for each pre and post sample within each region at each location.

I then tried to use this 'mean' data for the plot.

thank you

Matt

#treatments pre and post between regions across locations
loci = factor (dat$location, levels=c("Georgia","Maryland","North Carolina","South Carolina"))
regi = factor(dat$region, levels=c("Keys","North West","North East", "Texas"))
treat = factor(dat$treatment, levels=c("Pre","Post"),ordered=TRUE)

data <- tapply(dat$GB, list(treat,regi,loci), mean)
ses<- tapply(dat$GB, list(treat,regi,loci),function (x) sd(x)/sqrt(length(x)))
lower<-data-ses
upper<-data+ses

my.plot<-barplot(data,beside=TRUE, legend= F, main="",ylim=c(0,400),xlab="",ylab="",cex=1.3, 
                 cex.lab=1.3, cex.axis=1.3)
arrows(my.plot, data, my.plot, lower, angle=90,length=.1)
arrows(my.plot, data, my.plot, upper, angle=90, length=.1)
mtext(expression(ug~GB~(ml^-1)), side=2, line=3, cex=1.5)

DATA.

dat1 <-
structure(list(location = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L), .Label = c("Georgia", "Maryland", "North Carolina", 
"South Carolina"), class = "factor"), region = structure(c(1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 
4L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 
3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("Keys", "North East", 
"North West", "Texas"), class = "factor"), treatment = structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("post", "pre"), class = "factor"), 
    GB = c(354, 183, 182, 133, 44, 19, 70, 66, 102, 33, 106, 
    61, 101, 150, 187, 152, 148, 100, 637, 52, 43, 156, 38, 166, 
    91, 167, 104, 113, 370, 895, 198, 137, 168, 95, 331, 163, 
    90, 217, 91, 148, 208, 18, 49, 60, 167, 82, 31, 47, 10, 207, 
    70, 214, 183, 162, 94, 102, 101, 155, 85, 67, 60, 173, 148, 
    575, 96, 51, 86, 34, 67, 124, 155, 183, 101, 344, 241, 486, 
    191, 128, 14, 192, 298, 540, 236, 172, 87, 357, 221, 131, 
    55, 196, 85, 90, 530, 477, 447.253, 509, 64, 1204, 756, 635, 
    948, 740, 567, 549, 271, 173, 114, 1159, 113, 176, 187, 279, 
    182, 103, 230, 117, 143, 358, 458, 102)), .Names = c("location", 
"region", "treatment", "GB"), class = "data.frame", row.names = c(NA, 
-120L))

Answer 1

In order to find a solution to your problem, let's have a look at your code at first.

When I add location to the script as follows
data <- tapply(dat$GB, list(treat,regi,loci), mean)
I can get R to calculate the mean GB for pre and post for each region at each location.

How do you do that? My R gives me back an empty three dimensional array with 2 rows, 4 columns and 4 elements in the third dimension. If I run

all(is.na(data))

I get the result TRUE. So no happy days for me here. However, the help page of tapply (?tapply) helped me getting further: The list representing the indexes should contain one or more factors with the length of X, but lapply will convert suitable vectors to factors itself. So if I run

data <- tapply(dat$GB, list(dat$treatment,dat$region,dat$location), mean)

the above mentioned array is filled with the mean values. Same is true for the confidence interval borders.

But now starts the real trouble. barplot() still throws the error "'height' must be a vector or a matrix". Means it won't accpet a three dimensional array as input for the heights of the bars, i.e. as the data. I guess this is, because there is no 'standard' for displaying 3-D data, i.e. 4 informations per element in a barplot (it's value, it's treatment, it's sampling region and it's sampling subregion). Even with a proper 3-D barplot, as shown in this post, you still need to condense two informations on one axis (here only done rudimentarily by combining "location" and "region", see later plots for detailed solution and how exactly region and location have been compined to "ExactReg"):

library(latticeExtra)
cloud(GB ~ exactReg + treatment, data = dat, panel.3d.cloud=panel.3dbars, xbase=0.5, ybase=0.2)

produces: (It would be possible to add standard deviation or even confidence interval later, as can be seen here (I'm quite sure this is part of a SO post on R as well, but I can't find it).)

So we have to think of an alternative solution. Here it depends what aspect of your data you want to emphasize. The best information about the distribution of the data, and therefore very likely about the process forming this data, is imho provided by Hmisc::bpplot. But it focuses on robust statistics and therefore does not show mean values and standard deviations (but medians and quartiles). However, these can be added later, but adding both of them will end up in quite a crowdy plot, that is overloaded in my opinion. I therefore added only means, showing the concept. bpplots do also show outliers less clearly than box and whisker plots (hereafter known as boxplots). And another malus is that they don't have "at" and "add" arguments so far and are therefore very hard to be grouped nicely in one line, "pre" and "post" treatment measurements next to each other, which is best for comparing "pre" and "post". I would depict theem like this for now:

library(TeachingDemos)
set.seed(7)
ms.sun <- rbind(ms.polygon(n = 50, r = 1),
                c(NA, NA),
                ms.polygon(n = 50, r = 0.03),
                c(0, 0),
                c(0, 0.0075),
                c(0.0075, 0.0075),
                c(0.0075, 0),
                c(0.0075, -0.0075),
                c(0, -0.0075),
                c(-0.0075, -0.0075),
                c(-0.0075, 0),
                c(-0.0075, 0.0075),
                matrix(data = runif(80, min = -0.029, max = 0.029), ncol = 2),
                c(NA, NA))


library(Hmisc)
library(pBrackets)
x11(10,8)
par(mfrow = c(2,1), mar = c(0,4,2,1), oma = c(8,0,0,0))
      #2 plots in the same columns (one above the other) -- no margin at bottom -- more space for x-axis below lower plot
xpos <- bpplot(split(dat[which(dat$treatment == "pre"), "GB"], dat[which(dat$treatment == "pre"), "exactReg"]), xlab = "", name = FALSE, ylab = "pre treatment", main = "GB", plotopts = list("ylim" = c(0, 1300)))
      #we need the positions of the barplots on the x-axis later, that's what is saved in "xpos"
      #I want to give pre treatment and post treatment the same value range on the y-axis, for comparability. Therefore the "ylim" argument in plotopts (in both calls to bpplot)
  abline(v = c(seq(from = 0, by = 1.2, length.out = length(levels(dat$exactReg)))), lty = "dashed", col = gray(0.7))
      #add vertical lines in both plots to make assignment easy
  my.symbols(xpos, by(dat[which(dat$treatment == "pre"), "GB"], dat[which(dat$treatment == "pre"), "exactReg"], mean),
             symb = ms.sun, xsize = 0.3)
      #add mean values
  legend("topright", pch = c(NA), legend = c("mean value"), bg = "white")
  my.symbols(16.85487, 1238.545, symb = ms.sun, xsize = 0.3)
      #add a legend to clarify what this extra circle means

par(mar = c(1,4,0,1))
      #no margin on top of plot, same margins left and right!
bpplot(split(dat[which(dat$treatment == "post"), "GB"], dat[which(dat$treatment == "post"), "exactReg"]), main = "", ylab = "post treatment", plotopts = list("ylim" = c(0, 1300)))
  abline(v = c(seq(from = 0, by = 1.2, length.out = length(levels(dat$exactReg)))), lty = "dashed", col = gray(0.7))
  my.symbols(xpos, by(dat[which(dat$treatment == "post"), "GB"], dat[which(dat$treatment == "post"), "exactReg"], mean),
             symb = ms.sun, xsize = 0.3)
  #The Error occurs, because not every group has got values to calculate a mean from
  #I think, in this case there is no extra legend needed
  mtext("|", side = 1, line = -0.15, at = xpos)   #axis ticks
  text(x = xpos, y = par("usr")[3] - 1/2.5* diff(c(par("usr")[3], par("usr")[4])),
       labels = rep(levels(dat$region), length(levels(dat$location))), srt = 90, adj = c(0, NA), xpd = NA)
  brackets(xpos[4] + 0.5*diff(xpos[4:5]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           xpos[1] - 0.5*diff(xpos[1:2]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  brackets(xpos[8] + 0.5*diff(xpos[8:9]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           xpos[5] - 0.5*diff(xpos[4:5]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  brackets(xpos[12] + 0.5*diff(xpos[12:13]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           xpos[9] - 0.5*diff(xpos[8:9]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  brackets(xpos[16] + 0.5*diff(xpos[15:16]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           xpos[13] - 0.5*diff(xpos[12:13]), par("usr")[3] - 0.45*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  #mtext(dat$location[1], side = 1, line = 5, at = mean(xpos[2:3]), outer = T)   #I can't tell for now, why this doesn't work
  text(x = seq(from = mean(xpos[2:3]), by = length(levels(dat$region))*diff(xpos[1:2]), length.out = length(levels(dat$location))),
       y = par("usr")[3] - (1/1.8)* diff(c(par("usr")[3], par("usr")[4])), labels = levels(dat$location), xpd = NA)
  rm(xpos)

Giving:

NB: My solution works only like this, if you got at least one value for every region in every location of the "post" treated group as well. Right now the boxes in the "post treatment" plot don't coincide with the location stated. There might be a lot of fine work to be done still (like making the dot in the sun symbol for the means bigger, adding a bit of space between the curly brackets...), but the concept should be clear...

Boxplots on the other hand give less information about the distribution of the data, but show outliers more clearly and can be grouped ideally for comparing pre and post treatment data. They are also instruments of robust statistics and don't display means and standard deviations, but medians and quartiles. Again I added only means:

x11(10,5)
par(mfrow = c(1,1), mar = c(8,4,3,1))
boxplot(dat[which(dat$treatment == "pre"), "GB"] ~ dat[which(dat$treatment == "pre"), "exactReg"],
        at = seq(from = 1, by = 3, length.out = length(levels(dat$exactReg))), col = "brown1",
        xlim = c(1, 3*length(levels(dat$exactReg))-1), ylim = c(0, max(dat$GB)),
        names = character(length(levels(dat$exactReg))), xaxt = "n", ylab = "GB", main = "Header")
  my.symbols(seq(from = 1, by = 3, length.out = length(levels(dat$exactReg))),
             by(dat[which(dat$treatment == "pre"), "GB"], dat[which(dat$treatment == "pre"), "exactReg"], mean),
             symb = ms.sun, xsize = 0.7)
  boxplot(dat[which(dat$treatment == "post"), "GB"] ~ dat[which(dat$treatment == "post"), "exactReg"],
        at = seq(from = 2, by = 3, length.out = length(levels(dat$exactReg))), col = "lightblue",
        names = character(length(levels(dat$exactReg))), xaxt = "n", ylab = "", add = TRUE)
  my.symbols(seq(from = 2, by = 3, length.out = length(levels(dat$exactReg))),
             by(dat[which(dat$treatment == "post"), "GB"], dat[which(dat$treatment == "post"), "exactReg"], mean),
             symb = ms.sun, xsize = 0.7)
        #The Error occurs, because not every group has got values to calculate a mean from
  legend("topright", fill = c("brown1", "lightblue", NA), border = c("black", "black", NA), legend = c("prior to treatment", "after treatment", "mean value"))
  my.symbols(40.5, 964.8124, symb = ms.sun, xsize = 0.7)
  mtext("|", side = 1, line = -0.15, at = seq(from = 1.5, by = 3, length.out = length(levels(dat$exactReg))))   #axis ticks
  text(x = seq(from = 1.5, by = 3, length.out = length(levels(dat$exactReg))),
       y = par("usr")[3] - 1/2.5* diff(c(par("usr")[3], par("usr")[4])),
       labels = rep(levels(dat$region), length(levels(dat$location))), srt = 90, adj = c(0, NA), xpd = TRUE)
  brackets(11.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           0.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  brackets(23.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           12.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  brackets(35.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           24.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  brackets(47.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           36.5, par("usr")[3] - 0.41*diff(c(par("usr")[3], par("usr")[4])),
           h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
  mtext(levels(dat$location), side = 1, line = 6.5, at = seq(6, by = 12, length.out = length(levels(dat$location))) , xpd = TRUE)

This solution has the same problem with the missing values than the other: The boxes are not drawn at the position where they should be, but one after the other, i.e. as soon as one box is missing completely the following ones are one slot too far left and so on. Unfurtanately I can't do this with bpplots.

Finally, if you only want to show that your treatment did make a difference (not how), and if outliers aren't too big of an issue in your data, barplots indeed fit best for this task, since you can use classical statistics best.
I used your approach with tapply to make calculations with the data. It showed me clearly, that indeed data for the last four groups is missing, so the order is not messed up in the other plots. You can achieve a barplot similar to the boxplots:

data <- tapply(dat$GB, list(dat$exactReg, dat$treatment), mean)
ses <- tapply(dat$GB, list(dat$exactReg, dat$treatment), function(x){sd(x)/sqrt(length(x))})
lower <- data-ses
upper <- data+ses

x11(10,8)
par(mar = c(12,5,3,1), las = 3)
my.plot <- barplot(data[, "pre"], space = c(0, rep(2, 15)), main = "Header", ylim = c(0, 1100), xlab = "", ylab = "", names.arg = "", col = "lightblue")
arrows(my.plot, data[, "pre"], my.plot, lower[, "pre"], angle = 90, length = 0.08)
arrows(my.plot, data[, "pre"], my.plot, upper[, "pre"], angle = 90, length = 0.08)
my.plot <- barplot(data[, "post"], space = c(1, rep(2, 15)), main = "", ylim = c(0, 1100), xlab = "", ylab = "", names.arg = "", col = "brown1", add = TRUE)
arrows(my.plot, data[, "post"], my.plot, lower[, "post"], angle = 90, length = 0.08)
arrows(my.plot, data[, "post"], my.plot, upper[, "post"], angle = 90, length = 0.08)
mtext(expression(ug~GB~(ml^-1)), side=2, line=3, cex=1.5)
legend("topright", fill = c("lightblue", "brown1"), legend = c("prior to treatment", "after treatment"))
axis(1, at = seq(from = 1.5, by = 3, length.out = length(levels(dat$exactReg))), labels = FALSE, lwd.ticks = 2)
text(x = seq(from = 1.5, by = 3, length.out = length(levels(dat$exactReg))),
     y = par("usr")[3] - 1/5.1* diff(c(par("usr")[3], par("usr")[4])),
     labels = rep(levels(dat$region), length(levels(dat$location))), srt = 90, adj = c(0, NA), xpd = TRUE)
brackets(11.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         0.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
brackets(23.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         12.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
brackets(35.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         24.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
brackets(47.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         36.5, par("usr")[3] - 0.2*diff(c(par("usr")[3], par("usr")[4])),
         h = 0.04*diff(c(par("usr")[3], par("usr")[4])), xpd = NA)
mtext(levels(dat$location), side = 1, line = 6.2, at = seq(6, by = 12, length.out = length(levels(dat$location))) , xpd = TRUE)