Notebook 3.3 — Entropy Balancing vs. IPW vs. AIPW

Maya has a propensity score and three different ways to use it. In Ch 3.2 she matched on it and lost a third of her control group. This chapter’s three tools all start from the same place — the propensity score and the conditional-independence assumption (CIA): once we condition on the observed covariates \mathbf{X}, who enrolled is as good as random — but they push past matching’s limits.

We run all three on one simulated dataset where we know the true effect is \tau = -0.8 percentage points (enrolling in a free financial-literacy course lowers a student’s first-loan interest rate by 0.8pp). The recurring theme is stability: all three are consistent in theory, but they misbehave very differently when the propensity model is even slightly wrong — which, in real data, it always is.

The targets to watch for (these reproduce the chapter’s worked example): the naive difference is badly biased toward zero; IPW \approx -0.785, AIPW \approx -0.779, both close to the truth.

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")  # non-interactive backend: render to buffers, never pop a window
import matplotlib.pyplot as plt

import statsmodels
import statsmodels.formula.api as smf
from scipy.optimize import minimize  # used as a cross-check on our hand-rolled solver

# Notebook-wide RNG for any auxiliary randomness (bootstraps, jitter in plots).
rng = np.random.default_rng(42)

print("pandas", pd.__version__, "| statsmodels", statsmodels.__version__)

1. Simulate Maya’s data — a world where we know the truth

We build the chapter’s data-generating process exactly. Two standardized covariates, income and credit, drive both selection into treatment and the outcome — that is what makes them confounders. Enrollment follows a logit in those covariates (so higher-income, higher-credit students enroll more often — selection), and the true treatment effect is a constant -0.8.

Because this is a simulation, we can write down both potential outcomes for every student — Y_i(0) (rate if they do not enroll) and Y_i(1) = Y_i(0) - 0.8 (rate if they do) — and then reveal only the one that matches their actual treatment, exactly as the fundamental problem of causal inference forces on us in real data.

We seed this data cell with 33 to reproduce the chapter’s published numbers bit-for-bit; the notebook-wide rng (seed 42) is reserved for everything else.

# --- data-generating process (chapter seed 33 -> reproduces the chapter's numbers) ---
data_rng = np.random.default_rng(33)
N = 4000

income = data_rng.normal(0, 1, N)   # standardized household income
credit = data_rng.normal(0, 1, N)   # standardized starting credit score

# selection into treatment: enrollment depends on income & credit (the propensity TRUTH)
lin    = 0.9 * income + 0.7 * credit
p_true = 1.0 / (1.0 + np.exp(-lin))                 # true propensity score
D      = (data_rng.uniform(size=N) < p_true).astype(int)

# outcomes: TRUE effect is -0.8; the outcome ALSO depends on income & credit
tau_true = -0.8
y0 = 9.0 - 0.6 * income - 0.5 * credit + data_rng.normal(0, 1, N)  # rate if NOT enrolled
y1 = y0 + tau_true                                                 # rate if enrolled
Y  = np.where(D == 1, y1, y0)        # we observe only ONE potential outcome per student

df = pd.DataFrame({"Y": Y, "D": D, "income": income, "credit": credit})

# naive difference in means -- contaminated by selection
naive = df.loc[df.D == 1, "Y"].mean() - df.loc[df.D == 0, "Y"].mean()
print(f"N = {N},  treated = {int(df.D.sum())},  control = {int((1-df.D).sum())}")
print(f"true effect          : {tau_true:+.3f}")
print(f"naive diff-in-means  : {naive:+.3f}   <- biased toward zero by selection")

The naive difference comes out around -1.5 — almost twice the true effect. Why too big (in magnitude)? Enrollees have higher income and credit, and higher income/credit push the rate down on their own (-0.6 and -0.5 coefficients). So enrollees would have had lower rates even without the program. The naive comparison hands the program credit for a head start the students already had. That gap between -1.5 and -0.8 is exactly the selection bias every method below is built to remove.

2. IPW: build a pseudo-population by reweighting

2.1 Estimate the propensity score, then form weights

We estimate \hat e(\mathbf{X}) = \Pr(D=1 \mid \mathbf{X}) with a logit, then turn the scores into weights. We will build two flavors of weight and two estimators.

The raw Horvitz–Thompson weights are 1/\hat e for the treated and 1/(1-\hat e) for controls. The stabilized weights multiply by the marginal treatment share, \hat p/\hat e and (1-\hat p)/(1-\hat e) — same estimand, much tamer variance. We also clip (trim) the scores to [0.01, 0.99] so a single mis-estimated tiny probability cannot detonate the whole estimate.

