Regression

III 型平方和

  • April 19, 2013

我有一個帶有一個分類變量的線性回歸模型(男性和女性)和一個連續變量.

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

引用自:https://stats.stackexchange.com/questions/56625

comments powered by Disqus