Intro

This is a simple function that plots means as individual points (not bars) and approximate 95% confidence intervals using user-supplied means and standard errors or 95% CIs. It provides some feed back via the console and the plot itself regarding errors and how to make the plot more polished.

Base R has no default function to make plots like this, though there are numerous packages that have function that can do this (sciplot, ggplot, Hmisc, gplot, psych). I wrote this function as bare-bones alternative that my students could load via a script and which has only a few, simple arguements. I teach them to calculate means and standard deviations using the tapply() function.

I’ve found that loading packages can be difficult on computers behind campus firewalls, as can updating R if a package requires an up-to-date installation. At the end of this document I provide notes about and links to other ways to approach this problem with R published packages.

plot.means()

# means = your means. , contained in a vector, eg c(mean1, means2...) 
# SEs =  your standard errors, contained in a vector, 
# CI.hi = upper confidence intervals, in a vector, 
# CI.lo = lower confidence itnerval, in a vector
# categories = the names of the categories/groups, in the order that they appear
# x.axis.label = what should be plotted on the x axis
# y.axis.label = what should be plotted on the y axis
# adjust.y.max = allows you to adjust the y axis
# adjust.y.min
# adjust.x.spacing
#### The function STARTS here ####
plot.means <- function(means = NULL,
                           SEs = NULL,
                           CI.hi = NULL,
                           CI.lo = NULL,
                           categories = NULL,
                           x.axis.label = "Groups",
                           y.axis.label = "'y.axis.label' sets the axis label",
                       adjust.y.max = 0,
                       adjust.y.min = 0,
                       adjust.x.spacing = 5){
  
  
    # Error messages
    if(is.null(means)==T){
    stop("No means entered") } 
  
    if(is.null(SEs)==T & is.null(CI.lo) == T){
    stop("No standard errors entered") } 
  
    if(is.null(CI.lo) != T &  is.null(CI.hi)){
    stop("CI.lo entered but no CI.hi") } 
  
    if(is.null(CI.hi) != T &  is.null(CI.lo)){
    stop("CI.hi entered but no CI.lo") } 
  
  #check if both SE and and CI.is enter
  if(is.null(SEs) == F &  is.null(CI.hi) == F){
    stop("Both SEs and CI.hi entered; use eithe SEs or CIs") } 
  
    if(is.null(SEs) == F &  is.null(CI.lo) == F){
    stop("Both SEs and CI.lo entered; use eithe SEs or CIs") } 
  
  #Check if the number of means matches the number of SE
  n.means <- length(means)
  n.SEs    <- length(SEs)
  n.CI.lo  <- length(CI.lo)
  n.CI.hi  <- length(CI.hi)
  
  #CHeck standard errors against means
  if(n.means != n.SEs & is.null(CI.lo)==T){
  error.message <- paste("The number of means is",n.means,"but the number of standard errors is",n.SEs,sep = " ")
    stop(error.message) } 
  
  #Check CIs against CIs
    if(n.CI.lo != n.CI.hi){
  error.message <- paste("The number of CI.lo is",n.CI.lo,"but the number of standard errors is",n.CI.hi,sep = " ")
    stop(error.message) } 
  
  
  #assign arbitrary categories
  if(is.null(categories) == T) {
    categories <- paste("Group",1:n.means)
    print("Set categoris labls with 'categories=' ")
       }
  
  
  # calculate values for plotting limits 
  if(is.null(SEs)==F){
  y.max <- max(means+2*SEs) + adjust.y.max
  y.min <- min(means-2*SEs) - adjust.y.min
  }

  # calculate values for plotting limits 
  if(is.null(SEs)==T){
  y.max <- max(CI.hi) + adjust.y.max
  y.min <- min(CI.lo) - adjust.y.min
  }

  
  #determine where to plot points along x-axis
  x.values <- 1:n.means
  x.values <- x.values/adjust.x.spacing
  
  #set x axis min/max
  x.axis.min <- min(x.values)-0.05
  x.axis.max <- max(x.values)+0.05
  
  x.limits <- c(x.axis.min,x.axis.max)
  
  #Plot means
  plot(means ~ x.values,
       xlim = x.limits,
       ylim = c(y.min,y.max),
       xaxt = "n",
       xlab = "",
       ylab = "",
       cex = 1.25,
       pch = 16)
  
  #Add x labels
  axis(side = 1, 
       at = x.values,
       labels = categories
      )
  
  
  # Plot confidence intervals
  if(is.null(CI.hi) == FALSE & 
     is.null(CI.lo) == FALSE){
  #Plot upper error bar for CIs
  lwd. <- 2
  arrows(y0 = means,
         x0 = x.values,
         y1 = CI.hi,
         x1 = x.values,
         length = 0,
         lwd = lwd.)
  
  #Plot lower error bar
  arrows(y0 = means,
         x0 = x.values,
         y1 = CI.lo,
         x1 = x.values,
         length = 0,
         lwd = lwd.)
  }
    
    
  # Estimate CIs from SEs
  
  if(is.null(SEs) == FALSE & 
     is.null(CI.lo) == TRUE){
  lwd. <- 2
  arrows(y0 = means,
         x0 = x.values,
         y1 = means+2*SEs,
         x1 = x.values,
         length = 0,
         lwd = lwd.)
  
  #Plot lower error bar
  arrows(y0 = means,
         x0 = x.values,
         y1 = means-2*SEs,
         x1 = x.values,
         length = 0,
         lwd = lwd.) 
  }
  
  mtext(text = x.axis.label,side = 1,line = 2)
  mtext(text = y.axis.label,side = 2,line = 2)
  mtext(text = "Error bars = 95% CI",side = 3,line = 0,adj = 0)
  

  
}

