¿Por qué cambia el R-cuadrado con un modelo sin intercepción?

Question

Supongamos que tengo un predictor binario $x$ y ajusté un modelo OLS:

$$y = \alpha + \beta I(x=1) + \epsilon$$

Como alternativa, puedo ajustar el siguiente modelo

$$y = \beta I(x=1) + \beta I(x=0) + \epsilon$$

¿Por qué el segundo modelo me da un R-cuadrado más alto?

R ejemplo:

results <- numeric(1000)
for(i in 1:1000) {
  X <- rnorm(1000)
  b <- X > 0

  y <- rep(0, 1000)
  y[b == TRUE] <- 1

  y <- y + rnorm(1000)

  rsq0 <- summary(lm(y ~ b))$r.squared
  rsq1 <- summary(lm(y ~ b - 1))$r.squared

  results[i] <- rsq1 - rsq0
}

mean(results > 0) # This equals 1

Answer 1

Como escribe Roland se reduce a la definición de $R^2$ en los modelos con y sin intercepción. Más concretamente, la diferencia radica en el cálculo de la suma media de los cuadrados.

Podemos inspeccionar el código para summary() simplemente escribiendo summary.lm (sin () ). Aquí están las partes relevantes:

> summary.lm
function (object, correlation = FALSE, symbolic.cor = FALSE, 
    ...) 
{
    z <- object
    ...
    r <- z$residuals
    f <- z$fitted.values
    ...
        mss <- if (attr(z$terms, "intercept")) 
            sum((f - mean(f))^2)
        else sum(f^2)
        rss <- sum(r^2)
    ...
    ans <- z[c("call", "terms", if (!is.null(z$weights)) "weights")]
    ...
        ans$r.squared <- mss/(mss + rss)

    ans
}

El punto crucial es la distinción en attr(z$terms, "intercept")) que se evalúa en 1 (convertido a TRUE ) en su primer modelo, pero 0 (o FALSE ) en el segundo. En el primer caso, mss se calcula como la suma de las diferencias al cuadrado entre los ajustes y la media global, y en el segundo caso, es sólo la suma de los valores ajustados al cuadrado. De este modo, obtenemos valores diferentes para mss . El resto del cálculo es idéntico.

Podemos calcular un ejemplo y reconstruir las diferentes $R^2$ s a mano:

> nn <- 10
> set.seed(1)
> 
> X <- rnorm(nn)
> b <- X > 0
> y <- rep(0, nn)
> y[b == TRUE] <- 1
> y <- y + rnorm(nn)
>   
> model0 <- lm(y ~ b)
> model1 <- lm(y ~ b - 1)
> 
> mss0 <- sum((model0$fitted.values-mean(model0$fitted.values))^2)
> mss1 <- sum(model0$fitted.values^2)
> rss0 <- sum(model0$residuals^2)
> rss1 <- sum(model1$residuals^2)
> mss0/(mss0 + rss0)
[1] 0.1350045
> (rsq0 <- summary(model0)$r.squared)
[1] 0.1350045
> mss1/(mss1 + rss1)
[1] 0.4628321
> (rsq1 <- summary(model1)$r.squared)
[1] 0.4628321