Mixed ANOVA normality: which variables should be examined? (in universal and practical application with stats::aov)

Question

I have spent a lot of time reading book chapters, articles, online tutorials, etc., but with no clear answer (mostly because they only describe one-way ANOVA or other very specific applications). There have also been many similar questions on this site, but again no satisfactory answer for my purposes.

In essence, I'd like to know the clear and straightforward (non-technical), and completely generalizable (and practically implementable) answer for how to test/examine the (in)famous ANOVA normality assumption given any number of within-subject or between-subject factors (with any number of levels).

(Note: The only question here is which variables should be examined, not how they should be examined. By "testing/examining normality", I don't necessarily mean statistical hypothesis testing, it could also be based on density or Q-Q plots, etc., doesn't matter. The only problem would be if perhaps multivariate normality testing were needed, in which case again the question would be which variables should be included in it.)

At least this tutorial and this answer advises to examine the normality of every single cell, i.e. every possible combination of each level of each factor – but no references or detailed reasoning is given, and it seems quite extreme for complex designs. But most others (e.g. this or this or this answer or this book chapter or this video tutorial) suggests that only the residuals should be examined (regardless of within/between factors). Even if I assume that this is latter true, the question remains: which residuals should be examined?

In the following, I use the R function stats:aov output to illustrate in an example some potential answers.

I prepared an invented dataset for illustration. Each individual subject is denoted with "subject_id". There are two between-subject factors: "btwn_X" and "btwn_Y". There are also two within-subject factors: "wthn_X" and "wthn_Y".

# preparing some invented data    
dat_example = data.frame(
    subject = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
    btwn_X = c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2),
    btwn_Y = c(1, 2, 1, 2, 2, 1, 1, 1, 2, 1),
    measure_x1_yA = c(36.2, 45.2, 41, 24.6, 30.5, 28.2, 40.9, 45.1, 31, 16.9),
    measure_x2_yA = c(-14.1, 58.5, -25.5, 42.2, -13, 4.4, 55.5, -28.5, 25.6, -37.1),
    measure_x1_yB = c(83, 71, 111, 70, 92, 75, 110, 111, 110, 85),
    measure_x2_yB = c(8.024, -14.162, 3.1, -2.1, -1.5, 0.91, 11.53, 18.37, 0.3, -0.59),
    measure_x1_yC = c(27.4,-17.6,-32.7, 0.4, 37.2, 1.7, 18.2, 8.9, 1.9, 0.4),
    measure_x2_yC = c(7.7, -0.8, 2.2, 14.1, 22.1, -47.7, -4.8, 8.6, 6.2, 18.2)
)
dat_example$subject = as.factor(as.character(dat_example$subject))
dat_example$btwn_X = as.factor(as.character(dat_example$btwn_X))
dat_example$btwn_Y = as.factor(as.character(dat_example$btwn_Y))    
vars = c(
    'measure_x1_yA',
    'measure_x2_yA',
    'measure_x1_yB',
    'measure_x2_yB',
    'measure_x1_yC',
    'measure_x2_yC'
)
dat_l = stats::reshape(
    dat_example,
    direction = 'long',
    varying = vars,
    idvar = 'subject',
    timevar = "within_factor",
    v.names = "values",
    times = vars
)    
dat_l$wthn_X = sapply(strsplit(dat_l$within_factor, split = '_', fixed =
                                   TRUE), `[`, 2)
dat_l$wthn_Y = sapply(strsplit(dat_l$within_factor, split = '_', fixed =
                                   TRUE), `[`, 3)
dat_l$wthn_X = as.factor(as.character(dat_l$wthn_X))
dat_l$wthn_Y = as.factor(as.character(dat_l$wthn_Y))

# performing the ANOVA    
aov_BBWW = aov(values ~ btwn_X * btwn_Y * wthn_X * wthn_Y +
                   Error(subject / (wthn_X * wthn_Y)), data = dat_l)

(See also here an extended version with various within/between factor variations and lme4::lmer models.)

The aov object aov_BBWW returns the following:

Grand Mean: 23.6847

