# Function plot_means(A, Y, xlabel, ylabel, bars) plots means and error bars of Y at each level of A, where xlabel = x-axis label in quotes, ylabel = y-axis label in quotes, and bars = "SD" or "SE" or "CI", for error bars showing the standard deviation, or standard error, or 95% confidence intervals
plot_means <- function(x, y, xlabel, ylabel, BAR)
{
  if (class(xlabel) != "character") stop("request x-axis label in quotes")
  if (class(ylabel) != "character") stop("request y-axis label in quotes")
  if (class(BAR) != "character") stop("for bars request only the names SD, SE or CI, in quotes")
  x <- factor(x)
  xm <- tapply(as.numeric(x), x, mean)
  ym <- tapply(y, x, mean)
  SD <- function(x, y) {tapply(y, x, var)^0.5}
  SE <- function(x, y) {(tapply(y, x, var)/tapply(y, x, length))^0.5}
  CI <- function(x, y) {qt(0.975, (tapply(y, x, length) - 1))*(tapply(y, x, var)/tapply(y, x, length))^0.5}
  bars <- do.call(BAR, list(x, y)) # calculate appropriate error bars
  ylimits <- c(ym-bars, ym+bars)
  ylabelc <- paste(ylabel, " means and ", BAR)
  plot(xm, ym, type = "p", lwd = 6, xaxt = "n", xlim = c(min(xm), max(xm)), ylim = c(min(ylimits), max(ylimits)), las = 1, xlab = xlabel, ylab = toString(ylabelc)) # plot means
  axis(1, at = xm, labels = names(xm)) # plot x-axis labels at each level of x
  for (i in 1:length(xm)) {
    arrows(xm[i], ylimits[i], xm[i], ylimits[length(xm)+i], length = .05, angle = 90, code = 3) # add bars
  }
}