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
After a meeting with a Biostatistic Professor I can say: There is no other possibility besides "calculating all group combinations seperate".
If your contrasts are significant, you have to discuss this significance based on individual calculation of the group combinations and based on the profil plots