How to define contrast coefficient matrix?

Question

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!

Answer 1

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