# ANOVA in R

## 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:

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:

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:

The `detailed = TRUE`

line is useful for later calculating the mean squared error.

When you `print`

the above ANOVA, the output looks something like:

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:

~~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:

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: