Assuming your design is the following:
sex is a between-subjects IV (with two levels)stimulus is a within-subjects IV (with 3 assumed levels)condition is a within-subjects IV (with 2 levels)- all IVs are fully crossed
Then this is what you can do to run the full analysis, or to just test for a main effect of sex (generating some data first):
Nj <- 10 # number of subjects per sex
P <- 2 # number of levels for IV sex
Q <- 3 # number of levels for IV stimulus
R <- 2 # number of levels for IV condition
subject <- factor(rep(1:(P*Nj), times=Q*R)) # subject id
sex <- factor(rep(1:P, times=Q*R*Nj), labels=c("F", "M")) # IV sex
stimulus <- factor(rep(1:Q, each=P*R*Nj)) # IV stimulus
condition <- factor(rep(rep(1:R, each=P*Nj), times=Q), labels=c("EXP1", "EXP2"))
DV_t11 <- round(rnorm(P*Nj, 8, 2), 2) # responses for stimulus=1 and condition=1
DV_t21 <- round(rnorm(P*Nj, 13, 2), 2) # responses for stimulus=2 and condition=1
DV_t31 <- round(rnorm(P*Nj, 13, 2), 2)
DV_t12 <- round(rnorm(P*Nj, 10, 2), 2)
DV_t22 <- round(rnorm(P*Nj, 15, 2), 2)
DV_t32 <- round(rnorm(P*Nj, 15, 2), 2)
response <- c(DV_t11, DV_t12, DV_t21, DV_t22, DV_t31, DV_t32) # all responses
dfL <- data.frame(subject, sex, stimulus, condition, response) # long format
Now with the data set up, you can use aov(), but you won't get the $\hat{\epsilon}$ corrections for the within-effects.
> summary(aov(response ~ sex*stimulus*condition
+ + Error(subject/(stimulus*condition)), data=dfL))
Error: subject
Df Sum Sq Mean Sq F value Pr(>F)
sex 1 2.803 2.8030 0.51 0.4843 # ... snip ...
You can also use the Anova() function from the car package, which gives you the $\hat{\epsilon}$ corrections. However, it requires your data to be in wide format. You have to use multivariate notation for your model formula.
> sexW <- factor(rep(1:P, Nj), labels=c("F", "M")) # factor sex for wide format
> dfW <- data.frame(sexW, DV_t11, DV_t21, DV_t31, DV_t12, DV_t22, DV_t32) # wide format
> # between-model in multivariate notation
> fit <- lm(cbind(DV_t11, DV_t21, DV_t31, DV_t12, DV_t22, DV_t32) ~ sexW, data=dfW)
> # dataframe describing the columns of the data matrix
> intra <- expand.grid(stimulus=gl(Q, 1), condition=gl(R, 1))
> library(car) # for Anova()
> summary(Anova(fit, idata=intra, idesign=~stimulus*condition),
+ multivariate=FALSE, univariate=TRUE)
Univariate Type II Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 17934.1 1 98.930 18 3263.0403 < 2.2e-16 ***
sexW 2.8 1 98.930 18 0.5100 0.4843021 # ... snip ...
Using the ez package and the command suggested by @Mike Lawrence gives the same result:
> library(ez) # for ezANOVA()
> ezANOVA(data=dfL, wid=.(subject), dv=.(response),
+ within=.(stimulus, condition), between=.(sex), observed=.(sex))
$ANOVA
Effect DFn DFd F p p<.05 ges
2 sex 1 18 0.5099891 4.843021e-01 0.004660043 # ... snip ...
Finally, if the main effect for sex is really all you're interested in, it's equivalent to just average for each person across all the conditions created by the combinations of stimulus and condition, and then run a between-subjects ANOVA for the aggregated data.
# average per subject across all repeated measures
> mDf <- aggregate(response ~ subject + sex, data=dfL, FUN=mean)
> summary(aov(response ~ sex, data=mDf)) # ANOVA with just the between-effect
Df Sum Sq Mean Sq F value Pr(>F)
sex 1 0.4672 0.46716 0.51 0.4843
Residuals 18 16.4884 0.91602
Best Answer
You may not get a simple response to
residuals(npk.aovE)but that does not mean there are no residuals in that object. Dostrand see that within the levels there are still residuals. I would imagine you were most interested in the "Within" levelMy own training and practice has not been to use normality testing, instead to use QQ plots and parallel testing with robust methods.