# estimate the propensity score (a logit), as in Ch 3.2
ps = smf.logit("D ~ income + credit", data=df).fit(disp=0)
ehat_raw = ps.predict(df).to_numpy()          # raw fitted propensity scores
ehat     = ehat_raw.clip(0.01, 0.99)          # clip = trimming extreme scores

Dv = df.D.to_numpy()
Yv = df.Y.to_numpy()
p_marg = df.D.mean()                          # marginal P(D=1)

# raw (Horvitz-Thompson-style) weights -- used to show the explosion
w_raw  = np.where(Dv == 1, 1.0 / ehat, 1.0 / (1.0 - ehat))
# stabilized weights -- the modern default
w_stab = np.where(Dv == 1, p_marg / ehat, (1.0 - p_marg) / (1.0 - ehat))

print(f"fitted propensity range : [{ehat_raw.min():.4f}, {ehat_raw.max():.4f}]")
print(f"raw  weights : max = {w_raw.max():6.2f}")
print(f"stab weights : max = {w_stab.max():6.2f}")

2.2 Horvitz–Thompson vs. Hájek (the self-normalizing fix)

The Horvitz–Thompson estimator divides the weighted outcome sum by N:

\hat\tau_{\text{HT}} = \frac{1}{N}\sum_i \frac{D_i Y_i}{\hat e_i} - \frac{1}{N}\sum_i \frac{(1-D_i) Y_i}{1-\hat e_i}.

The defect: the weights do not add up to N in any given sample, so dividing by N is dividing by the wrong thing — the estimate can wander outside the range of the data. The Hájek fix divides each arm’s weighted outcome sum by that arm’s total weight (a self-normalizing weighted mean):

