In this chapter we will discuss how to conduct an Analysis of Variance (ANOVA) in R using the afex package. This chapter specifically focuses on ANOVA designs that are within subjects and mixed designs. For information about how to conduct between-subjects ANOVAs in R see Chapter 20. In this tutorial I will walk through the steps of how to run an ANOVA and the necessary follow-ups, first for a within subjects design and then a mixed design.

Before we begin, ensure that you have the necessary packages installed: (note: Use install.packages(“insert.package.name”) to install the packages if you are using them for the first time. See code below.)

#install.packages("afex")
#install.packages("emmeans")
#install.packages("ggplot2")
#install.packages("psych")
• afex: package to run ANOVAs

• emmeans: package to run follow-up comparisons

• ggplot2: package to visualize data (not required, but will be used in this tutorial)

• psych: package to get descriptive statistics that include SD and SE (not required, but will be used in this tutorial)

After you’ve installed the packages, make sure you also load them into R using the library command:

library(afex)
library(emmeans)
library(ggplot2)
library(psych)

# 1 Within Subjects ANOVAs

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.

NOTE: This section of the tutorial will be using the data file: Chapter_21_Within_Data.csv. If you want to follow along with this tutorial using this dataset, you will first need to download this file to your working directory (see step 1 below for how to set and find your working directory.

## 1.1 Data Background:

This fake data is a sample of 30 adults. All subjects were shown a set of words. Afterwards, for some of the words they were given a multiple choice test, and the rest of the words they simply restudied. Immediately afterwards their memory was tested for all words to determine if there were memory differences depending on whether the word was tested or restudied (immediate final test). Two days later, subjects came back and were given the same test on all words again (delayed final test). The dependent variable (DV) given in the file are a proportion of correctly remembered items for the immediate and delayed test for both conditions. Therefore, if we want to know if there are memory difference based on time delay and whether the word was tested or restudied, we need to conduct a within-subjects ANOVA.

The variables given in the data set:

Subject = Subject ID #

Within_Cond = Study Method (test or restudy)

Within_Time = Immediate or Delayed

DV = Memory Performance

Below we will walk through the steps to find our answer using the afex package in R.

## 1.2 Steps for Within Subject ANOVA

### 1.2.1 Set Working Directory

#setwd("path") - will set your working directory to the file path you specify
#getwd() - will tell you what your working directory currently is

### 1.2.2 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.

### 1.2.3 Load Data 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:

2. the separater [*note: by saving file as a .csv the separater by default will be a comma (“,”)].

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

### 1.2.4 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.

### 1.2.5 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. ### 1.2.6 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. ### 1.2.7 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 vars n mean sd median trimmed mad min max range skew kurtosis se 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 ### 1.2.8 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! ### 1.2.9 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: 1. A formula specifying the ANOVA Model (including specifying your error term). 2. 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)) Effect df MSE F ges p.value 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. ### 1.2.10 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. ### 1.2.11 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. ### 1.2.12 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)) Effect df MSE F ges p.value 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. ### 1.2.13 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)) Effect df MSE F pes p.value 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 ### 1.2.14 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)) Effect df MSE F ges p.value 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. ### 1.2.15 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 ### 1.2.16 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”. ### 1.2.17 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: #### 1.2.17.1 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. ### 1.2.18 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. ### 1.2.19 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. #### 1.2.19.1 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. #### 1.2.19.2 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. ### 1.2.20 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. # 2 Mixed Design ANOVAs (afex) Mixed designs are research designs in which there are both between subjects and within subjects variables being compared. In this section of the tutorial we are going to use the afex package to analyze a new data set that is a very simple 2 by 2 mixed design. ## 2.1 Data Background: This fake data is a sample of 15 younger adults and 15 older adults. Similar to our previous example, all subjects were initially shown a set of words. Afterwards, for some of the words they were given a multiple choice test, and the rest of the words they simply restudied. After a 2-day delay, their memory was tested for all words to determine if there were memory differences depending on whether the word was tested or restudied. The dependent variable (DV) given in the file are a proportion of correctly remembered items for the delayed memory test for both conditions. In this example, we want to know if there are memory differences between age group (our between subjects factor) and study method (our within subjects factor). We need to run an ANOVA to find our answer. The variables given in the data set: Subject = Subject ID # Within_Cond = Study Method (test or restudy) Btw_Cond = Age Group (Younger adult or Older adult) DV = Memory Performance Again, we will walk through the steps to find our answer using the afex package in R. NOTE: This section of the tutorial will be using the data file: Chapter_21_Mixed_Data.csv. If you want to follow along with this tutorial using this dataset, you will first need to download this file to your working directory. ## 2.2 Steps for Mixed Design ANOVAs ### 2.2.1 Set Working Directory. Just as a reminder, be sure to set your working directory to the file path where your data is saved (as a .csv file) before you can read the data into R. ### 2.2.2 Read in Data. First let’s read in our data, and call the head of the data to make sure R read in the data. Mixed_Data <- read.table("Chapter_21_Mixed_Data.csv", sep=",", header=T) head(Mixed_Data) ## Subject Btw_Cond Within_Cond DV ## 1 1 0 1 93 ## 2 1 0 2 66 ## 3 2 0 1 62 ## 4 2 0 2 55 ## 5 3 0 1 88 ## 6 3 0 2 64 ### 2.2.3 Check Structure of Data for Errors. Similar to our previous example, we want to check the structure of our dataset as R may have read our data differently than we wanted. str(Mixed_Data) ## 'data.frame': 60 obs. of 4 variables: ##$ Subject    : int  1 1 2 2 3 3 4 4 5 5 ...
##  $Btw_Cond : int 0 0 0 0 0 0 0 0 0 0 ... ##$ Within_Cond: int  1 2 1 2 1 2 1 2 1 2 ...
##  $DV : num 93 66 62 55 88 64 84 66 95 66 ... ### 2.2.4 Correct Structure of Dataset (if necessary). Here, again, we see that R has our subject IDs listed as integers. Let’s change them to factors and recheck the structure. Mixed_Data$Subject<-as.factor(Mixed_Data$Subject) str(Mixed_Data) ## 'data.frame': 60 obs. of 4 variables: ##$ Subject    : Factor w/ 30 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $Btw_Cond : int 0 0 0 0 0 0 0 0 0 0 ... ##$ Within_Cond: int  1 2 1 2 1 2 1 2 1 2 ...
##  $DV : num 93 66 62 55 88 64 84 66 95 66 ... ### 2.2.5 Name Variable(s). Now let’s add variable names to our between subject and within subject conditions so we can more easily interpret our data later. Mixed_Data$Btw_Cond <- factor(Mixed_Data$Btw_Cond, levels = c(0,1), labels = c("Younger_Adult", "Older_Adult")) Mixed_Data$Within_Cond <- factor(Mixed_Data$Within_Cond, levels = c(1,2), labels = c("Test", "Restudy")) And do one final check of our data. str(Mixed_Data) ## 'data.frame': 60 obs. of 4 variables: ##$ Subject    : Factor w/ 30 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $Btw_Cond : Factor w/ 2 levels "Younger_Adult",..: 1 1 1 1 1 1 1 1 1 1 ... ##$ Within_Cond: Factor w/ 2 levels "Test","Restudy": 1 2 1 2 1 2 1 2 1 2 ...
##  \$ DV         : num  93 66 62 55 88 64 84 66 95 66 ...

Here we can check that all of our changes were actually made and there were no errors.

### 2.2.6 Get Descriptives.

Next, we want to get descriptive statistics for our variables. Let’s use the psych package and describe command to create a descriptives table.

Mixed_Means.Table<-describe(Mixed_Data)

Mixed_Means.Table #call the object to see the means table
##              vars  n  mean    sd median trimmed   mad min max range skew
## Subject*        1 60 15.50  8.73   15.5   15.50 11.12   1  30    29 0.00
## Btw_Cond*       2 60  1.50  0.50    1.5    1.50  0.74   1   2     1 0.00
## Within_Cond*    3 60  1.50  0.50    1.5    1.50  0.74   1   2     1 0.00
## DV              4 60 64.73 14.31   64.5   64.23 17.79  41  95    54 0.27
##              kurtosis   se
## Subject*        -1.26 1.13
## Btw_Cond*       -2.03 0.07
## Within_Cond*    -2.03 0.07
## DV              -0.84 1.85
#use knitr function to call the object to see the means table.
knitr::kable(Mixed_Means.Table, digits = 2, caption = 'Mixed Means Table')
Mixed Means Table
vars n mean sd median trimmed mad min max range skew kurtosis se
Subject* 1 60 15.50 8.73 15.5 15.50 11.12 1 30 29 0.00 -1.26 1.13
Btw_Cond* 2 60 1.50 0.50 1.5 1.50 0.74 1 2 1 0.00 -2.03 0.07
Within_Cond* 3 60 1.50 0.50 1.5 1.50 0.74 1 2 1 0.00 -2.03 0.07
DV 4 60 64.73 14.31 64.5 64.23 17.79 41 95 54 0.27 -0.84 1.85

### 2.2.7 Plot Data (ggplot2).

Let’s make a plot of our data to help visualize. Using the ggplot command (see Ch. 8 & 9) we can make a bar graph of memory performance by study method, faceted by age group.

Mixed_Data.BarGraph<-ggplot(Mixed_Data, aes(Within_Cond, DV, fill=Within_Cond)) +
geom_bar(stat="summary", fun.y="mean", position = "dodge") +
scale_y_continuous(breaks = seq(0, 101, 10), limits =c(0,101)) +
facet_grid(.~Btw_Cond) +
xlab("Study Method") + ylab("Memory Performance") +
scale_fill_brewer(palette="Dark2") +
theme(legend.position="none")

Mixed_Data.BarGraph

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).

