why does the same repeated measures anova using ezANOVA() vs. aov() yield different distributions of model residuals?

Question

I am attempting to do a repeated measures anova using r with the aov() command from the {car} package. I wanted to be sure that I wrote my code for this approach correctly (see below), so I cross-checked the results against same repeated measures anova model using the ezANOVA() command from the {ez} package.

What I found was that the results from each anova were the same, but the residuals of each model using the respective packages were distributed differently. Using the aov() command, the residuals were non-normally distributed whereas using ezANOVA() they were normally distributed.

###############################
#aov() using the {car} package#
###############################

mod.aov <-aov(rating~(factor1*factor2) + Error(id/(factor1*factor2)), data=datalong)

shapiro.test(residuals(mod.aov))

Shapiro-Wilk normality test
data:  residuals(mod.aov)
W = 0.9961, p-value = 0.036 

pearson.test(residuals(mod.aov))

Pearson chi-square normality test
data:  residuals(gain.aov)
P = 50.4, p-value = 0.004086

##################################
#ezANOVA() using the {ez} package#
##################################

mod.ez <- ezANOVA(data = datalong,
                  dv = .(rating), 
                  wid = .(id),
                  within = .(factor1, factor2),
                  detailed = TRUE,
                  return_aov = TRUE)

shapiro.test(residuals(mod.ez$aov))

Shapiro-Wilk normality test
data:  residuals(mod.ez$aov)
W = 0.9876, p-value = 0.449

pearson.test(residuals(mod.aov$aov))

Pearson chi-square normality test
data:  residuals(mod.ez$aov)
P = 7.75, p-value = 0.6532

I do not understand whether the residuals of each model are different because of my code, or because of the differences between way the packages perform the analyses. Can anybody help? I can provide additional information if needed. Thanks.

Answer 1

First, there is no aov function in the car package, so I'm guessing you are referring to aov from the stats package.

Second, by comparing the QQ-plots of the residuals obtained via aov and ezANOVA, you are apparently not using the same data - the aov plot shows many more data points than the ezANOVA plot.

Third, I suspect that since you're not using identical data, this is the reason why you get different residuals. I just compared the residuals with a toy example, and both aov and ezANOVA yield identical results:

library(ez)

rating <- c(8, 9, 6, 5, 8, 7, 10, 12, 7, 5, 2, 3, 4, 5, 2, 6, 1, 2, 3, 1, 5, 6, 7, 8, 6, 5, 8, 9, 8, 7, 2, 1)
group <- factor(rep(1:4, each=8, len=32))
id <- factor(rep(1:8, times=4))
df <- data.frame(id, group, rating)

mod.ez <- ezANOVA(df, rating, id, within=group, return_aov=T)
mod.aov <- aov(rating ~ group + Error(id/group), df)

res.ez <- sort(proj(mod.ez$aov)[[3]][, "Residuals"])
res.aov <- sort(proj(mod.aov)[[3]][, "Residuals"])

res.ez - res.aov

plot(res.ez, res.aov)

So except for numerical errors, the residuals from the two methods are identical.