¿Cómo se define la matriz del coeficiente de contraste?

Question

Tengo estos datos

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 ...

Y ajusté el siguiente modelo

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

Mi matriz de diseño 'x' es

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

Estoy tratando de utilizar este diseño para reproducir la siguiente tabla:

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

La primera parte está bien

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)

La matriz del coeficiente de contraste para "pre" parece correcta.

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)

Pero no consigo entender cómo definir la matriz de coeficientes de contraste para x1, x2 y x1*x2. ¿Cuál es el problema con el resto de mi código? A continuación un ejemplo de cómo creo que debo calcular para 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) 

El resultado que obtengo para la segunda parte de la tabla ANOVA es

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 

Gracias.

Answer 1

La respuesta integral a su pregunta está contenida aquí https://stats.stackexchange.com/a/23198/11523

Para responder a su pregunta específica sobre la creación de matrices de contraste

Una vez que haya establecido options(contrasts=c("contr.sum","contr.poly"))

entonces puede utilizar

 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

para crear matrices de contraste válidas y para recrear sus tablas a mano..

Lea las respuestas (y los enlaces) de las preguntas enlazadas para decidir si las SS de tipo III son útiles o no.

Para evitar los tediosos cálculos manuales, utilice lo siguiente (de 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