¿Hay alguna manera de usar la matriz de covarianza para encontrar coeficientes de regresión múltiple?

Para la regresión lineal simple, el coeficiente de regresión es calculable directamente a partir de la matriz de varianza-covarianza $C$ , por

Si uno solo tiene la matriz de covarianza, ¿es posible calcular los coeficientes para un modelo con múltiples variables explicativas?

ETA: para dos variables explicativas, parece que y análogamente paraβ2. No veo de inmediato cómo extender esto a tres o más variables.

Sí, la matriz de covarianza de todas las variables, explicativa y de respuesta, contiene la información necesaria para encontrar todos los coeficientes, siempre que se incluya un término de intercepción (constante) en el modelo. (Aunque las covarianzas no proporcionan información sobre el término constante, se puede encontrar a partir de los medios de los datos).

Análisis

Dejar que los datos para las variables explicativas pueden organizar como vectores columna -dimensional x 1 , x 2 , ... , x p y la variable de respuesta sean el vector columna y , considerado como una realización de una variable aleatoria Y . Los ordinarios estimaciones de mínimos cuadrados beta de los coeficientes en el modelo$n$$x_1, x_2, \ldots, x_p$$y$$Y$$\hat\beta$

se obtienen ensamblando los vectores de columna X 0 = ( 1 , 1 , … , 1 ) ′ , X 1 , … , X p en una matriz n × p + 1 X y resolviendo el sistema de ecuaciones lineales$p+1$$X_0 = (1, 1, \ldots, 1)^\prime, X_1, \ldots, X_p$$n \times p+1$$X$

Es equivalente al sistema

La eliminación gaussiana resolverá este sistema. Se procede por contigua a la de la matriz $p+1\times p+1$y elp+1-vector$\frac{1}{n}X^\prime X$$p+1$en unamatrizp+1×p+2Ay reduciéndola por filas. $\frac{1}{n}X^\prime y$$p+1 \times p+2$$A$

El primer paso inspeccionará . Al encontrar que esto no es cero, procede a restar los múltiplos apropiados de la primera fila deAde las filas restantes para poner a cero las entradas restantes en su primera columna. Estos múltiplos serán$\frac{1}{n}(X^\prime X)_{11} = \frac{1}{n}X_0^\prime X_0 = 1$$A$y el número restado de la entradaAi+1,j+1=X ′ i Xjserá igual a ¯ X i ¯ X j. Esta es solo la fórmula para la covarianza deXiyXj. Además, el número que queda en la posicióni+1,p+2es igual a$\frac{1}{n}X_0^\prime X_i = \overline X_i$$A_{i+1,j+1} = X_i^\prime X_j$$\overline X_i \overline X_j$$X_i$$X_j$$i+1, p+2$ , la covarianza deXicony.$\frac{1}{n}X_i^\prime y - \overline{X_i}\overline{y}$$X_i$$y$

Por lo tanto, después del primer paso de eliminación gaussiana, el sistema se reduce a la resolución

y obviamente, dado que todos los coeficientes son covarianzas, esa solución se puede encontrar a partir de la matriz de covarianza de todas las variables.

(Cuando es invertible, la solución se puede escribir C - 1 ( Cov ( X i , y ) ) ' . Las fórmulas dadas en la pregunta son casos especiales de esto cuando p = 1 y p = 2. Escribir estas fórmulas explícitamente se vuelven más y más complejos a medida que crece P. Además, son inferiores para el cómputo numérico, que se realiza mejor resolviendo el sistema de ecuaciones en lugar de invertir la matriz C ).$C$$C^{-1}(\text{Cov}(X_i, y))^\prime$$p=1$$p=2$$p$$C$

$y$$X\hat{\beta}$

Ejemplo

Para ilustrar, el siguiente Rcódigo crea algunos datos, calcula sus covarianzas y obtiene las estimaciones de coeficientes de mínimos cuadrados únicamente a partir de esa información. Los compara con las estimaciones obtenidas del estimador de mínimos cuadrados lm.

#
# 1. Generate some data.
#
n <- 10        # Data set size
p <- 2         # Number of regressors
set.seed(17)
z <- matrix(rnorm(n*(p+1)), nrow=n, dimnames=list(NULL, paste0("x", 1:(p+1))))
y <- z[, p+1]
x <- z[, -(p+1), drop=FALSE]; 
#
# 2. Find the OLS coefficients from the covariances only.
#
a <- cov(x)
b <- cov(x,y)
beta.hat <- solve(a, b)[, 1]  # Coefficients from the covariance matrix
#
# 2a. Find the intercept from the means and coefficients.
#
y.bar <- mean(y)
x.bar <- colMeans(x)
intercept <- y.bar - x.bar %*% beta.hat  

El resultado muestra un acuerdo entre los dos métodos:

(rbind(`From covariances` = c(`(Intercept)`=intercept, beta.hat),
       `From data via OLS` = coef(lm(y ~ x))))

                  (Intercept)        x1        x2
From covariances     0.946155 -0.424551 -1.006675
From data via OLS    0.946155 -0.424551 -1.006675