¿Cuáles son los valores predichos devueltos por la función predict() en R cuando se utilizan los datos originales como entrada?

Question

Después de ejecutar una regresión de la forma reg <- lm(y ~ x1 + x2, data=example) en un conjunto de datos, puedo obtener valores predichos utilizando

predict(reg, example, interval="prediction", level=0.95)

Me pregunto a qué se refieren realmente los valores predichos cuando utilizo la regresión para predecir el conjunto de datos real. ¿No debería obtener los valores originales?

Answer 1

El modelo con el que se trabaja tiene la forma

$y_{i} = \mu + \beta_{1} x_{1i} + \beta_{2} x_{2i} + \epsilon_{i}$ $\hspace{0.75cm}$ (1)

donde $\epsilon_{i}$ es un término de error que se supone proviene de una distribución normal de media cero.

Ha ajustado el modelo y ha obtenido estimaciones: $\hat{\mu}$ , $\hat{\beta}_{1}$ y $\hat{\beta}_{2}$ .

Ahora, si se fijan los valores de las covariables dentro de su rango, digamos $x^{\star}_{1i}$ y $x^{\star}_{2i}$ un valor previsto para $y_{i}$ se puede obtener calculando

$y^{\star}_{i} = \hat{\mu} + \hat{\beta}_{1} x^{\star}_{1i} + \hat{\beta}_{2} x^{\star}_{2i}$ $\hspace{0.75cm}$ (2)

Si el modelo se ajusta perfectamente a los datos, los valores predichos son valores reales. Pero, en general, $y$ no pueden obtenerse exactamente como una simple combinación lineal de $x$ valores (" Todos los modelos son erróneos, pero algunos son útiles "). En otros términos, la varianza del término de error en (1) no es cero en general. Pero, básicamente, el modelo (1) es una buena aproximación si los residuos $y_{i} - y_{i}^{\star}$ (o una versión a escala de éstos) son "pequeños".

Editar

En sus comentarios, ha preguntado qué predict() realmente lo hace. He aquí un sencillo ejemplo ilustrativo.

    #generate a simple illustrative data set
> x <- runif(10)
> y <- 5 + 2.7 * x + rnorm(10, mean=0, sd=sqrt(0.15))
> 
>   #fit the model and store the coefficients
> regLin <- lm(y~x)
> coef <- coef(regLin)
> 
>   #use the predict() function
> y_star2 <- predict(regLin)
>   #use equation (2)
> y_star1 <- coef[1] + coef[2] * x
>   #compare
> cbind(y, y_star1, y_star2) 
          y  y_star1  y_star2
1  7.100217 6.813616 6.813616
2  6.186333 5.785473 5.785473
3  7.141016 7.492979 7.492979
4  5.121265 5.282990 5.282990
5  4.681924 4.849776 4.849776
6  6.102339 6.106751 6.106751
7  7.223215 7.156512 7.156512
8  5.158546 5.253380 5.253380
9  7.160201 7.198074 7.198074
10 5.555289 5.490793 5.490793