3

我在玩 R 的函数，在计算决定系数 R^2lm()时发现了一些奇怪的东西。所以假设我有一个玩具回归问题：

set.seed(0); N= 12;
x_ <- 1:N; y_ <- 2* x_+ rnorm(N, sd=2); 
lm1_ <- lm( y_ ~ x_ ); S1 <- summary(lm1_); 
lm2_ <- lm( y_ ~ -1 + model.matrix(lm1_) ); S2 <- summary(lm2_);

所以在这一点上，我们实际上将 lm1_ 和 lm2_ 设置为同一件事。使用相同的因变量，使用相同的模型矩阵，我们应该得到相同的精确拟合。让我们检查一下：

identical( fitted(lm2_),fitted(lm1_))       
[1] TRUE                                    #SURE!
identical( residuals(lm2_), residuals(lm1_))
[1] TRUE                                    #YEAP!
identical( lm1_$df.residual, lm2_$df.residual)
[1] TRUE                                    #SORTED!
identical( as.numeric(model.matrix(lm2_ )), as.numeric(model.matrix(lm1_ ))) 
[1] TRUE                                    #WHY NOT?
identical( S2$r.squared, S1$r.squared)        
[1] FALSE                                   #HUH?

发生了什么？让我们使用维基百科文章中的定义来确定系数（或如果您愿意，可以解释方差百分比）：

1 - ((sum( residuals(lm1_)^2))/( sum( (y_ -mean(y_))^2)))
[1] 0.8575376
1 - ((sum( residuals(lm2_)^2))/( sum( (y_ -mean(y_))^2)))
[1] 0.8575376
identical(1 - ((sum( residuals(lm1_)^2))/(sum( (y_ -mean(y_))^2))), S1$r.squared)
[1] TRUE
identical(1 - ((sum( residuals(lm2_)^2))/(sum( (y_ -mean(y_))^2))), S1$r.squared)
[1] TRUE
identical(1 - ((sum( residuals(lm2_)^2))/(sum( (y_ -mean(y_))^2))), S2$r.squared)
[1] FALSE
S2$r.squared
[1] 0.9700402 # ???

现在显然整个把戏都被这个函数弄糊涂了summary.lm()。R^2 计算为残差平方和的mss/ (rss+ mss)位置，并且是取决于截距的存在与否的平方均和（拟合）（来自 summary.lm() 的 cp）：rsssum(residuals(lmX_)^2)mss

mss <- if (attr(z$terms, "intercept")) 
            sum((f - mean(f))^2)
        else sum(f^2)

或者在我们的例子中：

sum(scale(fitted(lm2_),scale=F)^2) / 
(sum(scale(fitted(lm2_),scale=F)^2)+sum(residuals(lm2_)^2))
[1] 0.8575376
sum(fitted(lm2_)^2) / (sum(fitted(lm2_)^2) + sum(residuals(lm2_)^2))
[1] 0.9700402

很明显，R 没有将第二个模型矩阵识别为具有截距。但是，如果像我一样有人在移动模型矩阵，那将是一个重大问题。我知道显而易见的解决方法是采用没有截距的模型矩阵，但这并不能缓解这样一个事实，即这是一个有点容易出错的设计选择。更一般的定义难道不是只需稍作调整就可以为我们解决这个麻烦吗？我可以使用其他一些明显的解决方法吗？我不会遇到这些问题吗？

如果没有截距并且使用具有交互作用的因子变量，这个问题会变得非常严重。实际上，R 可能会尝试将其中一个级别视为截距。那么整个情况就悲剧了：

library(MASS)
lm_oats <- lm( Y ~ -1+  V*N , oats)
S_oats <- summary(lm_oats)
S_oats$r.squared
[1] 0.9640413 # :*D
1-  ( sum( residuals(lm_oats)^2)  / sum( ( oats$Y- mean(oats$Y))^2 )  )
[1] 0.4256653 # :*(

#For the record:
sum((fitted(lm_oats))^2 )/(sum( residuals(lm_oats)^2) +sum((fitted(lm_oats))^2)) 
[1] 0.9640413
sum((scale(scale=F,fitted(lm_oats)))^2 ) /( sum( residuals(lm_oats)^2) 
+sum((scale(scale=F,fitted(lm_oats)))^2 )) 
[1] 0.4256653

这是 R 的设计特征，可以以某种方式被规避吗？（我很感激我的问题非常开放，感谢所有一路下来的人！）

我正在使用 R 版本 3.0.1、Matrix_1.0-12、lattice_0.20-15、MASS_7.3-26。（Linux 32 位）

4

1 回答 1

4

R 邮件列表上关于 R^2 对被迫通过拦截的模型的不敏感/不同含义进行了长时间的讨论。R 尝试检测您何时拟合没有截距的模型（“语义上”，即通过查找模型公式中的-1或存在+0）并相应地修改 R^2 计算。如果你做了一些棘手的事情，你就会绕过它的尝试。

进一步来说：

  • summary.lm()测试terms拟合模型的元素是否具有intercept等于 1 的属性
  • 拟合模型的terms元素是通过调用构建的model.frame()
  • model.frame()来电terms.formula
  • terms.formula调用了一些非常繁琐的内部 C 代码C_termsform，但是您可以通过尝试terms(y~x-1)vs terms(y~x+0)vsterms(y~x)看到 R 只是在语义上进行检查（正如我上面建议的......）

底线：如果您想对围绕 R 的正常约定工作的模型矩阵做一些花哨的事情，您可能应该继续根据您自己的首选公式计算 R^2；您可以从中复制相关代码summary.lm()。如果您要深入研究，您可能希望使用它lm.fit()来提高计算效率。

于 2013-06-03T21:39:26.010 回答