Error estándar para predicciones individuales de regresión lineal - ¿Qué diablos es eso?

Question

Estoy tratando de leer el libro "Una inferencia estadística para la era informática", que está disponible aquí: https://web.stanford.edu/~hastie/CASI_files/PDF/casi.pdf

Al comienzo del libro (desde la página 4) hay un ejemplo de regresión lineal con estimaciones de "error estándar" de las predicciones, y me ha confundido completamente.

La configuración es la siguiente. Tienen un modelo de regresión lineal,

$$= \beta_0 + \beta_1x$$

y utilizando "mínimos cuadrados" deducen que los parámetros tienen los valores

$$\beta_0 = 2.86, \beta_1 = -0.079$$

Vale. Luego tienen una tabla de posibles valores de $x$, con los respectivos $y$ predichos y el "error estándar" para cada predicción, que es diferente para cada uno (esto es lo que no entiendo).

Para dar dos ejemplos, tienen:

$$x=20, y=1.29, stderror=0.21$$

$$x=30, y=0.5, stderror=0.15$$

Los errores estándar son diferentes para los diferentes ejemplos.

Ahora, unas páginas antes, discuten cómo calcular la media de una muestra y calcular el error estándar según la fórmula:

$$se = [\sum_{i=1}^n \frac{(x_i-\bar x)^2}{n(n-1)}]^{0.5}$$

Indican que en el caso de regresión lineal, calculan los errores estándar usando una "versión extendida" de la fórmula anterior, pero nunca dicen cuál es. No entiendo cómo han calculado el "error estándar" para los ejemplos individuales en el caso de regresión lineal, y por qué es diferente para cada ejemplo. ¿Cuál es la fórmula?

Al leer la página de wikipedia sobre el error estándar me hace pensar que el se es simplemente la desviación estándar de la muestra, pero esto no concuerda con lo que han hecho aquí.

Cualquier ayuda se agradece.

Answer 1

Puedes descargar su ejemplo desde https://web.stanford.edu/~hastie/CASI_files/DATA/kidney.txt y fácilmente replicar sus resultados.

> kidney <- read.table("kidney.txt", header=TRUE)
> str(kidney)
'data.frame':   157 obs. of  2 variables:
 $ age: int  18 19 19 20 21 21 21 22 22 22 ...
 $ tot: num  2.44 3.86 -1.22 2.3 0.98 -0.5 2.74 -0.12 -1.21 0.99 ...
> fit <- lm(tot ~ age, data=kidney)
> fit$coefficients
(Intercept)         age 
 2.86002680 -0.07858842 

En cuanto a los errores estándar, los errores estándar de los valores ajustados, $\text{se}(\hat{y})$, son diferentes de los errores estándar de los coeficientes, $\text{se}(\hat\beta)$.

La matriz modelo $X$ es:

> X <- model.matrix(fit)
> head(X)
  (Intercept) age
1           1  18
2           1  19
3           1  19
4           1  20
5           1  21
6           1  21

Colocando $S=(X^TX)^{-1}$, $\text{cov}(\hat\beta)=\sigma^2_yS$ (ver mi respuesta a esta pregunta). Dado un único valor ajustado, $\hat{y}_h$ y la fila correspondiente $h$ de $X$, por ejemplo $$y_1=2.44,\qquad x_1=\begin{bmatrix}1 \\ 18\end{bmatrix}$$ la varianza de $\hat{y}_h$ es: $$\text{var}(\hat{y}_h)=\text{var}(x_h^T\hat\beta)=x_h^T\text{cov}(\hat\beta)x_h=x_h^T(S\sigma^2_y)x_h =\sigma^2_y(x_h^TSx_h)$$ Estimas $\sigma^2_y$ por la media cuadrática residual, RMS, el error estándar de $\hat{y}_h$ es: $$\text{se}(\hat{y}_h)=\sqrt{RMS(x_h^TSx_h)}$$ y depende de $x_h$.