It again looks like we have some interesting effects here. Let’s run our ANOVA to find out if they are significant.

### 2.2.8 Run ANOVA.

Similar to our within subject example, we are going to use the aov_car command to run our ANOVA. But because we are using a mixed design, we will need to make some changes to our formula. Let’s first review what the aov_car command requires:

1. A formula specifying the ANOVA Model (including specifying your error term).

2. The data to analyze.

Below is an example using our “Mixed_Data” dataset. After the aov_car command, we enter the formula for our ANOVA model. In this case, we are looking at DV by (“~”) Btw_Cond and Within_Cond. Also we know that we must specify the error term we wish to use. In this example, the error term is represented by subject error divided by the within-subjects error (Within_Cond). IMPORTANTLY: If you have any between subjects variables, they do not need to be included in the error term. And finally, don’t forget to specify which dataset you are running the ANOVA on (Mixed_Data).

Let’s try it.

Mixed.aov.1<-aov_car(DV ~ Btw_Cond*Within_Cond + Error(Subject/Within_Cond), data=Mixed_Data)
## Contrasts set to contr.sum for the following variables: Btw_Cond
knitr::kable(nice(Mixed.aov.1)) ##call for formatted ANOVA table using knitr
Effect df MSE F ges p.value
Btw_Cond 1, 28 119.81 11.86 ** .17 .002
Within_Cond 1, 28 128.83 18.49 *** .25 .0002
Btw_Cond:Within_Cond 1, 28 128.83 10.29 ** .16 .003

