Number of contrasts in an orthogonal contrast test

Question

This is my first time asking a question here, so please let me know if anything is unclear.

So I was trying to learn how to set up contrasts from this page, and here is the code:

irrigation<-factor(c(rep("Control",10),rep("Irrigated 10 mm",10),rep("Irrigated 20 mm",10)))
biomass<-1:30
summary(aov(biomass~irrigation))
plot(irrigation, biomass)
contrasts(irrigation)
summary(aov(biomass~irrigation))
summary.lm(aov(biomass~irrigation))

contrastmatrix<-cbind(c(0,1,-1),c(-1,1,0))
contrastmatrix
contrasts(irrigation)<-contrastmatrix
contrasts(irrigation)
summary.lm(aov(biomass~irrigation))

This is fine and it produces the following outputs:

             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  15.5000     0.5528   28.04  < 2e-16 ***
irrigation1 -10.0000     0.7817  -12.79 5.67e-13 ***
irrigation2  10.0000     0.7817   12.79 5.67e-13 ***

However, I thought since these are orthogonal contrasts, I can just test the first contrast, as follows:

contrastmatrix<-cbind(c(0,1,-1))
contrastmatrix
contrasts(irrigation, 1)<-contrastmatrix
contrasts(irrigation)
summary.lm(aov(biomass~irrigation))

And here is the output:

              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   15.500      1.442   10.75 1.92e-11 ***
irrigation1   -5.000      1.767   -2.83  0.00851 ** 

However, for the same contrast, I am getting different coefficients. For "irrigation1", the coefficient was -10.0 before, but now it has changed to -5.0. Shouldn't they be the same since they are the same contrast?

Answer 1

Are you sure your contrasts are orthogonal? It's relatively easy to check.

contrastmatrix<-cbind(c(0,1,-1),c(-1,1,0))
crossprod(contrastmatrix)

which gives:

     [,1] [,2]
[1,]    2    1
[2,]    1    2

Nope, not orthogonal. Orthogonal contrasts would have 0 on the off-diagonal.

Answer 2

I came across a similar problem recently, and I got a comprehensive understanding from here. contrasts() command isn’t actually expecting contrast weights, it’s actually using the inverse of the matrix of desired contrast weights. You need to check the contrast matrix after lm/aov analysis if you specified your own contrast matrix:

#check your contrast matrix, making sure it fits your expectations.
#fit = aov(...)
attributes(fit$qr$qr)$contrasts

Also, if you add fewer than j-1 contrasts (j is total levels), your results will vary depending on whether you specify the contrasts variable or not.

Note that if you add fewer than j-1 contrasts to the contrasts argument in lm(), it will NOT fill out the remaining contrasts for you. Rather, any group differences other than those represented in your contrast will get lumped into the error term!

Also, note that you could not add more than j-1 contrasts in the contrasts matrix. In practice, when you want to specify your own contrasts, lsmeans/multcomp packages are recommended, which are much more flexible with more adjust options.

I wrote a blog summarizing all scenarios. Below is an example using lsmeans package.

irrigation<-factor(c(rep("Control",10),rep("Irrigated 10 mm",10),rep("Irrigated 20 mm",10)))
biomass<-1:30
fit <- aov(biomass~irrigation)
summary(fit)
plot(irrigation, biomass)
lsm <- lsmeans(fit, ~irrigation)
(ref1 <- lsmeans(fit, c("irrigation")))
#irrigation      lsmean    SE df lower.CL upper.CL
#Control            5.5 0.957 27     3.54     7.46
#Irrigated 10 mm   15.5 0.957 27    13.54    17.46
#Irrigated 20 mm   25.5 0.957 27    23.54    27.46
contrast1 <- list("compare1" = c(0,1,-1), #"Irrigated 10 mm" - "Irrigated 20 mm"
                  "compare2" = c(-1,1,0)) #"Irrigated 10 mm" - "Control"
contrast2 <- list("compare1" = c(0,1,-1)) #"Irrigated 10 mm" - "Irrigated 20 mm"

summary(contrast(ref1, contrast1), adjust = "none")
#contrast estimate   SE df t.ratio p.value
#compare1      -10 1.35 27  -7.385  <.0001
#compare2       10 1.35 27   7.385  <.0001

summary(contrast(ref1, contrast2), adjust = "none")
#contrast estimate   SE df t.ratio p.value
#compare1      -10 1.35 27  -7.385  <.0001



#DIY contrasts matrix1: correct way
contrastmatrix1<-cbind(c(0,1,-1),c(-1,1,0))
contrastmatrix1
mat.temp <- rbind(constant=1/3, t(contrastmatrix1))
mat.temp
mat <- solve(mat.temp)
mat
mat <- mat[ , -1]

contrasts(irrigation)<-mat
contrasts(irrigation)
fit_contrastmatrix1 <- aov(biomass~irrigation)
summary.lm(aov(biomass~irrigation))
#Coefficients:
#  Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  15.5000     0.5528  28.041  < 2e-16 ***
#  irrigation  -10.0000     1.3540  -7.385 6.05e-08 ***
#  irrigation   10.0000     1.3540   7.385 6.05e-08 ***
#Check your contrast matrix
attributes(fit_contrastmatrix1$qr$qr)$contrasts

#DIY contrasts matrix2: correct way
contrastmatrix2<-cbind(c(0,1,-1))
contrastmatrix2
contrasts(irrigation)<-contrastmatrix2
contrasts(irrigation)#Note: fullfill automatically!!
mat.temp <- rbind(constant=1/3, t(contrasts(irrigation)))
mat.temp
mat <- solve(mat.temp)
mat
mat <- mat[ , -1]
contrasts(irrigation)<-mat
contrasts(irrigation)
fit_contrastmatrix2 <- aov(biomass~irrigation)
summary.lm(fit_contrastmatrix2)
#Coefficients:
#  Estimate Std. Error t value Pr(>|t|)    
#(Intercept)  15.5000     0.5528  28.041  < 2e-16 ***
#  irrigation  -10.0000     1.3540  -7.385 6.05e-08 ***
#  irrigation   12.2474     0.9574  12.792 5.67e-13 ***
#Check your contrast matrix
attributes(fit_contrastmatrix2$qr$qr)$contrasts