Matriz de varianza-covarianza de los errores en regresión lineal

Question

¿Cómo se calcula en la práctica la matriz de errores de varianza/covarianza por paquetes de análisis estadístico?

Esta idea me resulta clara en teoría. Pero no en la práctica. Es decir, si tengo un vector de variables aleatorias $\textbf{X}=(X_{1}, X_{2}, \ldots, X_{n})^\top$, entiendo que la matriz de varianza/covarianza $\Sigma$ se obtendrá a partir del producto externo de los vectores de desviación de la media: $\Sigma=\mathrm{E}\left[(\textbf{X}-\mathrm{E}(\textbf{X}))(\textbf{X}-\mathrm{E}(\textbf{X}))^\top\right]$.

Pero cuando tengo una muestra, los errores de mis observaciones no son variables aleatorias. O mejor dicho, lo son, pero solo si tomo un número de muestras idénticas de la misma población. De lo contrario, están dados. Entonces, nuevamente mi pregunta es: ¿cómo puede un paquete estadístico producir una matriz de varianza/covarianza a partir de una lista de observaciones (es decir, una muestra) proporcionada por el investigador?

Answer 1

La matriz de covarianza para un modelo del tipo $y = X\beta + \epsilon$ generalmente se calcula como $$(X^t X)^{-1}\frac{\sigma^2}{d}$$ donde $\sigma^2$ es la suma residual de cuadrados, $\sigma^2=\sum_i (y_i - X_i\hat\beta)^2$ y $d$ es el número de grados de libertad (típicamente el número de observaciones menos el número de parámetros).

Para errores estándar robustos o agrupados, el producto $X^t X$ se modifica ligeramente. También puede haber otras formas de calcular la matriz de covarianza, por ejemplo, como sugiere la esperanza de los productos externos.

Answer 2

1. Estimación OLS de la varianza del error, $\sigma^2$:

$$s^2=\frac{\hat \varepsilon^\top\hat \varepsilon}{n-p}$$

Esto está incluido en Regresión Práctica y Anova usando R por Julian J. Faraway, página 21 .

Ejemplo de su cálculo en R, basado en un modelo lineal de millas por galón regresado sobre las especificaciones de múltiples modelos de automóviles incluidos en la base de datos mtcars: ols = lm(mpg ~ disp + drat + wt, mtcars). Estos son los cálculos manuales y la salida de la función lm():

> rdf = nrow(X) - ncol(X)                    # Grados de libertad residuales
> s.sq = as.vector((t(ols$residuals) %*% ols$residuals) / rdf) 
>                                            # s cuadrado (Estimación OLS de sigma cuadrado)
> (sigma = sqrt(s.sq))                       # Error estándar residual
[1] 2.950507
> summary(ols)

Call:
lm(formula = mpg ~ disp + drat + wt, data = mtcars)
...
Residual standard error: 2.951 on 28 degrees of freedom

2. Matriz de Varianza-Covarianza de los coeficientes estimados, $\hat \beta$:

$$\mathrm{Var}\left[\hat \beta \mid X \right] =\sigma^2 \left(X^\top X\right)^{-1}$$

estimada como en la página 8 de este documento en línea como

$$\hat{\mathrm{Var}}\left[\hat \beta \mid X \right] =s^2 \left(X^\top X\right)^{-1}$$

> X = model.matrix(ols)                             # Matriz del modelo X
> XtX = t(X) %*% X                                  # X transpuesta X
> Sigma = solve(XtX) * s.sq                         # Matriz de varianza-covarianza
> all.equal(Sigma, vcov(ols))                       # Igual que la fórmula integrada
[1] TRUE
> sqrt(diag(Sigma))                                 # Errores estándar calculados de los coeficientes
(Intercept)        disp          drat          wt 
7.099791769 0.009578313 1.455050731 1.217156605 
> summary(ols)[[4]][,2]                             # Salida de la función lm()
(Intercept)        disp        drat          wt 
7.099791769 0.009578313 1.455050731 1.217156605

Answer 3

Con la regresión lineal estamos ajustando un modelo $Y = \beta*X +\varepsilon$. $Y$ es la variable dependiente, las $X$ son las variables predictoras (explicativas). Utilizamos los datos proporcionados (el conjunto de entrenamiento o la muestra) para estimar los $\beta$ de la población. Las $X$ no se consideran variables aleatorias. Los $Y$ son aleatorios debido al componente de error.