Here we see our familar ANOVA table, that shows a main effect of for study method and age group, and an interaction!

Remember that the aov_car command automatically corrects for sphericity violations, and reports Generalized Eta Square. See within subjects example above for customizations of this analysis and output table if you want something other than these default corrections and output. Or call the aov_car help (?aov_car).

### 2.2.9 Follow-up Tests.

Now that we know we have some significant effects, let’s use the emmeans package to do some follow-up test to learn more.

First we need to name our object that will hold the marginal means table. Then we call the emmeans function and specify:

1.) Our ANOVA object we created when we ran the ANOVA using afex (Mixed.aov.1).

2.) The variable we want to see the marginal means for.

Let’s look at the marginal means for Study Method and then Age.

##Main effect of StudyMethod
Mixed_Fitted_StudyMethod<-emmeans(Mixed.aov.1, ~Within_Cond)
## NOTE: Results may be misleading due to involvement in interactions
Mixed_Fitted_StudyMethod
##  Within_Cond   emmean       SE    df lower.CL upper.CL
##  Test        71.03333 2.035674 55.93 66.95527  75.1114
##  Restudy     58.43333 2.035674 55.93 54.35527  62.5114
##
## Results are averaged over the levels of: Btw_Cond
## Confidence level used: 0.95
##Main effect of Age

