# How are partial regression slopes calculated in multiple regression?

I'm trying to understand how multiple regression statistically controls for the effects of other predictor variables when calculating partial regression slopes. In a multiple regression of Y~X1+X2, would the partial regression slope of X1 be given by [Y]~[residuals of X1~X2], or by [residuals of Y~X2] ~ [residuals of X1~X2]? Different pages on the internet tell me different things.

I've done some simulations to try and figure this out (see below), and it seems that both methods give the same estimates of slopes as multiple regression, but only the latter method has similar standard errors around those estimates. This makes me think the latter method is the one that multiple regression uses, but it would be nice to know for sure.

Similarly, if I wanted to plot Y against X1 so that I could visualise how strongly the two were related while also controlling for any confounding with X2, would I plot [Y]~[residuals of X1~X2], or [residuals of Y~X2] ~ [residuals of X1~X2]? These two plots in the code below look very different in terms of the strength of the relationship.

Jay

#1. simulate data, where x1 and x2 are correlated due to lurking variable,
#...and y is explained by both.
lurker <- rnorm(n=100)
x1 <- rnorm(n=100, mean=lurker*2, sd=1)
x2 <- rnorm(n=100, mean=lurker*5, sd=1)
y <- rnorm(n=100, mean=x1*2 + x2*5, sd=1)

#2. multiple regn model to estimate partial slopes:
summary(lm(y~x1+x2))      #partial slopes pretty close to simulated values

#3. calculate partial slopes manually, using either
#....(1) Y~[resids of X1~X2] OR (2) [resids of Y~X2]~[resids of X1~X2]

#3.a. based on (1) Y~[resids of X1~X2]
m.x1x2 <- lm(x1~x2)
resids.x1x2 <- m.x1x2$residuals summary(lm(y ~resids.x1x2)) #slope pretty close to true value, but conf intervals MUCH larger than those for MR estimate #3.b. based on (2) [resids of Y ~ X2]~[resids of X1~X2] m.y1x2 <- lm(y~x2) resids.y1x2 <- m.y1x2$residuals
summary(lm(resids.y1x2 ~resids.x1x2))  #also very close to true value, but conf intervals now similar scale to those from MR.

# plot the relationship between y and x1 after controlling for x2, based on the different methods:
op <- par(mfrow=c(2,1))
plot(y ~resids.x1x2)
plot(resids.y1x2 ~resids.x1x2)
par(op)


if I wanted to plot Y against X1 so that I could visualise how strongly the two were related while also controlling for any confounding with X2, would I plot [Y]~[residuals of X1~X2], or [residuals of Y~X2] ~ [residuals of X1~X2]?

The second one. You want to plot $e(Y|X_2) \sim e(X_1|X_2)$ to obtain the added variable plot for $X_1$ for the model already containing $X_2$. Intuitively think of it this way: the presence of $X_2$ in the model you already have is reducing the variance of the residuals (which is relative to the response). You want to determine if adding $X_1$ will reduce this variance further, so you want to consider $e(Y|X_2)$, which contains the variance remaining in the residuals after taking the predictor $X_2$ into consideration.

As you put it: you want to "partial out" the effects of the predictor(s) currently in your model on both the predictor(s) you are considering adding as wells as the response.

To avoid blackening the place, I will not use bold symbols -but the answer will be carried out in matrix form. Vectors are column vectors, a prime will denote the transpose.

Let a linear regression model

$$y = X_1b_1 + X_2b_2 + u_A \qquad [A]$$

The normal equations for the OLS estimator are

