Problem Statement
So I read in Elements of Statistical Learning, p. 54 that another way of doing OLS for multiple predictors is the following way. Let's assume two predictors $X_1, X_2$ and a target variable: $Y$.
Regress: $Y \sim X_1$, that gives you $\beta_1$
then take the residuals from that, $e_1$ and run: $X_2 \sim e_1$ to get your $\beta_2$
Intuitive Explanation
Does make sense intuitively as the first step gives us: "How much of $Y$ is explainable by $X_1$", and then the second step tells us: "How much of Y not explained by $X_1$, can be explained by $X_2$. i.e. How much residual variance in $Y$ is explainable by $X_2$."
Empirical Simulation
But then I had to ask: what happens when you flip things the other way? So I wrote a small bit of code to simulate this.
Trivial example: If $X_1$ and $X_2$ are standard gaussian, I get 1's all around...
But it gives the different values if I flip the order and they don't have the same variance...
Here is the code (Python3) that runs what I just said for a fictional $Y = X_1 + X_2$:
where $X_1 = N(0, 4)$ and $X_2 = N(0, 1)$
($X_1$ in the code is just $2*N(0,1)$)
import numpy as np
from sklearn import linear_model
# Generating Y ~ X1 + X2
N = int(1e4)
X1 = []
X2 = []
X1X2 = []
Y = []
for i in range(N):
x1 = np.random.normal()*2 # (multiplying by 2 makes variance=4)
x2 = np.random.normal()
y = x1 + x2
X1.append(x1)
X2.append(x2)
# stack them both for the full regression:
X1X2.append([x1,x2])
Y.append(y)
X1 = np.array(X1).reshape(-1,1)
X2 = np.array(X2).reshape(-1,1)
X1X2 = np.array(X1X2)
Y = np.array(Y)
# Running the regression starting with X1 first
reg = linear_model.LinearRegression()
reg.fit(X1, np.array(Y))
resid = Y - reg.predict(X1)
print(f"Y~X1 --> Beta_2 is: {reg.coef_}")
resid_fit = linear_model.LinearRegression()
resid_fit.fit(resid.reshape(-1,1), X2)
print(f"X2~e1 --> Beta_2 is: {resid_fit.coef_}")
print("-------------")
# Then flip:
reg = linear_model.LinearRegression()
reg.fit(X2, np.array(Y))
resid = Y - reg.predict(X2)
print(f"Y~X2 --> Beta_2 is: {reg.coef_}")
resid_fit = linear_model.LinearRegression()
resid_fit.fit(resid.reshape(-1,1), X1)
print(f"X1~e2 --> Beta_1 is: {resid_fit.coef_}")
print("-------------")
print("-------------")
# Then both:
reg = linear_model.LinearRegression()
reg.fit(X1X2, np.array(Y))
print(f"Full Regression Coeffs: ", reg.coef_)
If I start with $X_1$, I get the true coefficients: $\beta_1=2, \beta_2=1$, but if I start with $X_2$, I get: $ \beta_1=1, \beta_2=0.5$, half the real values?.
I ran it with $X_1 = N(0,4)$ and $X_2 = N(0, 9)$, and got the following output (rounded for clarity):
Y~X1 --> Beta_1 is: [2.01558849]
X2~e1 --> Beta_2 is: [[0.33333333]]
-------------
Y~X2 --> Beta_2 is: [3.00462912]
X1~e2 --> Beta_1 is: [[0.5]]
-------------
-------------
Full Regression Coeffs: [ 2.00000000e+00 3.00000000e+00 -3.33066907e-16]
So it looks like it's doing:
$$\hat{\beta_i} = \frac{\beta_i}{Var(X_i)}$$
for whichever $X_i$ goes second... why?
Question:
- What's going on here?
- Is there an intuitive explanation?
- How can I derive that final formula for $\beta_i$ I wrote above?