Stratum 1: subject

Terms:
                  btwn_X   btwn_Y btwn_X:btwn_Y Residuals
Sum of Squares    61.549  351.672        18.969  3221.628
Deg. of Freedom        1        1             1         6

Residual standard error: 23.17192
15 out of 18 effects not estimable
Estimated effects may be unbalanced

Stratum 2: subject:wthn_X

Terms:
                   wthn_X btwn_X:wthn_X btwn_Y:wthn_X btwn_X:btwn_Y:wthn_X Residuals
Sum of Squares  23432.120       612.948       712.387              773.779   513.165
Deg. of Freedom         1             1             1                    1         6

Residual standard error: 9.248106
8 out of 12 effects not estimable
Estimated effects may be unbalanced

Stratum 3: subject:wthn_Y

Terms:
                   wthn_Y btwn_X:wthn_Y btwn_Y:wthn_Y btwn_X:btwn_Y:wthn_Y Residuals
Sum of Squares  19262.400       982.159      1561.578             1836.188  5860.787
Deg. of Freedom         2             2             2                    2        12

Residual standard error: 22.09975
8 out of 16 effects not estimable
Estimated effects may be unbalanced

Stratum 4: subject:wthn_X:wthn_Y

Terms:
                wthn_X:wthn_Y btwn_X:wthn_X:wthn_Y btwn_Y:wthn_X:wthn_Y
Sum of Squares      20248.558              159.421              986.331
Deg. of Freedom             2                    2                    2
                btwn_X:btwn_Y:wthn_X:wthn_Y Residuals
Sum of Squares                      604.163  4789.399
Deg. of Freedom                           2        12

Residual standard error: 19.9779
Estimated effects may be unbalanced

I can access the following residuals (see here for more details):

aov_BBWW$subject$residuals
aov_BBWW$`subject:wthn_X`$residuals
aov_BBWW$`subject:wthn_Y`$residuals
aov_BBWW$`subject:wthn_X:wthn_Y`$residuals
aov_BBWW$`(Intercept)`$residuals

According to some of the sources cited above, these residuals should be used for normality testing, though it is not clear whether all or just one (and in that case which one).

---

EDIT:

After a lot of digging (and with the help of EdM's answer and comments), the most authoritative solution appears to be that in case of an ANOVA with only between-subject factors the correct variable is simply the residuals vector from the aov object (e.g. aov_BB$residuals), while in case there is any within-subject variable, I should do something like this:

aov_proj = proj(aov_BBWW)
aov_proj[[length(aov_proj)]][,"Residuals"]

Where the latter is the variable to be examined for normality and other related assumptions. Why this is so is beyond me, but several seemingly confident sources give this solution: this and this R mailing list replies, this and this and this CV answers (the latter two ironically not the accepted ones), this tutorial, and the MASS documentation. Most or perhaps all these sources originate from Venables and Ripley (2002), but I'd assume they would not all blindly copy something incorrect.

The question nonetheless is still open: I would be happy to receive further verification (or refutation) and explanation on the matter.

(Btw, if the above sources are to be trusted, the fitted values can apparently be accessed as: fitted(aov_BBWW[[length(aov_BBWW)]]) )

Answer 1

TL;DR: ANOVA pools information among all observations to get the best estimates of fixed effects, random effects, and error variance. If you want to examine normality of ANOVA residuals, doing so after all fixed and random effects are taken into account thus makes the most sense. Reliable ANOVA estimates don't require normality of residuals; the issue is the distribution of the test statistics. In repeated-measures ANOVA, issues like imbalance or mis-specification of correlation structures might be even more substantial obstacles to reliable statistical tests.

ANOVA is simply a particular type of a linear model, as described for example on this page of one of the sites that was linked from the question, and discussed extensively here. Like all linear models, ANOVA combines information from the combinations of predictor values to model the outcome values as a function of the predictors plus an error term. The error term is assumed to have a certain distribution shared among all cases, Gaussian with zero mean for standard ANOVA. Information about the distribution of the error terms is obtained by pooling across all the observations, smoothing out the vagaries that can happen just by chance within individual cells of the ANOVA design. A standard normal q-q diagnostic plot thus examines all the residual values, not those within individual cells.

Despite the usual assumption of Gaussian errors in an ANOVA model, the significance tests don't necessarily require that assumption to be met. Significance tests in ANOVA are tests on regression coefficients. It's thus the sampling distributions of those regression coefficients that must adequately meet assumptions when one performs a standard parametric test.

As @whuber put it in a crucially important comment:

What you really want to know is whether the assumed distributions of the ANOVA test statistics are sufficiently accurate to compute the p-values in which you are interested.

If the model assumptions are met and the shared error term has a Gaussian distribution then you know that tests on regression coefficients will be valid.* But strict normality of the error term isn't required for tests on the regression coefficients to be valid. Think about normally distributed error terms as sufficient but not always necessary for an adequately reliable significance test on linear model regression coefficients, including ANOVA.

That's not to say that it's useless to examine the distribution of residuals around model predictions that incorporate information from all cases. For example, the R lme4 package provides a normal q-q plot as one of its diagnostic plots; see page 33 of the vignette. What you will often find, however, is that substantial deviations from normality in such a plot of residuals mean that the model itself is poorly specified. That might be the most useful information from such a plot.

With a mixed ANOVA model having only fixed categorical predictors and including all interactions, you shouldn't have to worry about linearity in the fixed-effect predictors themselves. But there could be an incorrect handling of the outcome variable (e.g., if it's fundamentally log-normal rather than normal), omission of critical covariates associated both with outcome and with the included predictors, or mis-specification of the random-effects structure. Fix those problems exposed by the diagnostic plot rather than obsess about the normality per se.

