Refutation Tests and Specification Curves

63.1 Two other ways an estimate can be fragile (30 min)

Days 61–62 asked about one threat: a confounder outside the data. Two more are inside the analysis:

1. A broken pipeline. A bug, a leak of the outcome into the features, an estimator that finds effects anywhere. Refutation tests check for these.
2. Analyst choices. Which features, which learner, how to treat outliers, which estimand. Each choice is defensible; together they can move the answer a lot. A specification curve shows how much.

63.2 Refutation tests (60 min)

Concept. Change the data in ways that should change the answer in a known way, and check the estimator reacts correctly. You met DoWhy's versions in Module 2; here they are by hand with DoubleML, so you can apply them to any estimator.

Create day63/refutations.py (about a minute):

# day63/refutations.py   (python -m day63.refutations)
"""
Refutation tests (the idea behind DoWhy's refuters, Module 2, done by hand with DoubleML): change the data in
ways that SHOULD (or should NOT) change the answer, and check that the estimator reacts correctly.
  placebo treatment:   shuffle who got a coupon. The effect should vanish.
  random common cause: add a column of noise as a 'confounder'. The effect shouldn't move.
  data subsets:        re-estimate on random 80% subsets. The effect should be stable.
  placebo outcome:     an outcome the coupon can't affect (last year's basket size). Effect should be 0.
A refutation test can only show that something is wrong; passing doesn't prove the estimate right.
"""
import doubleml as dml
import numpy as np

from causalkit.coupons import X_COLS, coupons, public
from causalkit.learners import clf, reg

np.random.seed(63)  # DoubleML draws its cross-fitting folds from numpy's global RNG


def plr(df, x_cols=X_COLS, y="spend", d="coupon"):
    data = dml.DoubleMLData(df[x_cols + [d, y]], y, d, x_cols)
    fit = dml.DoubleMLPLR(data, ml_l=reg(), ml_m=clf(), n_folds=5).fit()
    return fit.coef[0], fit.se[0]


if __name__ == "__main__":
    df = public(coupons()).to_pandas()
    rng = np.random.default_rng(0)
    est, se = plr(df)
    print(f"original estimate:                  {est:6.2f} (SE {se:.2f})")

    placebo = [plr(df.assign(coupon=rng.permutation(df.coupon.to_numpy())))[0] for _ in range(5)]
    print(f"placebo treatment (5 shuffles):     {np.mean(placebo):6.2f}  (each: {', '.join(f'{p:.2f}' for p in placebo)})")

    noisy = df.assign(random_cause=rng.normal(size=len(df)))
    print(f"with a random common cause:         {plr(noisy, X_COLS + ['random_cause'])[0]:6.2f}")

    subsets = [plr(df.sample(frac=0.8, random_state=s))[0] for s in range(5)]
    print(f"80% subsets (5):                    {np.mean(subsets):6.2f}  (SD {np.std(subsets):.2f})")

    x_wo = [c for c in X_COLS if c != "basket"]
    p_est, p_se = plr(df, x_wo, y="basket")
    print(f"placebo outcome (basket size):      {p_est:6.2f}  (SE {p_se:.2f})")

python -m day63.refutations

original estimate:                   16.19 (SE 0.56)
placebo treatment (5 shuffles):       0.06  (each: -0.68, 0.16, 0.35, -0.15, 0.60)
with a random common cause:          16.13
80% subsets (5):                     16.18  (SD 0.28)
placebo outcome (basket size):        0.12  (SE 0.22)

All four pass. Shuffled coupons give effects scattered around zero, roughly within one standard error; noise doesn't move the estimate; subsets agree; and the coupon "changes" last year's basket size by $0.12, indistinguishable from zero.

63.3 What refutation tests can and can't do (40 min)

