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
This method is generally considered "old-fashioned" so while it may be possible, the syntax is difficult and I suspect fewer people know how to manipulate the anova commands to get what you want. The more common method is using glht with a likelihood-based model from nlme or lme4. (I'm certainly welcome to be proved wrong by other answers though.)
That said, if I needed to do this, I wouldn't bother with the anova commands; I'd just fit the equivalent model using lm, pick out the right error term for this contrast, and compute the F test myself (or equivalently, t test since there's only 1 df). This requires everything to be balanced and have sphericity, but if you don't have that, you should probably be using a likelihood-based model anyway. You might be able to somewhat correct for non-sphericity using the Greenhouse-Geiser or Huynh-Feldt corrections which (I believe) use the same F statistic but modify the df of the error term.
If you really want to use car, you might find the heplot vignettes helpful; they describe how the matrices in the car package are defined.
Using caracal's method (for the contrasts 1&2 - 3 and 1&2 - 4&5), I get
psiHat tStat F pVal
1 -3.0208333 -7.2204644 52.1351067 2.202677e-09
2 -0.2083333 -0.6098777 0.3719508 5.445988e-01
This is how I'd get those same p-values:
Reshape the data into long format and run lm to get all the SS terms.
library(reshape2)
d <- OBrienKaiser
d$id <- factor(1:nrow(d))
dd <- melt(d, id.vars=c(18,1:2), measure.vars=3:17)
dd$hour <- factor(as.numeric(gsub("[a-z.]*","",dd$variable)))
dd$phase <- factor(gsub("[0-9.]*","", dd$variable),
levels=c("pre","post","fup"))
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)
anova(m)
Make an alternate contrast matrix for the hour term.
foo <- matrix(0, nrow=nrow(dd), ncol=4)
foo[dd$hour %in% c(1,2) ,1] <- 0.5
foo[dd$hour %in% c(3) ,1] <- -1
foo[dd$hour %in% c(1,2) ,2] <- 0.5
foo[dd$hour %in% c(4,5) ,2] <- -0.5
foo[dd$hour %in% 1 ,3] <- 1
foo[dd$hour %in% 2 ,3] <- 0
foo[dd$hour %in% 4 ,4] <- 1
foo[dd$hour %in% 5 ,4] <- 0
Check that my contrasts give the same SS as the default contrasts (and the same as from the full model).
anova(lm(value ~ hour, data=dd))
anova(lm(value ~ foo, data=dd))
Get the SS and df for just the two contrasts I want.
anova(lm(value ~ foo[,1], data=dd))
anova(lm(value ~ foo[,2], data=dd))
Get the p-values.
> F <- 73.003/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 2.201148e-09
> F <- .5208/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 0.5445999
Optionally adjust for sphericity.
pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)
Best Answer
I'm not sure what you plan to test but looking at your first graph it seems pretty clear. There's generally an effect of colour but on the left side and at larger offsets it disappears.
I'm guessing you wanted to test all of the effects of colour to see where they were significant and where they weren't. If you found they were all significant or all not it wouldn't give your more information than your interaction (and the interaction is independent of such findings). If you found some were and some weren't it doesn't test your interaction because the difference between significant and not significant is not itself significant.
I suppose you could try to explore something else about the difference in patterns but given that offset seems like a continuous variable that should have some clear function if it is doing anything I'm not sure establishing anything else different about the waviness of the lines on the two sides would be where you'd want to go.
UPDATE after viewing comments
The explanation for your interaction would be that the effect of colour and side are consistent until offsets 4 and 5 where the effect of colour only exists for one of the sides. That's just a recasting of what I said before in line with your hypotheses.
Keep in mind what you'd have to test post hoc or in planned contrasts to really demonstrate this. Finding an effect just at large offsets is useless because colour itself is interacting with side; therefore you have to show it interacts with offsets to show they're having an effect as well. That's what your ANOVA is telling you. It's already the planned contrast you want. Is there anything else there that could be making that interaction occur? Do you need to explain anything else?
If you do the ANOVA at 4 and 5 you likely won't get an interaction with offset, just one between colour and side, which will be substantially less evidence for what you want to say, not more.
Keep in mind, interactions mean something. Look at your data and figure out what they mean before considering further statistical tests. If they're relatively clear, as in this case, then you're done.