9 votos

Tipo III sumas de cuadrados

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

Puedo configurar los contrastes de códigos en R con options(contrasts=c("contr.sum","contr.poly")). Y ahora he de Tipo III suma de cuadrados de $A$, $B$, y su interacción (a:B) utilizando drop1(model, .~., test="F").

Lo que estoy atascado con es cómo las sumas de cuadrados se calculan para $B$. Yo creo que es sum((predicted y of the full model - predicted y of the reduced model)^2). El modelo reducido se vería y~A+A:B. Pero cuando utilizo predict(y~A+A:B), R es devolver los valores de predicción que son las mismas que el modelo completo los valores de la predicción. Por lo tanto, las sumas de los cuadrados de 0.

(Para las sumas de los cuadrados de $A$, he utilizado un modelo reducido de y~B+A:B, que es el mismo que y~A:B.)

Aquí es el código de ejemplo al azar de los datos generados por:

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    

13voto

Chirag Shekhar Puntos 6

He encontrado diferencias en la estimación de los regresores entre R 2.15.1 y SAS 9.2, pero después de la actualización de la R a la versión 3.0.1 los resultados fueron los mismos. Así, en primer lugar me suggets para la actualización de R para la versión más reciente.

Estás usando el enfoque equivocado porque se trata de calcular la suma de cuadrados en contra de dos modelos diferentes, lo que implica dos diferentes diseño de matrices. Esto llevará a totalmente diferente de la estimación de los regresores utilizados por lm() para calcular los valores pronosticados (estás usando regresores con valores diferentes entre los dos modelos). SS3 se calcula con base en una prueba de hipótesis, suponiendo que todos los de acondicionamiento de regresores igual a cero, mientras que el condicionada de la variable es igual a 1. Para los cálculos se utiliza el mismo diseño de la matriz utilizada para estimar el modelo completo, como para el regresor a la estimada en el modelo completo. Recuerde que el SS3s no están llenos de aditivos. Esto significa que si usted suma la estimación de la SS3, no obtiene el modelo SS (SSM).

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

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

Espero que esta ayuda :)

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

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by: