计算R中的回归平方和

机器算法验证 r 回归 多重回归
2022-04-05 13:09:11

这是示例数据

        brainIQ <- 
      read.table (file= "https://onlinecourses.science.psu.edu/stat501/sites/
  onlinecourses.science.psu.edu.stat501/files/data/iqsize.txt",
     head = TRUE)

我正在尝试拟合多元线性回归。

mylm <- lm(PIQ ~  Brain + Height + Weight, data = brainIQ)
anova(mylm)

R 中的默认函数 anova 提供顺序平方和(类型 I)平方和。

Analysis of Variance Table

Response: PIQ
          Df  Sum Sq Mean Sq F value  Pr(>F)  
Brain      1  2697.1 2697.09  6.8835 0.01293 *
Height     1  2875.6 2875.65  7.3392 0.01049 *
Weight     1     0.0    0.00  0.0000 0.99775  
Residuals 34 13321.8  391.82                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我相信,因此SS是Brain,Height | 大脑,体重 | (大脑,重量)和残差分别。

使用包车我们也可以得到类型 II 的平方和。

library(car)
Anova(mylm, type="II")
Anova Table (Type II tests)

Response: PIQ
           Sum Sq Df F value    Pr(>F)    
Brain      5239.2  1 13.3716 0.0008556 ***
Height     1934.7  1  4.9378 0.0330338 *  
Weight        0.0  1  0.0000 0.9977495    
Residuals 13321.8 34                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这里的平方和如下:Brian | (身高、体重)、身高 | (大脑,体重),体重 | (大脑,身高)。

看起来很像 Mintab 输出:

在此处输入图像描述

我的问题是如何计算 R 上表中的回归行?

3个回答
SS(Regression) = SS(Total) - S(Residual)

您可以通过以下方式获得 SS(Total):

SSTotal <- var( brainIQ$PIQ ) * (nrow(brainIQ)-1)
    SSE     <- sum( mylm$resid^2 )
SSreg   <- SSTotal - SSE

“回归”行的自由度是回归的相应组件(在本例中为:大脑、身高和体重)的自由度之和。

然后得到其余的:

dfE   <- mylm$df.residual
dfReg <- nrow(brainIQ) - 1 - dfE
MSreg <- SSreg / dfReg
MSE   <- SSE / dfE
Fstat <- MSreg / MSE
pval  <- pf( Fstat , dfReg, dfE , lower.tail=FALSE )

如果您已经建立了线性模型,则可以用一条线计算回归平方和。使用您的模型:

sum((mylm$fitted.values - mean(mylm$fitted.values))^2)

这利用了响应的平均值等于拟合值的平均值这一事实。

更新:正如@whuber 指出的那样,上面的行假设模型在其设计矩阵的跨度中包含一个常数(“截距”)。下面的行不做这个假设,并且会更普遍地工作:

sum((mylm$fitted.values - mean(mylm$model$y))^2)

我找到了解决这个问题的矩阵方法:

在这里我们可以计算 SSR 和总平方和。

    # Y matrix 
    Y <- as.matrix(brainIQ$PIQ, ncol=1)
    n= nrow(Y)
    J = matrix(1, nrow=n, ncol=n)
    # Total sum of Square 
    SSTO = t(Y) %*% Y - (1/n)*t(Y)%*%J%*%Y 
    X <- as.matrix(cbind(1, brainIQ[,-1]))
    b = solve(t(X)%*%X)%*%t(X)%*%Y
# regression sum of square 
    SSR = t(b)%*%t(X)%*%Y - (1/n)%*%t(Y)%*%J%*%Y
> SSR
         [,1]
[1,] 5572.744

> SSTO
         [,1]
[1,] 18894.55

这里 - SSR(Brain, Height, Weight) = SSR(Brain) + SSR(Height|Brain) + SSR (Weight|Height, Brain) - 这实际上是来自 anova 的平方输出的顺序总和。

mylm <- lm(PIQ ~  Brain + Height + Weight, data = brainIQ)

# SSR(Brain, Height, Weight)
sum(anova(mylm)[-4,2])

# Total sum of square 
sum(anova(mylm)[,2])