# Multiple regression and partial correlations with more than 3 predictors

I need to estimate slopes of a multiple regression that has 5 independent variables:

$$y = \alpha + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 + \beta_5 x_5$$

The challenge is to estimate the slopes by estimating correlations and variances separately. For a single regression this is rather simple because you can just use the formula:

$$\beta = r_{xy} \frac{\sigma_y}{\sigma_x}$$

I did some research and found out that I need partial correlations in the case of a multiple regression to control for my other variables. Now I wonder how I can calculate the slopes. Obviously the estimation of the variances is rather simple, and the partial correlations I calculate with the corpcor package in R. I just don't understand how from that point I can calculate the slopes. For the case of 2 independent variables I found the following formula:

$$\beta_1 = \frac{r_{yx_1}-(r_{yx_2}r_{x_1x_2})}{(1-r_{x_1x_2}^2)} \frac{\sigma_y}{\sigma_{x_1}}$$

But I cannot find a good source for the case of more than 2 predictors. An explanation, a good reference or a package in R would be highly appreciated.

Context of my study: I try to regress log stock returns ($$y$$) on common factors in the finance literature ($$x_1$$ to $$x_5$$). I need to estimate correlations and volatilities separately because I am asked to use different data windows for the two. This is the tricky part for me because so far I mostly just did the typical OLS regression in R with the lm() function.

• Please check the last formula; it's not clear to me what $r_{23}^2$ might mean with only 2 predictors. Also, it would help if you could provide a reference (with a web link, if possible) for that formula.
– EdM
Dec 25 '19 at 21:24
• @EdM thank you for pointing this out. I updated the formula, it was a typo. Unfortunately, I did not find an accessible reference. I have the formula from a statistics course I took a while ago in university. Dec 26 '19 at 8:49
• This could be relevant: stats.stackexchange.com/q/76815/3277 Dec 26 '19 at 19:31

This can be done, but it would be terribly unwieldy and might be numerically unstable if the predictors are highly correlated.

It's simplest if we assume that each of the $$X_i$$ and $$Y$$ values is expressed as its difference from its respective mean. Then we can dispense with the intercept term of the regression equation and the column of 1s more generally included in the design matrix $$\mathbf{X}$$ (see this page) and write the general solution to the multiple regression in matrix form as:

$$\mathbf{\hat\beta = (X^TX)^{-1}X^T}Y,$$

where $$\mathbf{X}$$ is the matrix of (mean-centered) predictor values (number of rows $$n$$ equal to the number of observations, and number of columns $$p$$ equal to the number of predictors) and $$Y$$ is the vector of $$n$$ (mean-centered) outcome observations.

To put this into the context and notation of your question, note that the correlation coefficient $$r_{XY}$$ is $$\mathsf{Cov}(X,Y)/(\sigma_X \sigma_Y),$$ where $$\mathsf{Cov}(X,Y)$$ is the covariance between $$X$$ and $$Y$$. For a single predictor, $$X$$ is a vector with $$n$$ elements, and with $$X$$ and $$Y$$ centered around their mean values

$$X^TY = n \mathsf{Cov}(X,Y).$$

With a single predictor centered around its mean, $$\mathbf{X^TX}$$ is simply $$n \sigma_X^2$$, with inverse $$1/(n \sigma_X^2)$$. By the general matrix formula, $$\hat\beta = \mathsf{Cov}(X,Y)/\sigma_X^2$$, while your formula gives the same result:

$$\hat\beta = r_{XY} \frac{\sigma_Y}{\sigma_X} = \frac{\mathsf{Cov}(X,Y)}{\sigma_X \sigma_Y}\frac{\sigma_Y}{\sigma_X}=\frac{\mathsf{Cov}(X,Y)}{\sigma_X^2}.$$

As you increase the number of predictors $$p$$, $$\mathbf{X^T}Y$$ is just a simple generalization of the above: it's a vector of the $$p$$ individual covariances of the $$X_i$$ with $$Y$$, multiplied by $$n$$. If you wish, you can express each covariance in terms of the correlation coefficient and the square roots of individual variances:

$$\mathsf{Cov}(X_i,Y) = r_{X_iY}\sigma_{X_i}\sigma_Y.$$

