How are partial regression slopes calculated in multiple regression?

Question

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.

Thanks for your help,

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)

Answer 1

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.

Answer 2

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.