Plotting means and error bars (ggplot2)

Problem

You want to plot means and error bars for a dataset.

Solution

To make graphs with ggplot2, the data must be in a data frame, and in “long” (as opposed to wide) format. If your data needs to be restructured, see this page for more information.

Sample data

The examples below will the ToothGrowth dataset. Note that dose is a numeric column here; in some situations it may be useful to convert it to a factor.

  1. tg <- ToothGrowth
  2. head(tg)
  3. #>    len supp dose
  4. #> 1  4.2   VC  0.5
  5. #> 2 11.5   VC  0.5
  6. #> 3  7.3   VC  0.5
  7. #> 4  5.8   VC  0.5
  8. #> 5  6.4   VC  0.5
  9. #> 6 10.0   VC  0.5
  11. library(ggplot2)

First, it is necessary to summarize the data. This can be done in a number of ways, as described on this page. In this case, we’ll use the summarySE() function defined on that page, and also at the bottom of this page. (The code for the summarySE function must be entered before it is called here).

  1. # summarySE provides the standard deviation, standard error of the mean, and a (default 95%) confidence interval
  2. tgc <- summarySE(tg, measurevar="len", groupvars=c("supp","dose"))
  3. tgc
  4. #>   supp dose  N   len       sd        se       ci
  5. #> 1   OJ  0.5 10 13.23 4.459709 1.4102837 3.190283
  6. #> 2   OJ  1.0 10 22.70 3.910953 1.2367520 2.797727
  7. #> 3   OJ  2.0 10 26.06 2.655058 0.8396031 1.899314
  8. #> 4   VC  0.5 10  7.98 2.746634 0.8685620 1.964824
  9. #> 5   VC  1.0 10 16.77 2.515309 0.7954104 1.799343
  10. #> 6   VC  2.0 10 26.14 4.797731 1.5171757 3.432090

Line graphs

After the data is summarized, we can make the graph. These are basic line and point graph with error bars representing either the standard error of the mean, or 95% confidence interval.

  1. # Standard error of the mean
  2. ggplot(tgc, aes(x=dose, y=len, colour=supp)) + 
  3.     geom_errorbar(aes(ymin=len-se, ymax=len+se), width=.1) +
  4.     geom_line() +
  5.     geom_point()
  8. # The errorbars overlapped, so use position_dodge to move them horizontally
  9. pd <- position_dodge(0.1) # move them .05 to the left and right
  11. ggplot(tgc, aes(x=dose, y=len, colour=supp)) + 
  12.     geom_errorbar(aes(ymin=len-se, ymax=len+se), width=.1, position=pd) +
  13.     geom_line(position=pd) +
  14.     geom_point(position=pd)
  17. # Use 95% confidence interval instead of SEM
  18. ggplot(tgc, aes(x=dose, y=len, colour=supp)) + 
  19.     geom_errorbar(aes(ymin=len-ci, ymax=len+ci), width=.1, position=pd) +
  20.     geom_line(position=pd) +
  21.     geom_point(position=pd)
  23. # Black error bars - notice the mapping of 'group=supp' -- without it, the error
  24. # bars won't be dodged!
  25. ggplot(tgc, aes(x=dose, y=len, colour=supp, group=supp)) + 
  26.     geom_errorbar(aes(ymin=len-ci, ymax=len+ci), colour="black", width=.1, position=pd) +
  27.     geom_line(position=pd) +
  28.     geom_point(position=pd, size=3)

plot of chunk unnamed-chunk-4plot of chunk unnamed-chunk-4plot of chunk unnamed-chunk-4plot of chunk unnamed-chunk-4