The inner-product matrix product $$\mathbf{X^TX}$$ of mean-centered predictor values is a less-simple generalization of the variance of the single predictor; it's $$n$$ times the variance-covariance matrix of $$\mathbf{X}$$. That variance-covariance matrix has the individual variances of the $$X_i$$ along the diagonal, with each off-diagonal $$(i,j)$$ element equal to the corresponding covariance, $$\mathsf{Cov}(X_i,X_j)$$. You can express each of these covariances between predictors in terms of the corresponding correlation coefficients and individual (square roots of) variances, similarly to the covariance between an individual predictor and the outcome as above.

You must, however, take the inverse of the predictor variance-covariance matrix to calculate the vector of regression coefficients $$\mathbf{\beta}.$$ That's where expressing things in terms of correlation coefficients and individual (square roots of) variances gets complicated.

You could get a closed-form solution for $$\mathbf{X^TX}$$ for any number of predictors with elements expressed in that way. Brute force with the formula for a matrix inverse would do the trick, or you could start with $$\mathbf{(X^TX)^{-1}}=1/(n \sigma_x^2)$$ for the single-predictor case and add one predictor at a time following the formula on Page 19 of the Matrix Cookbook for updating $$\mathbf{(X^TX)^{-1}}$$ when you add one predictor vector to it. This gets messy very quickly, however.

Even if you got a closed-form solution in terms of correlation coefficients and individual variances, it would be unwise just to plug into that formula particularly if there were strong correlations among the predictors. In that case you might get numerically unstable results. Standard statistical software is designed to solve this matrix equation in a way that minimizes those problems.

Heuristically, think of the multiple-predictor formula

$$\mathbf{\hat\beta = (X^TX)^{-1}X^T}Y$$

as a generalization of the single-predictor formula

$$\hat\beta =\frac{\mathsf{Cov}(X,Y)}{\sigma_x^2}.$$

In both cases, you are multiplying the covariance(s) between the outcome and the predictor(s) by the inverse of the (co)variance(s) of the predictors. That's perhaps simpler than thinking in terms of correlation coefficients.* To do the calculation in practice though, let well vetted statistical software to the work for you.

In response to update and comment: As you are working with time-series data this approach is probably inappropriate for your application. That's particularly true if you are estimating the variances and the correlation coefficients from different data windows.

*If you want to think in terms of correlation coefficients, start with outcomes and predictors each not only centered around their means but also scaled to unit standard deviation/variance. Then the diagonal entries of $$\mathbf{X^TX}$$ are all 1, the off-diagonal entries are the pairwise correlation coefficients between the predictors, and $$\mathbf{X^T}Y$$ is the vector of correlation coefficients between each of the predictors and $$Y.$$ After applying the formula for $$\mathbf{\hat\beta}$$ re-scale the regression coefficients to take the initial scaling into account; that’s where the factors like $$\sigma_Y/\sigma_{X_i}$$ come from in your formula for the two-predictor case. You could use a symbolic mathematics manipulation program like Mathematica to do the algebra. I gave that a try, recovering your formula for the 2-predictor case. The result wasn't terribly complicated for 3 predictors, but with 4 or 5 predictors you end up with very many terms both in the numerator and denominator.

• Thank you for this detailed answer. I think I have to emphasize that the difficulty for me is that in my particular case I need to estimate the correlations and volatilities (standard deviations) of (log) stock returns separately to estimate the regression coefficients. I do this separately, because I am asked to use different data windows for the two. The estimation of correlations and volatilities is easy. For me the problem is how to get from here to the estimation of regression coefficients in a multiple regression and implement this in R. Any ideas how to implement this in practice? Dec 26 '19 at 9:07
• @Greenym interpreting correlations in time-series data is tricky. See this page for an introduction. Using standard deviations from one time window and correlations from another seems even more problematic. Even if you used the principles in this answer to get a formula, I worry that the results might not mean what you think. Try asking a separate question focusing on the time-series nature of your data and the specific task you are trying to accomplish. I’m not a time-series expert, but I suspect that there are better approaches.
– EdM
Dec 26 '19 at 11:03