# How to define contrast coefficient matrix?

I have this data

data <- as.data.frame(cbind(y,x1,x2,pre))
data$x1 <- as.factor(data$x1)
data$x2 <- as.factor(data$x2)

data$x1 <- C(data$x1, contr.treatment, base=3)
data$x2 <- C(data$x2, contr.treatment, base=2)

y x1 x2 pre
1  16  1  1  14
2  15  1  1  13
3  14  1  2  14
4  13  1  2  13
5  12  2  1  12
6  11  2  1  12
7  11  2  2  13
8  13  2  2  13
9  10  3  1  10
10 11  3  1  11
11 11  3  2  11
12  9  3  2  10

'data.frame':   12 obs. of  4 variables:
$y : num 16 15 14 13 12 11 11 13 10 11 ...$ x1 : Factor w/ 3 levels "1","2","3": 1 1 1 1 2 2 2 2 3 3 ...
..- attr(*, "contrasts")= num [1:3, 1:2] 1 0 0 0 1 0
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$: chr "1" "2" "3" .. .. ..$ : chr  "1" "2"
$x2 : Factor w/ 2 levels "1","2": 1 1 2 2 1 1 2 2 1 1 ... ..- attr(*, "contrasts")= num [1:2, 1] 1 0 .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr  "1" "2"
.. .. ..$: chr "1"$ pre: num  14 13 14 13 12 12 13 13 10 11 ...


And I fitted the following model

lm(y ~ x1 + x2 + x1*x2)


My design matrix 'x' is

x <- as.matrix(cbind(rep(1,12), data\$pre, c(1,1,1,1,0,0,0,0,0,0,0,0),
c(0,0,0,0,1,1,1,1,0,0,0,0), c(1,1,0,0,1,1,0,0,1,1,0,0),
c(1,1,0,0,0,0,0,0,0,0,0,0), c(0,0,0,0,1,1,0,0,0,0,0,0)))

x
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1   14    1    0    1    1    0
[2,]    1   13    1    0    1    1    0
[3,]    1   14    1    0    0    0    0
[4,]    1   13    1    0    0    0    0
[5,]    1   12    0    1    1    0    1
[6,]    1   12    0    1    1    0    1
[7,]    1   13    0    1    0    0    0
[8,]    1   13    0    1    0    0    0
[9,]    1   10    0    0    1    0    0
[10,]    1   11    0    0    1    0    0
[11,]    1   11    0    0    0    0    0
[12,]    1   10    0    0    0    0    0


I'm trying to use this design to reproduce the following table:

Source DF Squares Mean Square F Value Pr > F
Model 6 44.79166667 7.46527778 12.98 0.0064
Error 5 2.87500000 0.57500000
Corrected Total 11 47.66666667

Source DF Type III SS Mean Square F Value Pr > F
pre 1 3.12500000 3.12500000 5.43 0.0671
x1 2 4.58064516 2.29032258 3.98 0.0923
x2 1 3.01785714 3.01785714 5.25 0.0706
x1*x2 2 1.25000000 0.62500000 1.09 0.4055


The first part is fine

XtX <- t(x) %*% x
XtXinv <- solve(XtX)
betahat <- XtXinv %*% t(x) %*% y

H <- x %*% XtXinv %*% t(x)
IH <- (diag(1,12) - H)
yhat <- H %*% y
e <- IH %*% y
ybar <- mean(y)

MSS <- t(betahat) %*% t(x) %*% y - length(y)*(ybar^2)
ESS <- t(e) %*% e
TSS <- MSS + ESS

dfM <- sum(diag(H)) - 1
dfE <- sum(diag(IH))
dfT <- dfM + dfE

MSM <- MSS/dfM
MSE <- ESS/dfE

Ftest <- MSM / MSE
pr <- 1 - pf(Ftest, dfM, dfE)


The contrast coefficient matrix for 'pre' seems correct.

