1
$\begingroup$

I was wondering how to assess residual normality of a repeated measures ANOVA. In some threads, users refer to Venables and Ripley: Residuals in multistratum analyses: Projections and recommend to extract residuals using proj() function. Is there any argument speaking against application of statistical normality tests, like shapioro wilk test?

Or would a mixed model (lmer from lme4 package) be the better way?

Reproducible example:

library(MASS)
set.seed(123)

data<- data.frame(id =  factor(rep(1:10, each = 4)),
       cond1 = factor(rep(c("a", "b"), 20)),
       cond2 = factor(rep(rep(c("x", "y"), each = 2), 10)),
       Y = rnorm(40, 5, 2))

model<- aov(Y ~ cond1*cond2 +
      Error(id/(cond1*cond2)), data = data)

model.pr <- proj(model)

shapiro.test(model.pr[[5]][, "Residuals"])
Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ You might want to read the page on whether normality testing is essentially useless. Evaluating plots that document heteroscedasticity or major deviations from normality are important, but formal normality tests tend to be underpowered for small sample sizes, where they might be considered most important. Remember that for inference it's the distribution of the coefficient estimates that needs to be normal; that can be close enough to true even if the residuals aren't, particularly with large samples. $\endgroup$
    – EdM
    Commented Mar 22, 2023 at 19:27
  • $\begingroup$ Dear @EdM, thanks for your comment. Do you have any source by chance that explains why coefficient estimates are important (and not residuals) and how to assess distribution of those? $\endgroup$
    – a.henrietty
    Commented Mar 23, 2023 at 14:18
  • $\begingroup$ Schmidt and Finan have an extensive discussion in Linear Regression and the Normality Assumption, J. Clin Epidemiol. 98: 146-151, 2018. Repeated modeling on multiple bootstrap samples of the data (re-sampled by the individual for repeated measures, rather than by the individual observation) is probably a better way to get confidence intervals and the like anyway, in a way that doesn't depend on normality, but you could also use that to document normality of the distributions of coefficient estimates. $\endgroup$
    – EdM
    Commented Mar 23, 2023 at 14:53
  • $\begingroup$ @EdM thank you so much. Suppose one would be interested in the residuals. As my example does only feature factors with two levels, there are no within residuals - would the use of the last stratum (interaction of cond1:cond2) be valid? (Refering to Venables and Ripley) $\endgroup$
    – a.henrietty
    Commented Mar 23, 2023 at 16:10

1 Answer 1

Reset to default
1
$\begingroup$

You have to be careful with the residuals estimated by projection. Your example (without what you call "within residuals," given that there are no replicates within an id under the same condition) shows this nicely.

With your data and model:

length(model.pr[[5]][,"Residuals"])
# [1] 40

But there are patterns within those residuals. Each individual has 4 residual values, in pairs of pairs: all 4 have the same absolute values, with 2 positive and 2 negative. For example, for id=1:

model.pr[[5]][1:4,"Residuals"]
#          1          2          3          4 
# -0.6504411  0.6504411  0.6504411 -0.6504411

That doesn't seem like an appropriate set of data to submit to a Shapiro-Wilk test, even if you think that normality testing isn't essentially useless. A truly normal data set of 40 shouldn't have only 10 unique absolute values.

In Section 10.2 of Venables and Ripley, which discusses these matters, the residuals from projections are only used for a qqnorm() plot, not for a formal normality test. For a plot of residuals versus fitted, they further recommend using the fitted() and studres() from the model itself, at the last stratum. In your model there are only 10 such values, one for each individual:

length(studres(model[[5]]))
# [1] 10

The recommended diagnostic plots for your model would be based on:

plot(fitted(model[[5]]),studres(model[[5]]))
qqnorm(model.pr[[5]][,"Residuals"])

This page discusses some other issues in evaluating this type of model. Note that your data lead to a singular fit with a 0 estimate for random-intercept variance when modeled with lmer().

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for this elaborated answer. It was very helpful! $\endgroup$
    – a.henrietty
    Commented Apr 24, 2023 at 8:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.