Existe-t-il un moyen d'utiliser la matrice de covariance pour trouver des coefficients de régression multiple?

23

Pour une régression linéaire simple, le coefficient de régression peut être calculé directement à partir de la matrice de variance-covariance C , par

Cd,eCe,e
dest l'indice de la variable dépendante eteest l'indice de la variable explicative.

Si l'on n'a que la matrice de covariance, est-il possible de calculer les coefficients d'un modèle à multiples variables explicatives?

ETA: Pour deux variables explicatives, il apparaît que et de façon analogue pourβ2. Je ne vois pas immédiatement comment étendre cela à trois variables ou plus.

β1=Cov(y,x1)var(x2)Cov(y,x2)Cov(x1,x2)vuner(X1)vuner(X2)-Cov(X1,X2)2
β2
David
la source
3
Le vecteur de coefficient β est une solution à X ' Y = ( X ' X ) - 1 β . Certaines manipulations algébriques révèlent qu'il s'agit en fait de la même formule que celle que vous donnez dans le cas des 2 coefficients. Joliment présenté ici: stat.purdue.edu/~jennings/stat514/stat512notes/topic3.pdf . Je ne sais pas si cela aide du tout. Mais j'oserais deviner que cela est impossible en général sur la base de cette formule. β^XY=(XX)1β
shadowtalker
1
@David Avez-vous compris comment étendre cela à un nombre arbitraire de variables explicatives (au-delà de 2)? J'ai besoin de l'expression.
Jane Wayne
1
@JaneWayne Je ne suis pas sûr de comprendre votre question: whuber a donné la solution ci-dessous sous forme matricielle, C1(Cov(Xi,y))
David
1
yup je l'ai étudié et il a raison.
Jane Wayne

Réponses:

36

Oui, la matrice de covariance de toutes les variables - explicatives et réponses - contient les informations nécessaires pour trouver tous les coefficients, à condition qu'un terme d'interception (constant) soit inclus dans le modèle. (Bien que les covariances ne fournissent aucune information sur le terme constant, elles peuvent être trouvées à partir des moyennes des données.)


Une analyse

Que les données pour les variables explicatives être disposés comme vecteurs colonnes de dimension x 1 , x 2 , ... , x p , et la variable de réponse soit le vecteur colonne y , considéré comme une réalisation d'une variable aléatoire Y . Les estimations des moindres carrés ordinaires ß des coefficients dans le modèlenx1,x2,,xpyYβ^

E(Y)=α+Xβ

sont obtenus en assemblant les vecteurs de colonnes X 0 = ( 1 , 1 , , 1 ) , X 1 , , X p en un tableau n × p + 1 X et en résolvant le système d'équations linéairesp+1X0=(1,1,,1),X1,,Xpn×p+1X

XXβ^=Xy.

Il est équivalent au système

1nXXβ^=1nXy.

L'élimination gaussienne résoudra ce système. Il procède en joignant le matrice + 1 1p+1×p+1et lep1nXXvecteur + 1 1p+1dans untableaup+1×p+2Aet en le réduisant en ligne. 1nXyp+1×p+2A

La première étape inspectera 1n(XX)11=1nX0X0=1A1nX0Xi=X¯iAi+1,j+1=XiXj will equal X¯iX¯j. This is just the formula for the covariance of Xi and Xj. Moreover, the number left in the i+1,p+2 position equals 1nXiyXi¯y¯, the covariance of Xi with y.

Thus, after the first step of Gaussian elimination the system is reduced to solving

Cβ^=(Cov(Xi,y))

and obviously--since all the coefficients are covariances--that solution can be found from the covariance matrix of all the variables.

(When C is invertible the solution can be written C1(Cov(Xi,y)). The formulas given in the question are special cases of this when p=1 and p=2. Writing out such formulas explicitly will become more and more complex as p grows. Moreover, they are inferior for numerical computation, which is best carried out by solving the system of equations rather than by inverting the matrix C.)

The constant term will be the difference between the mean of y and the mean values predicted from the estimates, Xβ^.


Example

To illustrate, the following R code creates some data, computes their covariances, and obtains the least squares coefficient estimates solely from that information. It compares them to the estimates obtained from the least-squares estimator 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  

The output shows agreement between the two methods:

(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
whuber
la source
1
Thanks, @whuber! This is exactly what I was looking for, and my atrophied brain was unable to get to. As an aside, the motivation for the question is that for various reasons we essentially do not have the full X available, but have cov(z) from previous calculations.
David
7
Answers like this raise the bar of this Cross Validated
jpmuc
@whuber In your example, you computed the intercept from y and x and beta.hat. The y and x are part of the original data. Is it possible to derive the intercept from the covariance matrix and means alone? Could you please provide the notation?
Jane Wayne
@Jane Given only the means X¯, apply β^ to them:
X¯β^=Xβ^¯.
I have changed the code to reflect this.
whuber
very helpful +1 for the code
Michael