NOTE: This post only contains information on repeated measures ANOVAs, and not how to conduct a comparable analysis using a linear mixed model. For that, be on the lookout for an upcoming post!

When I was studying psychology as an undergraduate, one of my biggest frustrations with R was the lack of quality support for repeated measures ANOVAs.They’re a pretty common thing to run into in much psychological research, and having to wade through incomplete and often contradictory advice for conducting them was (and still is) a pain, to put it mildly.

Thankfully, though, they’re not too tricky to set up once you figure out what you’re doing.

To get started, let’s construct a phony data set where we’re measuring participant stress on a 100-point scale. Higher numbers mean the participant is more stressed out. For our experimental manipulation, let’s say that participants are exposed to a series of several images presented with various background music playing. The images can depict scenes that are happy or angry. The background music can be a Disney soundtrack or music from a horror movie. Each participant sees multiple images and listens to multiple music samples. (Your variables can have more than 2 factors, and you can include more than 2 IVs. We’re just keeping it simple for the purposes of explanation!)

First, here’s the code we’ll use to generate our phony data:

set.seed(5250)

myData <- data.frame(PID = rep(seq(from = 1,
                               to = 50, by = 1), 20),
                     stress = sample(x = 1:100,
                                     size = 1000,
                                     replace = TRUE),
                     image = sample(c("Happy", "Angry"),
                                    size = 1000,
                                    replace = TRUE),
                     music = sample(c("Disney", "Horror"),
                                    size = 1000,
                                    replace = TRUE)
)

myData <- within(myData, {
  PID   <- factor(PID)
  image <- factor(image)
  music <- factor(music)
})

myData <- myData[order(myData$PID), ]
head(myData)

PID stress image  music
  1     90 Happy Disney
  1     70 Angry Horror
  1     61 Angry Horror
  1     87 Happy Horror
  1     79 Happy Disney
  1     95 Happy Horror

So we see that we have one row per observation per participant. If your dataset is in wide form rather than long, I’d suggest checking out our article on converting between wide and long since everything from this point out assumes that your data look like what’s shown above!

Extracting Condition Means

Before we can run our ANOVA, we need to find the mean stress value for each participant for each combination of conditions. We’ll do that with:

myData.mean <- aggregate(myData$stress,
                      by = list(myData$PID, myData$music,
                              myData$image),
                      FUN = 'mean')

colnames(myData.mean) <- c("PID","music","image","stress")

myData.mean <- myData.mean[order(myData.mean$PID), ]
head(myData.mean)

PID  music   image   stress
  1 Disney   Angry 39.33333
  1 Horror   Angry 65.50000
  1 Disney   Happy 68.00000
  1 Horror   Happy 69.57143
  1 Disney Neutral 40.00000
  1 Horror Neutral 52.66667

So now we’ve gone from one row per participant per observation to one row per participant per condition. At this point we’re ready to actually construct our ANOVA!

Building the ANOVA

Now, our actual ANOVA is going to look something like this:

stress.aov <- with(myData.mean,
                   aov(stress ~ music * image +
                       Error(PID / (music * image)))
)

But what’s all that mean? What’s with that funky Error() term we threw in there? Pretty simple: what we’re saying is that we want to look at how stress changes as a function of the music and image that participants were shown. (Thus the stress ~ music * image) The asterisk specifies that we want to look at the interaction between the two IVs as well. But since this was a repeated measures design, we need to specify an error term that accounts for natural variation from participant to participant. (E.g., I might react a little differently to scary music than you do because I love zombie movies and you hate them!) We do this with the Error() function: specifically, we are saying that we want to control for that between-participant variation over all of our within-subjects variables.

Now that we’ve specified our model, we can go ahead and look at the results:

summary(stress.aov)

Error: PID
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 49   8344   170.3               

Error: PID:music
          Df Sum Sq Mean Sq F value Pr(>F)
music      1      1    0.78   0.003  0.954
Residuals 49  11524  235.19               

Error: PID:image
          Df Sum Sq Mean Sq F value Pr(>F)