A finished graph with error bars representing the standard error of the mean might look like this. The points are drawn last so that the white fill goes on top of the lines and error bars.

  1. ggplot(tgc, aes(x=dose, y=len, colour=supp, group=supp)) + 
  2.     geom_errorbar(aes(ymin=len-se, ymax=len+se), colour="black", width=.1, position=pd) +
  3.     geom_line(position=pd) +
  4.     geom_point(position=pd, size=3, shape=21, fill="white") + # 21 is filled circle
  5.     xlab("Dose (mg)") +
  6.     ylab("Tooth length") +
  7.     scale_colour_hue(name="Supplement type",    # Legend label, use darker colors
  8.                      breaks=c("OJ", "VC"),
  9.                      labels=c("Orange juice", "Ascorbic acid"),
  10.                      l=40) +                    # Use darker colors, lightness=40
  11.     ggtitle("The Effect of Vitamin C on\nTooth Growth in Guinea Pigs") +
  12.     expand_limits(y=0) +                        # Expand y range
  13.     scale_y_continuous(breaks=0:20*4) +         # Set tick every 4
  14.     theme_bw() +
  15.     theme(legend.justification=c(1,0),
  16.           legend.position=c(1,0))               # Position legend in bottom right

plot of chunk unnamed-chunk-5

Bar graphs

The procedure is similar for bar graphs. Note that tgc$size must be a factor. If it is a numeric vector, then it will not work.

  1. # Use dose as a factor rather than numeric
  2. tgc2 <- tgc
  3. tgc2$dose <- factor(tgc2$dose)
  5. # Error bars represent standard error of the mean
  6. ggplot(tgc2, aes(x=dose, y=len, fill=supp)) + 
  7.     geom_bar(position=position_dodge(), stat="identity") +
  8.     geom_errorbar(aes(ymin=len-se, ymax=len+se),
  9.                   width=.2,                    # Width of the error bars
  10.                   position=position_dodge(.9))
  13. # Use 95% confidence intervals instead of SEM
  14. ggplot(tgc2, aes(x=dose, y=len, fill=supp)) + 
  15.     geom_bar(position=position_dodge(), stat="identity") +
  16.     geom_errorbar(aes(ymin=len-ci, ymax=len+ci),
  17.                   width=.2,                    # Width of the error bars
  18.                   position=position_dodge(.9))

plot of chunk unnamed-chunk-6plot of chunk unnamed-chunk-6

A finished graph might look like this.

  1. ggplot(tgc2, aes(x=dose, y=len, fill=supp)) + 
  2.     geom_bar(position=position_dodge(), stat="identity",
  3.              colour="black", # Use black outlines,
  4.              size=.3) +      # Thinner lines
  5.     geom_errorbar(aes(ymin=len-se, ymax=len+se),
  6.                   size=.3,    # Thinner lines
  7.                   width=.2,
  8.                   position=position_dodge(.9)) +
  9.     xlab("Dose (mg)") +
  10.     ylab("Tooth length") +
  11.     scale_fill_hue(name="Supplement type", # Legend label, use darker colors
  12.                    breaks=c("OJ", "VC"),
  13.                    labels=c("Orange juice", "Ascorbic acid")) +
  14.     ggtitle("The Effect of Vitamin C on\nTooth Growth in Guinea Pigs") +
  15.     scale_y_continuous(breaks=0:20*4) +
  16.     theme_bw()

plot of chunk unnamed-chunk-7

Error bars for within-subjects variables

When all variables are between-subjects, it is straightforward to plot standard error or confidence intervals. However, when there are within-subjects variables (repeated measures), plotting the standard error or regular confidence intervals may be misleading for making inferences about differences between conditions.

The method below is from Morey (2008), which is a correction to Cousineau (2005), which in turn is meant to be a simpler method of that in Loftus and Masson (1994). See these papers for a more detailed treatment of the issues involved in error bars with within-subjects variables.

One within-subjects variable

