nb2.3 — FWL Residualization, Visualized

Companion to Ch 2.3, “The Frisch–Waugh–Lovell Theorem.”

In the chapter we proved a single, slightly surprising fact. When a paper reports a coefficient “controlling for the other variables,” the word controlling is not a filter and not an experiment — it is a subtraction. You residualize the other regressors out of both the outcome and the regressor you care about, then run a plain two-variable regression on what is left. That is the Frisch–Waugh–Lovell (FWL) theorem, and the slope it produces is identical, to machine precision, to the multiple-regression coefficient.

This notebook makes that identity concrete on Leah’s six firms. We will:

  1. Run the full multiple regression of patents on R&D and firm size and recover the chapter’s targets \hat\beta_1 \approx 1.838, \hat\beta_2 \approx 3.676, \hat\beta_0 \approx 1.919.
  2. Do FWL by hand — residualize y on size, residualize x on size, regress leftover on leftover — and watch the slope land exactly on \hat\beta_1.
  3. Visualize it as a partial-regression (added-variable) plot: the scatter of residualized patents against residualized R&D, with the FWL slope drawn through it.
  4. Show that demeaning = the within transformation = fixed effects: regressing on group dummies gives the identical slope to regressing on group-demeaned data.

Then a Your Turn that absorbs a high-cardinality fixed effect by demeaning, and times it against the brute-force dummy regression — the trick that makes modern fixed-effects software fast.

Setup

We fix a seed so the notebook reproduces identically (the rule from CONVENTIONS.md), and force matplotlib’s non-interactive Agg backend so everything runs headless — on a laptop, a server, or in CI — without needing a display.

import matplotlib
matplotlib.use("Agg")  # headless, non-interactive backend

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)  # one seeded generator for the whole notebook

print("numpy", np.__version__, "| pandas", pd.__version__, "| statsmodels", sm.__version__)

Leah’s six firms

Leah studies innovation. Her hunch: firms that spend more on R&D produce more patents. The obvious objection is that big firms do everything more — they spend more on R&D and file more patents simply because they are big. So a raw correlation might be measuring nothing but size. The honest question is the partial one: does R&D predict patents after netting out firm size?

Here is the tiny six-firm dataset from the chapter, small enough to check by hand. y is patents per year, x is R&D spending ($millions), z is firm size (log assets). Notice R&D and size move together across these firms — exactly the entanglement Leah worries about.

firms = pd.DataFrame({
    "firm":    [1, 2, 3, 4, 5, 6],
    "patents": [10, 14, 22, 24, 33, 39],   # y
    "rnd":     [2,  3,  5,  4,  7,  8],     # x  (R&D, $M)
    "size":    [1,  2,  3,  4,  5,  6],     # z  (log assets)
}).astype({"patents": float, "rnd": float, "size": float})

y = firms["patents"].to_numpy()
x = firms["rnd"].to_numpy()
z = firms["size"].to_numpy()

print(firms.to_string(index=False))
print(f"\nmeans:  y-bar={y.mean():.3f}   x-bar={x.mean():.3f}   z-bar={z.mean():.3f}")
print(f"corr(R&D, size) = {np.corrcoef(x, z)[0, 1]:.3f}   <- R&D and size are entangled")

Step 0: the target — the full multiple regression

Before any FWL trickery, run the regression Leah would actually report:

y_i = \beta_0 + \beta_1\,x_i + \beta_2\,z_i + \varepsilon_i.

The coefficient \hat\beta_1 is “the effect of R&D, controlling for size.” This is the number FWL must reproduce. The chapter’s targets are \hat\beta_1 \approx 1.838, \hat\beta_2 \approx 3.676, \hat\beta_0 \approx 1.919.

X_full = sm.add_constant(np.column_stack([x, z]))   # columns: [const, rnd, size]
full_fit = sm.OLS(y, X_full).fit()
b0, b1, b2 = full_fit.params

print(f"beta_0 (const) = {b0:.4f}    (target ~ 1.919)")
print(f"beta_1 (R&D)   = {b1:.4f}    (target ~ 1.838)  <- the one we will reproduce")
print(f"beta_2 (size)  = {b2:.4f}    (target ~ 3.676)")

# Confirm we hit the chapter's numbers.
assert np.isclose(b1, 1.83783784, atol=1e-6)
assert np.isclose(b2, 3.67567568, atol=1e-6)
assert np.isclose(b0, 1.91891892, atol=1e-6)
print("\nMatches the chapter's targets.")

