[R] contrasts

Craig Brozinsky cjbroz at ucdavis.edu
Wed Nov 9 01:48:52 CET 2005


Hi,

I'm having difficulty specifying contrasts for a within subjects  
factor with 3 levels. I can do it correctly for my factors with 2- 
levels, but i'm not getting the correct results for a 3-level factor.

My design has 1 between subjects factor (gp) and 3 within subjects  
factors. Within factor "w" has 2 levels, within factor "x" has two  
factors, and within
factor "y" has 3 factors. i first ran an ombibus F-test using aov. i  
then compared all of my contrasts to that output to determine if i  
was calculating the
contrasts correctly.

i first calculated the within factor w's main effect and interaction  
with the aov and then using a contrast:
 > summary(aov(effect ~ gp * w * x * y + Error(subj/(w+x+y)),  
data=both.uni))
Error: subj:w
           Df Sum Sq Mean Sq F value Pr(>F)
w          1 1.3889  1.3889  1.1364 0.3465
gp:w       1 0.2222  0.2222  0.1818 0.6918
Residuals  4 4.8889  1.2222
...

with contrasts i got the same results using the following code:
contr <- matrix(c(
1, 1, 1, 1, 1, 1,-1, -1, -1,-1,-1,-1,), ncol=1)

taov <- aov(cbind (
w1x1y, w1x1y2, w1x1y3, w1x2y, w1x2y2, w1x2y3,w2x1y, w2x1y2, w2x1y3,  
w2x2y, w2x2y2, w2x2y3)  %*% contr ~ gp, data = both.mul)
summary(taov,intercept=T)

now, things get a little more complicated for the three-level factor,  
'y.' If my memory serves me correctly, i can achieve this by writing  
a contrast that compares level 1 with level3 and level 2 with level  
3. is that correct? if so the code should look something like:

contr <- matrix(c(
1, 0,
0, 1,
-1, -1,
1, 0,
0, 1,
-1, -1,
1, 0,
0, 1,
-1, -1,
1, 0,
0, 1,
-1, -1, ), nrow=12, byrow=T)

tmp<-manova(cbind(
w1x1y, w1x1y2, w1x1y3, w1x2y, w1x2y2, w1x2y3, w2x1y, w2x1y2,  
w2x1y3,w2x2y, w2x2y2, w2x2y3)  %*% contr ~ gp, data = both.mul)
summary(tmp, test="Wilks", intercept=T)

when i run this code, it gives me the incorrect degrees of freedom  
and a p value different from what i obtained in my omnibus anova.  
does anyone have an idea of where i went wrong?

if access to my dataset is useful, i list it below.

# set up matrix for control subjects
cont.mat <- c 
(3,1,5,2,1,3,2,3,5,2,2,4,1,2,4,1,2,6,0,1,4,1,1,4,3,3,4,3,1,5,3,2,6,4,4,5 
,)
nsubs_c <- length(cont.mat)/12
con.mat <- matrix(cont.mat,  nsubs_c, 12, T)
con.mul <- cbind.data.frame(subj=1:nsubs_c, conds=factor(rep(1,rep 
(nsubs_c,1))), con.mat)

dimnames(con.mul)[[2]] <-
c("subj","gp",
"w1x1y", "w1x1y2", "w1x1y3",
"w1x2y", "w1x2y2", "w1x2y3",
"w2x1y", "w2x1y2", "w2x1y3",
"w2x2y", "w2x2y2", "w2x2y3",)

con.uni <- data.frame(effect = as.vector(con.mat),
subj = factor(paste("s", rep(1:nsubs_c, 12), sep="")),
gp = factor( paste("gp", rep(c(1), nsubs_c*12), sep="")),
w = factor(paste("w", rep(c(1,2), c(nsubs_c*3*2,nsubs_c*3*2)), sep="")),
x = factor( paste("x", rep(rep(c(1, 2), c(nsubs_c*3, nsubs_c*3)),2),  
sep="")),
y = factor(paste("y", rep(rep(c(1:3), rep(c(nsubs_c),3)),4), sep="")),
row.names = NULL)

# set up matrix for patients
pati.mat <- c 
(3,1,4,0,2,3,2,2,5,1,1,5,5,0,3,1,1,2,3,1,4,0,1,4,1,1,4,1,0,3,0,2,4,2,1,4 
,)

nsubs_p <- length(pati.mat)/12
pat.mat <- matrix(pati.mat,  nsubs_p, 12, T)
pat.mul <- cbind.data.frame(subj=1:nsubs_p, conds=factor(rep(2,rep 
(nsubs_p,1))), pat.mat)

dimnames(pat.mul)[[2]] <-
c("subj","gp",
"w1x1y", "w1x1y2", "w1x1y3",
"w1x2y", "w1x2y2", "w1x2y3",
"w2x1y", "w2x1y2", "w2x1y3",
"w2x2y", "w2x2y2", "w2x2y3",)

pat.uni <- data.frame(effect = as.vector(pat.mat),
subj = factor(paste("s", rep((nsubs_c+1):(nsubs_c+nsubs_p), 12),  
sep="")),
gp = factor( paste("gp", rep(c(2), nsubs_p*12), sep="")),
w = factor(paste("w", rep(c(1,2), c(nsubs_p*3*2,nsubs_p*3*2)), sep="")),
x = factor( paste("x", rep(rep(c(1, 2), c(nsubs_p*3, nsubs_p*3)),2),  
sep="")),
y = factor(paste("y", rep(rep(c(1:3), rep(c(nsubs_p),3)),4), sep="")),
row.names = NULL)

both.uni <- rbind(con.uni, pat.uni) # univariate file containing  
everone's data
both.mul <- rbind(con.mul, pat.mul) # multivariate file containing  
everone's data

############ 1: give me the whole anova ##########

summary(aov(effect ~ gp * w * x * y + Error(subj/(w+x+y)),  
data=both.uni))

############ 2: give me the w main effect and interaction with gp  
interaction ############
contr <- matrix(c(
-1, -1, -1,
-1,-1,-1,
1, 1, 1,
1, 1, 1,
), ncol=1)

taov <- aov(cbind (
w1x1y, w1x1y2, w1x1y3,
w1x2y, w1x2y2, w1x2y3,
w2x1y, w2x1y2, w2x1y3,
w2x2y, w2x2y2, w2x2y3)  %*% contr ~ gp, data = both.mul)
summary(taov,intercept=T)

############ 3: give me the y main effect and interaction with gp  
interaction ############

contr <- matrix(c(
1, 0,
0, 1,
-1, -1,
1, 0,
0, 1,
-1, -1,
1, 0,
0, 1,
-1, -1,
1, 0,
0, 1,
-1, -1, ), nrow=12, byrow=T)

tmp<-manova(cbind(
w1x1y, w1x1y2, w1x1y3,
w1x2y, w1x2y2, w1x2y3,
w2x1y, w2x1y2, w2x1y3,
w2x2y, w2x2y2, w2x2y3)  %*% contr ~ gp, data = both.mul)

summary(tmp, test="Wilks", intercept=T)




More information about the R-help mailing list