Doubly Robust Estimation (AIPW)

36.1 Two models, one estimator (45 min)

Concept. G-computation relies on the outcome model; IPW relies on the propensity model. Augmented IPW (AIPW) uses both. For each customer:

ψ = μ₁(X) − μ₀(X)                       ← g-computation's prediction of their effect
    + D · (Y − μ₁(X)) / e(X)             ← IPW-weighted correction: joiners' prediction errors
    − (1 − D) · (Y − μ₀(X)) / (1 − e(X)) ← IPW-weighted correction: non-joiners' prediction errors

The ATE estimate is the average of ψ. Read it as "g-computation, plus a weighted average of g-computation's mistakes".

• If the outcome model is right, its prediction errors average to zero, so the correction does nothing harmful, even with a wrong propensity score.
• If the propensity model is right, the weighted correction removes whatever bias the outcome model has.

Either one being right is enough: that's double robustness. A bonus: the variance of ψ gives the standard error directly, SE = sd(ψ) / √n. No bootstrap needed. (ψ is the efficient influence function; you'll meet it again in Module 6.)

36.2 AIPW and its standard errors (50 min)

Type aipw and its two helpers into estimators.py. Two details:

• Cross-fitting: customers are split into 5 folds, and each customer's predictions come from models fitted on the other four. With simple regressions it hardly matters; with machine learning it matters a lot (36.5).
• Clipping: propensity scores are clipped to [0.01, 0.99] so no weight exceeds 100.

Create day36/aipw.py:

# day36/aipw.py   (python -m day36.aipw)
"""
AIPW = g-computation + an IPW correction of its mistakes. Consistent if EITHER the outcome model
OR the propensity model is right ('doubly robust'), with standard errors from a simple formula.
"""
from causalkit.estimators import aipw, g_computation, ipw
from causalkit.loyalty import public, truth, world

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

g, w, a = g_computation(df), ipw(df), aipw(df)
print(f"{'':<16}{'ATE':>8}{'ATT':>8}")
print(f"{'g-computation':<16}{g['ATE']:8.1f}{g['ATT']:8.1f}")
print(f"{'IPW':<16}{w['ATE']:8.1f}{w['ATT']:8.1f}")
print(f"{'AIPW':<16}{a['ATE']:8.1f}{a['ATT']:8.1f}   SEs {a['ATE_se']:.1f} / {a['ATT_se']:.1f}")
print(f"{'Truth':<16}{t['ATE']:8.1f}{t['ATT']:8.1f}")

python -m day36.aipw

                     ATE     ATT
g-computation      156.8   183.9
IPW                176.1   223.5
AIPW               157.8   184.6   SEs 2.1 / 2.7
Truth              158.4   186.0

AIPW sits with g-computation, close to the truth, not with IPW. When the outcome model is good, AIPW inherits its precision; the IPW part only corrects what the model gets wrong.

36.3 Double robustness, demonstrated (60 min)

Break each model on purpose, over 40 datasets. The "wrong" outcome model is Day 31's linear one; the "wrong" propensity model forgets engaged. Create day36/double_robustness.py:

# day36/double_robustness.py   (python -m day36.double_robustness)
"""
Double robustness, demonstrated: break the outcome model, the propensity model, or both,
across 40 simulated datasets. Wrong outcome model = the plain linear model of Day 31 (it misses
the curvature in pre_spend and tenure); wrong propensity model = forgets 'engaged' (a real driver of joining).
"""
import numpy as np
import polars as pl

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

WRONG_OUTCOME = "spend ~ " + " + ".join(COVARIATES)
WRONG_PS = "joined ~ np.log(pre_spend) + pre_orders + tenure_years + age + mobile"
NO_ENGAGED_OUTCOME = "spend ~ " + " + ".join(c for c in COVARIATES if c != "engaged")

cases = {
    "g-computation, wrong outcome model": lambda df: g_computation(df, WRONG_OUTCOME)["ATE"],
    "IPW, wrong propensity model": lambda df: ipw(df, propensity(df, WRONG_PS))["ATE"],
    "AIPW, both right": lambda df: aipw(df)["ATE"],
    "AIPW, outcome wrong, propensity right": lambda df: aipw(df, outcome_formula=WRONG_OUTCOME)["ATE"],
    "AIPW, outcome right, propensity wrong": lambda df: aipw(df, ps_formula=WRONG_PS)["ATE"],
    "AIPW, both models miss 'engaged'": lambda df: aipw(df, outcome_formula=NO_ENGAGED_OUTCOME,
                                                         ps_formula=WRONG_PS)["ATE"],
}
errors = {k: [] for k in cases}
for s in range(40):
    full = world(seed=3000 + s)
    df = public(full).to_pandas()
    ate = truth(full)["ATE"]
    for k, f in cases.items():
        errors[k].append(f(df) - ate)

with pl.Config(fmt_str_lengths=45, tbl_width_chars=100):
    print(pl.DataFrame([{"estimator": k, "average error": round(np.mean(v), 1), "SD": round(np.std(v), 1)}
                        for k, v in errors.items()]))

python -m day36.double_robustness

