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

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

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

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

El modelo con el que está trabajando toma la forma

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

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

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

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

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

Si su modelo se ajusta perfectamente a sus datos, los valores pronosticados son valores reales. Pero, en general, los valores de no pueden obtenerse exactamente como una combinación lineal simple de valores de (" Todos los modelos están equivocados, 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 (o una versión a escala de estos) son "pequeños".$y$$x$$y_{i} - y_{i}^{\star}$

Editar

En tus comentarios, preguntaste qué predict()hace realmente. Aquí hay un ejemplo ilustrativo simple.

    #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