Soma de quadrados do tipo III

9

Eu tenho um modelo de regressão linear com uma variável categórica (masculino e feminino) e uma variável contínua .UMAB

Eu configurei códigos de contraste em R com options(contrasts=c("contr.sum","contr.poly")). E agora eu tenho somas de quadrados do Tipo III para , e sua interação (A: B) usando .UMABdrop1(model, .~., test="F")

O que eu estou preso com é como somas de quadrados é calculada para . EuB acho que sim sum((predicted y of the full model - predicted y of the reduced model)^2). O modelo reduzido seria semelhante y~A+A:B. Mas quando eu uso predict(y~A+A:B), R está retornando valores previstos iguais aos valores previstos do modelo completo. Portanto, a soma dos quadrados seria 0.

(Para as somas dos quadrados de , usei um modelo reduzido de , que é o mesmo que .)UMAy~B+A:By~A:B

Aqui está um código de exemplo para dados gerados aleatoriamente:

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    
Jo Lewis
fonte
11
Essa é uma boa pergunta e tenho algumas idéias sobre como uma resposta pode ser. Mas sem um exemplo reproduzível, não estou investindo meu tempo. OP, entregue!
Henrik
11
O que faz você querer testes do tipo III ("Senado dos EUA") em oposição aos testes do tipo II ("Câmara dos Deputados dos EUA")? (analogias devido a Paul Gallo, Novartis)
Frank Harrell
o código ajuda?
Jo Lewis
Pergunta relacionada: Como interpretar ANOVA e MANOVA tipo I (seqüencial)?
gung - Restabelece Monica

Respostas:

3

Eu encontrei diferenças na estimativa de regressores entre o R 2.15.1 e o SAS 9.2, mas após a atualização do R para a versão 3.0.1, os resultados foram os mesmos. Portanto, primeiro sugiro que você atualize o R para a versão mais recente.

Você está usando a abordagem errada porque está calculando a soma do quadrado em relação a dois modelos diferentes, o que implica duas matrizes de design diferentes. Isso leva a uma estimativa totalmente diferente nos regressores usados ​​por lm () para calcular os valores previstos (você está usando regressores com valores diferentes entre os dois modelos). O SS3 é calculado com base em um teste de hipotese, assumindo que todos os regressores condicionantes são iguais a zero, enquanto o regressor condicionado é igual a 1. Para os cálculos, você usa a mesma matriz de projeto usada para estimar o modelo completo, assim como para o regressor estimado na totalidade. modelo. Lembre-se de que os SS3s não são totalmente aditivos. Isso significa que, se você soma o SS3 estimado, não obtém o modelo SS (SSM).

Aqui, sugiro uma implementação R da matemática que implementa o algoritmo GLS usado para estimar SS3 e regressores.

Os valores gerados por esse código são exatamente os mesmos gerados usando o SAS 9.2 como para os resultados que você forneceu no seu código, enquanto o SS3 (B | A, AB) é 0,166486 em vez de 0,15075. Por esse motivo, sugiro novamente atualizar sua versão R para a mais recente disponível.

Espero que isto ajude :)

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