Sommes de carrés de type III

J'ai un modèle de régression linéaire avec une variable catégorique (mâle et femelle) et une variable continue .$A$$B$

J'ai mis en place des codes de contrastes en R avec options(contrasts=c("contr.sum","contr.poly")). Et maintenant, j'ai des sommes de carrés de type III pour , et leur interaction (A: B) en utilisant .$A$$B$drop1(model, .~., test="F")

Ce que je suis coincé est de savoir comment la somme des carrés est calculé pour . Je$B$ pense que oui sum((predicted y of the full model - predicted y of the reduced model)^2). Le modèle réduit ressemblerait y~A+A:B. Mais lorsque j'utilise predict(y~A+A:B), R renvoie des valeurs prédites qui sont les mêmes que les valeurs prédites du modèle complet. Par conséquent, les sommes des carrés seraient égales à 0.

(Pour les sommes de carrés de , j'ai utilisé un modèle réduit de , qui est le même que .)$A$y~B+A:By~A:B

Voici un exemple de code pour des données générées aléatoirement:

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    

J'ai trouvé des différences dans l'estimation des régresseurs entre R 2.15.1 et SAS 9.2, mais après la mise à jour de la version R vers 3.0.1, les résultats étaient les mêmes. Donc, je vous suggère d'abord de mettre à jour R vers la dernière version.

Vous utilisez la mauvaise approche car vous calculez la somme des carrés par rapport à deux modèles différents, ce qui implique deux matrices de conception différentes. Cela vous amène à une estimation totalement différente dans les régresseurs utilisés par lm () pour calculer les valeurs prédites (vous utilisez des régresseurs avec des valeurs différentes entre les deux modèles). SS3 est calculé sur la base d'un test d'hypotèse, en supposant que tous les régresseurs de conditionnement sont égaux à zéro, tandis que le régresseur conditionné est égal à 1. Pour les calculs, vous utilisez la même matrice de conception utilisée pour estimer le modèle complet, que pour le régresseur estimé dans son intégralité. modèle. N'oubliez pas que les SS3 ne sont pas complètement additifs. Cela signifie que si vous additionnez le SS3 estimé, vous n'obtenez pas le modèle SS (SSM).

Ici, je suggère une implémentation R des mathématiques qui implémente l'algorithme GLS utilisé pour estimer SS3 et les régresseurs.

Les valeurs générées par ce code sont exactement les mêmes générées à l'aide de SAS 9.2 que pour les résultats que vous avez donnés dans votre code, tandis que le SS3 (B | A, AB) est 0.167486 au lieu de 0.15075. Pour cette raison, je suggère à nouveau de mettre à jour votre version R vers la dernière version disponible.

J'espère que cela t'aides :)

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