为什么摘要用“无截距”模型公式高估 R 平方答案

【问题标题】：Why does summary overestimate the R-squared with a "no-intercept" model formula为什么摘要用“无截距”模型公式高估 R 平方
【发布时间】：2013-12-18 11:22:48
【问题描述】：

我想制作一个没有截距系数的简单线性模型 (lm())，因此我将 -1 放入我的模型公式中，如下例所示。问题是summary(myModel) 的 R 平方回报似乎被高估了。 lm()、summary() 和 -1 是 R 中非常经典的函数/功能。因此我有点惊讶，我想知道这是一个错误还是这种行为是否有任何原因。

这是一个例子：

x <- rnorm(1000, 3, 1)
mydf <- data.frame(x=x, y=1+x+rnorm(1000, 0, 1))
plot(y ~ x, mydf, xlim=c(-2, 10), ylim=c(-2, 10))

mylm1 <- lm(y ~ x, mydf)
mylm2 <- lm(y ~ x - 1, mydf)

abline(mylm1, col="blue") ; abline(mylm2, col="red")
abline(h=0, lty=2) ; abline(v=0, lty=2)

r2.1 <- 1 - var(residuals(mylm1))/var(mydf$y)
r2.2 <- 1 - var(residuals(mylm2))/var(mydf$y)
r2 <- c(paste0("Intercept - r2: ", format(summary(mylm1)$r.squared, digits=4)),
        paste0("Intercept - manual r2: ", format(r2.1, digits=4)),
        paste0("No intercept - r2: ", format(summary(mylm2)$r.squared, digits=4)),
        paste0("No intercept - manual r2: ", format(r2.2, digits=4)))
legend('bottomright', legend=r2, col=c(4,4,2,2), lty=1, cex=0.6)

【问题讨论】：

这真是一道统计理解题；我将其标记为要迁移到Cross Validated--即 stats.SE（请不要交叉发布，但是）。在此期间，您需要阅读这篇简历：removal of statistically significant intercept term boosts r2 in linear model，尤其是@cardinal 的出色回答。
R 在这种情况下所做的事情对许多人来说并不明显，以至于它需要自己的常见问题解答entry。您可能需要检查并确认您使用的是相同的方法进行比较。
@Gavin 我很高兴我不是唯一一个发现这种 R 行为很奇怪的人，感谢您的链接！我自己发现了困难的方式，但至少我了解 what 会发生（不幸的是不是 why :-)），请参阅我的答案
请注意，不仅 R^2 很奇怪，而且 F-statistic 及其 p 值也很奇怪。
感谢大家分享您的经验和知识。特别要感谢 gung、Gavin Simpson、Dirk Eddelbuettel 和 Tomas，感谢您的 cmets 和回答，这节省了我的时间。我已经习惯了在这个网站上找到我的 R 编程问题的答案，以至于我忘记查看交叉验证和 R 常见问题解答网站。无论如何，我现在绝对清楚这个问题。

标签： r summary intercept lm

【解决方案1】：

哦，对了，我也掉进了这个陷阱！很好的问题！！这是因为

和

在模型 with 截距的情况下（您的 mylm1），y̅ 是平均值（y_i） - 这就是您所期望的，这就是 SStot 你基本上想要适当的 R²
而在模型没有截距的情况下，y̅被取为0 - 所以SS_tot会非常高，所以R^{2 sup> 将非常接近 1！ SS_res 会根据更差的拟合而有所不同（在没有截距的情况下会稍微高一点），但不会太大。}

代码：

attach(mylm1) # in general be careful with attach, here only for code clarity

y_fit <- mylm1$fitted.values
SSE <- sum((y_fit - y)^2)
SST <- sum((y - mean(y))^2)
1-SSE/SST  # R^2 with intercept

y_fit2 <- mylm2$fitted.values
SSE2 <- sum((y_fit2 - y)^2) # SSE2 only slightly higher than SSE..
SST2 <- sum((y - 0)^2)  # !!! the key difference is here !!!
1-SSE2/SST2 # R^2 without intercept

注意：我不清楚为什么在没有截距的模型中 y̅ 为 0 而不是均值 (y_i)，但事实就是如此。通过调查和破解上述代码，我自己发现了一条艰难的道路。

【讨论】：

FAQ 的最后一句话解释了为什么没有截距的模型与零比较而不是均值：“...如果您先验地知道当 x 时 E[Y]=0 =0，那么您应该与拟合线比较的“空”模型，即 x 不解释任何方差的模型，是处处 E[Y]=0 的模型。”
这是加文的链接，很高兴知道官方的答案。不过，这个论点是有争议的：它们不应该隐含地假设 E[Y]=0。当 E[Y]0 并且您不希望截距时存在合法情况（例如，接收所有分类协变量的“合理”系数，请参阅example 1 in this answer）。
我想，所以，是的。但这是编码与意图不匹配的情况。也就是说，即使公式的编码中没有截距项，您也希望与具有截距项的模型进行比较。

【解决方案2】：

你的公式是错误的。以下是我在计算摘要时在RcppArmadillo 中的fastLm.R 中所做的：

## cf src/library/stats/R/lm.R and case with no weights and an intercept 
f <- object$fitted.values    
r <- object$residuals         
mss <- if (object$intercept) sum((f - mean(f))^2) else sum(f^2)     
rss <- sum(r^2)   

r.squared <- mss/(mss + rss) 
df.int <- if (object$intercept) 1L else 0L    

n <- length(f)     
rdf <- object$df     
adj.r.squared <- 1 - (1 - r.squared) * ((n - df.int)/rdf)

您需要在两个地方跟踪是否存在拦截。

【讨论】：