Search code examples
rggplot2graphline

Adding Error bars from ANOVA + Tukey's posthoc in a Line Graph in R (ggplot2)

I'm trying to add error bars from an ANOVA analysis with Tukey's post-hoc to my line graph in R using ggplot2. Thanks in advance, and I'm happy to post any more info that will help me solve this problem!

Here is my code so far:

#### ANOVA ####
aov_CW <- aov(CW ~ Subject, data = VisualAcuity)
summary(aov_CW)
model.tables(aov_CW, "means")

#### Tukey HSD post-hoc Test ####
TukeyHSD(aov_CW, conf.level = 0.95)

#### Line Plot ####

VisualAcuity$Subject <- as.factor(VisualAcuity$Subject)
VisualAcuity$Day <- as.factor(VisualAcuity$Day)
ggplot(data=VisualAcuity, aes(x=Day, y=CW, group = Subject)) +
  geom_line(size=1, aes(color = Subject)) +
  geom_point(size=2, shape=21, aes(color = Subject, fill = Subject)) +
  ylim(0, max(1)) + ylab ('Visual Acuity (CW)') +
  geom_errorbar(aes(ymin = CW - se, ymax = CW + se))

Sample data for question

Solution

  • The values produced by TukeyHSD are for pairwise comparisons, and thus can't be plotted on the graph you provided code fore.

    However, here is an approach to add standard error sd/sqrt(n) bars with dplyr. Obviously you can't add error bars for 1 observation, so I expanded your dataset with some random numbers. 

    library(dplyr)
library(ggplot2)
VisualAcuity$Subject <- as.factor(VisualAcuity$Subject)
VisualAcuity$Day <- as.factor(VisualAcuity$Day)

VisualAcuity %>% 
  dplyr::group_by(Day,Subject) %>% 
  dplyr::summarize(Mean = mean(CW),std_err = sd(CW)/sqrt(n()), n = n()) %>%
ggplot(aes(x=Day, y=Mean, color = Subject, group = Subject)) +
  geom_line(size=1,) +
  geom_point(size=2, shape=21, aes(fill = Subject)) +
  ylim(0, max(1)) + ylab ('Visual Acuity (CW)') +
  geom_errorbar(aes(ymin = Mean - std_err, ymax = Mean + std_err))

    enter image description here

    Data

    VisualAcuity <- structure(list(Subject = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 
3L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 
4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("1", 
"2", "3", "4"), class = "factor"), Day = structure(c(1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L), .Label = c("1", "2", "3"), class = "factor"), CW = c(0.528, 
0.56, 0.486, 0.436, 0, 0.525, 0.622, 0.6, 0.522, 0.453, 0.5, 
0.494, 0.566, 0.606, 0.56, 0.509, 0.054, 0.593, 0.645, 0.668, 
0.56, 0.51, 0.508, 0.533, 0.518, 0.502, 0.413, 0.412, 0.042, 
0.431, 0.582, 0.508, 0.435, 0.368, 0.417, 0.485)), row.names = c(NA, 
-36L), class = "data.frame")