RStudio Workshop Series: Repeated Measures ANOVA

Welcome everyone! Today we are doing RStudio workshop series lesson #9: repeated measures ANOVA, and then the non-parametric version which is the Friedman's test. So a repeated measures ANOVA, or “ANalysis Of VAriance”, is used to determine whether the means of three or more continuous variables (or measures at three or more time points on the same continuous variable) differs for one group. So fancy way of saying: if I was a participant in some kind of study, I might have a score at timepoint 1, timepoint 2, timepoint 3, and we're measuring the same three things across those three timepoints. So it might be happiness at time 1, 2, and 3. But because I have three scores and I'm a participant in the study, this is a repeated measures ANOVA. 
It's a parametric test. So what we mean when we say it's a parametric test, is anytime we say parametric, we mean “normality”: highest amount of data in the middle, lowest amount of data at the ends. That typical bell-shaped curve. Our data have to look a certain way in order for this test to be appropriate. 
I've left you a couple guides in the code if you wanted to follow along. So we've got my RStudio LibGuide, where we've got some help on this available. We've got the RStudio documentation, where there's some help available. Today, we are using the ezANOVA() function; there are other functions to run this test, but we're just going to cover one of them today. And I've left you a guide that I have not written, but it's pretty well done, it's called Laerd statistics. So they also have a guide on how to run this test.  
If you are new here, the first thing we need to do is we have to tell the computer: “Where is the file located that we are looking to work with today?”. So on line 23, this is a function. We know it's a function because it's some text next to round brackets. So it's the setwd() function. The easiest way to get this for your own directory, because right now this is pointing to MY directory, is: you can click the “Session” button, you can click “Set Working Directory”, and you can click “Choose Directory”. What this will do, is if you click Session > Set Working Directory > Choos1e Directory, it will open a File Explorer where you can go through the files on your computer to select where YOUR file is located today. So if I do this, it will open up my directory, and it says Documents > Classes and Workshops > Lindsay's workshops > MICRO WORKSHOPS > RStudio. On my computer, it doesn't look like there's anything in this file folder structure. But that's okay, I know this is where my file is located, so I can say “Open”. You might have your file on your Desktop or in your Downloads, or in a different folder structure; just go to where you know where the file is, and say “Open”. What this does, is if it worked, you'll get some blue text down in the bottom-left in the Console that starts with setwd(). So that set working directory function. And you can copy this and replace what is on line 23. If you run line 23 as written, it won't work on your computer because it's for my file folder structure, and you probably don't have the same file folder structure. So copy this and paste it on line 23 and the next time you run this code, it will point to where your file is located. The other way we can tell this worked is down in the bottom-right in the Files window, we can actually see all of the files that I have in this folder structure on my computer, so we know it worked because we're pointing at the location that has the file we want to open. 
So we've set our working directory; we have told the computer we would like to work out of this location. The next thing we need to do is actually open the file we're trying to open, and give it a name. So my file is called “RStudio micro 9 RM ANOVA.R” … oh, I lied. We're actually opening “Fake_Data.csv”. That's the R code, silly me. We're opening: Fake_Data.csv. So that's what we're going to open. We want to give it a name; you generally want to give your stuff a meaningful title so that you know which dataset you are working with. I have just called mine Fake_Data. So you give it a name. We're going to open a .csv file, so we use the function read.csv(). 
Inside those brackets and those quotation marks, we're going to give it the exact name as the file as written. So in my case, it's: Fake_Data.csv. You need the file extension. If you had a different file type, you would need a different package to run this and a different function to run this. But today it's a .csv file, so we can just highlight this line or keep the cursor anywhere on this line [Fake_Data = read.csv(“Fake_Data.csv”)], and click the Run button, kind of in the top middle up here. We get some blue text in the bottom-left in the Console, so it means it worked. And in the top-right in our Environment, we now have a dataset called Fake_Data. We can look at this if we want, we can just click where it says Fake_Data and it will open it like it’s in an Excel sheet. So we know that this worked! We have data open in RStudio that we can work with. 
All right. Every single test we do has certain assumptions that must be met in order for us to be able to trust the results of the test. So the first thing we're going to do, is we're going to check a few assumptions for this test, and then we're going to run the test. So assumption #1: our dependent variable that we want to use must be continuous. So to run a repeated measures ANOVA, we need a continuous dependent variable. RStudio has a few different data types that mean continuous, so that could be: <int> for integer, <dbl> double, or <num> numeric. Those all count as continuous data types within RStudio. If you haven't come to a workshop before, you will need to potentially Run line 30. This is installing the package tidyverse (so taking it from the Internet, and putting it on your computer so you can use it). I have tidyverse on my computer, I don't need to Run line 30. But I do need to Run line 31: I'm going to Run library(tidyverse). This is telling the computer: “I have this package on my computer, I would like to use this package right now”. So we can Run line 31…I can click Run, we get some blue text. We get a little bit of other stuff in the Console down here as well, but that's okay, there's nothing saying it didn't work, so we're good to go. So what I want to do right now, is look at my Fake_Data dataset to be able to say: “Do I have the correct data type to run this test?”. So on line 32, we're going to take a glimpse() of our Fake_Data dataset [glimpse(Fake_Data)]. I can click Run and it spits out some information down here in the Console. So it says you took a glimpse() of Fake_Data. It has 30 rows and 7 columns (which is good because that matches what's in our Global Environment). And it gives the dollar sign ($) indicator to say: “Which column am I looking at right now?”. We've got a $Gender column which is labeled as <int>. We've got $Fake_Data1 for <dbl>, and this is actually one of the things we're going to be looking at today, so this is the correct data type because it's listed as double, which is RStudio’s way of saying “continuous”. We're also going to use Fake_Data2 (which is listed as <dbl>), and Fake_Data3 (which is listed as <dbl>). So we pass this assumption: our three variables that we're looking for here are listed as continuous (double) format in RStudio. So we pass assumption one, we can move to assumption two.  
Assumption #2: our independent variable must be categorical with three or more related or matched groups. So in this case, if we View() our Fake_Data (which we've already done: you can either click the data, click the tab that we already have, or click View(Fake_Data) here)...if we click Run, we're going to say each row is a unique participant. And our independent variable is going to be going across here: so Fake_Data1, 2, and 3. These might be three separate timepoints, so maybe this is the score for Fake_Data at timepoint 1, the score for Fake_Data at timepoint 2, and the score for Fake_Data at timepoint 3. We have three separate groups, and each person or each observation (each row) has a score for each of those three groups. So we've passed that second assumption; we have three different related or matched or paired (there's different ways to say it)…we've paired data for each person. So we pass that. 
Assumption #3. I said at the beginning that this test is a parametric test, which means we must pass the normality assumption. So assumption 3, the dependent variable is approximately normally distributed for each group or condition. That means if we have three timepoints, we need to check normality for each timepoint (so for Fake_Data1, 2, and 3). If you have a smaller dataset, we normally say 50 or fewer observations, we can run a statistic called the Shapiro-Wilk statistic to check for normality. If we have a larger data set, we normally use a different statistic, it's called the Kolmogorov-Smirnov statistic. Today we have smaller groups, so we're going to run Shapiro-Wilk. Some fields check normality in a different way, they'll do something like a histogram, and then they'll say: “Does it look like a bell-shaped curve?”. You could do that as well, it depends what's common in your field. But today we will do Shapiro-Wilk. To run this, it's: shapiro.test(Fake_Data$Fake_Data1) to check that column (Fake_Data1). And then you do the same thing for Fake_Data2 and Fake_Data3. I can Run all three of these at the same time if I want; I can highlight all three and click Run, and then we'll get three pieces of output down in the Console. So it tells me you ran Shapiro-Wilk test on the Fake_Data dataset. And which column? It's the Fake_Data1 column. This is an assumption check, so here we are expecting p to be greater than (>) .05 to say we passed this assumption. If p was less than (<) .05 for even one of the three groups, if you fail even one, you failed for the entire test. So let's look and see. So it tells you which data you used, and my p-value is listed right here. So if my first group is .4551, I have passed normality for my first group. My second group though is .002933, which is less than (<) .05; group 2 has failed normality. And Group 3, my p-value is .009086, which is also less than (<) .05. Two of our three groups have failed normality. Because even one failed, it means we've failed normality for this test. And you might be like: “Oh, shoot. Now what do I do Lindsay?”. If you failed, normality for even one of your groups and you're doing a repeated measures ANOVA, you can switch to the non-parametric version of the test if you so choose. And the non-parametric version for this test is the Friedman's test, which we'll cover in just a couple minutes. Or if you're watching the re-run, we'll cover in the next video. 
Today, we're going to pretend everything was fine, and we passed normality, and we're allowed to use the repeated measures ANOVA. We didn't actually pass this, so we want to switch to the other test, but today just for practice, we're going to run it anyway. 
All right, assumption #4: there should be no significant outliers for each group or condition. So again, if I have three timepoints, I need to be checking for outliers in each of those three timepoints. Different fields or areas of study have different opinions on how to check for outliers and remove outliers, so I haven't left you any code here, but I've left you some text that you can read. So there's different methods. You might do the mean and standard deviation for the groups, and then if it's too far away from that mean you would remove it. You might be doing something like a histogram or a boxplot; if some of the data looks really far away, you would remove the data that way. So if you're on a research team, check in with your research lead and ask: “How do we find and identify and then remove outliers if there are outliers present?”. So an expert tip here as well: if you have outliers and you remove them, check your normality again because sometimes removing outliers can fix your normality assumption and then you can use the repeated measures ANOVA. That's our expert tip.  
Okay, assumption #5: this is called sphericity. This is any time you have a within factor for an ANOVA, you need to check the assumption of sphericity. We actually check this when we run the test. So ezANOVA() has this sphericity check built into it so we can check this in just a minute so we can pretend we passed all of our assumptions today? How do you actually run the repeated measures ANOVA and what are we doing? We would like to check whether the means of Fake_Data1, Fake_Data2, and Fake_Data3 are different. Maybe this is a happiness, or I give people treats at certain timepoints and see if their happiness goes up. Maybe this is an insulin score, I'm looking to see what their insulin is if I give them treats versus if I give them some drug treatment. But we're going to see: “Are the means different here?”. We're going to be using the ezANOVA() function which requires the “ez” package. If you don't have this package already, you will need to highlight and Run line 74, which is: install.packages(“ez”). I already have this on my computer, so I will just Run line 75 which is loading the library saying: “I already have this, I would like to use this right now”. So I can Run this [library(ez)]. I do get a little red text in the Console, it says: “Warning message: package ‘ez’ was built under R version 4.3.3”. You can always check your RStudio version by clicking the “Help” button and then saying “About RStudio” (I'm actually on 2023.06.1). This isn't an error code even though it's red text; it's literally just a notice saying it was built under a different version than what we have. 

