[R] contrast matrix for aov

[Prof Brian Ripley]

On Wed, 9 Mar 2005, Darren Weber wrote:

> How do we specify a contrast interaction matrix for an ANOVA model?
>
> We have a two-factor, repeated measures design, with

Where does `repeated measures' come into this?  You appear to have 
repeated a 2x2 experiment in each of 8 blocks (subjects).  Such a design 
is usually analysed with fixed effects.  (Perhaps you averaged over 
repeats in the first few lines of your code?)

> Cue Direction (2) x  Brain Hemisphere(2)
>
> Each of these has 2 levels, 'left' and 'right', so it's a simple 2x2 design 
> matrix.  We have 8 subjects in each cell (a balanced design) and we want to 
> specify the interaction contrast so that:
>
> CueLeft>CueRght for the Right Hemisphere
> CueRght>CueLeft for the Left Hemisphere.
>
> Here is a copy of the relevant commands for R:
>
> ########################################
> lh_cueL <- rowMeans( LHroi.cueL[,t1:t2] )
> lh_cueR <- rowMeans( LHroi.cueR[,t1:t2] )
> rh_cueL <- rowMeans( RHroi.cueL[,t1:t2] )
> rh_cueR <- rowMeans( RHroi.cueR[,t1:t2] )
> roiValues <- c( lh_cueL, lh_cueR, rh_cueL, rh_cueR )
>
> cuelabels <- c("CueLeft", "CueRight")
> hemlabels <- c("LH", "RH")
>
> roiDataframe <- data.frame( roi=roiValues, Subject=gl(8,1,32,subjectlabels), 
> Hemisphere=gl(2,16,32,hemlabels), Cue=gl(2,8,32,cuelabels) )
>
> roi.aov <- aov(roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)), 
> data=roiDataframe)

I think the error model should be Error(Subject).  In what sense are `Cue' 
and `Cue:Hemisphere' random effects nested inside `Subject'?

Let me fake some `data':

set.seed(1); roiValues <- rnorm(32)
subjectlabels <- paste("V"1:8, sep = "")
options(contrasts = c("contr.helmert", "contr.poly"))
roi.aov <- aov(roi ~ Cue*Hemisphere + Error(Subject), data=roiDataframe)

> roi.aov

Call:
aov(formula = roi ~ Cue * Hemisphere + Error(Subject), data = roiDataframe)

Grand Mean: 0.1165512

Stratum 1: Subject

Terms:
                 Residuals
Sum of Squares   4.200946
Deg. of Freedom         7

Residual standard error: 0.7746839

Stratum 2: Within

Terms:
                       Cue Hemisphere Cue:Hemisphere Residuals
Sum of Squares   0.216453   0.019712       0.057860 21.896872
Deg. of Freedom         1          1              1        21

Residual standard error: 1.021131
Estimated effects are balanced

Note that all the action is in one stratum, and the SSQs are the same 
as

aov(roi ~ Subject + Cue * Hemisphere, data = roiDataframe)

(and also the same as for your fit).

> print(summary(roi.aov))

It auto-prints, so you don't need print().

> ########################################
>
>
> I've tried to create a contrast matrix like this:
>
> cm <- contrasts(roiDataframe$Cue)
>
> which gives the main effect contrasts for the Cue factor.  I really want to 
> specify the interaction contrasts, so I tried this:
>
> ########################################
> # c( lh_cueL, lh_cueR, rh_cueL, rh_cueR )
> # CueRight>CueLeft for the Left Hemisphere.
> # CueLeft>CueRight for the Right Hemisphere
>
> cm <- c(-1, 1, 1, -1)
> dim(cm) <- c(2,2)

(That is up to sign what Helmert contrasts give you.)

> roi.aov <- aov( roi ~ (Cue*Hemisphere) + Error(Subject/(Cue*Hemisphere)),
> contrasts=cm, data=roiDataframe)
> print(summary(roi.aov))
> ########################################
>
> but the results of these two aov commands are identical.  Is it the case that 
> the 2x2 design matrix is always going to give the same F values for the 
> interaction regardless of the contrast direction?

Yes, as however you code the design (via `contrasts') you are fitting the 
same subspaces.  Not sure what you mean by `contrast direction', though.

However, you have not specified `contrasts' correctly:

contrasts: A list of contrasts to be used for some of the factors in
           the formula.

and cm is not a list, and an interaction is not a factor.

> OR, is there some way to get a summary output for the contrasts that is 
> not available from the print method?

For more than two levels, yes: see `split' under ?summary.aov.
Also, see se.contrasts which allows you to find the standard error for any 
contrast.

For the fixed-effects model you can use summary.lm:

> fit <- aov(roi ~ Subject + Cue * Hemisphere, data = roiDataframe)
> summary(fit)
                Df  Sum Sq Mean Sq F value Pr(>F)
Subject         7  4.2009  0.6001  0.5756 0.7677
Cue             1  0.2165  0.2165  0.2076 0.6533
Hemisphere      1  0.0197  0.0197  0.0189 0.8920
Cue:Hemisphere  1  0.0579  0.0579  0.0555 0.8161
Residuals      21 21.8969  1.0427
> summary.lm(fit)

Call:
aov(formula = roi ~ Subject + Cue * Hemisphere, data = roiDataframe)

Residuals:
     Min      1Q  Median      3Q     Max
-1.7893 -0.4197  0.1723  0.5868  1.3033

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
[...]
Cue1             -0.08224    0.18051  -0.456    0.653
Hemisphere1       0.02482    0.18051   0.137    0.892
Cue1:Hemisphere1 -0.04252    0.18051  -0.236    0.816

where the F values are the squares of the t values.


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595