Testing the normality assumption for repeated measures anova? (in R) So assuming that there is a point in testing the normality assumption for anova (see 1 and 2)
How can it be tested in R?
I would expect to do something like:
## From Venables and Ripley (2002) p.165.
utils::data(npk, package="MASS")
npk.aovE <- aov(yield ~  N*P*K + Error(block), npk)
residuals(npk.aovE)
qqnorm(residuals(npk.aov))

Which doesn't work, since "residuals" don't have a method (nor predict, for that matter) for 
the case of repeated measures anova.
So what should be done in this case?
Can the residuals simply be extracted from the same fit model without the Error term?  I am not familiar enough with the literature to know if this is valid or not, thanks in advance for any suggestion.
 A: You may not get a simple response to residuals(npk.aovE) but that does not mean there are no residuals in that object. Do str and see that within the levels there are still residuals. I would imagine you were most interested in the "Within" level
> residuals(npk.aovE$Within)
          7           8           9          10          11          12 
 4.68058815  2.84725482  1.56432584 -5.46900749 -1.16900749 -3.90234083 
         13          14          15          16          17          18 
 5.08903669  1.28903669  0.35570336 -3.27762998 -4.19422371  1.80577629 
         19          20          21          22          23          24 
-3.12755705  0.03910962  2.60396981  1.13730314  2.77063648  4.63730314 

My own training and practice has not been to use normality testing, instead to use QQ plots and parallel testing with robust methods.
A: Another option would be to use the lme function of the nlme package (and then pass the obtained model to anova). You can use residuals on its output.
A: Venables and Ripley explain how to do residual diagnostics for a repeated-measures design later in their book (p. 284), in the section on random and mixed effects.
The residuals function (or resid) is implemented for the aov results for each stratum:
from their example:
oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)
To get the fitted values or residuals:
"Thus fitted(oats.aov[[4]]) and resid(oats.aov[[4]]) are vectors of length 54 representing the fitted values and residuals from the last stratum."
Importantly, they add:
"It is not possible to associate them uniquely with the plots of the original experiment."
For diagnostics, they use a projection:
plot(fitted(oats.aov[[4]]), studres(oats.aov[[4]]))
abline(h=0, lty=2)
oats.pr <- proj(oats.aov)
qqnorm(oats.pr[[4]][, "Residuals"], ylab = "Stratum 4 residuals")
qqline(oats.pr[[4]][, "Residuals"])

They also show that the model can be done using lme, as another user posted.
A: Here are two options, with aov and with lme (I think the 2nd is preferred):
require(MASS) ## for oats data set
require(nlme) ## for lme()
require(multcomp) ## for multiple comparison stuff

Aov.mod <- aov(Y ~ N * V + Error(B/V), data = oats)
the_residuals <- aov.out.pr[[3]][, "Residuals"]

Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)
the_residuals <- residuals(Lme.mod)

The original example came without the interaction (Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)) but it seems to be working with it (and producing different results, so it is doing something). 
And that's it...
but for completeness:
1 - The summaries of the model
summary(Aov.mod)
anova(Lme.mod)

2 - The Tukey test with repeated measures anova (3 hours looking for this!!).
summary(Lme.mod)
summary(glht(Lme.mod, linfct=mcp(V="Tukey")))

3 - The normality and homoscedasticity plots
par(mfrow=c(1,2)) #add room for the rotated labels
aov.out.pr <- proj(aov.mod)                                            
#oats$resi <- aov.out.pr[[3]][, "Residuals"]
oats$resi <- residuals(Lme.mod)
qqnorm(oats$resi, main="Normal Q-Q") # A quantile normal plot - good for checking normality
qqline(oats$resi)
boxplot(resi ~ interaction(N,V), main="Homoscedasticity", 
        xlab = "Code Categories", ylab = "Residuals", border = "white", 
        data=oats)
points(resi ~ interaction(N,V), pch = 1, 
       main="Homoscedasticity",  data=oats)


A: I think that the normality assumption can be assessed for each of the repeated measures, before performing the analysis. I would reshape the dataframe so that each column corresponds to a repeated measure, and then perform a shapiro.test to each one of those columns. 
apply(cast(melt(npk,measure.vars="yield"), ...~N+P+K)[-c(1:2)],2,function(x) shapiro.test(x)$p.value)