L <- matrix(c(0,1,0,0,0,0,0), 1, 7, byrow=T)
Lb <- L %*% betahat
LXtXinvLt <- round(L %*% XtXinv %*% t(L), digits=4)
SSpre <- t(Lb) %*% solve(LXtXinvLt) %*% (Lb)
MSpre <- SSpre / 1
Fpre <- MSpre / MSE
PRpre <- 1 - pf(Fpre, 1, 12-7)


But I can't understand how to define the contrast coefficient matrix for x1, x2, and x1*x2. What's the problem with the rest of my code? Below an example for how I think I should calculate for x1

L <- matrix(c(0,0,1,1,0,0,0), 1, 7, byrow=T)
Lb <- L %*% betahat
LXtXinvLt <- round(L %*% XtXinv %*% t(L), digits=4)
SSX1 <- t(Lb) %*% solve(LXtXinvLt) %*% (Lb)
MSX1 <- SSX1 / 1
FX1 <- MSX1 / MSE
PRX1 <- 1 - pf(FX1, 1, 12-7)


The result I get for the second part of the ANOVA table is

pre     1    3.1250 6.0000  3.6822  5.4348 0.06711
x1      2    1.0439 3.9189 -3.4291  0.9078 0.46098
x2      1    0.2500 3.1250 -4.1457  0.4348 0.53881
x1:x2   2    1.2500 4.1250 -2.8141  1.0870 0.40554


Thanks!

• Read this question and the answers stats.stackexchange.com/questions/23197/… – mnel Oct 22 '12 at 3:50
• Thanks for your reply, but I don't think that answer helps me. My problem is that I don't know how to define a contrast for a categorical variable and then how to use it in matrix form. I just need help with the last part of my code. – user1172558 Oct 22 '12 at 4:42
• You are attempting to recreate a type III SS table from SAS in R. You should read and understand the answers in the linked question. – mnel Oct 22 '12 at 4:44
• You should not be specifying contr.treatments within your definition of x1 and x2 – mnel Oct 22 '12 at 22:17

The holistic answer to your question is contained here https://stats.stackexchange.com/a/23198/11523

To answer your specific question regarding creating contrast matrices

Once you have set options(contrasts=c("contr.sum","contr.poly"))

then you can use

 contrasts(x1)
##   [,1] [,2]
## 1    1    0
## 2    0    1
## 3   -1   -1
contrasts(x2)
##   [,1]
## 1    1
## 2   -1
contrasts(interaction(x1,x2))
##     [,1] [,2] [,3] [,4] [,5]
## 1.1    1    0    0    0    0
## 2.1    0    1    0    0    0
## 3.1    0    0    1    0    0
## 1.2    0    0    0    1    0
## 2.2    0    0    0    0    1
## 3.2   -1   -1   -1   -1   -1


to create valid contrast matrices and to recreate your tables by hand..

Read the answers (an links) from the linked questions to decide whether type III SS are useful or not.

To avoid the tedious hand-calculations, use the following (from https://stats.stackexchange.com/a/23198/11523)

model <- lm(y ~ pre + x1*x2)

drop1(.model, .~., test = 'F')
Single term deletions

Model:
y ~ pre + x1 * x2
Df Sum of Sq    RSS     AIC F value  Pr(>F)
<none>              2.8750 -3.1462
pre     1    3.1250 6.0000  3.6822  5.4348 0.06711 .
x1      2    4.5806 7.4556  4.2888  3.9832 0.09234 .
x2      1    3.0179 5.8929  3.4660  5.2484 0.07057 .
x1:x2   2    1.2500 4.1250 -2.8141  1.0870 0.40554

• I don't know why, but I can't get this result. Even using drop1() options(contrasts=c("contr.sum","contr.poly")) m2 <- lm(y ~ pre + x1*x2, data) summary(m2) anova(m2) drop1(m2, .~., test='F') # Type III SS The only result which is still correct is the one for 'pre'. – user1172558 Oct 22 '12 at 5:27
• Are you sure that x1 and x2 are coded as factors? Should pre be a factor? – mnel Oct 22 '12 at 5:30
• Yes. x1 and x2 are factors, but not pre. I posted the str(data) above. – user1172558 Oct 22 '12 at 5:35