为什么摘要会使用“非拦截”模型公式高估R平方

SES*_*man 4 r summary intercept lm

我想创建一个lm()没有截距系数的简单线性模型(),所以我将-1模型公式输入如下例所示。问题在于R平方的回报summary(myModel)似乎被高估了。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)
Run Code Online (Sandbox Code Playgroud)

在此处输入图片说明

TMS*_*TMS 6

哦,是的,我也掉进了这个陷阱!很好的问题!!这是因为

在此处输入图片说明

在此处输入图片说明

  • 具有截距(your mylm1)的模型的情况下,y̅是均值(y i)-这是您所期望的,这是您基本上想要获得适当R 2的SS 总和
  • 而在没有截距的模型中,y̅取为0-因此SS tot将非常高,因此R 2将非常接近1!SS的分辨率会根据最差的拟合度而有所不同(如果没有截距,则较高),但差异不大。

码:

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
Run Code Online (Sandbox Code Playgroud)

注意:对我来说尚不清楚为什么在没有截距的模型中y̅为0而不是mean(y i),但是事实就是如此。我本人通过调查和破解以上代码发现了困难的方法。

  • [FAQ](http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-does- summary_0028_0029-report-strange-results-for-the-R_005e2-estimate-when-I-fit-a-linear-model-with-no-intercept_003f):“...如果你知道E[Y]=的先验当 x=0 时为 0,那么您应该与拟合线进行比较的“空”模型,即 x 不解释任何方差的模型,就是 E[Y]=0 处处的模型。” (2认同)

Dir*_*tel 2

你的公式是错误的。以下是我在RcppArmadillofastLm.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)      
Run Code Online (Sandbox Code Playgroud)

有两个地方需要跟踪是否有拦截。