The short regression: what omitting size does

For contrast, run the short regression — patents on R&D alone, with size left out:

y_i = \gamma_0 + \gamma_1\,x_i + u_i.

Because nothing is residualized, the raw R&D variable still carries the size signal inside it, so patents respond to both the genuine R&D effect and the size effect riding along with R&D. The chapter reports the short slope as about 4.65 — more than double \hat\beta_1 = 1.838. That gap is the omitted-variable bias from leaving size out; Ch 2.5 turns it into a formula.

short_fit = sm.OLS(y, sm.add_constant(x)).fit()
g1 = short_fit.params[1]
print(f"short slope (R&D only)         = {g1:.4f}   (target ~ 4.646)")
print(f"full beta_1 (R&D | size)       = {b1:.4f}")
print(f"gap = short - full             = {g1 - b1:.4f}   <- the omitted-variable bias from dropping size")
assert np.isclose(g1, 4.64596273, atol=1e-6)

Steps 1–2: residualize both sides against “everything else”

Now FWL by hand. “Everything else” here is the \mathbf{X}_2 block: the intercept plus firm size. We sweep that block out of both variables.

  • Step 1 — residualize the key regressor: regress R&D on size (with intercept) and keep the residual \tilde x_i = x_i - \hat x_i. This is the part of R&D that size cannot explain — the firm-specific R&D intensity that has nothing to do with mere bigness.
  • Step 2 — residualize the outcome: regress patents on size and keep \tilde y_i = y_i - \hat y_i, the part of patenting that size cannot explain.

Both residual vectors average to zero (residuals from a regression with an intercept always do), so in Step 3 we will not need an intercept.

Z = sm.add_constant(z)   # the X2 block: intercept + size

# statsmodels .resid gives y - fitted directly; equivalent to the residual maker M2 applied to each.
x_tilde = sm.OLS(x, Z).fit().resid   # R&D with size netted out
y_tilde = sm.OLS(y, Z).fit().resid   # patents with size netted out

resid_tbl = pd.DataFrame({"firm": firms["firm"], "x_tilde": x_tilde, "y_tilde": y_tilde})
print(resid_tbl.to_string(index=False, float_format=lambda v: f"{v:+.3f}"))
print(f"\nsum(x_tilde) = {x_tilde.sum():.2e}   sum(y_tilde) = {y_tilde.sum():.2e}   (both ~ 0)")
assert np.isclose(x_tilde.sum(), 0.0, atol=1e-10)
assert np.isclose(y_tilde.sum(), 0.0, atol=1e-10)

Step 3: regress leftover on leftover — the punchline

Run a simple, no-intercept regression of the residualized outcome on the residualized regressor,

\tilde y_i = \beta_1\,\tilde x_i + (\text{error}), \qquad \hat\beta_1 = \frac{\sum_i \tilde x_i\,\tilde y_i}{\sum_i \tilde x_i^2}.

FWL guarantees this slope equals the multiple-regression \hat\beta_1 exactly — not approximately, not “close,” but to machine precision, on these six firms and on any dataset. We check the difference is below 10^{-10}.

beta_fwl = (x_tilde @ y_tilde) / (x_tilde @ x_tilde)   # no-intercept slope = sum(xt*yt)/sum(xt^2)

print(f"full multiple-regression beta_1 = {b1:.12f}")
print(f"FWL residualize-then-regress    = {beta_fwl:.12f}")
print(f"absolute difference             = {abs(b1 - beta_fwl):.3e}")

assert np.isclose(b1, beta_fwl, atol=1e-10), "FWL slope must equal the multiple-regression coef"
assert np.isclose(beta_fwl, 1.83783784, atol=1e-6), "and it must match the chapter's 1.838"
print("\nIdentical to ~1e-10, and equal to the chapter's 1.838. FWL holds.")

One side is enough for the point estimate (but not for the rest)

The chapter’s §3 algebra notes that residualizing only the regressor still nails the point estimate, because \tilde x is already orthogonal to the size block, so the size-explained part of y contributes nothing to the slope. We confirm it below — (x_tilde @ y) / (x_tilde @ x_tilde) gives the same 1.838. The catch: you must clean both sides to get the correct residuals, standard error, and R^2. So “residualize both” is the safe default to memorize.

beta_one_side = (x_tilde @ y) / (x_tilde @ x_tilde)   # raw y, residualized x only
print(f"clean both sides : {beta_fwl:.10f}")
print(f"clean x only     : {beta_one_side:.10f}")
assert np.isclose(beta_one_side, b1, atol=1e-10)
print("Same point estimate -- but only 'clean both' gives the right residuals/SE/R^2.")

