Suma de cuadrados de tipo III

Question

Tengo un modelo de regresión lineal con una variable categórica $A$ (hombre y mujer) y una variable continua $B$ .

Configuré los códigos de contraste en R con options(contrasts=c("contr.sum","contr.poly")) . Y ahora tengo sumas de cuadrados de Tipo III para $A$ , $B$ y su interacción (A:B) usando drop1(model, .~., test="F") .

Lo que me queda es cómo se calculan las sumas de los cuadrados para $B$ . I creo que es sum((predicted y of the full model - predicted y of the reduced model)^2) . El modelo reducido se vería como y~A+A:B . Pero cuando uso predict(y~A+A:B) R devuelve valores predichos que son los mismos que los del modelo completo. Por lo tanto, las sumas de los cuadrados serían 0.

(Por las sumas de los cuadrados de $A$ usé un modelo reducido de y~B+A:B que es lo mismo que y~A:B .)

He aquí un ejemplo de código para los datos generados al azar:

A<-as.factor(rep(c("male","female"), each=5))
set.seed(1)
B<-runif(10)
set.seed(5)
y<-runif(10)

model<-lm(y~A+B+A:B)

options(contrasts = c("contr.sum","contr.poly"))

#type3 sums of squares
drop1(model, .~., test="F")
#or same result:
library(car)
Anova(lm(y~A+B+A:B),type="III")

#full model
predFull<-predict(model)

#Calculate sum of squares
#SS(A|B,AB)
predA<-predict(lm(y~B+A:B))
sum((predFull-predA)^2) 

#SS(B|A,AB) (???)
predB<-predict(lm(y~A+A:B))
sum((predFull-predB)^2) 
#Sums of squares should be 0.15075 (according to anova table)
#but calculated to be 2.5e-31

#SS(AB|A,B)
predAB<-predict(lm(y~A+B))
sum((predFull-predAB)^2)

#Anova Table (Type III tests)
#Response: y
#             Sum Sq Df F value Pr(>F)
#(Intercept) 0.16074  1  1.3598 0.2878
#A           0.00148  1  0.0125 0.9145
#B           0.15075  1  1.2753 0.3019
#A:B         0.01628  1  0.1377 0.7233
#Residuals   0.70926  6

Answer 1

He encontrado diferencias en la estimación de los regresores entre R 2.15.1 y SAS 9.2, pero después de actualizar R a la versión 3.0.1 los resultados fueron los mismos. Así que, primero le sugiero que actualice R a la última versión.

Estás usando el enfoque equivocado porque estás calculando la suma de los cuadrados contra dos modelos diferentes, lo que implica dos matrices de diseño diferentes. Esto te lleva a una estimación totalmente diferente en los regresores utilizados por lm() para calcular los valores predichos (estás utilizando regresores con valores diferentes entre los dos modelos). SS3 se calcula en base a una prueba de hipótesis, asumiendo que todos los regresores condicionados son iguales a cero, mientras que el regresor condicionado es igual a 1. Para los cálculos se utiliza la misma matriz de diseño utilizada para estimar el modelo completo, como para el regresor estimado en el modelo completo. Recuerda que los SS3 no son aditivos completos. Esto significa que si sumas los SS3 estimados, no obtienes el modelo SS (SSM).

Aquí sugiero una implementación R de las matemáticas que implementa el algoritmo GLS utilizado para estimar SS3 y los regresores.

Los valores generados por este código son exactamente los mismos generados usando SAS 9.2 que los resultados que diste en tu código, mientras que el SS3(B|A,AB) es 0.167486 en lugar de 0.15075. Por esta razón sugiero de nuevo que actualice su versión R a la última disponible.

Espero que esto ayude :)

A<-as.factor(rep(c("male","female"), each=5))
set.seed(1)
B<-runif(10)
set.seed(5)
y<-runif(10)

# Create a dummy vector of 0s and 1s
dummy <- as.numeric(A=="male")

# Create the design matrix
R <- cbind(rep(1, length(y)), dummy, B, dummy*B)

# Estimate the regressors
bhat <- solve(t(R) %*% R) %*% t(R) %*% y
yhat <- R %*% bhat
ehat <- y - yhat

# Sum of Squares Total
# SST <- t(y)%*%y - length(y)*mean(y)**2
# Sum of Squares Error
# SSE <- t(ehat) %*% ehat
# Sum of Squares Model
# SSM <- SST - SSE

# used for ginv()
library(MASS)

# Returns the Sum of Squares of the hypotesis test contained in the C matrix
SSH_estimate <- function(C)
{
    teta <- C%*%bhat
    M <- C %*% ginv(t(R)%*%R) %*% t(C)
    SSH <- t(teta) %*% ginv(M) %*% teta
    SSH
}

# SS(A|B,AB)
# 0.001481682
SSH_estimate(matrix(c(0, 1, 0, 0), nrow=1, ncol=4))
# SS(B|A,AB)
# 0.167486
SSH_estimate(matrix(c(0, 0, 1, 0), nrow=1, ncol=4))
# SS(AB|A,B)
# 0.01627824
SSH_estimate(matrix(c(0, 0, 0, 1), nrow=1, ncol=4))