To evaluate the model all the diagnostic plots should be examined: not only the q-q plot for normality of residuals but also the fitted vs. residual plot and the scale-location plot and the various profile plots (see page 36 of the vignette) for mixed models and their random effects. Examine undue influence of particular observations, e.g. with the influence.ME package in R. This process, rather than a simple examination of normality, is critical to evaluating and improving the quality of the model specification.

If the model is properly specified then the normality assumption on the sampling distribution of the regression coefficients can be reasonably reliable. With enough data the Central Limit Theorem can help with that despite non-normal residuals, although how much data is "enough" depends on the particular case. See this answer, for example. If you don't want to rely on that assumption, bootstrapping provides a way to get non-parametric confidence intervals. But that should be done only when the model itself is adequately specified.

---

As an edit to the question notes, some diagnostic plots can be generated from repeated-measures data analyzed by aov, which according to its manual page fits "an analysis of variance model by a call to lm for each stratum." Each stratum is a portioning of the means of the observations by progressively complex models, starting with the overall mean. As Venables and Ripley say on page 283 with respect to a simpler split-plot design:

Multistratum models may be fitted using aov, and are specified by a model formula of the form

response ~ mean.formula + Error (strata.formula)

In our example the strata.formula is B/V, specifying strata 2 and 3; the fourth stratum is included automatically as the "within" stratum, the residual stratum from the strata formula.

For more complicated models, the last stratum is thus the automatically included "within" stratum. Continuing on page 284: "It is not possible to associate [fitted values and residuals from the last stratum] uniquely with the plots of the original experiment." You need the residuals from "the projections of the original data vector onto the subspaces defined by each line in the analysis of variance tables." The residuals can be examined for every stratum, but only the final stratum takes all aspects of the model into account. This answer shows the code for the Venables and Ripley example in which the fourth stratum is the "within" stratum.

Before proceeding with aov, however, pay attention to the following quote from its help page:

Note

aov is designed for balanced designs, and the results can be hard to interpret without balance: beware that missing values in the response(s) will likely lose the balance. If there are two or more error strata, the methods used are statistically inefficient without balance, and it may be better to use lme in package nlme.

---

*This is more complicated with mixed models, for which there is dispute about the number of degrees of freedom to use in the test. But that dispute won't be resolved by examining the distribution of residuals. Tests on mixed models can also involve assumptions about the covariance structure of correlated observations.