机器算法验证 - III 型平方和 - 吾爱随笔录

III 型平方和

机器算法验证回归方差分析平方和

2022-03-16 05:29:43

我有一个带有一个分类变量的线性回归模型 $A$ （男性和女性）和一个连续变量 $B$ .

我在 R 中用options(contrasts=c("contr.sum","contr.poly")). 现在我有 III 型平方和 $A$ , $B$ , 以及它们的交互 (A:B) 使用drop1(model, .~., test="F").

我坚持的是如何计算平方和 $B$ . 我认为是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。

（对于平方和 $A$ ，我使用了的简化模型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

1个回答

我发现 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))

其它你可能感兴趣的问题

上一篇用于描述当线条处于较大梯度时看起来彼此更接近的两条线条的术语下一篇从真实数据中创建“演示”数据：伪装而不毁容