Here is a data set (from Morey 2008) with one within-subjects variable: pre/post-test.

  1. dfw <- read.table(header=TRUE, text='
  2.  subject pretest posttest
  3.        1    59.4     64.5
  4.        2    46.4     52.4
  5.        3    46.0     49.7
  6.        4    49.0     48.7
  7.        5    32.5     37.4
  8.        6    45.2     49.5
  9.        7    60.3     59.9
  10.        8    54.3     54.1
  11.        9    45.4     49.6
  12.       10    38.9     48.5
  13.  ')
  15. # Treat subject ID as a factor
  16. dfw$subject <- factor(dfw$subject)

The first step is to convert it to long format. See this page for more information about the conversion.

  1. # Convert to long format
  2. library(reshape2)
  3. dfw_long <- melt(dfw,
  4.                  id.vars = "subject",
  5.                  measure.vars = c("pretest","posttest"),
  6.                  variable.name = "condition")
  8. dfw_long
  9. #>    subject condition value
  10. #> 1        1   pretest  59.4
  11. #> 2        2   pretest  46.4
  12. #> 3        3   pretest  46.0
  13. #> 4        4   pretest  49.0
  14. #> 5        5   pretest  32.5
  15. #> 6        6   pretest  45.2
  16. #> 7        7   pretest  60.3
  17. #> 8        8   pretest  54.3
  18. #> 9        9   pretest  45.4
  19. #> 10      10   pretest  38.9
  20. #> 11       1  posttest  64.5
  21. #> 12       2  posttest  52.4
  22. #> 13       3  posttest  49.7
  23. #> 14       4  posttest  48.7
  24. #> 15       5  posttest  37.4
  25. #> 16       6  posttest  49.5
  26. #> 17       7  posttest  59.9
  27. #> 18       8  posttest  54.1
  28. #> 19       9  posttest  49.6
  29. #> 20      10  posttest  48.5

Collapse the data using summarySEwithin (defined at the bottom of this page; both of the helper functions below must be entered before the function is called here).

  1. dfwc <- summarySEwithin(dfw_long, measurevar="value", withinvars="condition",
  2.                         idvar="subject", na.rm=FALSE, conf.interval=.95)
  4. dfwc
  5. #>   condition  N value value_norm       sd        se       ci
  6. #> 1  posttest 10 51.43      51.43 2.262361 0.7154214 1.618396
  7. #> 2   pretest 10 47.74      47.74 2.262361 0.7154214 1.618396
  9. library(ggplot2)
  10. # Make the graph with the 95% confidence interval
  11. ggplot(dfwc, aes(x=condition, y=value, group=1)) +
  12.     geom_line() +
  13.     geom_errorbar(width=.1, aes(ymin=value-ci, ymax=value+ci)) +
  14.     geom_point(shape=21, size=3, fill="white") +
  15.     ylim(40,60)

plot of chunk unnamed-chunk-10

The value and value_norm columns represent the un-normed and normed means. See the section below on normed means for more information.

Understanding within-subjects error bars

This section explains how the within-subjects error bar values are calculated. The steps here are for explanation purposes only; they are not necessary for making the error bars.

The graph of individual data shows that there is a consistent trend for the within-subjects variable condition, but this would not necessarily be revealed by taking the regular standard errors (or confidence intervals) for each group. The method in Morey (2008) and Cousineau (2005) essentially normalizes the data to remove the between-subject variability and calculates the variance from this normalized data.

  1. # Use a consistent y range
  2. ymax <- max(dfw_long$value)
  3. ymin <- min(dfw_long$value)
  5. # Plot the individuals
  6. ggplot(dfw_long, aes(x=condition, y=value, colour=subject, group=subject)) +
  7.     geom_line() + geom_point(shape=21, fill="white") + 
  8.     ylim(ymin,ymax)
  11. # Create the normed version of the data
  12. dfwNorm.long <- normDataWithin(data=dfw_long, idvar="subject", measurevar="value")
  14. # Plot the normed individuals
  15. ggplot(dfwNorm.long, aes(x=condition, y=value_norm, colour=subject, group=subject)) +
  16.     geom_line() + geom_point(shape=21, fill="white") + 
  17.     ylim(ymin,ymax)

plot of chunk unnamed-chunk-11plot of chunk unnamed-chunk-11

The differences in the error bars for the regular (between-subject) method and the within-subject method are shown here. The regular error bars are in red, and the within-subject error bars are in black.

  1. # Instead of summarySEwithin, use summarySE, which treats condition as though it were a between-subjects variable
  2. dfwc_between <- summarySE(data=dfw_long, measurevar="value", groupvars="condition", na.rm=FALSE, conf.interval=.95)
  3. dfwc_between
  4. #>   condition  N value       sd       se       ci
  5. #> 1   pretest 10 47.74 8.598992 2.719240 6.151348
  6. #> 2  posttest 10 51.43 7.253972 2.293907 5.189179
  8. # Show the between-S CI's in red, and the within-S CI's in black
  9. ggplot(dfwc_between, aes(x=condition, y=value, group=1)) +
  10.     geom_line() +
  11.     geom_errorbar(width=.1, aes(ymin=value-ci, ymax=value+ci), colour="red") +
  12.     geom_errorbar(width=.1, aes(ymin=value-ci, ymax=value+ci), data=dfwc) +
  13.     geom_point(shape=21, size=3, fill="white") +
  14.     ylim(ymin,ymax)

plot of chunk unnamed-chunk-12

Two within-subjects variables

If there is more than one within-subjects variable, the same function, summarySEwithin, can be used. This data set is taken from Hays (1994), and used for making this type of within-subject error bar in Rouder and Morey (2005).

  1. data <- read.table(header=TRUE, text='
  2.  Subject RoundMono SquareMono RoundColor SquareColor
  3.        1        41         40         41          37
  4.        2        57         56         56          53
  5.        3        52         53         53          50
  6.        4        49         47         47          47
  7.        5        47         48         48          47
  8.        6        37         34         35          36
  9.        7        47         50         47          46
  10.        8        41         40         38          40
  11.        9        48         47         49          45
  12.       10        37         35         36          35
  13.       11        32         31         31          33
  14.       12        47         42         42          42
  15. ')

The data must first be converted to long format. In this case, the column names indicate two variables, shape (round/square) and color scheme (monochromatic/colored).

  1. # Convert it to long format
  2. library(reshape2)
  3. data_long <- melt(data=data, id.var="Subject",
  4.                   measure.vars=c("RoundMono", "SquareMono", "RoundColor", "SquareColor"),
  5.                   variable.name="Condition")
  6. names(data_long)[names(data_long)=="value"] <- "Time"
  8. # Split Condition column into Shape and ColorScheme
  9. data_long$Shape <- NA
  10. data_long$Shape[grepl("^Round",  data_long$Condition)] <- "Round"
  11. data_long$Shape[grepl("^Square", data_long$Condition)] <- "Square"
  12. data_long$Shape <- factor(data_long$Shape)
  14. data_long$ColorScheme <- NA
  15. data_long$ColorScheme[grepl("Mono$",  data_long$Condition)] <- "Monochromatic"
  16. data_long$ColorScheme[grepl("Color$", data_long$Condition)] <- "Colored"
  17. data_long$ColorScheme <- factor(data_long$ColorScheme, levels=c("Monochromatic","Colored"))
  19. # Remove the Condition column now
  20. data_long$Condition <- NULL
  22. # Look at first few rows 
  23. head(data_long)
  24. #>   Subject Time Shape   ColorScheme
  25. #> 1       1   41 Round Monochromatic
  26. #> 2       2   57 Round Monochromatic
  27. #> 3       3   52 Round Monochromatic
  28. #> 4       4   49 Round Monochromatic
  29. #> 5       5   47 Round Monochromatic
  30. #> 6       6   37 Round Monochromatic

Now it can be summarized and graphed.

  1. datac <- summarySEwithin(data_long, measurevar="Time", withinvars=c("Shape","ColorScheme"), idvar="Subject")
  2. datac
  3. #>    Shape   ColorScheme  N     Time Time_norm       sd        se        ci
  4. #> 1  Round       Colored 12 43.58333  43.58333 1.212311 0.3499639 0.7702654
  5. #> 2  Round Monochromatic 12 44.58333  44.58333 1.331438 0.3843531 0.8459554
  6. #> 3 Square       Colored 12 42.58333  42.58333 1.461630 0.4219364 0.9286757
  7. #> 4 Square Monochromatic 12 43.58333  43.58333 1.261312 0.3641095 0.8013997
  9. library(ggplot2)
  10. ggplot(datac, aes(x=Shape, y=Time, fill=ColorScheme)) +
  11.     geom_bar(position=position_dodge(.9), colour="black", stat="identity") +
  12.     geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=Time-ci, ymax=Time+ci)) +
  13.     coord_cartesian(ylim=c(40,46)) +
  14.     scale_fill_manual(values=c("#CCCCCC","#FFFFFF")) +
  15.     scale_y_continuous(breaks=seq(1:100)) +
  16.     theme_bw() +
  17.     geom_hline(yintercept=38)

