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
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 "
# 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
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).
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: