Day 75: Cross-Fitting and the doubleml Package · Applied Decision Intelligence & Data Architecture

Written lesson

Time	Block
0:00 – 0:23	75.1 Sample splitting and cross-fitting
0:23 – 0:48	75.2 DML by hand: a reusable implementation
0:48 – 0:55	Break
0:55 – 1:20	75.3 The `doubleml` package
1:20 – 1:43	75.4 Choosing and diagnosing learners
1:43 – 1:50	Break
1:50 – 2:08	75.5 Repeated cross-fitting and inference
2:08 – 2:30	75.6 Mini-project: coverage of DML confidence intervals

75.1 Sample Splitting and Cross-Fitting (23 min)

To avoid overfitting bias, the nuisance predictions used for observation i must come from models that didn't see i.

Simple sample splitting fits nuisances on half the data and estimates θ on the other half. It works, but wastes half the sample.

Cross-fitting gets the benefit without the waste:

data split into K folds:  [ 1 ][ 2 ][ 3 ][ 4 ][ 5 ]

fold 1: fit ℓ̂, m̂ on folds 2–5 → predict on fold 1 → residuals for fold 1
fold 2: fit ℓ̂, m̂ on folds 1,3–5 → predict on fold 2 → residuals for fold 2
…
then: θ̂ from ALL out-of-fold residuals together

Every observation's residual is out-of-sample, and every observation is used. That's the "cross" in cross-fitting. It's the same idea as K-fold cross-validation, but used to construct residuals rather than to choose hyperparameters.

Forecasting vs causal cross-fitting. In Module 7, training on the future to predict the past was forbidden: it was leakage, because we were simulating real-time forecasts. In DML, we aren't forecasting. We estimate a structural relationship that's assumed stable over the sample, so learning nuisance functions from other periods, including later ones, is legitimate. What still matters with time series is dependence between neighbouring observations, which we handle on Day 76.

75.2 DML by Hand: A Reusable Implementation (25 min)

Writing it yourself makes the method transparent, and gives you full control later (blocked folds, HAC errors). Create causal/dml.py:

"""Double/debiased ML for the partially linear model, written out explicitly."""
from dataclasses import dataclass

import numpy as np
import statsmodels.api as sm
from sklearn.base import clone
from sklearn.model_selection import KFold


@dataclass
class DMLResult:
    theta: float
    se: float
    y_res: np.ndarray
    d_res: np.ndarray
    r2_l: float           # out-of-fold R² of E[Y|X]
    r2_m: float           # out-of-fold R² of E[D|X] (how predictable policy is)

    def ci(self, z=1.96):
        return self.theta - z * self.se, self.theta + z * self.se


def cross_fit_residuals(y, d, X, learner_l, learner_m, splits):
    """Out-of-fold residuals Y − ℓ̂(X) and D − m̂(X) for the given (train, test) splits."""
    y_res = np.full(len(y), np.nan)
    d_res = np.full(len(d), np.nan)
    for train, test in splits:
        y_res[test] = y[test] - clone(learner_l).fit(X[train], y[train]).predict(X[test])
        d_res[test] = d[test] - clone(learner_m).fit(X[train], d[train]).predict(X[test])
    return y_res, d_res


def plr_from_residuals(y_res, d_res, hac_lags: int | None = None):
    """θ̂ = Σ d̃ỹ / Σ d̃²; SE either heteroskedasticity-robust or HAC (Newey–West)."""
    ok = ~np.isnan(y_res) & ~np.isnan(d_res)
    model = sm.OLS(y_res[ok], d_res[ok])
    if hac_lags is None:
        fit = model.fit(cov_type="HC1")
    else:
        fit = model.fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags})
    return float(fit.params[0]), float(fit.bse[0])


def r2(y, residual):
    ok = ~np.isnan(residual)
    return 1 - np.var(residual[ok]) / np.var(y[ok])


def dml_plr(y, d, X, learner_l, learner_m, splits=None, n_folds=5, seed=0,
            hac_lags=None) -> DMLResult:
    y, d, X = np.asarray(y, float), np.asarray(d, float), np.asarray(X, float)
    if splits is None:
        splits = list(KFold(n_folds, shuffle=True, random_state=seed).split(X))
    y_res, d_res = cross_fit_residuals(y, d, X, learner_l, learner_m, splits)
    theta, se = plr_from_residuals(y_res, d_res, hac_lags)
    return DMLResult(theta, se, y_res, d_res, r2(y, y_res), r2(d, d_res))

The final step is the FWL regression of residuals on residuals (74.2), so its standard error is the familiar regression standard error on the residualised data: heteroskedasticity-robust (HC1) by default, or HAC for time series (Day 76). Asymptotically this matches the variance of the orthogonal score.

from sklearn.ensemble import RandomForestRegressor

from causal.dml import dml_plr
from causal.simulate import nonlinear_world

df = nonlinear_world(T=1000, seed=1)
X = df[[f"x{i}" for i in range(5)]]
rf = RandomForestRegressor(200, min_samples_leaf=5, random_state=0, n_jobs=2)
res = dml_plr(df["housing_growth"], df["rate_change"], X, rf, rf)
print(f"θ̂ = {res.theta:.2f} (se {res.se:.2f}), 95% CI {tuple(round(v, 2) for v in res.ci())}")
print(f"out-of-fold R²: E[Y|X] {res.r2_l:.2f}, E[D|X] {res.r2_m:.2f}")

75.3 The `doubleml` Package (25 min)

The DoubleML package (Bach, Chernozhukov, Kurz & Spindler) implements these estimators with a consistent API, many model types (PLR, interactive models, IV variants), repeated cross-fitting, bootstrap inference and sensitivity analysis.

import doubleml as dml
from sklearn.base import clone
from sklearn.ensemble import RandomForestRegressor

from causal.simulate import nonlinear_world

df = nonlinear_world(T=1000, seed=1)
data = dml.DoubleMLData(df, y_col="housing_growth", d_cols="rate_change",
                        x_cols=[f"x{i}" for i in range(5)])

rf = RandomForestRegressor(n_estimators=200, min_samples_leaf=5, random_state=0, n_jobs=2)
plr = dml.DoubleMLPLR(data, ml_l=clone(rf), ml_m=clone(rf),
                      n_folds=5, n_rep=1, score="partialling out")
plr.fit(store_predictions=True)

print(plr.summary)
print(plr.confint(level=0.95))

Argument	Meaning
`ml_l`	learner for ℓ₀(X) = E[Y \| X]
`ml_m`	learner for m₀(X) = E[D \| X] (the reaction function)
`n_folds`	number of cross-fitting folds
`n_rep`	number of repetitions with different random splits (75.5)
`score="partialling out"`	the orthogonal score from 74.5

plr.summary reports θ̂, its standard error, t, p and the 95% interval.

Cross-check with your implementation

from causal.dml import dml_plr

mine = dml_plr(df["housing_growth"], df["rate_change"], df[[f"x{i}" for i in range(5)]], rf, rf)
print(f"by hand:   {mine.theta:.3f} (se {mine.se:.3f})")
print(f"doubleml:  {plr.coef[0]:.3f} (se {plr.se[0]:.3f})")

The numbers won't be identical (different random folds), but they should be close. Use the same splits to make them match almost exactly:

from sklearn.model_selection import KFold

splits = list(KFold(5, shuffle=True, random_state=42).split(df))
plr_same = dml.DoubleMLPLR(data, ml_l=clone(rf), ml_m=clone(rf), score="partialling out")
plr_same.set_sample_splitting(all_smpls=[splits])       # one repetition, our splits
plr_same.fit()
mine_same = dml_plr(df["housing_growth"], df["rate_change"],
                    df[[f"x{i}" for i in range(5)]], rf, rf, splits=splits)
print(plr_same.coef[0], mine_same.theta)

Two implementations that agree is strong evidence that both are correct. The same principle drove the Day 30 and Day 40 regression tests: verify new machinery against an independent version.

75.4 Choosing and Diagnosing Learners (23 min)

DML is only as good as its nuisance estimates. Compare several learners, and report their quality.

from sklearn.linear_model import LassoCV, LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from xgboost import XGBRegressor

LEARNERS = {
    "ols": LinearRegression(),
    "lasso_poly2": make_pipeline(PolynomialFeatures(2, include_bias=False), StandardScaler(),
                                 LassoCV(cv=5, max_iter=5000)),
    "random_forest": RandomForestRegressor(300, min_samples_leaf=5, random_state=0, n_jobs=2),
    "xgboost": XGBRegressor(n_estimators=300, learning_rate=0.05, max_depth=3,
                            subsample=0.8, colsample_bytree=0.8, random_state=0, n_jobs=2),
}

import pandas as pd
rows = []
for name, learner in LEARNERS.items():
    r = dml_plr(df["housing_growth"], df["rate_change"], df[[f"x{i}" for i in range(5)]],
                learner, learner)
    rows.append({"learner": name, "theta": r.theta, "se": r.se,
                 "r2_E[Y|X]": r.r2_l, "r2_E[D|X]": r.r2_m})
print(pd.DataFrame(rows).round(3))

Interpreting the diagnostics

