Tipe III jumlah kuadrat

Saya memiliki model regresi linear dengan satu variabel kategoris (laki-laki & perempuan) dan satu variabel kontinu .$A$$B$

Saya mengatur kode kontras dalam R dengan options(contrasts=c("contr.sum","contr.poly")). Dan sekarang saya memiliki tipe III jumlah kuadrat untuk , , dan interaksinya (A: B) menggunakan .$A$$B$drop1(model, .~., test="F")

Apa yang saya terjebak dengan adalah bagaimana jumlah kuadrat dihitung untuk . Saya$B$ pikir begitu sum((predicted y of the full model - predicted y of the reduced model)^2). Model yang dikurangi akan terlihat seperti y~A+A:B. Tetapi ketika saya gunakan predict(y~A+A:B), R mengembalikan nilai prediksi yang sama dengan nilai prediksi model penuh. Oleh karena itu, jumlah kuadrat akan menjadi 0.

(Untuk jumlah kuadrat , saya menggunakan model tereduksi , yang sama dengan .)$A$y~B+A:By~A:B

Berikut adalah contoh kode untuk data yang dihasilkan secara acak:

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    

Saya telah menemukan perbedaan dalam estimasi regressor antara R 2.15.1 dan SAS 9.2, tetapi setelah memperbarui versi R ke 3.0.1 hasilnya sama. Jadi, pertama saya menyarankan Anda untuk memperbarui R ke versi terbaru.

Anda menggunakan pendekatan yang salah karena Anda menghitung jumlah kuadrat terhadap dua model yang berbeda, yang menyiratkan dua matriks desain yang berbeda. Ini membawa Anda ke estimasi yang sama sekali berbeda dalam regressor yang digunakan oleh lm () untuk menghitung nilai yang diprediksi (Anda menggunakan regressor dengan nilai yang berbeda antara kedua model). SS3 dihitung berdasarkan uji hipotesis, dengan asumsi bahwa semua regresi pengkondisian sama dengan nol, sedangkan regresi terkondisi sama dengan 1. Untuk perhitungan, Anda menggunakan matriks desain yang sama yang digunakan untuk memperkirakan model lengkap, seperti untuk regresi yang diperkirakan secara penuh. model. Ingat bahwa SS3 tidak aditif penuh. Ini berarti bahwa jika Anda menjumlahkan perkiraan SS3, Anda tidak mendapatkan model SS (SSM).

Di sini saya menyarankan implementasi R dari matematika yang mengimplementasikan algoritma GLS yang digunakan untuk memperkirakan SS3 dan regressor.

Nilai yang dihasilkan oleh kode ini sama persis dengan yang dihasilkan menggunakan SAS 9.2 seperti untuk hasil yang Anda berikan dalam kode Anda, sedangkan SS3 (B | A, AB) adalah 0,167486 bukannya 0,15075. Untuk alasan ini saya sarankan lagi untuk memperbarui versi R Anda ke versi terbaru yang tersedia.

Semoga ini membantu :)

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