│ estimator                             ┆ average error ┆ SD   │
│ g-computation, wrong outcome model    ┆ -24.2         ┆ 4.9  │
│ IPW, wrong propensity model           ┆ 25.1          ┆ 14.0 │
│ AIPW, both right                      ┆ 0.6           ┆ 1.9  │
│ AIPW, outcome wrong, propensity right ┆ -0.5          ┆ 7.1  │
│ AIPW, outcome right, propensity wrong ┆ 0.5           ┆ 1.9  │
│ AIPW, both models miss 'engaged'      ┆ 25.8          ┆ 7.1  │

Read it row by row:

1. Each single-model method is biased by about $25 when its model is wrong.
2. AIPW is unbiased if either model is right. When the outcome model is wrong, AIPW pays in variance (SD 7.1 instead of 1.9), because the IPW correction now does real work.
3. The last row is the warning. Double robustness protects against a wrong functional form, not against a missing confounder. If engaged isn't in the data, no model can use it, and AIPW is as biased as everything else.

36.4 Do the confidence intervals work? (30 min)

Try it. Using the same 40 seeds (world(seed=3000 + s)), count how often AIPW's 95% interval (ATE ± 1.96 × SE) contains the true ATE. Create day36/coverage.py.

# day36/coverage.py   (python -m day36.coverage)
import numpy as np

from causalkit.estimators import aipw
from causalkit.loyalty import public, truth, world

hits, n = 0, 40
for s in range(n):
    full = world(seed=3000 + s)
    a = aipw(public(full).to_pandas())
    hits += abs(a["ATE"] - truth(full)["ATE"]) <= 1.96 * a["ATE_se"]
print(f"95% CIs that contain the true ATE: {hits} of {n}")

95% CIs that contain the true ATE: 38 of 40

36.5 Day 36 mini-project: machine learning and cross-fitting (65 min)

On Day 31, gradient boosting in g-computation overshot the ATT by $24. Does AIPW with boosted models do better, and does cross-fitting matter? Create day36/aipw_ml.py:

# day36/aipw_ml.py   (python -m day36.aipw_ml)
"""
Mini-project: machine learning for the nuisance models, with no formulas to get wrong.
Day 31 showed that plugging gradient boosting into g-computation overshoots. AIPW with cross-fitting
corrects it, and the influence-function SE gives a confidence interval. (This is the core of
Double Machine Learning, Module 6.)
"""
import polars as pl

from causalkit.estimators import aipw
from causalkit.loyalty import public, truth, world
from day31.specification_table import g_comp_boosting

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

rows = []
plain = g_comp_boosting(df)
rows.append({"estimator": "boosting g-computation (Day 31)", "ATE": round(plain["ATE"], 1),
             "ATT": round(plain["ATT"], 1), "ATT 95% CI": "none"})
for label, folds in [("AIPW + boosting, no cross-fitting", 1), ("AIPW + boosting, 5-fold cross-fitting", 5)]:
    a = aipw(df, outcome="boosting", ps_model="boosting", folds=folds)
    rows.append({"estimator": label, "ATE": round(a["ATE"], 1), "ATT": round(a["ATT"], 1),
                 "ATT 95% CI": f"[{a['ATT'] - 1.96 * a['ATT_se']:.1f}, {a['ATT'] + 1.96 * a['ATT_se']:.1f}]"})
rows.append({"estimator": "TRUTH", "ATE": round(t["ATE"], 1), "ATT": round(t["ATT"], 1), "ATT 95% CI": ""})
with pl.Config(fmt_str_lengths=45, tbl_width_chars=110):
    print(pl.DataFrame(rows))

python -m day36.aipw_ml

│ estimator                             ┆ ATE   ┆ ATT   ┆ ATT 95% CI     │
│ boosting g-computation (Day 31)       ┆ 170.0 ┆ 209.7 ┆ none           │
│ AIPW + boosting, no cross-fitting     ┆ 167.5 ┆ 205.2 ┆ [196.1, 214.2] │
│ AIPW + boosting, 5-fold cross-fitting ┆ 164.2 ┆ 198.9 ┆ [180.4, 217.5] │
│ TRUTH                                 ┆ 158.4 ┆ 186.0 ┆                │

• Without cross-fitting, each customer's prediction error comes from a model that has seen that customer, so errors look smaller than they are. The correction is too small, and the CI is too narrow: it misses the truth.
• With cross-fitting, the estimate moves towards the truth, and the CI is honest (and wider): it contains $186.
• With no formulas to write, boosting still leaves some bias at n = 20,000. With hand-crafted features (the flexible model), AIPW was better. Machine learning is a tool for when you don't know the functional form, not a replacement for thinking about it.

This combination, ML models + an orthogonal score like ψ + cross-fitting, is Double/Debiased Machine Learning, the core of Module 6.

36.6 Save your work with Git (20 min)

git add .
git commit -m "Day 36: AIPW, double robustness, coverage, cross-fitting with ML"
git push

Coming next (Days 37–40). Everything so far assumed that every confounder is measured. What if loyal fans, people who simply like NorthMart, both join and spend more, and nobody measures "liking NorthMart"? Days 37–38 use the random $20 invitation as an instrument; Day 39 races every estimator in both worlds; Day 40 is the capstone.

Day 36 checklist

• I can write the AIPW score and explain each of its three terms
• I can show by simulation that AIPW is consistent if either model is right
• I can explain why double robustness doesn't protect against a missing confounder
• I use cross-fitting whenever the nuisance models are machine-learning models
• I check a standard-error formula by its coverage in simulation