#### The function ENDS here ####
#### The function ENDS here ####
#### The function ENDS here ####
#### The function ENDS here ####
#### The function ENDS here ####

Test the function with 3 groups

The following code is used to test the function. DO NOT LOAD THIS CODE for general use

Calcualte mean and SE from Fisher’s iris data

# THIS IS NOT PART OF THE FUNCTION

data(iris)

iris.3mean <- tapply(iris$Sepal.Length,
                     iris$Species,
                     mean)

iris.3sd <- tapply(iris$Sepal.Length,
                     iris$Species,
                     sd)

iris.3n <- tapply(iris$Sepal.Length,
                     iris$Species,
                     length)

iris.3se <- iris.3sd/sqrt(iris.3n)

Plot Fisher’s iris data with plot.means()

par(mfrow = c(1,1))
plot.means(means = iris.3mean,
          SEs = iris.3se*3)

## [1] "Set categoris labls with 'categories=' "

Test confidence interval function

Make fake confidence intervals

CI.hi <-iris.3mean+10
CI.lo <-iris.3mean-10

plot.means(means = iris.3mean,
          CI.lo = CI.lo,
          CI.hi = CI.hi)

## [1] "Set categoris labls with 'categories=' "

Compare to the output of the sciplot package

This requires the sciplot package to be downloded

library(sciplot)
ci.fun.use <- function(x) c(mean(x)-1.96*se(x), mean(x)+1.96*se(x))
par(mfrow = c(1,2))
plot.means(means = iris.3mean,
          SEs = iris.3se)

## [1] "Set categoris labls with 'categories=' "

lineplot.CI(x.factor = Species,
            response = Sepal.Length,
            ci.fun = ci.fun.use,
            data = iris)

Test the function with 4 groups

library(MASS)
data(crabs)
crab.means <- coef(lm(FL ~ -1 + sp:sex, data = crabs))
crab.categories <- names(crab.means)
crab.se.fake <- rep(sd(crabs$FL)/sqrt(dim(crabs)[1]),4)

par(mfrow = c(1,1))
plot.means(means = crab.means,
          SEs = crab.se.fake,
          categories = crab.categories)

Other approaches

These approach vary in if/how they generate the intial plot of the means before adding the error bars.

gplot::plotmeans http://svitsrv25.epfl.ch/R-doc/library/gplots/html/plotmeans.html

Using ggplot w/geom_errorbar ggplot is the standard now for advanced graphics. http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/

https://www.r-bloggers.com/line-plot-for-two-way-designs-using-ggplot2/

Using sciplot:lineplot.CI See examples above

Hmisc::errbar http://svitsrv25.epfl.ch/R-doc/library/Hmisc/html/errbar.html

Using psych::error.bars https://www.personality-project.org/r/html/error.bars.html

If you’re into “dynamite plots” https://www.r-bloggers.com/building-barplots-with-error-bars/ http://monkeysuncle.stanford.edu/?p=485

plot.means(): plot means and 95%CIs

brouwern@gmail.com