1
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Read this question and the answers stats.stackexchange.com/questions/23197/… $\endgroup$ – mnel Oct 22 '12 at 3:50
  • $\begingroup$ 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. $\endgroup$ – user1172558 Oct 22 '12 at 4:42
  • $\begingroup$ 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. $\endgroup$ – mnel Oct 22 '12 at 4:44
  • $\begingroup$ You should not be specifying contr.treatments within your definition of x1 and x2 $\endgroup$ – mnel Oct 22 '12 at 22:17
1
$\begingroup$

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  
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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'. $\endgroup$ – user1172558 Oct 22 '12 at 5:27
  • 1
    $\begingroup$ Are you sure that x1 and x2 are coded as factors? Should pre be a factor? $\endgroup$ – mnel Oct 22 '12 at 5:30
  • $\begingroup$ Yes. x1 and x2 are factors, but not pre. I posted the str(data) above. $\endgroup$ – user1172558 Oct 22 '12 at 5:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.