Visualize: the partial-regression (added-variable) plot

Because FWL collapses a multivariate coefficient into a bivariate slope, you can plot it. The scatter of \tilde y against \tilde x, with the FWL slope drawn through the origin, is the honest two-dimensional picture of a coefficient that otherwise lives in higher dimensions. It is the single best way to see whether a control coefficient is carried by the bulk of the data or by one or two influential firms. The fitted line’s slope is exactly \hat\beta_1 = 1.838.

fig, ax = plt.subplots(figsize=(6.5, 5))
ax.axhline(0, color="0.8", lw=0.8, zorder=0)
ax.axvline(0, color="0.8", lw=0.8, zorder=0)
ax.scatter(x_tilde, y_tilde, s=70, color="#1f77b4", zorder=3, label="firms (residualized)")

# FWL line through the origin (both residual vectors are mean-zero).
xs = np.linspace(x_tilde.min() - 0.3, x_tilde.max() + 0.3, 100)
ax.plot(xs, beta_fwl * xs, color="#d62728", lw=2,
        label=f"FWL slope = {beta_fwl:.3f}", zorder=2)

for i, fid in enumerate(firms["firm"]):
    ax.annotate(f"firm {fid}", (x_tilde[i], y_tilde[i]),
                textcoords="offset points", xytext=(6, 4), fontsize=8)

ax.set_xlabel(r"$\tilde{x}$ = R&D residualized on size  (R&D beyond what size predicts)")
ax.set_ylabel(r"$\tilde{y}$ = patents residualized on size")
ax.set_title("Partial-regression plot: a multivariate coef as a 2-D slope")
ax.legend(loc="upper left", frameon=False)
fig.tight_layout()
fig.savefig("fwl_partial_regression.png", dpi=110)
plt.close(fig)
print("Saved fwl_partial_regression.png  (slope of the red line = multiple-regression beta_1 = 1.838)")

Demeaning = the within transformation = fixed effects

The most consequential special case of FWL. Take the \mathbf{X}_2 block to be a full set of group dummies. Regressing any variable on a complete set of group dummies fits each observation its own group’s mean, so the residual maker subtracts the group-specific mean:

(\mathbf{M}_2\mathbf{y})_i = y_i - \bar y_{g(i)}.

That operation — subtract each group’s own mean — is the within transformation. FWL then gives:

Dummies ⇔ demeaning. Regressing y on x plus a full set of group dummies gives the identical slope on x as regressing the within-demeaned y on the within-demeaned x.

Group Leah’s firms by industry: firms 1–3 are biotech, firms 4–6 are hardware. We compute the slope two ways and confirm they agree to the last digit. (Note the slope here, \approx 3.725, differs from §2’s 1.838 — controlling for industry asks a different question than controlling for size, and FWL makes that visible as a different residualization.)

firms2 = firms.copy()
firms2["industry"] = ["biotech", "biotech", "biotech", "hardware", "hardware", "hardware"]

# --- Way 1: within transformation (subtract each industry's own mean), no-intercept slope ---
grp_mean = firms2.groupby("industry")[["patents", "rnd"]].transform("mean")
x_within = firms2["rnd"].to_numpy()     - grp_mean["rnd"].to_numpy()
y_within = firms2["patents"].to_numpy() - grp_mean["patents"].to_numpy()
slope_within = (x_within @ y_within) / (x_within @ x_within)

# --- Way 2: regress y on R&D + an industry dummy (hardware = 1) ---
hardware = (firms2["industry"] == "hardware").astype(float).to_numpy()
X_dummy = sm.add_constant(np.column_stack([firms2["rnd"].to_numpy(), hardware]))
slope_dummy = sm.OLS(firms2["patents"].to_numpy(), X_dummy).fit().params[1]

print(f"within-demeaning slope  = {slope_within:.12f}")
print(f"industry-dummy slope    = {slope_dummy:.12f}")
print(f"absolute difference     = {abs(slope_within - slope_dummy):.3e}")
assert np.isclose(slope_within, slope_dummy, atol=1e-10)
assert np.isclose(slope_within, 3.725, atol=1e-3)
print("\nIdentical to ~1e-10. Demeaning IS fixed effects -- this is what powers Weeks 3-4.")

Your Turn — absorb a high-cardinality fixed effect by demeaning

