Within subjects designs occur when each subject is exposed to both the manipulation and control conditions. This design has strengths and weaknesses compared to a between subjects design that are discussed in more detail here: https://web.mst.edu/~psyworld/within_subjects.htm. For this tutorial, we will be using fake data sets to show how to run a within subjects ANOVA in R.
Steps for Within Subject ANOVA
Set Working Directory
Before loading your dataset, you will need to set your working directory to where your dataset is saved. Set your working directory using the setwd command:
#setwd("path") - will set your working directory to the file path you specify
#getwd() - will tell you what your working directory currently is
Save Dataset as a .csv File.
If your file is in a different format, you will need to save it as a comma-separated values (CSV) file before loading the data into R.
Load Data into R.
Next, load your data set into R.
For this you can use the read.table command. This command requires us to name our data as a variable. In this example we will call our dataset Within_Data.
The read.table command requires three inputs:
- name of the file you are uploading in quotes "".
- the separater [*note: by saving file as a .csv the separater by default will be a comma (“,”)].
- Does your data file have headers? (type “T” if your data file includes headers, or “F” if your data file does not include headers.).
Here’s our code for loading our example data.
Within_Data = read.table("Chapter_21_Within_Data.csv", sep=",", header=T) #loads dataset into object called "Within_Data"
If using R-studio, you can confirm the data was loaded by seeing it appear in your “global environment”. You can also confirm the data was loaded correctly by calling the dataset by the variable name you gave it (in our case “Within_Data”). Here, for brevity I’ve used the head and tail commands to only show the first 5 rows and last 5 rows of our data. To show the entire dataset, simply type the name of the dataset “Within_Data” without the head or tail command.
head(Within_Data, n=5) #shows first 5 rows of our data set
## Subject Within_Cond Within_Time DV
## 1 1 1 1 96
## 2 2 1 1 81
## 3 3 1 1 73
## 4 4 1 1 90
## 5 5 1 1 72
tail(Within_Data, n=5) #shows last 5 rows of our data set
## Subject Within_Cond Within_Time DV
## 116 26 2 2 53
## 117 27 2 2 74
## 118 28 2 2 77
## 119 29 2 2 82
## 120 30 2 2 52
Check the Structure of the Dataset.
Next we want to check the structure of the dataset to ensure that R read our data correctly. Use the str command to look at the structure of your data.
str(Within_Data)
## 'data.frame': 120 obs. of 4 variables:
## $ Subject : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Within_Cond: int 1 1 1 1 1 1 1 1 1 1 ...
## $ Within_Time: int 1 1 1 1 1 1 1 1 1 1 ...
## $ DV : num 96 81 73 90 72 83 83 96 86 72 ...
By calling the structure of our dataset, we can see that our “subject” numbers were coded as integers but they are actually factors. Be sure to check for any other discrepancies in your dataset, between what you want and how R reads the data, that are specific to your dataset.
Correct Structure of Dataset (if necessary).
In this example we see that our subject IDs are listed as integers. We need to manually tell R that these are factors. We do this using the as.factor command. (note: “$” specifies where you want a command to act on your dataset, in this case we want to specify Subject) See Chapter 2 on indexing for more info.
Within_Data$Subject<-as.factor(Within_Data$Subject) #Changes variable "Subject" from an integer to a factor.
Do a final check of the structure of the dataset after making changes to see if it has done what we expect.
str(Within_Data)
## 'data.frame': 120 obs. of 4 variables:
## $ Subject : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Within_Cond: int 1 1 1 1 1 1 1 1 1 1 ...
## $ Within_Time: int 1 1 1 1 1 1 1 1 1 1 ...
## $ DV : num 96 81 73 90 72 83 83 96 86 72 ...
Now we see that R recognizes our subject ID variable as a factor like we want.
Label Independent Variable(s).
Next, we will label our independent variables so we know what condition and time point each memory score corresponds to later on. In our example, our independent variables are Within_Cond (test v. restudy) and Within_Time (immediate v. delayed) In our data they are listed as 1’s and 2’s. Let’s make these columnns include the labels for both variables.
Within_Data$Within_Cond <- factor(Within_Data$Within_Cond,
levels = c(1,2),
labels = c("Test", "Restudy")) #converts 1's and 2's in the data file to correspond to factor labels of the Within Subjects variable of study method.
Within_Data$Within_Time <-factor(Within_Data$Within_Time,
levels = c(1,2),
labels = c("Immediate", "Delayed")) #converts 1's and 2's in the data file to correspond to factor labels of the Within Subjects variable of time.
Now lets look at the head and structure of our data again.
head(Within_Data, n=5)
## Subject Within_Cond Within_Time DV
## 1 1 Test Immediate 96
## 2 2 Test Immediate 81
## 3 3 Test Immediate 73
## 4 4 Test Immediate 90
## 5 5 Test Immediate 72
str(Within_Data)
## 'data.frame': 120 obs. of 4 variables:
## $ Subject : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Within_Cond: Factor w/ 2 levels "Test","Restudy": 1 1 1 1 1 1 1 1 1 1 ...
## $ Within_Time: Factor w/ 2 levels "Immediate","Delayed": 1 1 1 1 1 1 1 1 1 1 ...
## $ DV : num 96 81 73 90 72 83 83 96 86 72 ...
Now we have labels for our IVs.
Check Descpritives.
Before we analyze our data, we should check the descriptive statistics of our variables. Let’s use the psych package and describe command to look at a table of descriptive statistics of our variables. Be sure to place the table into an object so we can view it.
#use psych package to call for descriptives of your dataset and place into an object.
Within_Means.Table<-describe(Within_Data)
Within_Means.Table #call the object to see the means table
## vars n mean sd median trimmed mad min max range skew
## Subject* 1 120 15.50 8.69 15.5 15.50 11.12 1 30 29 0.00
## Within_Cond* 2 120 1.50 0.50 1.5 1.50 0.74 1 2 1 0.00
## Within_Time* 3 120 1.50 0.50 1.5 1.50 0.74 1 2 1 0.00
## DV 4 120 79.19 11.77 81.0 80.31 10.38 47 99 52 -0.82
## kurtosis se
## Subject* -1.23 0.79
## Within_Cond* -2.02 0.05
## Within_Time* -2.02 0.05
## DV 0.76 1.07
We can use the knitr function to create a table that is formatted and easier to read. More on this later when we run the ANOVA.
#use knitr function to call the object to see the means table.
knitr::kable(Within_Means.Table, digits = 2, caption = 'Within Means Table') #Here we've also specified that we want to round to "2 digits", and we want our table to have a title of "Within Means Table". These are not necessary, but can be included.
Within Means Table
| Subject* |
1 |
120 |
15.50 |
8.69 |
15.5 |
15.50 |
11.12 |
1 |
30 |
29 |
0.00 |
-1.23 |
0.79 |
| Within_Cond* |
2 |
120 |
1.50 |
0.50 |
1.5 |
1.50 |
0.74 |
1 |
2 |
1 |
0.00 |
-2.02 |
0.05 |
| Within_Time* |
3 |
120 |
1.50 |
0.50 |
1.5 |
1.50 |
0.74 |
1 |
2 |
1 |
0.00 |
-2.02 |
0.05 |
| DV |
4 |
120 |
79.19 |
11.77 |
81.0 |
80.31 |
10.38 |
47 |
99 |
52 |
-0.82 |
0.76 |
1.07 |
Plot Data Using ggplot2 to Visualize.
Next, it can be useful to plot your data to see what our data look like visually. This can help in interpretation.
Here, we call for a simple bar graph of memory performance by study method, faceted by time. For more info on how to make and customize bar graphs using ggplot2 see Chapters 10 & 11.
##GGplot command is placed into an object "Within_Data.BarGraph" so we can view the graph after the command runs.
Within_Data.BarGraph<-ggplot(Within_Data, aes(Within_Cond, DV, fill=Within_Cond)) +
geom_bar(stat="summary", fun.y="mean") +
scale_y_continuous(breaks = seq(0, 101, 10), limits =c(0,101)) +
facet_grid(.~Within_Time) +
xlab("Study Method") + ylab("Memory Performance") +
scale_fill_brewer(palette="Dark2") +
theme(legend.position="none")
Within_Data.BarGraph ### call for object to see plot
NOTE: This plot does not include error bars. To include error bars that are corrected for within subjects design, doing so requires a few more steps. See the following for a step-by-step guide if you wish to include error bars on your plot: http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2).
Here we can see that we likely have some effects. Let’s run our ANOVA to find out if they’re significant!
Run ANOVA.
Now that we’ve plotted our data and have it organized how we want, let’s run our ANOVA to test for differences in memory performance between our condition and time. There is a built-in R function to run ANOVAs (aov). For more info on this built in function see Chapter 20 or type “?aov” in the R console.
For this tutorial we are using the AFEX package. The call for the ANOVA function we will use in afex is “aov_car” (aov_ez & aov4 are alternatives that can be used). We can get more info on the aov_car function (and alternatives) by typing “?aov_car” in the console.
?aov_car
## starting httpd help server ... done
To run the ANOVA, aov_car requires two things:
A formula specifying the ANOVA Model (including specifying your error term).
The data to analyze.
Below is an example using our dataset. After the aov_car command, we enter the formula for our ANOVA model. In this case, we are looking at DV by (“~”) Within_Cond and Within_Time. Also in your formula you must specify the error term you wish to use. In our example, the error term is represented by subject error divided by the within-subjects error (Within_Cond x Within_Time). If you have any between subjects variables, they do not need to be included in the error term (as we will see when we get to our Mixed Design example later). Lastly, we specify the data we want to use as Within_Data. Also, note that in the output given by R, the system uses the error term you specify to calculate these values, this is different than SPSS).
It should also be noted that in order to see the results of the ANOVA, you must place your aov_car command into an object. Then call that object to see the results of the ANOVA. Here we’ve named our object Within.aov.1.
Let’s try it.
Within.aov.1 <- aov_car(DV ~ Within_Cond*Within_Time + Error(Subject/Within_Cond*Within_Time), data=Within_Data)
Within.aov.1 #call object to see ANOVA table
## Anova Table (Type 3 tests)
##
## Response: DV
## Effect df MSE F ges p.value
## 1 Within_Cond 1, 29 109.51 12.24 ** .10 .002
## 2 Within_Time 1, 29 89.63 13.21 ** .09 .001
## 3 Within_Cond:Within_Time 1, 29 115.65 20.17 *** .17 .0001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
Interesting. Here we see that our data shows two main effects and an interaction!
However, this table might be hard to read, and isn’t great for publishing. Let’s try calling for a formatted ANOVA table using the knitr command.
knitr::kable(nice(Within.aov.1))
| Within_Cond |
1, 29 |
109.51 |
12.24 ** |
.10 |
.002 |
| Within_Time |
1, 29 |
89.63 |
13.21 ** |
.09 |
.001 |
| Within_Cond:Within_Time |
1, 29 |
115.65 |
20.17 *** |
.17 |
.0001 |
That looks better and easier to read. We will use the knitr function throughout the rest of the tutorial when we call our ANOVA table.
AFEX Defaults
Before we move to follow-up test, there are some things we should note about the aov_car function in the afex package.
1.) First we should note, that unlike the built in R function (aov), the afex package defaults to using Type III Sum of Squares.
2.) The afex package also automatically detects and corrects for violations of sphericity in your data. If a correction is made to your DF, R will report which correction was used below the ANOVA table. NOTE: that in a 2x2 design there is no sphericity, therefore these corrections in our example have no effect
3.) Additionally, when reporting effect sizes, the afex package defaults to reporting Generalized Eta Squared (GES). This is in contrast to Partial Eta Squared (PES) that is also commonly reported in some fields. For more information on GES v. PES, see Olejnik & Algina (2003); Bakeman (2005) to help decide which one is best for your analysis and reporting.
AFEX Customizations
Although these are the defaults for the afex package. There are some commands you can use to customize your analysis and output to meet your specifications. Let’s look at some common customizations that might be useful in conducting your ANOVA.
When you have an observed IV.
If you have an observed IV, include the observed line and specify which variable you are designating as observed to have afex correct the effect size.
Within.aov.2<-aov_car(DV ~ Within_Cond*Within_Time + Error(Subject/Within_Cond*Within_Time), data=Within_Data,
observed = "Within_Time") #this corrects GES if a variable is observed and not manipulated
knitr::kable(nice(Within.aov.2))
| Within_Cond |
1, 29 |
109.51 |
12.24 ** |
.08 |
.002 |
| Within_Time |
1, 29 |
89.63 |
13.21 ** |
.08 |
.001 |
| Within_Cond:Within_Time |
1, 29 |
115.65 |
20.17 *** |
.15 |
.0001 |
| *Note that our Generalize |
d Eta Sq |
uared val |
ues are diff |
erent |
than our original ANOVA table above after making this correction. |
Reporting Partial Eta Squared instead of Generalized Eta Squared.
Include the anova_table line and specify es = “pes” to have your ANOVA report PES instead of the default (GES).
Within.aov.3<-aov_car(DV ~ Within_Cond*Within_Time + Error(Subject/Within_Cond*Within_Time), data=Within_Data,
anova_table = list(es = "pes")) #This reports PES not GES
knitr::kable(nice(Within.aov.3))
| Within_Cond |
1, 29 |
109.51 |
12.24 ** |
.30 |
.002 |
| Within_Time |
1, 29 |
89.63 |
13.21 ** |
.31 |
.001 |
| Within_Cond:Within_Time |
1, 29 |
115.65 |
20.17 *** |
.41 |
.0001 |
Preventing Sphericity Corrections.
The afex package will check your data for sphericity assumptions and if necessary correct for any violations. However, you may want to turn off these corrections in some instances. Include the following anova_table line and specify correction = “none” to prevent these default corrections.
Within.aov.4<-aov_car(DV ~ Within_Cond*Within_Time + Error(Subject/Within_Cond*Within_Time), data=Within_Data,
anova_table = list(correction = "none")) #This turns off GG correction (sphericity corrections disabled)
knitr::kable(nice(Within.aov.4))
| Within_Cond |
1, 29 |
109.51 |
12.24 ** |
.10 |
.002 |
| Within_Time |
1, 29 |
89.63 |
13.21 ** |
.09 |
.001 |
| Within_Cond:Within_Time |
1, 29 |
115.65 |
20.17 *** |
.17 |
.0001 |
NOTE: in a 2x2 design there is no sphericity, therefore these corrections in our example have no effect, but if using a more complex design, this will remove any corrections that were made.
Viewing Univariate Output.
Sometimes it can be useful to see the univariate output, that is not given by default in the afex package. However, we can call for this output to be generataed with our ANOVA table by including the return command and specifying “univariate”.
Within.aov.5<-aov_car(DV ~ Within_Cond*Within_Time + Error(Subject/Within_Cond*Within_Time), data=Within_Data,
return="univariate") ##Changes the output to Univariate
Within.aov.5 #since we've requested a specific output, we will just call the object, instead of using knitr for this ANOVA table.
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 752558 1 2497.3 29 8738.970 < 2.2e-16
## Within_Cond 1340 1 3175.7 29 12.237 0.0015321
## Within_Time 1184 1 2599.3 29 13.214 0.0010667
## Within_Cond:Within_Time 2332 1 3353.7 29 20.165 0.0001042
##
## (Intercept) ***
## Within_Cond **
## Within_Time **
## Within_Cond:Within_Time ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Additional Customization.
These are just a few examples of customizations to the output given by the afex package. For a full list and more info on how to use other customizations, utilize the R help file on the aov_car command by typing “?aov_car”.
Follow-up Tests (emmeans).
Now that we know we have some significant effects, we should follow up these effects with pairwise comparisons or contrasts. We will use the emmeans package to conduct our follow-ups. Using emmeans we will need:
1.) The name of the object our ANOVA table is saved as.
2.) The name of the variable we want to compare.
Let’s start with the main effects:
Main effect of Study Method
First we create an object name, then we call the emmeans function and include our ANOVA table object from before (Within.aov.1) and then specify that we want the marginal means of the Study Method (~Within_Cond). Then we call the object we created to see the means table.
Within_Fitted_StudyMethod<-emmeans(Within.aov.1, ~Within_Cond)
## NOTE: Results may be misleading due to involvement in interactions
Within_Fitted_StudyMethod #call object to see marginal means table.
## Within_Cond emmean SE df lower.CL upper.CL
## Test 82.53333 1.276791 57.18 79.97678 85.08989
## Restudy 75.85000 1.276791 57.18 73.29345 78.40655
##
## Results are averaged over the levels of: Within_Time
## Confidence level used: 0.95
Notice that emmeans always gives the warning that we should take into account any interactions when interpreting these data. Which in our case is applicable.
Also notice our funny degrees freedom in this table. emmeans automatically corrects for certain statistical violations that are relatively opaque to most users. For more information on these corrections and whether emmeans is correct for your data, see the “introduction” section of the emmeans vignette: https://cran.r-project.org/web/packages/emmeans/vignettes/using-emmeans.pdf
Now, let’s quickly see the means for our time variable before moving to our interaction.
Main Effect of Time
Within_Fitted_TimePoint<-emmeans(Within.aov.1, ~Within_Time)
## NOTE: Results may be misleading due to involvement in interactions
Within_Fitted_TimePoint
## Within_Time emmean SE df lower.CL upper.CL
## Immediate 82.33333 1.210192 57.98 79.91085 84.75582
## Delayed 76.05000 1.210192 57.98 73.62752 78.47248
##
## Results are averaged over the levels of: Within_Cond
## Confidence level used: 0.95
Since there are only two levels for each of our variables, our ANOVA tells us these means are significantly different in both conditions. Let’s look at our interaction to see an example of how to do pairwise comparisons if you’re comparing more than 2 levels.
Interaction
Most importantly, our ANOVA showed an interaction between study method and time. Let’s use pairwise comparisons to interpret our interaction. Again, we start by getting the marginal means, this time for both conditions by calling emmeans and specifying we want means for our study method condition (Within_cond) by time (Within_Time). emmeans uses the “|” symbol to look at an interaction between two variables.
Within_Fitted_Interaction <- emmeans(Within.aov.1, ~ Within_Cond|Within_Time)
Within_Fitted_Interaction #call the object we created
## Within_Time = Immediate:
## Within_Cond emmean SE df lower.CL upper.CL
## Test 81.26667 1.8278 114.19 77.64588 84.88746
## Restudy 83.40000 1.8278 114.19 79.77921 87.02079
##
## Within_Time = Delayed:
## Within_Cond emmean SE df lower.CL upper.CL
## Test 83.80000 1.8278 114.19 80.17921 87.42079
## Restudy 68.30000 1.8278 114.19 64.67921 71.92079
##
## Confidence level used: 0.95
Now, to test the significance of the mean differences, we will use the pairs function. For this, we simply call the pairs function and then specify the object that holds our marginal means. In this case (Within_Fitted_Interaction)
pairs(Within_Fitted_Interaction) ## pairwise comparison with no correction
## Within_Time = Immediate:
## contrast estimate SE df t.ratio p.value
## Test - Restudy -2.133333 2.739553 57.96 -0.779 0.4393
##
## Within_Time = Delayed:
## contrast estimate SE df t.ratio p.value
## Test - Restudy 15.500000 2.739553 57.96 5.658 <.0001
Interesting. From our pairwise comparisons, we can see that in the immediate test there was no difference between study methods, but for the delayed test, the multiple choice test provided better memory performance than restudying the words.
A Note About Follow-ups
If your follow-up tests are not a priori, you may want to correct your significance values using Tukey’s, Bonferroni, etc. emmeans allows you to change the correction method used by simply adding the adjust command and then specifying the correction you want. Let’s try a bonferroni correction.
pairs(Within_Fitted_Interaction, adjust="bon") ## pairwise comparison with Bonferroni correction
## Within_Time = Immediate:
## contrast estimate SE df t.ratio p.value
## Test - Restudy -2.133333 2.739553 57.96 -0.779 0.4393
##
## Within_Time = Delayed:
## contrast estimate SE df t.ratio p.value
## Test - Restudy 15.500000 2.739553 57.96 5.658 <.0001
Note that in a 2x2 design, these corrections make no changes to our results as we only have one test “per family”. But if using a more complex design, this is how you would correct for unplanned follow-up tests.
A Final Note:
emmeans also allows you to create your own contrasts (dummy coding, effect coding, etc.). For more information and examples for how to create your own contrasts, as well as other post hoc tests in R using afex, see: https://cran.r-project.org/web/packages/afex/vignettes/anova_posthoc.html. This site also provides examples of more complex designs with greater than 2 levels per variable.
Interpret Results.
Lastly, now that we’ve ran our ANOVA and follow-ups, we can interpret our results.
In our example, we tested whether self-testing or restudying previously learned words would improve memory on an immediate and delayed test. Our ANOVA showed an interaction, which our Bonferroni corrected follow-up tests showed that testing improves memory over restudying only on a delayed test, not on an immediate test.
Interesting. Now, let’s turn to a similar example, but using a mixed design.