Diagnostic	Tells you	Warning sign
out-of-fold R² of E[Y \| X]	how well controls explain the outcome	very low: noisy θ̂ (wide CI)
out-of-fold R² of E[D \| X]	how predictable policy is from the information set	very high (close to 1): almost no residual variation in D, so θ̂ rests on tiny surprises; the continuous-treatment version of poor overlap (Day 72.5)
θ̂ stable across good learners	conclusions not driven by one model	large disagreement between learners with similar R²: investigate

For our real question, the R² of E[D | X] has a direct economic meaning: how much of Bank of Canada rate setting is explained by its observable reaction function. If it's 0.9, then 90% of the variation in rate changes is "systematic policy", and θ is identified only from the remaining 10%: the timing and size of moves beyond what the economic situation dictated.

Try it 75.4

Learners also have hyperparameters (Day 64). Using GridSearchCV inside the learner (e.g. GridSearchCV(RandomForestRegressor(), {"min_samples_leaf": [2, 5, 20]}, cv=3) as the learner object), does tuning change θ̂ or its CI on the simulation? Why is tuning inside each cross-fitting fold legitimate here?

Answer

Tuning inside the learner happens within each training fold, so the held-out fold remains unseen and there's no overfitting bias. It often improves nuisance R² modestly, with little change in θ̂, which is reassuring. If θ̂ changed substantially, it would signal that nuisance quality was poor before tuning.

75.5 Repeated Cross-Fitting and Inference (18 min)

Different random fold assignments give slightly different θ̂. To make results not depend on one lucky split, repeat cross-fitting R times and aggregate. DoubleML reports the median θ̂, with a standard error that adds the across-split variability:

θ̂ = median(θ̂ᵣ), SE = √ median(SEᵣ² + (θ̂ᵣ − θ̂)²)

plr5 = dml.DoubleMLPLR(data, ml_l=clone(rf), ml_m=clone(rf), n_folds=5, n_rep=5)
plr5.fit()
print(plr5.summary)
print("estimates across repetitions:", plr5.all_coef.round(3))

Add the same aggregation to causal/dml.py:

def aggregate_reps(thetas, ses):
    thetas, ses = np.asarray(thetas), np.asarray(ses)
    theta = float(np.median(thetas))
    se = float(np.sqrt(np.median(ses ** 2 + (thetas - theta) ** 2)))
    return theta, se

Setting	Typical choice
`n_folds`	5 (2–10); more folds = nuisances trained on more data
`n_rep`	5–20 for final results; 1 while developing

75.6 Mini-Project: Coverage of DML Confidence Intervals (23 min)

A 95% confidence interval should contain the true θ in about 95% of repeated samples. That's coverage, the same idea as forecast-interval coverage on Day 63. Measure it:

Add to causal/montecarlo.py:

from causal.dml import dml_plr


def coverage(n_sims=100, T=500, true_effect=-5.0, learner=RF):
    rows = []
    for s in range(n_sims):
        df = nonlinear_world(T=T, true_effect=true_effect, seed=1000 + s)
        X = df[[f"x{i}" for i in range(5)]]
        r = dml_plr(df["housing_growth"], df["rate_change"], X, learner, learner, seed=s)
        lo, hi = r.ci()
        rows.append({"theta": r.theta, "se": r.se, "covered": lo <= true_effect <= hi})
    res = pd.DataFrame(rows)
    return {"mean_theta": res.theta.mean(), "sd_theta": res.theta.std(),
            "mean_se": res.se.mean(), "coverage_95": res.covered.mean()}

Report, for T = 250, 500, 1000:

the mean θ̂ (bias);
the standard deviation of θ̂ across simulations vs the average reported SE (they should be close if the SEs are valid);
the coverage of the 95% CI.

Repeat with LinearRegression as the learner: coverage should collapse, because linear nuisances leave non-linear confounding: the CI is narrow around the wrong value. That's the most important lesson of the day: valid inference requires good nuisance estimates.

Day 75 checklist

I can explain sample splitting vs cross-fitting, and why learning from "other periods" is fine for causal (not forecasting) work
My causal/dml.py implements cross-fitted PLR with robust SEs and out-of-fold R² diagnostics
I can run DoubleMLPLR, read its summary, and reproduce its estimate by hand with identical splits
I compare learners and interpret nuisance R², including the reaction-function R² as an overlap signal
I use repeated cross-fitting and median aggregation
My Monte Carlo shows near-nominal coverage with good learners, and failure with misspecified ones

git add causal notebooks && git commit -m "Day 75: cross-fitting, DML by hand, doubleml, learners, coverage" && git push

Cross-Fitting and the `doubleml` Package

What you'll learn today

Today's plan