Regression Adjustment Done Right

31.1 The question and the naive answer (40 min)

Create day31/naive_and_ols.py:

# day31/naive_and_ols.py   (python -m day31.naive_and_ols)
"""
Who joined NorthMart's loyalty programme, and what did it do to their spending?
In this world, everything that drove joining is in the data (selection on observables holds).
So regression CAN work. Does it, and what exactly does it estimate?
"""
import polars as pl
import statsmodels.formula.api as smf

from causalkit.loyalty import public, truth, world

full = world()
df = public(full)
t = truth(full)
print(df.head(5))

naive = df.filter(pl.col("joined") == 1)["spend"].mean() - df.filter(pl.col("joined") == 0)["spend"].mean()
pdf = df.to_pandas()
specs = {
    "naive difference": "spend ~ joined",
    "linear controls": "spend ~ joined + pre_spend + pre_orders + engaged + tenure_years + age + mobile",
    "flexible controls": "spend ~ joined + pre_spend + I(pre_spend**2) + pre_orders + engaged"
                         " + tenure_years + I(tenure_years**0.5) + age + mobile",
}
print(f"\nJoiners spent ${naive:,.0f} more than non-joiners next year.")
for name, f in specs.items():
    fit = smf.ols(f, pdf).fit(cov_type="HC2")
    print(f"{name:<18} coefficient on joined {fit.params['joined']:+8.1f}  (SE {fit.bse['joined']:.1f})")
print(f"\nTruth: ATE {t['ATE']:.1f}   ATT {t['ATT']:.1f}   ATC {t['ATC']:.1f}")

python -m day31.naive_and_ols

Joiners spent $593 more than non-joiners next year.
naive difference   coefficient on joined   +593.0  (SE 9.9)
linear controls    coefficient on joined   +118.8  (SE 3.7)
flexible controls  coefficient on joined   +151.7  (SE 2.0)

Truth: ATE 158.4   ATT 186.0   ATC 132.7

Three lessons in four lines:

1. The naive gap ($593) is mostly selection: joiners were already big spenders.
2. Adding controls linearly removes most of it, but under-shoots ($119). Spending next year rises faster than linearly with spending last year, and a straight line can't capture that.
3. With pre_spend² and sqrt(tenure) the coefficient is close to the ATE. But which effect is it? The truth has three: ATE (everyone), ATT (joiners) and ATC (non-joiners). The effect is larger for big spenders, who are more likely to join, so ATT > ATE > ATC.

31.2 Why one regression coefficient isn't enough (40 min)

Concept. When the effect varies across customers, the coefficient on joined in spend ~ joined + X is a weighted average of customer-level effects, with weights chosen by the regression, not by you. It's usually closer to the effect for customers whose joining is hardest to predict, which is rarely the population a decision is about.

Recall the estimand card from Day 12. Before estimating, decide:

The estimator should then target that number directly. That's what g-computation does.

31.3 G-computation: predict both worlds for everyone (60 min)

Concept. Fit an outcome model that includes the treatment and the confounders. Then, for every customer, predict spending twice: once with joined = 1, once with joined = 0. The difference is each customer's predicted effect. Average it over everyone for the ATE, over joiners for the ATT, over non-joiners for the ATC. This is also called regression adjustment or the standardisation estimator, and it's Day 11's potential-outcomes table filled in by a model.

To let the effect vary, the outcome model interacts joined with the covariates. In causalkit/estimators.py, type in COVARIATES, OUTCOME_FORMULA, g_computation and bootstrap.

Create day31/g_computation.py:

# day31/g_computation.py   (python -m day31.g_computation)
"""
G-computation (regression adjustment done right): model the outcome separately for joiners and
non-joiners, predict what EVERY customer would spend in both cases, then average the difference
over whichever population the estimand names.
"""
from causalkit.estimators import bootstrap, g_computation
from causalkit.loyalty import public, truth, world

full = world()
df = public(full).to_pandas()
t = truth(full)

est = g_computation(df)
se = bootstrap(g_computation, df, reps=200)
for k in ["ATE", "ATT", "ATC"]:
    print(f"{k}: estimate {est[k]:7.1f}  (bootstrap SE {se[k]:.1f})   truth {t[k]:7.1f}")

python -m day31.g_computation

ATE: estimate   156.8  (bootstrap SE 1.9)   truth   158.4
ATT: estimate   183.9  (bootstrap SE 2.2)   truth   186.0
ATC: estimate   131.6  (bootstrap SE 2.1)   truth   132.7

All three estimands, each within about one standard error of the truth, from one fitted model.

31.4 Bootstrap standard errors (40 min)

