3.1 OLS Mechanics: the Line of Best Fit

Why regression is still the workhorse

Causal machine learning uses random forests and boosting, but its core move, partialling out covariates, is the logic of linear regression. Instrumental variables, difference-in-differences and regression discontinuity are all regressions in disguise. Understanding exactly what OLS does is the best investment you can make before Part II.

!pip install -q numpy scipy polars statsmodels matplotlib seaborn
import numpy as np, polars as pl
import statsmodels.formula.api as smf
from scipy.optimize import minimize
import matplotlib.pyplot as plt, seaborn as sns

rng = np.random.default_rng(11)
sns.set_theme(style="darkgrid")

Simulating a Mincer earnings world

The Mincer equation relates earnings to schooling and experience. We simulate 1,000 workers in a world where each extra year of education causes a $2,500 increase in annual wages. For now there is no confounding: education is assigned independently of everything else.

N = 1_000
TRUE_RETURN = 2_500
education = np.clip(np.round(rng.normal(13, 2, N)), 8, 20)
experience = rng.uniform(0, 30, N)
wage = 20_000 + TRUE_RETURN * education + 600 * experience + rng.normal(0, 6_000, N)

labor = pl.DataFrame({"education": education, "experience": experience, "wage": wage})
print(labor.describe())

Three ways to get the same line

OLS picks the intercept and slope that minimise the sum of squared residuals (SSR), the squared vertical distances between each point and the line.

1. Numerical optimisation. Let the computer search for the minimum, as you would without calculus.

x, y = labor["education"].to_numpy(), labor["wage"].to_numpy()

def ssr(beta):
    return np.sum((y - (beta[0] + beta[1] * x)) ** 2)

opt = minimize(ssr, x0=[0.0, 0.0], method="BFGS")
print("Optimiser:   intercept = {:,.0f}, slope = {:,.1f}".format(*opt.x))

2. The matrix formula β̂ = (X′X)⁻¹X′y, which is what every regression package computes.

X = np.column_stack([np.ones(N), x])
beta_hat = np.linalg.solve(X.T @ X, X.T @ y)
print("Matrix:      intercept = {:,.0f}, slope = {:,.1f}".format(*beta_hat))

3. statsmodels, which adds standard errors, t-statistics and p-values.

df = labor.to_pandas()
simple = smf.ols("wage ~ education", data=df).fit()
print("statsmodels: intercept = {:,.0f}, slope = {:,.1f}".format(*simple.params))

All three agree. The slope is close to $2,500 even though we left experience out. That works here only because education and experience are unrelated in this world. Lesson 3.2 breaks that.

Reading the output: coefficients, residuals, R²

print(simple.summary().tables[1])
print(f"R-squared: {simple.rsquared:.3f}")

plt.figure(figsize=(10, 5))
sns.scatterplot(x="education", y="wage", data=df, alpha=0.4, label="Workers")
grid = np.linspace(8, 20, 50)
plt.plot(grid, simple.params["Intercept"] + simple.params["education"] * grid, color="red", lw=3, label="OLS fit")
plt.title("Wage vs years of education"); plt.legend(); plt.show()

• Slope: the average difference in wages between workers who differ by one year of schooling.
• Intercept: the predicted wage at zero years of schooling, an extrapolation with no real meaning here.
• Residual: actual minus predicted wage, everything the model does not explain.
• R²: the share of wage variation explained. A low R² does not mean the slope is wrong. Causal work cares about an unbiased slope, not a high R².

Two mechanical facts of OLS you will use again: residuals average exactly zero, and they are uncorrelated with every regressor in the model.

print(f"Mean residual: {simple.resid.mean():.2e} | corr(resid, education): {np.corrcoef(simple.resid, x)[0, 1]:.2e}")

Logs and percentage effects

Economists usually model log wages: a coefficient b then means "about 100·b percent higher wages per year of schooling", which compares better across countries and decades.

df["log_wage"] = np.log(df["wage"])
log_model = smf.ols("log_wage ~ education + experience", data=df).fit()
b = log_model.params["education"]
print(f"Coefficient {b:.4f}  ->  about {100 * b:.1f}% per year (exact: {100 * (np.exp(b) - 1):.1f}%)")

mincer = smf.ols("log_wage ~ education + experience + I(experience**2)", data=df).fit()
b1, b2 = mincer.params["experience"], mincer.params["I(experience ** 2)"]
print(mincer.params.round(5))
print(f"Implied peak at {-b1 / (2 * b2):.0f} years of experience")
# Our DGP is linear in experience on the dollar scale, so the quadratic term is small and the
# implied 'peak' is unreliable. Real wage data usually show a clear concave profile.

Checklist and recap

• I can explain OLS as minimising squared residuals and reproduce it three ways.
• I can interpret a slope, an intercept and R², and know R² is not the goal.
• I can read log-outcome coefficients as approximate percentage effects.