If it [the error code] says it's not available for this version of R and you're on the 2026 version, it means you're probably on the newest version that exists. It might not let you use the ez package, which is not fantastic. What it means is if you wanted to use this specific function to run an ANOVA, you either need to download an older version of RStudio or you could do the fun thing: there's so many packages to do whatever you want. You don't have to use the ezANOVA() function; anova_test is another really common one that I recommend for a lot of folks as well. It's not in our slides today or in our code today, but if you are interested, there's another test that you could run instead that does the same test. 

The other thing I can mention while we're here actually is this version of the anova_test() works if you have unequal groups. So ezANOVA() only generally works if you have equal groups, no missing data, all that fun stuff. anova_test tends to be a little bit more flexible, but it's a little harder to find the sphericity test, so I show normally the other version of the test [ezANOVA()]. There are other versions as well, so like it gives you directly within the documentation, Anova() and aov(). Those are some options as well, but this one [anova_test()] is normally quite flexible and the one a lot of people are using if they're not using ezANOVA(). So you've got some options.  

So what we can do is we can actually review the ezANOVA package or function to see what this is doing. So we can say ?ezANOVA(). We can Run this and it will open our Help window to tell you a little bit about the package. 
So the function is ezANOVA(), the squiggle brackets “{ez}” tell you it's from this specific package, and it says: Compute ANOVA. Description: This function provides easy analysis of data from factorial experiments, including purely within-Ss [capital S lowercase s means within-“subjects”] designs (a.k.a. “repeated measures”), purely between-Ss (so between-subject) designs, and mixed within-and-between-Ss [subject] designs, yielding ANOVA results, generalized effect sizes and assumption checks”. And it says: “Usage”, and it's got a whole whack of different arguments you can give it, and it gives you information about the arguments that you can give it, and if you scroll all the way down it also gives you some examples. Might sound a little confusing, let me show you how this works. All right, this requires the data to be in long format going DOWN the page in order to run the test. Some software likes stuff going across the page so our fake dataset goes across the page in that we have Fake_Data1, Fake_Data2, Fake_Data3. If you're using this specific test, it wants it to go DOWN the page: so you'll have a “Participant” column (so like participant one, one, and one), and then they have a “Score” for Fake_Data1, Fake_Data2, Fake_Data3. I've left some code to do this, so if you want you can Run line 79 [Fake_Data_Long = gather(Fake_Data, Fake_Data_group, Fake_Data_score, c("Fake_Data1","Fake_Data2","Fake_Data3"))]; we're going to switch stuff from wide to long. We can Run line 81 [Fake_Data_Long = Fake_Data_Long[,-c(1:4)]]. We can Run line 83 [Fake_Data_Long$Participant = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30)]. And then on 85 we can View() the data to make sure all of that worked [View(Fake_Data_Long)]. When we View() the data, the new data we have is: we have which “group” they are in (so whether it's Fake_Data1, scroll, Fake_Data2, scroll a bit more, Fake_Data3), the “score” for that unique person, and then the “Participant” column. So each person has a score for [Fake_Data]1, and if we scroll down they have a score for [Fake_Data]2, and if we scroll down they have a score for [Fake_Data]3. So we just took it from going across the page to going down the page. So that worked. I can go into detail about that code if we want; all you need to know is it just took it from wide to long so that you can actually use it today. We can take a glimpse() of our data just to make sure it's in the right data type still. And we see the [Fake_Data_]group is listed as character. The [Fake_Data_]score is listed as double or in this case, this means continuous data. And Participant is listed as double as well, which is not great, but that's okay. The thing that we care about is that “group”, which we're going to use in a minute, should probably be listed as “factor” (i.e., which group are you in, 1, 2, or 3?). That's a categorical variable. We want it to be treated as factor within RStudio. 
And it would be great if Participant was treated as factor as well. So we can do that on the next two lines. So on line 88, we're going to take our long data set, we're going to take the “group” variable, and we're going to use the factor() function: as.factor() to say “treat it as a factor, change that data type so it's the right data type”. [Fake_Data_Long$Fake_Data_group = as.factor(Fake_Data_Long$Fake_Data_group)]. We're also going to do that in our “long” data set with the Participant column; we're going to use the as.factor() function to make sure it's treated as the correct data type. [Fake_Data_Long$Participant = as.factor(Fake_Data_Long$Participant)]. I can Run both of these. I can take another glimpse() of my data [glimpse(Fake_Data_Long)], and now we've got: “group” is factor, “score” is double or continuous, and “Participant” is factor. So we've got the correct data types. 
All right, the format we need when we are using the ezANOVA() function is we're going to say “dv” (or dependent variable) is equal to in this case it's our “score” column; what was the score that the participant had. “wid” is your ID [identification] column; our ID column is going to be the “Participant”. This tells us which column or which columns match up with our different groups. So if I was participant one in the study, I need to have a score for each of those things; the ID column says as we're going down the page, which three timepoints go for me, for example. And we have a “within” argument here; that's our grouping column (you're in Fake_Data1, Fake_Data2, Fake_Data3; those are our three groups, so you have to say which is your within-subject variable). You have to give it the dataset; we need the long form dataset, so it's going be Fake_Data_Long. And we say “detailed = TRUE” because that gives us our sphericity stuff; that gives us one of our assumption checks. So if I highlight this and click Run, we get some output in the Console down here. So we say we've got: ezANOVA(dv = Fake_Data_score, wid = Participant, within = Fake_Data_group, data = Fake_Data_Long, detailed = TRUE).  
There are three main pieces that we get output (three chunks to this output). The 1st is “ $ANOVA ”: this is the result of your ANOVA test. We'll cover that in a second. The second chunk we get is “ $’Mauchly’s Test for Sphericity’ ”: this is your sphericity check, so this is an assumption check we're expecting p to be greater than (>) .05 to say we passed this assumption. And here, our p-value is this score here, it's .6179884; we have passed this sphericity assumption. If you failed this sphericity assumption, if p is less than (<) .05, you need to look at the third chunk where it says “ $‘Sphericity Correction’ ”: you can apply a correction if you have failed sphericity that adjusts a few pieces of your test. We don't need to look at this right now because we pass sphericity, but if you have failed sphericity, you could either add the Greenhouse-Geisser [GGe] correction or the Huynd-Feldt [HFe] correction to your test. 
So we pass sphericity, we're going to look at this first chunk, the ANOVA chunk. We generally don't read the “(Intercept)” line for most things unless you're doing like fancy regression modelling. We're going to look where it says “Fake_Data_group” because this tells us about the main effect of “group”, which is saying: “Is there a difference somewhere in the scores of Fake_Data1, Fake_Data2, and Fake_Data3?”. We're looking at the means and saying: “Are they different somewhere?”. And when we look across, we've got our degree of freedom for the numerator [DFn], our degree of freedom for the denominator [DFd], we've got some sums of squares [SSn], some more sums of squares [SSd], our F-statistic [F], our p-value [p], a star (*) denoting that our p-value is statistically significant [p<.05], and here it's listing generalized eta-squared [ges]. We generally would want partial eta-squared [pes], and we could adjust that easily if we wanted to, but don't worry about that today, today we're just doing an intro. 
So the thing we generally care about the most is actually going to be this thing right here: this is our p-value. We have “7.613627e-56”. That means significant. The “e-56” means this decimal place needs to move 56 places to the left, so we're going to have .000…so many zeros…and then a “7”. Because our p-value is less than (<) .05, we can say there's a difference in our fake data values SOMEWHERE. We don't know where the differences are between groups 1, 2, and 3, but there are differences somewhere. 
If we would like to know where those differences are, we're going to have to follow them up. So we're going to have to do some post-hoc testing in a moment. I did say this is “ges” [generalized eta-squared], you can ask for “pes” [partial eta-squared]. I've left you a calculation for “pes” here as well, if you care. We're going to take the numerator and we're going to take the denominator sums of squares and do a little calculation to get that value if you want to report that for something like a paper. So this would be your partial eta-squared value; that's your effect size for an ANOVA. So it's an intro class, but I still gave you the hard and fast stuff anyway.  
Okay. If you have missing data or if you have unequal groups, this function might not work for you. You might have to switch to a different test like Anova_test() or aov() or ANOVA(). This is an easy test in that it gives you the sphericity output very quickly and easily, and it gives you corrections for that very quickly and easily. But if you've got a “not perfect” dataset (if you got some missing data, some unequal groups, something weird going on with your data), you might need a different test and that's okay.  
I've left you how you actually write the results here. So we would say capital and italicized F(2, 58) = 2277.921…so that's this value here, that's your F-statistic. I wrote out your p-value with all of the zeros for you, if you would like to count them (it's a lot of zeros). Some folks do scientific notation instead, that might vary by your field. And then we've got our partial eta-squared value. That would be how you would write this in a paper, for example.  
The other thing you might want is the mean of each group to say: “Well, what was the mean of Fake_Data1, 2, and 3?”. So we can Run all of these at the same time, it's just the mean() function. [mean(Fake_Data1)]. The mean of Fake_Data1 is this one here: 87.10212. The mean of Fake_Data2 is 92.4917. And the mean of Fake_Data3 is 112.4268. Importantly, the means with the p-value do not tell you which groups are different. We still don't know which groups are statistically significantly different from each other, we just know that they are different. So the next thing we want to do is visualize the data and then follow up and say: “Well, where are those differences?”. So this is not a graphing class, we can go into more detail if folks have questions about this later, but we're going to Run lines 112 to 115 to make a quick and simple graph of our data. And I can click the Zoom button from the Plots window to see it a little bit bigger. Here is our graph, we have the means on our y-axis, and on our x-axis we have Fake_Data1, Fake_Data2, and Fake_Data3. The graph is not enough to tell us where the differences are. They look different, but do we know (for example) if these two groups are different? Not yet. Let's go figure that out.  
So we found a significant effect for our repeated measure ANOVA. We don't know where the differences are, we need to follow that up with some t-tests. 
We need to do some post hoc comparisons to say which groups are different, so that's on line 120. Let me make this a little bigger for you. It says: pairwise.t.test(x = Fake_Data_Long$Fake_Data_score, g = Fake_Data_Long$Fake_Data_group, paired = TRUE, p.adjust.method = “bonf”). What this means, is we're going to compare the scores based on the three groups. It's paired t-tests because each person has all three scores. And we're going to apply a correction; we're going to adjust our p-values because we're doing 3 different tests here. 
So I can Run this…and what this gives me is the p-values for the three tests I just ran. So if we compare Fake_Data1 to Fake_Data2, our p-value is 2.4e-14. So we take that decimal place and we move it 14 times to the left. That's statistically significant; there's a difference between Fake_Data1 and Fake_Data2. We also get a p-value for Fake_Data1 versus Fake_Data3; this says it's 2e-16, this is also significant. 
You're going to take your decimal place and move it 16 times to the left. 
There's a difference between Fake_Data1 and Fake_Data3. And then we also have Fake_Data2 and Fake_Data3; that's the same, 2e-16. We have a difference between Fake_Data2 and Fake_Data3. So if we look at our graph…this might look like a small difference, but there is a statistically significant difference between Fake_Data1 and 2…between 2 and 3…and between 1 and 3. That's what those p-values are telling you: all three groups are statistically significantly different from each other. So it's just the p-values, but it says they're all significant. And that is how you run a repeated measures ANOVA.