Concept. G-computation's estimate is an average of predictions, so the regression's standard errors don't apply to it directly. The bootstrap gets a standard error by brute force: resample customers with replacement, re-run the whole estimator, repeat 200 times, and take the standard deviation of the 200 estimates. It works for any estimator, as long as you resample the units that were independently sampled (here, customers) and repeat every step, including model fitting.

Read bootstrap in estimators.py: it takes any estimator function and its arguments, so we'll reuse it for other estimators.

Try it. invited was randomised, so it isn't a confounder. Add + invited to the outcome formula. What happens to the estimates and their standard errors? Create day31/add_invited.py.

# day31/add_invited.py   (python -m day31.add_invited)
from causalkit.estimators import OUTCOME_FORMULA, bootstrap, g_computation
from causalkit.loyalty import public, world

df = public(world()).to_pandas()
for name, f in [("without invited", OUTCOME_FORMULA), ("with invited", OUTCOME_FORMULA + " + invited")]:
    est = g_computation(df, f)
    se = bootstrap(g_computation, df, reps=100, formula=f)
    print(f"{name:<16} ATE {est['ATE']:.1f} (SE {se['ATE']:.2f})   ATT {est['ATT']:.1f} (SE {se['ATT']:.2f})")

without invited  ATE 156.8 (SE 2.01)   ATT 183.9 (SE 2.18)
with invited     ATE 156.8 (SE 2.10)   ATT 184.0 (SE 2.30)

31.5 Day 31 mini-project: a specification table (70 min)

The estimate is only as good as the outcome model. Let's see how much the model matters. Create day31/specification_table.py:

# day31/specification_table.py   (python -m day31.specification_table)
"""
Mini-project: the outcome model decides the answer. G-computation with three outcome models:
a linear one, the flexible one, and a machine-learning model (gradient boosting) that finds the
shape on its own.
"""
import numpy as np
import polars as pl
from sklearn.ensemble import HistGradientBoostingRegressor

from causalkit.estimators import COVARIATES, g_computation
from causalkit.loyalty import public, truth, world

full = world()
df = public(full).to_pandas()
t = truth(full)

LINEAR = "spend ~ " + " + ".join(COVARIATES)


def g_comp_boosting(df):
    X = df[COVARIATES]
    preds = {}
    for g in [0, 1]:
        m = HistGradientBoostingRegressor(max_iter=300, learning_rate=0.05, random_state=0)
        m.fit(X[df.joined == g], df.spend[df.joined == g])
        preds[g] = m.predict(X)
    diff = preds[1] - preds[0]
    d = (df.joined == 1).to_numpy()
    return {"ATE": diff.mean(), "ATT": diff[d].mean(), "ATC": diff[~d].mean()}


if __name__ == "__main__":
    rows = []
    for name, est in [("linear outcome model", g_computation(df, LINEAR)),
                      ("flexible outcome model", g_computation(df)),
                      ("gradient boosting", g_comp_boosting(df))]:
        rows.append({"outcome model": name, **{k: round(v, 1) for k, v in est.items()}})
    rows.append({"outcome model": "TRUTH", **{k: round(t[k], 1) for k in ["ATE", "ATT", "ATC"]}})
    with pl.Config(fmt_str_lengths=30):
        print(pl.DataFrame(rows))

python -m day31.specification_table

│ outcome model          ┆ ATE   ┆ ATT   ┆ ATC   │
│ linear outcome model   ┆ 139.5 ┆ 195.0 ┆ 88.0  │
│ flexible outcome model ┆ 156.8 ┆ 183.9 ┆ 131.6 │
│ gradient boosting      ┆ 170.0 ┆ 209.7 ┆ 133.1 │
│ TRUTH                  ┆ 158.4 ┆ 186.0 ┆ 132.7 │

• The linear model gets the ATC badly wrong ($88 vs $133): it extrapolates a straight line into regions where it doesn't fit.
• The flexible model is right because we wrote it with the true functional form. In real data you never know that form.
• Gradient boosting learns the shape from the data, but overstates the ATT by $24. Machine-learning models shrink towards the average to avoid overfitting, and that shrinkage biases a difference of two predictions. Day 36 shows how to fix this.

The specification table is a habit: if your conclusion changes across reasonable outcome models, the model is doing the work, not the data. Write two sentences in day31/NOTES.md on which row you'd report and why.

31.6 Save your work with Git (20 min)

git add .
git commit -m "Day 31: loyalty simulator, g-computation, bootstrap, specification table"
git push

Day 31 checklist

• I can explain why the naive difference and a single regression coefficient may not answer the decision question
• I can choose between ATE, ATT and ATC from the decision being made
• I can estimate all three by g-computation, with bootstrap standard errors
• I can explain why controlling for a cause of treatment only (like invited) doesn't help
• I check how much my answer depends on the outcome model