Solved – How to specify specific contrasts for repeated measures ANOVA using car

anovacontrastsrrepeated measuressums-of-squares

I am trying to run a repeated measures Anova in R followed by some specific
contrasts on that dataset. I think the correct approach would be to use
Anova() from the car package.

Lets illustrate my question with the example taken from ?Anova using the
OBrienKaiser data (Note: I ommited the gender factor from the example):
We have a design with one between subjects factor, treatment (3 levels: control,
A, B), and 2 repeated-measures (within subjects) factors,
phase (3 levels: pretest, posttest, followup) and hour (5 levels: 1 to 5).

The standard ANOVA table is given by (in difference to example(Anova) I switched
to Type 3 Sums of Squares, that is what my field wants):

require(car)
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser)
av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = 3)
summary(av.ok, multivariate=FALSE)

Now, imagine that the highest order interaction would have been significant
(which is not the case) and we would like to explore it further with the
following contrasts:
Is there a difference between hours 1&2 versus hours 3 (contrast 1) and between
hours 1&2 versus hours 4&5 (contrast 2) in the treatment conditions (A&B
together)?
In other words, how do I specify these contrasts:

((treatment %in% c("A", "B")) & (hour %in% 1:2)) versus ((treatment %in% c("A", "B")) & (hour %in% 3))
((treatment %in% c("A", "B")) & (hour %in% 1:2)) versus ((treatment %in% c("A", "B")) & (hour %in% 4:5))

My idea would be to run another ANOVA ommitting the non-needed treatment
condition (control):

mod2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser, subset = treatment != "control")
av2 <- Anova(mod2, idata=idata, idesign=~phase*hour, type = 3)
summary(av2, multivariate=FALSE)

However, I still have no idea how to set up the appropriate
within-subject contrast matrix comparing hours 1&2 with 3 and 1&2 with 4&5.
And I am not sure if omitting the non-needed treatment group is indeed a good
idea as it changes the overall error term.

Before going for Anova() I was also thinking going for lme. However, there are
small differences in F and p values between textbook ANOVA and what is returned
from anove(lme) due to possible negative variances in standard ANOVA (which are not allowed in lme). Relatedly, somebody pointed me to gls which allows for fitting repeated measures ANOVA, however, it has no contrast argument.

To clarify: I want an F or t test (using type III sums of squares) that answers whether or not the desired contrasts are significant or not.

Update:

I already asked a very similar question on R-help, there was no answer.

A similar questions was posed on R-help some time ago. However, the answers did also not solve the problem.

Update (2015):

As this question still generates some activity, specifying theses and basically all other contrasts can now be done relatively easy with the afex package in combination with the lsmeans package as described in the afex vignette.

Best Answer

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)

Related Solutions

Repeated Measures ANOVA – Why Do lme and aov Return Different Results in R? A Comprehensive Analysis

They are different because the lme model is forcing the variance component of id to be greater than zero. Looking at the raw anova table for all terms, we see that the mean squared error for id is less than that for the residuals.

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131

When we compute the variance components, this means that the variance due to id will be negative. My memory of expected mean squares memory is shaky, but the calculation is something like

(0.15052-0.16131)/3 = -0.003597.

This sounds odd but can happen. What it means is that the averages for each id are closer to each other than you would expect to each other given the amount of residual variation in the model.

In contrast, using lme forces this variance to be greater than zero.

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

This reports standard deviations, squaring to get the variance yields 9.553e-10 for the id variance and 0.1586164 for the residual variance.

Now, you should know that using aov for repeated measures is only appropriate if you believe that the correlation between all pairs of repeated measures is identical; this is called compound symmetry. (Technically, sphericity is required but this is sufficient for now.) One reason to use lme over aov is that it can handle different kinds of correlation structures.

In this particular data set, the estimate for this correlation is negative; this helps explain how the mean squared error for id was less than the residual squared error. A negative correlation means that if an individual's first measurement was below average, on average, their second would be above average, making the total averages for the individuals less variable than we would expect if there was a zero correlation or a positive correlation.

Using lme with a random effect is equivalent to fitting a compound symmetry model where that correlation is forced to be non-negative; we can fit a model where the correlation is allowed to be negative using gls:

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

This ANOVA table agrees with the table from the aov fit and from the lm fit.

OK, so what? Well, if you believe that the variance from id and the correlation between observations should be non-negative, the lme fit is actually more appropriate than the fit using aov or lm as its estimate of the residual variance is slightly better. However, if you believe the correlation between observations could be negative, aov or lm or gls is better.

You may also be interested in exploring the correlation structure further; to look at a general correlation structure, you'd do something like

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

Here I only limit the output to the correlation structure. The values 1 to 4 represent the four levels of factor; we see that factor 1 and factor 4 have a fairly strong negative correlation:

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

One way to choose between these models is with a likelihood ratio test; this shows that the random effects model and the general correlation structure model aren't statistically significantly different; when that happens the simpler model is usually preferred.

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965

Solved – a valid post-hoc analysis for a three-way repeated measures ANOVA

I think statisticians will tell you that there is always a problem with any post hoc analysis because seeing the data may influence what you look at and you could be biased becuase you are hunting for significant results. The FDA in clinical trial studies requires that the statistical plan be completely spelled out in the protocol. in a linear model you certainly could prespecify the contrasts that you would like to look at in the event that the ANOVA or ANCOVA finds an overall difference. Such prespecified contrasts would be fine to look at as long as the usual treatment for multiplicity is also part of it.

Best Answer

Related Solutions

Repeated Measures ANOVA – Why Do lme and aov Return Different Results in R? A Comprehensive Analysis

Solved – a valid post-hoc analysis for a three-way repeated measures ANOVA

Related Question