image      1     61   61.11   0.296  0.589
Residuals 49  10127  206.66               

Error: PID:music:image
            Df Sum Sq Mean Sq F value Pr(>F)
music:image  1    564   563.8   2.626  0.112
Residuals   49  10520   214.7  

We see that there is no main effect of either music:

F(1, 49) = 0.003; p-value = 0.954

or image:

F(1, 49) = 0.296; p-value = 0.589

on participant stress. Likewise, we see that there is not a significant interaction effect between the two independent variables:

F(1, 49) = 2.626; p-value = 0.112

What do I do with my Between-Subjects Effects?

This has all been fine and good, but what if you have an independent variable that’s between-subjects? To continue our previous example, let’s say that some participants could only come in during the day and some could only come in at night. Our data might instead look like this:

set.seed(5250)

myData <- data.frame(PID = rep(seq(from = 1,
                               to = 50, by = 1), 20),
                     stress = sample(x = 1:100,
                                     size = 1000,
                                     replace = TRUE),
                     image = sample(c("Happy", "Angry"),
                                    size = 1000,
                                    replace = TRUE),
                     music = sample(c("Disney", "Horror"),
                                    size = 1000,
                                    replace = TRUE),
                     time = rep(sample(c("Day", "Night"),
                                       size = 50,
                                       replace = TRUE), 2))

head(myData)


PID stress image  music  time
  1     66 Happy Disney   Day
  2     21 Happy Disney Night
  3     25 Angry Horror   Day
  4     61 Happy Disney   Day
  5     11 Angry Disney Night
  6     85 Angry Horror   Day

From there, the steps we take look pretty similar to before:

myData <- within(myData, {
  PID   <- factor(PID)
  image <- factor(image)
  music <- factor(music)
  time  <- factor(time)
})

myData <- myData[order(myData$PID), ]
head(myData)

myData.mean <- aggregate(myData$stress,
                         by = list(myData$PID, myData$music,
                                 myData$image, myData$time),
                         FUN = 'mean')

colnames(myData.mean) <- c("PID", "music", "image",
                           "time", "stress")
myData.mean <- myData.mean[order(myData.mean$PID), ]

stress.aov <- with(myData.mean, aov(stress ~ time * music *
                                    image + Error(PID /
                                    (music * image))))
summary(stress.aov)

Error: PID
          Df Sum Sq Mean Sq F value Pr(>F)
time       1      6     5.7   0.033  0.857
Residuals 48   8338   173.7               

Error: PID:music
           Df Sum Sq Mean Sq F value Pr(>F)
music       1      1    0.78   0.003  0.955
time:music  1     22   21.96   0.092  0.763
Residuals  48  11502  239.63               

Error: PID:image
           Df Sum Sq Mean Sq F value Pr(>F)
image       1     61   61.11   0.292  0.591
time:image  1     81   80.77   0.386  0.537
Residuals  48  10046  209.29               

Error: PID:music:image
                 Df Sum Sq Mean Sq F value Pr(>F)
music:image       1    564   563.8   2.578  0.115
time:music:image  1     24    23.7   0.109  0.743
Residuals        48  10496   218.7

The only big difference is that we don’t include out between-subjects factor (time) in the Error() term. In any case, we see that there are no significant main effects (of time, music, or image) nor any significant interactions (between time and music, time and image, music and image, or music and time and image).

Dealing with “Error() model is singular”

Sometimes you might be unlucky enough to get this error when you try to specify your aov() object. It’s not the end of the world, it just means that you don’t have an observation for every between-subjects condition for every participant. This can happen due to a bug in your programming, a participant being noncompliant, data trimming after the fact, or a whole host of other reasons. The moral of the story, though, is that you need to find the participant that is missing data and drop him or her from this analysis for the error to go away. Or, if the idea of dropping a participant entirely rubs you the wrong way, you could look into conducting the analysis as a linear mixed model. We don’t have a tutorial for that (yet!), but keep your eyes peeled: as soon as it’s written, we’ll update this post and link you to it!