III 型平方和
我有一個帶有一個分類變量的線性回歸模型(男性和女性)和一個連續變量.
我在 R 中用
options(contrasts=c("contr.sum","contr.poly"))
. 現在我有 III 型平方和,, 以及它們的交互 (A:B) 使用drop1(model, .~., test="F")
.我堅持的是如何計算平方和. 我認為是
sum((predicted y of the full model - predicted y of the reduced model)^2)
。簡化後的模型看起來像y~A+A:B
。但是當我使用 時predict(y~A+A:B)
,R 會返回與完整模型預測值相同的預測值。因此,平方和為 0。(對於平方和,我使用了 的簡化模型
y~B+A:B
,與 相同y~A:B
。)以下是隨機生成數據的示例代碼:
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
我發現 R 2.15.1 和 SAS 9.2 之間的回歸量估計存在差異,但在將 R 更新到 3.0.1 版本後,結果是相同的。因此,首先我建議您將 R 更新到最新版本。
您使用了錯誤的方法,因為您正在計算兩個不同模型的平方和,這意味著兩個不同的設計矩陣。這會導致您對 lm() 用於計算預測值的回歸量進行完全不同的估計(您在兩個模型之間使用具有不同值的回歸量)。SS3 是基於假設檢驗計算的,假設所有條件回歸量都等於 0,而條件回歸量等於 1。對於計算,您使用用於估計完整模型的相同設計矩陣,就像在完整模型中估計的回歸量一樣模型。請記住,SS3 不是完全添加劑。這意味著如果對估計的 SS3 求和,則不會獲得模型 SS (SSM)。
在這裡,我建議使用 R 實現用於估計 SS3 和回歸量的 GLS 算法的數學實現。
此代碼生成的值與使用 SAS 9.2 生成的值與您在代碼中給出的結果完全相同,而 SS3(B|A,AB) 是 0.167486 而不是 0.15075。出於這個原因,我再次建議將您的 R 版本更新到可用的最新版本。
希望這可以幫助 :)
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))