aov() troubles

Doing analysis of variance (ANOVA) – specifically the repeated measures kind – in R is a frustrating task that took me many hours to figure out. Here are some examples of the problem.

R has the aov() function, which can be used to perform a regular one-way ANOVA like so:

aov(myDV ~ firstGroup * secondGroup,
    data=myData)

The problems happen when you try to do a repeated measures ANOVA.

These posts explain that you can do a repeated measures ANOVA in R like so:

aov(myDV ~ firstWithinVar * secondWithinVar +
	Error(subject/(firstWithinVar * secondWithinVar)),
	data=myData)

The key difference is that you set the Error term to factor out between-subject variance.

However, I can report that this does not give you the same result as when you run the analysis in SPSS. While you will get the same Sums of Squares (SS) for each level, the Mean Squared Error (MSE) and F values are wildly different from what SPSS reports.

So wildly different that a highly significant three-way interaction that was reported as F = 24 and p < .001 in SPSS was not at all significant as reported by R’s aov() function. Knowing and visually examining the data involved leads me to easily trust SPSS’s analysis over R’s.

That is disturbing.

The posts I linked to earlier attempt various explanations as to why this happens. My understanding is that it has to do with Type I vs. Type III sums of squares calculations.

I can’t claim that I fully understand the reasons yet, but I can say that I spent many, many hours trying to get aov() to give the analysis I want and I failed.

I came across this post that explains some options you can set and incantations you can cast to get aov() to give you the analysis you want, but the instructions did not work for me. That is no criticism of the post itself - who knows why it didn’t work and what could have changed since 2008.

ezANOVA

The right and good way to perform repeated measures ANOVA in R is using the ez package, and its ezANOVA function.

In an imaginary data-frame myData, imagine I have two within-subjects variables, block and check. My dependent measure is response time (RT) measured in milliseconds (ms). You can run the repeated measures ANOVA for that model like so:

demoAnova <- ezANOVA(myData, # specify data frame
                     dv = RT, # specify dependent variable 
                     wid = subject, # specify the subject variable
                     within = .(block, check), # specify within-subject variables
                     detailed = TRUE # get a detailed table that includes SS
                     )

The detailed = TRUE line is useful for later calculating the mean squared error.

When you print the above ANOVA, the output looks something like:

$ANOVA
       Effect DFn DFd          SSn        SSd          F            p p<.05          ges
1 (Intercept)   1  24 5.450216e+07 2462038.77 531.288094 6.989868e-18     * 0.9515381648
2       block   1  24 1.287267e+06  209098.11 147.750811 9.573487e-12     * 0.3168219488
3       check   1  24 2.307161e+05   89495.63  61.871034 4.252016e-08     * 0.0767388245
4 block:check   1  24 8.199221e+02   15162.40   1.297824 2.658545e-01       0.0002952956

These results will match what you get from running a repeated measures ANOVA on the data in SPSS.

The table reports that you have a significant main effect of block, a significant main effect of check, but no interaction between the two, since the p value for the block:check term is > .05.

Mean Squared Error

The results include the F value, p value, and the degrees of freedom, which you need to report your analyses in any scientific/statistical context. However, it is also customary to report the MSE (Mean Squared Error) value. The results table doesn’t provide that, but there is a way to calculate it.

The MSE is calculated by dividing the sums of squares (SS) for the error term (denominator) of the highest order interaction by the degrees of freedom (df) of the error (denominator) for that interaction.

If you run the ANOVA without setting detailed = TRUE, the results stored in the ANOVA object do not include the SS values. Including that line gives you the SSn, which is the SS for the effect (numerator), and SSd, which is the SS for the error (denominator) for all main effects and interactions.

To calculate MSE for an ANOVA of any size:

anovaLength = length(demoAnova$ANOVA)
myMSE = demoAnova$ANOVA$SSd[anovaLength-1]/demoAnova$ANOVA$DFd[anovaLength-1]

The first line gets the number of elements in the ANOVA object inside demoAnova. The second line divides the SSd for the last element of the ANOVA object by the corresponding DFd, and stores the value in myMSE.

Update (2/2/2015)

The above MSE calculation is wrong. That is how you calculate MSE for a between-subjects design, in which the MSE you use is the one for the highest order interaction. For repeated-measures designs, each main effect or interaction gets their own MSE.

The MSE for each level is the SS of the error (SSd) divided by the DF of the error (DFd).

Therefore, we can calculate the MSE for each level like so:

demoAnova$ANOVA$MSE = demoAnova$ANOVA$SSd/demoAnova$ANOVA$DFd

This creates a new column titled MSE with the MSE values for each level.

As a bonus, you can also calculate effect sizes (\(\eta^2\)). You can calculate the effect size for each main effect and interaction by dividing the \(SS_{effect}\) by the \(SS_{total}\).

\[\eta^2 = \frac{SS_{effect}}{SS_{effect} + SS_{error}}\]

Like so:

demoAnova$ANOVA$eta <- demoAnova$ANOVA$SSn/(demoAnova$ANOVA$SSn + demoAnova$ANOVA$SSd)