plot of chunk unnamed-chunk-15

Note about normed means

The summarySEWithin function returns both normed and un-normed means. The un-normed means are simply the mean of each group. The normed means are calculated so that means of each between-subject group are the same. These values can diverge when there are between-subject variables.

For example:

  1. dat <- read.table(header=TRUE, text='
  2. id trial gender dv
  3.  A     0   male  2
  4.  A     1   male  4
  5.  B     0   male  6
  6.  B     1   male  8
  7.  C     0 female 22
  8.  C     1 female 24
  9.  D     0 female 26
  10.  D     1 female 28
  11. ')
  13. # normed and un-normed means are different
  14. summarySEwithin(dat, measurevar="dv", withinvars="trial", betweenvars="gender",
  15.                 idvar="id")
  16. #> Automatically converting the following non-factors to factors: trial
  17. #>   gender trial N dv dv_norm sd se ci
  18. #> 1 female     0 2 24      14  0  0  0
  19. #> 2 female     1 2 26      16  0  0  0
  20. #> 3   male     0 2  4      14  0  0  0
  21. #> 4   male     1 2  6      16  0  0  0

Helper functions

The summarySE function is also defined on this page. If you only are working with between-subjects variables, that is the only function you will need in your code. If you have within-subjects variables and want to adjust the error bars so that inter-subject variability is removed as in Loftus and Masson (1994), then the other two functions, normDataWithin and summarySEwithin must also be added to your code; summarySEwithin will then be the function that you call.

  1. ## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
  2. ##   data: a data frame.
  3. ##   measurevar: the name of a column that contains the variable to be summariezed
  4. ##   groupvars: a vector containing names of columns that contain grouping variables
  5. ##   na.rm: a boolean that indicates whether to ignore NA's
  6. ##   conf.interval: the percent range of the confidence interval (default is 95%)
  7. summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
  8.                       conf.interval=.95, .drop=TRUE) {
  9.     library(plyr)
  11.     # New version of length which can handle NA's: if na.rm==T, don't count them
  12.     length2 <- function (x, na.rm=FALSE) {
  13.         if (na.rm) sum(!is.na(x))
  14.         else       length(x)
  15.     }
  17.     # This does the summary. For each group's data frame, return a vector with
  18.     # N, mean, and sd
  19.     datac <- ddply(data, groupvars, .drop=.drop,
  20.       .fun = function(xx, col) {
  21.         c(N    = length2(xx[[col]], na.rm=na.rm),
  22.           mean = mean   (xx[[col]], na.rm=na.rm),
  23.           sd   = sd     (xx[[col]], na.rm=na.rm)
  24.         )
  25.       },
  26.       measurevar
  27.     )
  29.     # Rename the "mean" column    
  30.     datac <- rename(datac, c("mean" = measurevar))
  32.     datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  34.     # Confidence interval multiplier for standard error
  35.     # Calculate t-statistic for confidence interval: 
  36.     # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  37.     ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  38.     datac$ci <- datac$se * ciMult
  40.     return(datac)
  41. }
  1. ## Norms the data within specified groups in a data frame; it normalizes each
  2. ## subject (identified by idvar) so that they have the same mean, within each group
  3. ## specified by betweenvars.
  4. ##   data: a data frame.
  5. ##   idvar: the name of a column that identifies each subject (or matched subjects)
  6. ##   measurevar: the name of a column that contains the variable to be summariezed
  7. ##   betweenvars: a vector containing names of columns that are between-subjects variables
  8. ##   na.rm: a boolean that indicates whether to ignore NA's
  9. normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
  10.                            na.rm=FALSE, .drop=TRUE) {
  11.     library(plyr)
  13.     # Measure var on left, idvar + between vars on right of formula.
  14.     data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
  15.      .fun = function(xx, col, na.rm) {
  16.         c(subjMean = mean(xx[,col], na.rm=na.rm))
  17.       },
  18.       measurevar,
  19.       na.rm
  20.     )
  22.     # Put the subject means with original data
  23.     data <- merge(data, data.subjMean)
  25.     # Get the normalized data in a new column
  26.     measureNormedVar <- paste(measurevar, "_norm", sep="")
  27.     data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
  28.                                mean(data[,measurevar], na.rm=na.rm)
  30.     # Remove this subject mean column
  31.     data$subjMean <- NULL
  33.     return(data)
  34. }
  1. ## Summarizes data, handling within-subjects variables by removing inter-subject variability.
  2. ## It will still work if there are no within-S variables.
  3. ## Gives count, un-normed mean, normed mean (with same between-group mean),
  4. ##   standard deviation, standard error of the mean, and confidence interval.
  5. ## If there are within-subject variables, calculate adjusted values using method from Morey (2008).
  6. ##   data: a data frame.
  7. ##   measurevar: the name of a column that contains the variable to be summariezed
  8. ##   betweenvars: a vector containing names of columns that are between-subjects variables
  9. ##   withinvars: a vector containing names of columns that are within-subjects variables
  10. ##   idvar: the name of a column that identifies each subject (or matched subjects)
  11. ##   na.rm: a boolean that indicates whether to ignore NA's
  12. ##   conf.interval: the percent range of the confidence interval (default is 95%)
  13. summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
  14.                             idvar=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
  16.   # Ensure that the betweenvars and withinvars are factors
  17.   factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
  18.     FUN=is.factor, FUN.VALUE=logical(1))
  20.   if (!all(factorvars)) {
  21.     nonfactorvars <- names(factorvars)[!factorvars]
  22.     message("Automatically converting the following non-factors to factors: ",
  23.             paste(nonfactorvars, collapse = ", "))
  24.     data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
  25.   }
  27.   # Get the means from the un-normed data
  28.   datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
  29.                      na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
  31.   # Drop all the unused columns (these will be calculated with normed data)
  32.   datac$sd <- NULL
  33.   datac$se <- NULL
  34.   datac$ci <- NULL
  36.   # Norm each subject's data
  37.   ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)
  39.   # This is the name of the new column
  40.   measurevar_n <- paste(measurevar, "_norm", sep="")
  42.   # Collapse the normed data - now we can treat between and within vars the same
  43.   ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
  44.                       na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
  46.   # Apply correction from Morey (2008) to the standard error and confidence interval
  47.   #  Get the product of the number of conditions of within-S variables
  48.   nWithinGroups    <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
  49.                            FUN.VALUE=numeric(1)))
  50.   correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )
  52.   # Apply the correction factor
  53.   ndatac$sd <- ndatac$sd * correctionFactor
  54.   ndatac$se <- ndatac$se * correctionFactor
  55.   ndatac$ci <- ndatac$ci * correctionFactor
  57.   # Combine the un-normed means with the normed results
  58.   merge(datac, ndatac)
  59. }