The equivalence above is not just elegant; it is fast. The dummy way builds and inverts an \mathbf{X}'\mathbf{X} with one column per group — ruinous when there are thousands of groups. The demeaning way never forms the dummy matrix at all: it just subtracts group means, a cost linear in the data. That is exactly how packages like pyfixest “absorb” fixed effects.

Below we build a panel with 800 firms (a high-cardinality firm fixed effect), where each firm has its own baseline patenting level \alpha_i and the true within-firm R&D slope is \beta = 2.0. We then estimate \beta two ways and time them:

  1. Demeaning — subtract each firm’s own mean from y and x, then one tiny no-intercept slope.
  2. Brute-force dummies — regress on R&D plus 800 firm-dummy columns.

Both give the same \hat\beta. Watch the timing gap.

import time

n_firms, n_years = 800, 8
firm_ids = np.repeat(np.arange(n_firms), n_years)
alpha = rng.normal(50, 15, n_firms)[firm_ids]          # firm fixed effect (baseline patents)
x_panel = rng.normal(10, 3, n_firms * n_years)         # within-firm R&D
BETA_TRUE = 2.0
y_panel = alpha + BETA_TRUE * x_panel + rng.normal(0, 5, n_firms * n_years)

panel = pd.DataFrame({"firm": firm_ids, "rnd": x_panel, "patents": y_panel})

# --- Method 1: demeaning (absorb the firm FE) ---
t0 = time.perf_counter()
gm = panel.groupby("firm")[["patents", "rnd"]].transform("mean")
xw = panel["rnd"].to_numpy()     - gm["rnd"].to_numpy()
yw = panel["patents"].to_numpy() - gm["patents"].to_numpy()
beta_demean = (xw @ yw) / (xw @ xw)
t_demean = time.perf_counter() - t0

# --- Method 2: brute-force firm dummies ---
t0 = time.perf_counter()
D = pd.get_dummies(panel["firm"], prefix="f", drop_first=True).to_numpy(dtype=float)
X_big = np.column_stack([np.ones(len(panel)), panel["rnd"].to_numpy(), D])
beta_dummies = np.linalg.lstsq(X_big, panel["patents"].to_numpy(), rcond=None)[0][1]
t_dummies = time.perf_counter() - t0

print(f"true beta                       = {BETA_TRUE}")
print(f"demeaning  beta = {beta_demean:.6f}   in {t_demean*1e3:7.2f} ms   (design: 1 column)")
print(f"dummies    beta = {beta_dummies:.6f}   in {t_dummies*1e3:7.2f} ms   (design: {X_big.shape[1]} columns)")
print(f"\nsame estimate?  {np.isclose(beta_demean, beta_dummies, atol=1e-6)}")
print(f"demeaning was ~{t_dummies / t_demean:.0f}x faster, and never built the {D.shape[1]}-column dummy matrix.")
assert np.isclose(beta_demean, beta_dummies, atol=1e-6)

# YOUR TURN:
#   1. Bump n_firms to 5000. How does the timing gap scale? Does the dummy method still finish?
#   2. The demeaning point estimate is right, but its naive SE would use the wrong degrees of
#      freedom (Ch 2.3 sec.7). What is the correct df here -- N - K counting all the absorbed
#      dummies? Compute it for this panel.
#   3. Add a SECOND fixed effect (year) and absorb both by iterating the demeaning
#      (alternating projections) until it converges. Confirm it matches a two-way dummy regression.

What to carry forward

Three things from this notebook do real work for the rest of the camp.

The reading rule. A multiple-regression coefficient is always a two-variable slope in disguise — the simple slope of the outcome on the regressor, after both have been residualized on the controls. “Controlling for size” is not a disclaimer; it names the exact subtraction that produced the number. The partial-regression plot is that slope made visible.

The equivalence that powers panels. Regressing on a full set of group dummies is identical to subtracting group means and regressing on the leftovers — to machine precision, as we verified twice. Demeaning is residualizing on dummies. When Week 3 writes y_{it} = \alpha_i + \beta x_{it} + \varepsilon_{it} and estimates \beta by subtracting each firm’s time-average, you will already know why the within estimator equals the dummy estimator, and why it is cheap: you saw demeaning beat the dummy matrix by an order of magnitude on 800 firms.

The threat-naming discipline. FWL draws a hard line: a coefficient controls for exactly the columns you residualized out, and nothing else. Everything outside that block still sits inside your residuals, uncontrolled. That is why every spec must name its controls and its identifying assumption — the list of what you residualized is the boundary of what “holding constant” covers.