How do I get regression coefficients from a variance covariance matrix in R?

Question

I want to work out a multiple regression example all the way through using matrix algebra to calculate the regression coefficients.

#create vectors -- these will be our columns
y <- c(3,3,2,4,4,5,2,3,5,3)
x1 <- c(2,2,4,3,4,4,5,3,3,5)
x2 <- c(3,3,4,4,3,3,4,2,4,4)

#create matrix from vectors
M <- cbind(y,x1,x2)
k <- ncol(M) #number of variables
n <- nrow(M) #number of subjects

#create means for each column
M_mean <- matrix(data=1, nrow=n) %*% cbind(mean(y),mean(x1),mean(x2)); M_mean

#creates a difference matrix which gives deviation scores
D <- M - M_mean; D

#creates the covariance matrix, the sum of squares are in the diagonal and the sum of cross products are in the off diagonals.  
C <-  t(D) %*% D; C

I can see what the final values should be (-.19, -.01) and what the matrices before this calculation look like.

E<-matrix(c(10.5,3,3,4.4),nrow=2,ncol=2)
F<-matrix(c(-2,-.6),nrow=2,ncol=1)

But I'm not sure how to create these from the variance-covariance matrix to get the coefficients using matrix algebra.

Hope you can help.

I'm not sure what you are doing there, but it's not OLS regression: http://stats.stackexchange.com/a/80889/11849 — Roland, Nov 23 '16 at 11:16
The inverse of the matrix Sum of squares.xx times the matrix Sum of squares.xy should give me regression coefficients. — adkane, Nov 23 '16 at 11:28

score 3 · Accepted Answer · edited May 23 '17 at 12:30

I can see that you are doing centred regression:

The answer by sandipan is not quite what you want, as it goes through the usual normal equation to estimate:

There is already a thread on the latter: Solving normal equation gives different coefficients from using lm? Here I focus on the former.

Actually you are already quite close. You have obtained the mixed covariance C:

#      y   x1   x2
#y  10.4 -2.0 -0.6
#x1 -2.0 10.5  3.0
#x2 -0.6  3.0  4.4

From your definition of E and F, you know you need sub-matrices to proceed. In fact, you can do matrix subsetting rather than manually imputing:

E <- C[2:3, 2:3]

#     x1  x2
#x1 10.5 3.0
#x2  3.0 4.4

F <- C[2:3, 1, drop = FALSE]  ## note the `drop = FALSE`

#      y
#x1 -2.0
#x2 -0.6

Then the estimate is just , and you can do in R (read ?solve):

c(solve(E, F))  ## use `c` to collapse matrix into a vector
# [1] -0.188172043 -0.008064516

Other suggestions

you can find column means by colMeans, instead of a matrix multiplication (read ?colMeans);
you can perform centring by using sweep (read ?sweep);
use crossprod(D) than t(D) %*% D (read ?crossprod).

Here is a session I would do:

y <- c(3,3,2,4,4,5,2,3,5,3)
x1 <- c(2,2,4,3,4,4,5,3,3,5)
x2 <- c(3,3,4,4,3,3,4,2,4,4)

M <- cbind(y,x1,x2)
M_mean <- colMeans(M)
D <- sweep(M, 2, M_mean)
C <- crossprod(D)

E <- C[2:3, 2:3]
F <- C[2:3, 1, drop = FALSE]
c(solve(E, F))
# [1] -0.188172043 -0.008064516

I should revise this answer when I have time. At the moment it is not well structured. Simply put, if we have a response vector `y` and a matrix `X` of independent variables (excluding intercept), slope for those variables is `b <- c(solve(cov(X), cov(X, y)))`. And `mean(y) - c(crossprod(b, colMeans(X)))` is the intercept. — Zheyuan Li, Aug 21 '18 at 02:12

score 1 · Answer 2 · answered Nov 23 '16 at 11:44

Probably you want something like this:

X <- cbind(1, x1, x2)
C <- t(X) %*% X # No need of centering the columns with means
D <- t(X) %*% y

coef <- t(solve(C) %*% D)
coef
#                      x1           x2
# [1,] 4.086022 -0.188172 -0.008064516

lm(y~x1+x2)$coef # coefficients with R lm()
# (Intercept)           x1           x2 
# 4.086021505 -0.188172043 -0.008064516

How do I get regression coefficients from a variance covariance matrix in R?

2 Answers2