Day 61's world with the hidden browsing signal would pass every one of these tests too: shuffling, noise and subsets can't detect a confounder that isn't in the data. Refutation tests catch mistakes in the analysis, not mistakes in the assumptions. The one exception is the placebo outcome: if the confounder affects the placebo outcome too, the test can catch it. That's why good placebo outcomes (things the treatment can't touch but the confounders would) are worth searching for: last year's spending, behaviour before the campaign was announced, a product category the coupon didn't cover.

Try it. Replace spend with pure noise and re-estimate (DoWhy calls this the "dummy outcome" refuter). Create day63/dummy_outcome.py.

# day63/dummy_outcome.py   (python -m day63.dummy_outcome)
import numpy as np

from causalkit.coupons import X_COLS, coupons, public
from day63.refutations import plr

df = public(coupons()).to_pandas()
df["noise_outcome"] = np.random.default_rng(1).normal(150, 60, len(df))
est, se = plr(df, X_COLS, y="noise_outcome")
print(f"Effect of the coupon on pure noise: {est:.2f} (SE {se:.2f})")

Effect of the coupon on pure noise: -0.04 (SE 0.94)

63.4 The garden of forking paths (40 min)

Concept. Gelman & Loken's phrase: an analysis isn't one path but a garden of them. Each fork is a reasonable choice (include these features? cap outliers? which learner?), and an analyst who tries a few and reports the one that "works" can find almost anything, without any intent to mislead.

The specification curve (Simonsohn, Simmons & Nelson, 2020) is the antidote: list every defensible combination of choices, run them all, and show every estimate, sorted, with a panel underneath marking which choices produced which estimate. Readers see at once whether the conclusion depends on the choices, and which choice matters.

63.5 Day 63 mini-project: a specification curve (75 min)

Four choices, two options each, 16 specifications. Create day63/spec_curve.py (under a minute):

# day63/spec_curve.py   (python -m day63.spec_curve)
"""
Mini-project: a specification curve (Simonsohn, Simmons & Nelson, 2020). An analyst makes many defensible
choices; report ALL combinations, sorted, so readers see how much the answer depends on them.
Choices here (2 x 2 x 2 x 2 = 16 specifications):
  features:  all 24 | without the 10 noise columns
  learners:  lasso (linear) | LightGBM
  model:     PLR (variance-weighted effect) | IRM (ATE)
  outcome:   spend | spend capped at its 99th percentile (common practice with heavy tails)
"""
import itertools

import doubleml as dml
import matplotlib.pyplot as plt
import numpy as np
import polars as pl
from sklearn.linear_model import LassoCV, LogisticRegressionCV
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from causalkit.coupons import BEHAVIOUR, X_COLS, coupons, public, truth
from causalkit.learners import clf, reg

np.random.seed(63)
full = coupons()
df = public(full).to_pandas()
df["spend_capped"] = df.spend.clip(upper=df.spend.quantile(0.99))
choices = {"features": ["all 24", "no noise"], "learners": ["lasso", "LightGBM"], "model": ["PLR", "IRM"],
           "outcome": ["spend", "capped"]}
rows = []
for feats, learn, model, outcome in itertools.product(*choices.values()):
    x = X_COLS if feats == "all 24" else BEHAVIOUR + ["lapsing"]
    y = "spend" if outcome == "spend" else "spend_capped"
    data = dml.DoubleMLData(df[x + ["coupon", y]], y, "coupon", x)
    if learn == "lasso":
        g, m = make_pipeline(StandardScaler(), LassoCV(cv=3)), make_pipeline(StandardScaler(), LogisticRegressionCV(cv=3))
    else:
        g, m = reg(), clf()
    fit = (dml.DoubleMLPLR(data, g, m, n_folds=5) if model == "PLR"
           else dml.DoubleMLIRM(data, g, m, n_folds=5, trimming_threshold=0.02)).fit()
    rows.append({"features": feats, "learners": learn, "model": model, "outcome": outcome,
                 "estimate": fit.coef[0], "lo": fit.confint().iloc[0, 0], "hi": fit.confint().iloc[0, 1]})
spec = pl.DataFrame(rows).sort("estimate")
with pl.Config(tbl_rows=16, tbl_width_chars=110):
    print(spec.with_columns(pl.col("estimate", "lo", "hi").round(2)))
print(f"Range of estimates: {spec['estimate'].min():.2f} to {spec['estimate'].max():.2f}; truth ATE {truth(full)['ATE']:.2f}")
print("Median estimate by choice:")
for c in choices:
    print("   " + c + ": " + ", ".join(f"{k} {v:.2f}" for k, v in
                                     spec.group_by(c).agg(pl.col("estimate").median()).sort(c).iter_rows()))

fig, (top, bottom) = plt.subplots(2, 1, figsize=(9, 6), sharex=True, gridspec_kw={"height_ratios": [2, 1.4]})
xs = np.arange(spec.height)
top.errorbar(xs, spec["estimate"], yerr=[spec["estimate"] - spec["lo"], spec["hi"] - spec["estimate"]],
             fmt="o", color="#2563eb", ms=4)
top.axhline(truth(full)["ATE"], color="#dc2626", ls="--", lw=1, label="true ATE")
top.set_ylabel("coupon effect ($)")
top.legend()
labels = [(c, v) for c, vals in choices.items() for v in vals]
for i, (c, v) in enumerate(labels):
    on = (spec[c] == v).to_numpy()
    bottom.scatter(xs[on], np.full(on.sum(), i), s=14, color="#334155")
bottom.set_yticks(range(len(labels)), [f"{c}: {v}" for c, v in labels], fontsize=8)
bottom.set_xlabel("specification (sorted by estimate)")
fig.tight_layout()
fig.savefig("figures/day63_spec_curve.png", dpi=120)
print("Saved figures/day63_spec_curve.png")

python -m day63.spec_curve

│ features ┆ learners ┆ model ┆ outcome ┆ estimate ┆ lo    ┆ hi    │
│ no noise ┆ lasso    ┆ IRM   ┆ capped  ┆ 11.32    ┆ 9.18  ┆ 13.46 │
│ all 24   ┆ lasso    ┆ IRM   ┆ capped  ┆ 11.51    ┆ 9.44  ┆ 13.57 │
...
│ all 24   ┆ LightGBM ┆ IRM   ┆ capped  ┆ 16.03    ┆ 14.86 ┆ 17.2  │
│ no noise ┆ LightGBM ┆ PLR   ┆ spend   ┆ 16.06    ┆ 14.95 ┆ 17.16 │
Range of estimates: 11.32 to 16.06; truth ATE 16.51
Median estimate by choice:
   features: all 24 14.95, no noise 14.94
   learners: LightGBM 15.94, lasso 13.00
   model: IRM 14.15, PLR 14.97
   outcome: capped 14.75, spend 14.94

Open figures/day63_spec_curve.png. The bottom panel makes it obvious: every low estimate uses lasso, and every LightGBM estimate is between $15.8 and $16.1. The feature set and outcome capping hardly matter. One choice does all the work, and it's one that Day 53's nuisance-RMSE comparison would settle (LightGBM predicts better).

What would a fishing analyst report? With a business case that needs coupons to look weak, $11.32 (lasso, IRM, capped); with one that needs them to look strong, $16.06. Both come with confidence intervals that look perfectly respectable. Only the full curve tells the truth about the analysis.

Write in day63/NOTES.md: which choices would you fix in advance in a pre-analysis plan (Day 21), and why?

63.6 Save your work with Git (20 min)

git add .
git commit -m "Day 63: refutation tests by hand, specification curve"
git push

Coming next (Days 64–66). The question changes from "how sure are we about the effect?" to "what do we do about it?" Uplift modelling, policy learning, and the cost-benefit analysis of a whole programme.

Day 63 checklist

• I can run placebo-treatment, random-cause, subset, placebo-outcome and dummy-outcome refutations
• I know that refutation tests catch analysis mistakes, not missing confounders (except good placebo outcomes)
• I can build and read a specification curve with a choice panel
• I can identify which analyst choice drives an answer, and how to settle it