Cuando solo hay una variable independiente, $$S=(X^TX)^{-1}=\frac{1}{n\sum(x_i-\bar{x})^2} \begin{bmatrix}\sum x_i^2 & -\sum x_i \\ -\sum x_i & n \end{bmatrix}$$ y $$\begin{align*} x_h^T(X^TX)^{-1}x_h &=\frac{\sum x_i^2-2x_hn\bar{x}+nx_h^2}{n\sum(x_i-\bar{x})^2}=\frac{\sum x_i^2 -n\bar{x}^2+n(x_h-\bar{x})^2}{n\sum(x_i-\bar{x})^2}\\ &=\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2} \end{align*}$$ (Recuerda que $\sum(x_i-\bar{x})^2=\sum x_i^2-n\bar{x}^2$).

La "versión extendida de la fórmula (1.2)" (que es simplemente el error estándar de una media) es: $$\text{se}(\hat{y}_h)=\left[RMS\left(\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)\right]^{\frac12}$$ Por cierto, así es como se calculan los intervalos de confianza.

Ver Kutner, Nachtsheim, Neter & Li, Applied Linear Statistical Models, McGraw-Hill, 2005, §2.4, o Seber & Lee, Linear Regression Analysis, John Wiley & Sons, 2003, §6.1.3.

En R:

> S <- solve(t(X) %*% X)
> RMS <- summary(fit)$sigma^2
> x_h <- matrix(c(1, 20), ncol=1)             # primer error estándar en la Tabla 1.1
> y_h_se <- sqrt(RMS * (t(x_h) %*% S %*% x_h)); y_h_se
          [,1]
[1,] 0.2066481
> x_h <- matrix(c(1, 80), ncol=1)             # último error estándar en la Tabla 1.1
> y_h_se <- sqrt(RMS * (t(x_h) %*% S %*% x_h)); y_h_se
         [,1]
[1,] 0.420226

EDITAR

Si estás interesado en el error estándar de $\hat{y}_{h(new)}=\hat\alpha+\hat\beta x_{h(new)}$, cuando $x_{h(new)}$ es una nueva observación, no sabes cuál sería $\hat{y}_h$ en una regresión en $n+1$ puntos. Diferentes muestras darían diferentes predicciones, por lo que debes tener en cuenta la desviación de $\hat{y}_{h(new)}$ alrededor de $\hat{y}_h=\hat\alpha+\hat\beta x_h$: $$\text{var}[y_{h(new)}-\hat{y}_h]=\text{var}(y_{h(new)})+\text{var}(\hat{y}_h)$$ Entonces, la varianza de tu predicción tiene dos componentes: la varianza de $y$, que estimas por RMS, y la varianza de la distribución muestral de $\hat{y}_h$, $RMS(x_h^TSx_h)$:

$$RMS + RMS\left(\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)$$ La "versión extendida de la fórmula (1.2)" se convierte en: $$\text{se}(\hat{y}_{h(new)})=\left[RMS\left(1+\frac1n+\frac{(x_{h(new)}-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)\right]^{\frac12}$$ Ver Kutner, Nachtsheim, Neter & Li, Applied Linear Statistical Models, McGraw-Hill, 2005, §2.5, o https://online.stat.psu.edu/stat501/lesson/3/3.3.

Answer 2

El valor predicho en $X=x$ es $\hat\mu=\hat\beta_0+\hat\beta_1x$. Esta es una función de una constante conocida, $x$, y variables aleatorias $(\hat\beta_0, \hat\beta_1)$. El error estándar de $\hat\mu$ es su desviación estándar, que es una función de la desviación estándar de $(\hat\beta_0, \hat\beta_1)$

Específicamente, la varianza de $x\hat\beta$ es $$x^2\mathrm{var}[\hat\beta_1]+2x\mathrm{cov}[\hat\beta_1,\hat\beta_0]+ \mathrm{var}[\hat\beta_0].$$

Esto depende de $x$, por lo que es diferente para cada observación. Dado que conocemos $x$ y tenemos un buen estimador de la matriz de varianza-covarianza de $\hat\beta$, podemos estimarlo.

La razón de la fórmula simplificada que mencionas es que el álgebra lineal se simplifica si la media de $X$ es cero, de modo que $\hat\beta_0$ y $\hat\beta_1$ son incorrelacionados. Puedes lograr eso transformando $x$ a $x-\bar x$.