\begin{align} \left(X_1'X_1\right)b_1+\left(X_1'X_2\right)b_2=& X_1'y \qquad [1]\\ \\ \left(X_2'X_1\right)b_1+\left(X_2'X_2\right)b_2=& X_2'y \qquad [2]\\ \end{align} Solving $[2]$ for $b_2$ we have $$[2]\rightarrow b_2= \left(X_2'X_2\right)^{-1}X_2'y-\left(X_2'X_2\right)^{-1}\left(X_2'X_1\right)b_1$$

Inserting this into $[1]$ we obtain $$\left(X_1'X_1\right)b_1+\left(X_1'X_2\right)\left(X_2'X_2\right)^{-1}X_2'y-\left(X_1'X_2\right)\left(X_2'X_2\right)^{-1}\left(X_2'X_1\right)b_1= X_1'y$$

Collecting terms w.r.t $b_1$ and $y$, $$X_1'\left[I-X_2\left(X_2'X_2\right)^{-1}X_2'\right]X_1b_1= X_1'\left[I-X_2\left(X_2'X_2\right)^{-1}X_2'\right]y$$ $$\Rightarrow X_1'M_2X_1b_1 = X_1'M_2y \qquad [3]$$

where $M_2$ is the "annihilator" or "residual maker"matrix related to $X_2$, namely the matrix that produces the residuals when a variable is regressed on $X_2$, by pre-multiplying this variable. This matrix is symmetric and idempotent, $M_2=M_2',\; M_2= M_2M_2$. So we can write

$$(M_2X_1)'(M_2X_1)b_1 = (M_2X_1)'y$$ $$\Rightarrow R_{1\sim2}'R_{1\sim2}b_1=R_{1\sim2}'y \Rightarrow \hat b_1 = \left(R_{1\sim2}'R_{1\sim2}\right)^{-1}R_{1\sim2}'y\qquad [4]$$

where $R_{1\sim2}$ denotes the residual vector from regressing $X_1$ on $X_2$.

This last formula is exactly the OLS formula from the regression model $$y= R_{1\sim2}d_1+u_B \qquad [B]$$

So eq. $[4]$ tells us that the coefficient estimate for $X_1$ that we will obtain in a multiple regression setting, will be exactly the same with what we will obtain if we regress the dependent variable on the residuals from the regression of $X_1$ on $X_2$.

Now consider the second case, regressing the residuals on the residuals. This is the model

$$R_{y\sim2} = R_{1\sim2}c_1+u_C \Rightarrow (M_2y)= (M_2X_1)c_1 +u_C \qquad [C]$$

The OLS estimator of $c$ is $$\hat c_1 = \left[(M_2X_1)'(M_2X_1)\right]^{-1}(M_2X_1)'(M_2y) \qquad [5]$$

By the properties of $M_2$ we have $$(M_2X_1)'(M_2y) = X_1'M_2'M_2y=X_1'M_2M_2y=X_1'M_2y=X_1'M_2'y=(M_2X_1)'y$$ Noting that $M_2X_1 = R_{1\sim2}$ eq. $[5]$ becomes

$$\hat c_1= \left(R_{1\sim2}'R_{1\sim2}\right)^{-1}R_{1\sim2}'y \qquad [6]$$

which is identical to eq. $[4]$, and so $\hat c_1 = \hat d_1 =\hat b_1$. In other words the three models give mathematically identical results.

Let's now consider the issue of the estimator variance. Models $[B]$ and $[C]$ have the same regressor matrix so the question is what happens with the estimated error variances, $\sigma^2_B$ and $\sigma^2_C$. We will denote $M(r)_{1\sim2}$ the annihilator matrix of the regressor $R_{1\sim2}$. It has analogous properties as $M_2$ For model $[B]$ we have

$$u'_Bu_B = \left(M(r)_{1\sim2}y\right)'\left(M(r)_{1\sim2}y\right) = y'M(r)_{1\sim2}y \qquad [7]$$

while for model $[C]$ we have

$$u'_Cu_C = \left(M(r)_{1\sim2}(M_2y)\right)'\left(M(r)_{1\sim2}(M_2y)\right) = y'M_2M(r)_{1\sim2}M_2y \qquad [8]$$

Are the RHS of eq. $[7]$ and $[8]$ equal? I 'll leave that to the reader.

• Thanks for that in depth answer...which unfortunately has gone way over my head! At this stage I'm after a more conceptual cf mathematical understanding - ie is the partial slope of X1 in a model also containing X2 found by partialling out effects of X2 on X1 AND partialling out the effects of X2 on Y, or just by partialling out the effects of X2 on X1? Cheers
– jay
Commented Nov 25, 2013 at 3:38
• "The three models give mathematically identical results" means that either way you will obtain the exact same result. Can you get the intuition as to why? Commented Nov 25, 2013 at 8:42
• I've read through it multiple times and I'm afraid no, don't get the intuition. I appreciate you taking the time to write out such a detailed explanation...it's just not in a language that I understand.
– jay
Commented Nov 28, 2013 at 7:14
• No problem Jay... try this: informally, residuals are "what is left unexplained". The X2~X1 residuals is what of X1 is not explained by X2. Now consider the Y~X2 residuals. These reflect what part of Y is not explained by X2 (and therefore is explained by X1 and by the error term). So even if you regress the original Y on the X2~X1 residuals, these will only explain that part of Y that is left unexplained by X2. So by using Y~X2 as the dependent variable, you "remove" only that part of Y that is explained by X2- leaving unaffected that part of Y that is explained by the residuals X2~X1. Commented Nov 28, 2013 at 14:57
• Yes that makes sense thanks. We got there in the end!
– jay
Commented Dec 1, 2013 at 5:09