R中决定系数的计算问题

Question

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

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)，其中rss 是残差平方和sum(residuals(lmX_)^2)，mss 是平方均和（拟合），取决于截距的存在与否（cp来自summary.lm())：

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 位）

Answer 1

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

更具体地说：

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

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