Mixed_Fitted_AgeGroup<-emmeans(Mixed.aov.1, ~Btw_Cond)
## NOTE: Results may be misleading due to involvement in interactions
Mixed_Fitted_AgeGroup
##  Btw_Cond        emmean       SE df lower.CL upper.CL
##  Younger_Adult 69.60000 1.998412 28 65.50644 73.69356
##  Older_Adult   59.86667 1.998412 28 55.77311 63.96023
##
## Results are averaged over the levels of: Within_Cond
## Confidence level used: 0.95

Since these variables only have two levels, we know from our ANOVA that they are significant differences. It seems that testing improves memory over restudying, and that younger adults performed better than older adults.

However, there is an interaction. Let’s look at our marginal means for the interaction, and do a pairwise comparison to help interpret.

##Study Method by Age Interaction

Mixed_Fitted_Interaction<-emmeans(Mixed.aov.1, ~Within_Cond|Btw_Cond)
Mixed_Fitted_Interaction
## Btw_Cond = Younger_Adult:
##  Within_Cond   emmean       SE    df lower.CL upper.CL
##  Test        80.60000 2.878878 55.93 74.83275 86.36725
##  Restudy     58.60000 2.878878 55.93 52.83275 64.36725
##
##  Within_Cond   emmean       SE    df lower.CL upper.CL
##  Test        61.46667 2.878878 55.93 55.69941 67.23392
##  Restudy     58.26667 2.878878 55.93 52.49941 64.03392
##
## Confidence level used: 0.95
##pairwise comparison with no correction.

pairs(Mixed_Fitted_Interaction)
## Btw_Cond = Younger_Adult:
##  contrast       estimate       SE df t.ratio p.value
##  Test - Restudy     22.0 4.144532 28   5.308  <.0001
##
##  contrast       estimate       SE df t.ratio p.value
##  Test - Restudy      3.2 4.144532 28   0.772  0.4465

Great! Now we can see that for younger adults, testing improved memory on a delayed test more so than restudying. However, for older adults, memory was equal for words that were tested and words that were restudied. Interesting.

### 2.2.10 Interpret Results.

Now that we’ve done our follow-up tests we can interpret our results.

In this example, we tested whether self-testing or restudying previously learned words would improve memory on a 2-day delayed test for younger and older adults. Our ANOVA showed an interaction, which our Bonferroni corrected follow-up tests showed that testing improves memory over restudying on a delayed final test, but only for younger adults. For older adults, memory was equal for tested items and restudied items.

If you’d like, you can recall our graph to see if this graph makes more sense now that we have our results.

Mixed_Data.BarGraph

# 3 Additional Resources for ANOVAs in R.

This tutorial outlines the basic steps to running a simple 2x2 ANOVA for within subjects and mixed designs. For more complex examples and tutorials please utilize the following resources:

-A overview of the theory and code to do RM and mixed ANOVAs with suggestions on how to do follow ups

-A tutorial on using afex aov_ez function as well as different constrasts/follow-ups with the emmeans package.

-A tutorial for using the emmeans package that explains contrasts.

# 4 References

Hall, R (1998), Psychology World https://web.mst.edu/~psyworld/within_subjects.htm.