\hat\tau_{\text{H\'ajek}} = \frac{\sum_i (D_i/\hat e_i)\,Y_i}{\sum_i D_i/\hat e_i} - \frac{\sum_i ((1-D_i)/(1-\hat e_i))\,Y_i}{\sum_i (1-D_i)/(1-\hat e_i)}.

Use Hájek, not raw Horvitz–Thompson. Below we compute both so you can see HT misbehave.

# --- Horvitz-Thompson: divide weighted outcome sum by N (the fragile version) ---
ht = np.mean(Dv * Yv / ehat) - np.mean((1 - Dv) * Yv / (1.0 - ehat))

# --- Hajek: self-normalizing weighted mean per arm (the stable version) ---
def hajek_mean(w, mask):
    ww = w[mask]
    return np.sum(ww * Yv[mask]) / np.sum(ww)

ipw_stab = hajek_mean(w_stab, Dv == 1) - hajek_mean(w_stab, Dv == 0)
ipw_raw  = hajek_mean(w_raw,  Dv == 1) - hajek_mean(w_raw,  Dv == 0)  # Hajek is scale-free

print(f"true effect           : {tau_true:+.3f}")
print(f"Horvitz-Thompson (/N) : {ht:+.3f}   <- numerically wild")
print(f"IPW Hajek (stabilized): {ipw_stab:+.3f}   <- our IPW estimate")
print(f"(Hajek is invariant to weight scaling, so raw vs stab Hajek agree: "
      f"{ipw_raw:+.3f})")

2.3 Diagnose the weights: max weight and effective sample size

The IPW estimate (\approx -0.785) lands near the truth, but the real discipline is to interrogate the weights. The headline diagnostic is the effective sample size

N_{\text{eff}} = \frac{\big(\sum_i w_i\big)^2}{\sum_i w_i^2},

which tells you how many “real” observations your weighted sample is worth. If a few units carry most of the weight, N_{\text{eff}} collapses far below N and your data simply do not contain a credible counterfactual. Always look at your weights.

def eff_n(w):
    return w.sum()**2 / np.sum(w**2)

for name, w in [("raw (HT-style)", w_raw), ("stabilized", w_stab)]:
    print(f"{name:16s}: max = {w.max():6.2f}   "
          f"99th pct = {np.percentile(w, 99):5.2f}   "
          f"eff. N = {eff_n(w):6.0f} of {N}")

2.4 When IPW explodes — thin overlap

The danger is not visible on this moderate-selection data. To see the explosion, crank up the selection strength so some students have covariates that make treatment nearly certain or nearly impossible — the thin-overlap regime. Then a treated unit with \hat e \to 0 (or a control with \hat e \to 1) gets an astronomically large weight, and a handful of such units swamp everything. We show that stabilizing + trimming tames the max weight and rescues the effective sample size — and we plot the weight distributions side by side.

# strong-selection world -> thin overlap -> exploding weights
strong_rng = np.random.default_rng(101)
Ns = 4000
inc_s = strong_rng.normal(0, 1, Ns)
cre_s = strong_rng.normal(0, 1, Ns)
lin_s = 2.5 * inc_s + 2.0 * cre_s                     # much steeper selection
p_s   = 1.0 / (1.0 + np.exp(-lin_s))
D_s   = (strong_rng.uniform(size=Ns) < p_s).astype(int)

ps_s   = smf.logit("D ~ income + credit",
                   data=pd.DataFrame({"D": D_s, "income": inc_s, "credit": cre_s})).fit(disp=0)
ehat_s = ps_s.predict(pd.DataFrame({"income": inc_s, "credit": cre_s})).to_numpy()

# unstabilized (no trim) vs stabilized (with trim)
w_explode = np.where(D_s == 1, 1.0 / ehat_s.clip(1e-6, 1),
                     1.0 / (1.0 - ehat_s).clip(1e-6, 1))
ec = ehat_s.clip(0.01, 0.99)
pm_s = D_s.mean()
w_tame = np.where(D_s == 1, pm_s / ec, (1.0 - pm_s) / (1.0 - ec))

print(f"thin-overlap propensity range: [{ehat_s.min():.4f}, {ehat_s.max():.4f}]")
print(f"UNSTABILIZED : max w = {w_explode.max():7.1f}   eff. N = {eff_n(w_explode):5.0f} of {Ns}")
print(f"stab + trim  : max w = {w_tame.max():7.1f}   eff. N = {eff_n(w_tame):5.0f} of {Ns}")

fig, ax = plt.subplots(1, 2, figsize=(11, 3.6))
ax[0].hist(w_explode, bins=60, color="firebrick", alpha=0.8)
ax[0].set_title(f"Unstabilized weights (max {w_explode.max():.0f})")
ax[0].set_xlabel("weight"); ax[0].set_ylabel("count")
ax[1].hist(w_tame, bins=60, color="steelblue", alpha=0.8)
ax[1].set_title(f"Stabilized + trimmed (max {w_tame.max():.0f})")
ax[1].set_xlabel("weight")
fig.suptitle("IPW weights explode under thin overlap; stabilizing + trimming helps")
fig.tight_layout()
plt.close(fig)
print("[figure rendered to buffer]")

3. Entropy balancing: solve for the weights directly

3.1 The idea

We never actually cared about the propensity score — balance was the goal, the logit just an indirect, fish-for-a-spec route to it. Entropy balancing (Hainmueller 2012) solves for balance directly. To estimate the ATT (effect for the kind of student who enrolled), reweight the control units so their weighted covariate means exactly equal the treated group’s means, and among all weight vectors that achieve that, pick the one closest to uniform. “Closest to uniform” is made precise by maximizing entropy:

\min_{\{w_i\}} \sum_{i\in\text{control}} w_i \log\frac{w_i}{q_i} \quad\text{s.t.}\quad \sum_i w_i X_{ik} = \bar X_k^{\text{treated}}\ \ \forall k,\quad \sum_i w_i = 1,\quad w_i \ge 0,

with base weights q_i = 1/N_c (uniform). The objective is the Kullback–Leibler divergence of the weights from uniform; minimizing it keeps the weights flat so no single unit dominates — exactly the explosion IPW could not control.

3.2 The closed-form weights, and a tiny convex solve

The Lagrangian gives weights that are exponential in the covariates, w_i \propto q_i\exp\!\big(-\sum_k \lambda_k X_{ik}\big), where the \lambda_k are Lagrange multipliers — one per balance condition, not one per unit. So the whole thing reduces to a small, convex optimization in K unknowns. We solve the convex dual with Newton’s method (the standard, robust route for entropy balancing): the dual objective is \log\sum_i q_i\exp(-\mathbf{Z}_i\boldsymbol\lambda) where \mathbf{Z}_i = \mathbf{X}_i - \text{target}, its gradient is the weighted mean of \mathbf{Z} (which we want to drive to zero), and its Hessian is the weighted covariance of \mathbf{Z}. Newton converges to machine precision in a handful of steps.

def entropy_balance(X_control, targets, base=None, tol=1e-12, maxit=100):
    """Max-entropy ATT weights on controls s.t. weighted column means of X_control
    equal `targets`, kept as close to uniform as the constraints allow.

    Solves the convex dual D(lambda) = log sum_i q_i exp(-Z_i . lambda), where
    Z = X_control - targets, by Newton's method.  Returns (w, max_abs_imbalance, iters).
    """
    Z = np.asarray(X_control, float) - np.asarray(targets, float)
    Nc = Z.shape[0]
    q = np.ones(Nc) / Nc if base is None else base / base.sum()
    lmb = np.zeros(Z.shape[1])
    for it in range(1, maxit + 1):
        a = -Z @ lmb
        a -= a.max()                       # stabilize the exponential
        ew = q * np.exp(a)
        w = ew / ew.sum()                  # normalized weights (sum to 1)
        grad = Z.T @ w                     # weighted mean of Z -> want this = 0
        if np.max(np.abs(grad)) < tol:
            break
        Zc = Z - grad
        H = (w[:, None] * Zc).T @ Zc       # weighted covariance = dual Hessian (PSD)
        lmb = lmb + np.linalg.solve(H, grad)   # Newton step
    a = -Z @ lmb; a -= a.max(); ew = q * np.exp(a); w = ew / ew.sum()
    return w, float(np.max(np.abs(Z.T @ w))), it

# split covariates by arm
Xc = df.loc[df.D == 0, ["income", "credit"]].to_numpy()
Xt = df.loc[df.D == 1, ["income", "credit"]].to_numpy()
yc = df.loc[df.D == 0, "Y"].to_numpy()
yt_mean = df.loc[df.D == 1, "Y"].mean()      # ATT uses raw treated mean
Nc = Xc.shape[0]

target_means = Xt.mean(axis=0)               # treated covariate means = balance targets
w_eb, imbalance, iters = entropy_balance(Xc, target_means)
print(f"entropy balance (means) converged in {iters} Newton steps")
print(f"max |weighted-control-mean  -  treated-mean| = {imbalance:.2e}")

3.3 Verify balance is exact, then form the ATT

Entropy balancing’s signature property: balance is exact, not approximate. The covariate means do not come out “close” — they come out equal, to machine precision, by construction. We assert it. Then the entropy-balanced ATT is the (unweighted, because ATT) treated mean minus the reweighted control mean.

reweighted_means = w_eb @ Xc
balance_table = pd.DataFrame({
    "treated_mean":            target_means,
    "control_mean_raw":        Xc.mean(axis=0),
    "control_mean_reweighted": reweighted_means,
    "abs_diff_after":          np.abs(reweighted_means - target_means),
}, index=["income", "credit"])
print(balance_table.to_string(float_format=lambda v: f"{v: .8f}"))

# HARD CHECK: reweighted control means must equal treated means to ~1e-6
assert np.allclose(reweighted_means, target_means, atol=1e-6), "entropy balance failed!"

att_eb = yt_mean - np.sum(w_eb * yc)
print(f"\nentropy-balanced ATT (means matched): {att_eb:+.3f}")
print(f"control weights: max (x uniform) = {w_eb.max()*Nc:5.1f}   "
      f"eff. N = {1.0/np.sum(w_eb**2):5.0f} of {Nc}")

3.4 Match means and variances — and watch the answer barely move

A strength of entropy balancing is that you choose which moments to match. We now add the variances: for each covariate we append a second constraint feature (X_{ik}-\bar X_k^{\text{treated}})^2 whose target is the treated-group variance. The balance table stays mechanically perfect, and — the key reassurance — the ATT barely moves. The IPW wobble (where adding an interaction shifted the estimate by half a point) never happens here, because the estimate no longer rides on a fragile functional form.

target_vars = Xt.var(axis=0)                       # treated variances (population, ddof=0)

# constraint features: [income, credit, (income-mu)^2, (credit-mu)^2]
C_mv = np.hstack([Xc, (Xc - target_means)**2])
tgt_mv = np.concatenate([target_means, target_vars])

w_eb2, imbalance2, iters2 = entropy_balance(C_mv, tgt_mv)

ach_means = w_eb2 @ Xc
ach_vars  = (w_eb2[:, None] * (Xc - target_means)**2).sum(axis=0)
print(f"means+variances converged in {iters2} steps, max imbalance = {imbalance2:.2e}")
assert np.allclose(ach_means, target_means, atol=1e-6), "mean balance failed!"
assert np.allclose(ach_vars,  target_vars,  atol=1e-6), "variance balance failed!"

att_eb2 = yt_mean - np.sum(w_eb2 * yc)
print(f"\nATT, means only      : {att_eb:+.3f}")
print(f"ATT, means+variances : {att_eb2:+.3f}   <- barely moves: robust to the spec")

4. AIPW / doubly robust: two shots at the truth

4.1 The estimator

There is a second, independent route to the counterfactual: a model of the outcome. Fit \hat\mu_1(\mathbf{X}) = \hat{\mathbb E}[Y\mid\mathbf{X},D=1] and \hat\mu_0(\mathbf{X}), impute every unit’s missing potential outcome, and average the difference (regression imputation). AIPW combines that outcome model with the propensity model, per arm:

\hat\mu_d^{\text{DR}} = \frac{1}{N}\sum_i\Big[\hat\mu_d(\mathbf{X}_i) + \frac{\mathbb 1\{D_i=d\}\,(Y_i-\hat\mu_d(\mathbf{X}_i))}{\hat e_d(\mathbf{X}_i)}\Big], \qquad \hat\tau_{\text{DR}} = \hat\mu_1^{\text{DR}} - \hat\mu_0^{\text{DR}},

with \hat e_1 = \hat e and \hat e_0 = 1-\hat e. Read it as the outcome model’s prediction, patched by an inverse-probability-weighted correction to the regression’s residuals. The headline guarantee — doubly robust — is that the estimate is consistent if either the propensity model or the outcome model is correct; you need only one of the two to be right.

def aipw(ps_formula, out_formula, data):
    """Augmented IPW (ATE). Returns the doubly-robust tau-hat."""
    d = data
    Dd = d.D.to_numpy(); Yd = d.Y.to_numpy()
    e = smf.logit(ps_formula, data=d).fit(disp=0).predict(d).to_numpy().clip(0.01, 0.99)
    m1 = smf.ols(out_formula, data=d[d.D == 1]).fit().predict(d).to_numpy()
    m0 = smf.ols(out_formula, data=d[d.D == 0]).fit().predict(d).to_numpy()
    dr1 = np.mean(m1 + Dd * (Yd - m1) / e)
    dr0 = np.mean(m0 + (1 - Dd) * (Yd - m0) / (1.0 - e))
    return dr1 - dr0

aipw_both_ok = aipw("D ~ income + credit", "Y ~ income + credit", df)
print(f"true effect            : {tau_true:+.3f}")
print(f"AIPW (both models good): {aipw_both_ok:+.3f}   <- our doubly-robust estimate")

4.2 Demonstrate double robustness — break one model at a time

Now the surprising part. We misspecify a model by dropping the confounder income from it — a genuinely wrong model, because income drives both selection and the outcome. We run four combinations. For contrast we also run the single-method estimators (plain IPW with the broken propensity model; plain regression imputation with the broken outcome model), which have only one shot and therefore fail when that shot is broken.

The prediction: AIPW survives either single failure (one good model carries it) but finally drifts when both are broken.

def ipw_only(ps_formula, data):
    d = data; Dd = d.D.to_numpy(); Yd = d.Y.to_numpy()
    e = smf.logit(ps_formula, data=d).fit(disp=0).predict(d).to_numpy().clip(0.01, 0.99)
    w = np.where(Dd == 1, d.D.mean() / e, (1 - d.D.mean()) / (1.0 - e))
    h = lambda m: np.sum(w[m] * Yd[m]) / np.sum(w[m])
    return h(Dd == 1) - h(Dd == 0)

def reg_only(out_formula, data):
    d = data
    m1 = smf.ols(out_formula, data=d[d.D == 1]).fit().predict(d).to_numpy()
    m0 = smf.ols(out_formula, data=d[d.D == 0]).fit().predict(d).to_numpy()
    return np.mean(m1 - m0)

PS_GOOD, PS_BAD   = "D ~ income + credit", "D ~ credit"      # BAD drops confounder income
OUT_GOOD, OUT_BAD = "Y ~ income + credit", "Y ~ credit"

rows = [
    ("AIPW",     "good", "good", aipw(PS_GOOD, OUT_GOOD, df)),
    ("AIPW",     "good", "BAD ", aipw(PS_GOOD, OUT_BAD,  df)),
    ("AIPW",     "BAD ", "good", aipw(PS_BAD,  OUT_GOOD, df)),
    ("AIPW",     "BAD ", "BAD ", aipw(PS_BAD,  OUT_BAD,  df)),
    ("IPW only", "BAD ", "  - ", ipw_only(PS_BAD, df)),
    ("reg only", "  - ", "BAD ", reg_only(OUT_BAD, df)),
]
print(f"true effect = {tau_true:+.3f}\n")
print(f"{'method':9s} {'ps':5s} {'outcome':8s} {'estimate':>9s}   survives?")
for meth, ps_s, out_s, est in rows:
    ok = "yes" if abs(est - tau_true) < 0.1 else "NO  <-- fails"
    print(f"{meth:9s} {ps_s:5s} {out_s:8s} {est:+9.3f}   {ok}")

5. Compare all three to the truth

We line up every estimator against the known -0.8. Two honest notes on what is being compared: IPW and AIPW here target the ATE (effect over everyone), while entropy balancing targets the ATT (effect for the enrolled) — they coincide here because the true effect is a constant -0.8, so all valid estimators chase the same number. And on a correctly specified logit DGP like this one, entropy balancing is mathematically equivalent to a particular IPW, so its max control weight is comparable to IPW’s; its guarantee is exact balance in every sample, which IPW only promises.

summary = pd.DataFrame({
    "estimate": [
        tau_true, naive, ht, ipw_stab, att_eb, att_eb2, aipw_both_ok,
    ],
    "estimand": ["truth", "ATE", "ATE", "ATE", "ATT", "ATT", "ATE"],
}, index=[
    "TRUE effect", "naive diff", "IPW (Horvitz-Thompson)", "IPW (Hajek, stabilized)",
    "Entropy balance (means)", "Entropy balance (means+var)", "AIPW (doubly robust)",
])
summary["abs_error"] = (summary["estimate"] - tau_true).abs()
print(summary.to_string(float_format=lambda v: f"{v: .4f}"))

fig, ax = plt.subplots(figsize=(8, 4))
labels = summary.index.tolist()
vals = summary["estimate"].to_numpy()
colors = ["black" if l == "TRUE effect" else
          ("firebrick" if l == "naive diff" else "steelblue") for l in labels]
ax.barh(range(len(vals)), vals, color=colors, alpha=0.85)
ax.axvline(tau_true, color="black", ls="--", lw=1, label="truth = -0.8")
ax.set_yticks(range(len(labels))); ax.set_yticklabels(labels)
ax.invert_yaxis(); ax.set_xlabel("estimated treatment effect (pp)")
ax.set_title("Naive is biased; IPW, entropy balancing, and AIPW all land near -0.8")
ax.legend()
fig.tight_layout(); plt.close(fig)
print("\n[comparison figure rendered to buffer]")

Your Turn — break BOTH models and watch AIPW fail too

You saw AIPW survive a single broken model. The doubly-robust guarantee is “either, not both” — so the natural stress test is to break both and confirm the safety net tears. The aipw(...) and single-method helpers are already written; this is one line plus a comparison.

  1. Break both at once. Call aipw(PS_BAD, OUT_BAD, df) (both formulas drop the confounder income). Compare it to the single-method estimators ipw_only(PS_BAD, df) and reg_only(OUT_BAD, df). You should find all three land near -1.25, far from -0.8 — with both models blind to income, there is no good shot left to carry the estimate, and AIPW is no better than the broken pieces it is built from.

  2. Why -1.25 and not -1.5? It is between the naive -1.5 and the truth -0.8. Both broken models still adjust for credit (which they keep), removing part of the selection bias but not the part that runs through income. Confirm by re-running with OUT_BAD = "Y ~ 1" and PS_BAD = "D ~ 1" (drop all covariates): now AIPW should collapse all the way back to the naive difference. The lesson: AIPW is insurance against getting one model’s form wrong, never against omitting a confounder from both.

  3. Entropy balancing under a missing moment. Entropy balancing is not magic either. Re-run entropy_balance but balance on credit only (pass Xc[:, [1]] and target_means[[1]]). Compute the ATT and the leftover income imbalance. You will see the income means are not balanced and the ATT drifts — entropy balancing guarantees balance only on the moments you chose. State your moments and defend them.

The unifying moral of Ch 3.3: every method here rides on the conditional-independence assumption. None can rescue you from an unobserved confounder — a variable driving both selection and outcome that you never measured. When you cannot defend CIA, you need a different source of variation entirely: an instrument (Ch 3.4).