Reasoning & Evaluation

Double Machine Learning: How to Use Flexible Models for Causal Estimates Without Inheriting Their Bias

Plug a random forest into a causal regression and the confidence interval tightens around the wrong number, because regularisation bias shrinks more slowly than the standard error. Double machine learning fixes this with two devices, a Neyman-orthogonal score and cross-fitting, that make the causal estimate insensitive to first-order learner error. It fixes estimation, not identification, and its guarantees are asymptotic in ways that bite at a few thousand rows.

Regress net financial assets on 401(k) eligibility in the 1991 Survey of Income and Program Participation and you get $19,559, with a standard error of $1,413. Control for income, age, education and six other covariates using random forests inside an orthogonal estimating equation, with 5-fold cross-fitting and 100 random splits, and the estimate falls to $9,247 in the partially linear model and $8,105 in the fully interactive one (Chernozhukov et al., 2018, Econometrics Journal 21(1)). Swap in lasso and it becomes $8,187; boosting, $9,110. Each carries a valid confidence interval.

The reason is not that the learners are accurate. It is that the estimating equation was built so that learner errors cancel to first order, and the data were split so that each unit's residual comes from a model that never saw it. Those two devices, Neyman orthogonality and cross-fitting, are what "double machine learning" (DML) means.

Why this matters: Every analytics team has gradient boosting on the shelf and a stakeholder asking what the campaign actually did. Dropping a flexible model into a regression and reading off a coefficient yields a number whose bias, measured in its own standard errors, grows with the dataset. Knowing which failure DML fixes, and which it leaves untouched, separates a defensible estimate from a precise wrong one.

TL;DR

  • Naive plug-in estimation is inconsistent at the \(\sqrt{n}\) scale. If a learner's error decays like \(n^{-\varphi}\) with \(\varphi < 1/2\) (typical of regularised learners), the causal coefficient's bias is of order \(\sqrt{n}\,n^{-\varphi}\) standard errors and diverges.
  • Orthogonality turns a first-order bias into a second-order one. The remainder becomes a product of nuisance errors, so two learners each converging faster than \(n^{-1/4}\) suffice for \(\sqrt{n}\) inference.
  • In a hand-computable example with 10% learner shrinkage at \(n=1{,}000\), the naive estimate sits 4.5 standard errors from the truth (95% interval coverage: 0.6%); the orthogonal one sits 0.63 away (90% coverage). At a million rows: 14.1 versus 0.2.
  • Cross-fitting removes a separate bias, from a unit's noise leaking into its own in-sample prediction, at \(K\) times the model fits.
  • In the 401(k) data, seven nuisance learners land within about one standard error of each other ($8,187 to $9,247), the practical signature of an orthogonal score working.
  • DML does nothing about identification. It is "very sensitive" to a few bad controls (Hünermund, Louw and Caspi, 2023).

At a Glance

flowchart LR
    W["Data W = (Y, D, X)"] --> S["Split into K folds"]
    S --> L["Learn outcome model l(X)"]
    S --> M["Learn treatment model m(X)"]
    L --> R["Out-of-fold residuals"]
    M --> R
    R --> O["Orthogonal score"]
    O --> T["theta with valid SE"]
    L -.->|"first-order error cancels"| O
    M -.->|"first-order error cancels"| O

    classDef blue fill:#1e40af,stroke:#3b82f6,stroke-width:1px,color:#fff
    classDef purple fill:#6d28d9,stroke:#a78bfa,stroke-width:1px,color:#fff
    classDef teal fill:#0e7490,stroke:#22d3ee,stroke-width:1px,color:#fff
    classDef emerald fill:#047857,stroke:#34d399,stroke-width:1px,color:#fff
    classDef slate fill:#334155,stroke:#64748b,stroke-width:1px,color:#e2e8f0

    class W blue
    class S slate
    class L,M purple
    class R,O teal
    class T emerald

Out-of-fold residuals from two learners feed a moment equation whose derivative with respect to the learners is zero at the truth.

Before Double Machine Learning

The residual-on-residual trick is ninety years old: Frisch and Waugh showed in 1933 that a multiple-regression coefficient equals the coefficient from regressing the outcome residual on that regressor's residual. The open question was what happens when the purging is nonparametric.

Peter Robinson answered it for kernel smoothers: in \(Y = \theta D + g(X) + U\), regressing \(Y - \hat{\mathbb{E}}[Y \mid X]\) on \(D - \hat{\mathbb{E}}[D \mid X]\) gives a \(\sqrt{n}\)-consistent, asymptotically normal \(\hat\theta\) even though the kernel estimates converge more slowly (Robinson, 1988, Econometrica 56(4)). Kernels need smoothness growing with dimension, so the result stayed low-dimensional. The property underneath is older still: Neyman's 1959 \(C(\alpha)\) tests used statistics insensitive to small nuisance errors, hence "Neyman orthogonality".

Biostatistics built the same machinery under another name. The augmented inverse-probability-weighted (AIPW) estimators of Robins, Rotnitzky and Zhao (1994, JASA 89(427)) are consistent if either the outcome or the propensity model is right. Van der Laan and Rubin's (2006) targeted maximum likelihood (TMLE) reached the same efficiency bound by nudging an initial fit until the efficient influence-function equation is solved.

Then high-dimensional data reached economics. Belloni, Chernozhukov and Hansen (2014, Review of Economic Studies 81(2)) showed that selecting controls from the outcome equation alone drops confounders that weakly predict \(Y\) but strongly predict \(D\), producing exactly the regularisation bias this post is about. Selecting on both equations is orthogonalisation in lasso clothing. The 2016 preprint that became the 2018 paper generalised it to any sufficiently fast learner and added cross-fitting to remove the Donsker conditions that had tied semiparametric theory to simple estimators.

timeline
    title From partial regression to double machine learning
    1933 : Frisch and Waugh show partial regression is residual on residual
    1959 : Neyman builds C(alpha) tests insensitive to nuisance error
    1988 : Robinson proves root-n partially linear regression with kernels
    1994 : Robins, Rotnitzky and Zhao introduce augmented IPW estimators
    2006 : van der Laan and Rubin introduce targeted maximum likelihood
    2014 : Belloni, Chernozhukov and Hansen publish post-double-selection lasso
    2018 : Chernozhukov et al. publish DML with cross-fitting
    2022 : DoubleML published in JMLR
         : Automatic debiased ML published in Econometrica

[IMAGE: Lineage chart with an "Econometrics" column (Frisch-Waugh, Robinson, Belloni-Chernozhukov-Hansen, DML) and a "Biostatistics" column (AIPW, TMLE), arrows converging on "efficient influence function". Caption: "Two communities built the same estimator."]

How Double Machine Learning Actually Works

The model and the target

The 2018 paper's lead example is the partially linear regression model:

\[ Y = \theta_0 D + g_0(X) + U, \qquad \mathbb{E}[U \mid X, D] = 0 \]
\[ D = m_0(X) + V, \qquad \mathbb{E}[V \mid X] = 0 \]

\(X\) is a possibly high-dimensional vector of controls; \(g_0\) and \(m_0\) are nuisance parameters, needed only as means to \(\theta_0\). The second equation, how confounders drive treatment, turns out to be the key to removing bias. \(\theta_0\) is causal only if \(X\) contains every common cause of \(D\) and \(Y\) and nothing caused by \(D\), an assumption no estimator below tests.

Why the naive plug-in fails

The obvious approach learns \(g_0\) on an auxiliary sample, then regresses \(Y - \hat g(X)\) on \(D\) in a main sample of size \(n\). Substituting the model and scaling by \(\sqrt n\):

\[ \sqrt{n}(\hat\theta - \theta_0) = \underbrace{\Big(\tfrac{1}{n}\textstyle\sum D_i^2\Big)^{-1}\tfrac{1}{\sqrt n}\textstyle\sum D_i U_i}_{a} + \underbrace{\Big(\tfrac{1}{n}\textstyle\sum D_i^2\Big)^{-1}\tfrac{1}{\sqrt n}\textstyle\sum D_i\big(g_0(X_i) - \hat g(X_i)\big)}_{b} \]

Term \(a\) is an average of mean-zero terms and converges to a normal. Term \(b\) is the problem. Since \(D_i = m_0(X_i) + V_i\) with \(V_i\) mean-zero given \(X_i\), to first order

\[ b \approx \big(\mathbb{E} D^2\big)^{-1} \tfrac{1}{\sqrt n}\textstyle\sum_i m_0(X_i)\big(g_0(X_i) - \hat g(X_i)\big) \]

a sum of \(n\) terms that do not have mean zero, divided by \(\sqrt n\). Every regularised learner trades variance for bias, and in high dimensions its error shrinks like \(n^{-\varphi_g}\) with \(\varphi_g < 1/2\). So \(b\) is of order \(\sqrt{n}\,n^{-\varphi_g}\) and diverges; the paper states it as \(|\sqrt n(\hat\theta - \theta_0)| \to \infty\) (Chernozhukov et al., 2018, section 1).

The learner's error reaches \(\theta\) through \(m_0(X)\). If treatment did not depend on \(X\) and were centred (\(m_0 = 0\)), term \(b\) would be mean-zero; the paper notes the naive score lacks orthogonality "unless \(m_0(X) = 0\)". Confounding is what lets the learner's bias leak into the causal estimate.

[IMAGE: Overlaid histograms of the studentised estimate (n = 500, p = 20), after Figure 1 of Chernozhukov et al. (2018): naive plug-in (amber) displaced from zero, DML (teal) centred under a normal curve. Caption: "Regularisation bias moves the distribution; it does not widen it."]

Orthogonalisation: residual on residual

Partial \(X\) out of both sides instead. With \(\ell_0(X) = \mathbb{E}[Y \mid X] = \theta_0 m_0(X) + g_0(X)\),

\[ Y - \ell_0(X) = \theta_0\big(D - m_0(X)\big) + U \]

a regression through the origin of one residual on another. The Robinson-style partialling-out score is

\[ \psi(W;\theta,\eta) = \Big(Y - \ell(X) - \theta\big(D - m(X)\big)\Big)\big(D - m(X)\big), \qquad \eta = (\ell, m) \]

and solving \(\frac{1}{n}\sum_i \psi(W_i; \hat\theta, \hat\eta) = 0\) gives

\[ \hat\theta = \frac{\sum_i \hat V_i \hat U_i}{\sum_i \hat V_i^2}, \qquad \hat V_i = D_i - \hat m(X_i),\quad \hat U_i = Y_i - \hat\ell(X_i) \]

"Double" refers to the two prediction problems, one for \(\ell\) and one for \(m\).

Neyman orthogonality, derived

Perturb the nuisances and differentiate the expected score at the truth. For the naive score \((Y - D\theta - g(X))D\), moving \(g\) to \(g_0 + r h\) gives

\[ \partial_r\, \mathbb{E}\big[\psi^{\text{naive}}(W;\theta_0, g_0 + r h)\big]\Big|_{r=0} = -\mathbb{E}\big[m_0(X)\, h(X)\big] \]

which is non-zero for almost any direction \(h\) once confounders drive treatment. A small error in \(g\) moves \(\hat\theta\) proportionally.

For the orthogonal score, perturb \(\ell\) to \(\ell_0 + r h_\ell\) and \(m\) to \(m_0 + r h_m\). The two factors become \(U + r\theta_0 h_m - r h_\ell\) and \(V - r h_m\), with expected product

\[ r\big(\theta_0\,\mathbb{E}[h_m V] - \mathbb{E}[h_\ell V] - \mathbb{E}[U h_m]\big) + r^2\big(\mathbb{E}[h_\ell h_m] - \theta_0\,\mathbb{E}[h_m^2]\big) \]

Every first-order term vanishes, because \(h_m\) and \(h_\ell\) are functions of \(X\) while \(\mathbb{E}[V \mid X] = \mathbb{E}[U \mid X] = 0\). A zero derivative at \(r = 0\) is the definition of Neyman orthogonality. What survives is the \(r^2\) term, a cross-product of the two nuisance errors plus the squared treatment-model error, and that is exactly the rate condition the paper imposes on this score (Assumption 4.1(ii)):

\[ \|\hat m - m_0\|_{2} \times \Big(\|\hat m - m_0\|_{2} + \|\hat \ell - \ell_0\|_{2}\Big) \le \delta_N N^{-1/2}, \qquad \delta_N \to 0 \]

Learners converging faster than \(n^{-1/4}\) satisfy it. For the paper's other orthogonal score, \((Y - D\theta - g(X))(D - m(X))\), the condition is the pure product \(\|\hat m - m_0\|_2 \cdot \|\hat g - g_0\|_2\); partialling-out penalises a poor treatment model twice.

Any moment-condition target has an orthogonal version, formed by adding the nuisance's influence-function correction; for linear functionals of a regression that correction uses the Riesz representer, which Chernozhukov, Newey and Singh (2022, Econometrica 90(3)) showed can itself be learned automatically.

Cross-fitting: the second, separate fix

Orthogonality removes regularisation bias, not overfitting bias, which arises when one observation both trains \(\hat\eta\) and is scored by it. The remainder contains terms like

\[ \tfrac{1}{\sqrt n}\textstyle\sum_i V_i\big(\hat g(X_i) - g_0(X_i)\big) \]

If \(\hat g\) never saw unit \(i\), this is conditionally a mean-zero sum with variance equal to the learner's mean squared error, and vanishes by Chebyshev's inequality. If it did, its error at \(X_i\) absorbed part of that unit's noise and nothing centres the term. The paper's contrived learner \(\hat g(X_i) = g_0(X_i) + (Y_i - g_0(X_i))/N^{1/2-\epsilon}\) converges at the nearly parametric rate \(N^{-1/2+\epsilon}\), yet makes the full-sample remainder grow like \(N^{\epsilon}\).

Splitting cures this but wastes data; cross-fitting recovers it. Partition the \(N\) units into \(K\) folds; for each fold, fit \(\hat\ell_k\) and \(\hat m_k\) on the other \(K-1\) and compute residuals for the held-out units; pool all out-of-fold residuals and solve the moment equation once ("DML2") rather than per fold ("DML1"). The authors "found moderate values of K, such as 4 or 5, to work better than K = 2" and "generally recommend DML2", because the pooled Jacobian \(\sum_i \hat V_i^2\) is more stable than \(K\) separate ones.

The variance comes from the same residuals:

\[ \hat\sigma^2 = \Big(\tfrac{1}{N}\textstyle\sum_i \hat V_i^2\Big)^{-2}\, \tfrac{1}{N}\textstyle\sum_i \hat V_i^2 \big(\hat U_i - \hat\theta \hat V_i\big)^2, \qquad \text{SE} = \hat\sigma/\sqrt{N} \]

with no term for uncertainty in \(\hat\ell\) or \(\hat m\); orthogonality licenses leaving it out.

Repeated splits and the median

Because the answer depends on the random split, the paper repeats cross-fitting \(S\) times, reports the median \(\tilde\theta\), and inflates the variance for split-to-split dispersion:

\[ \hat\sigma^2_{\text{median}} = \operatorname{median}_s\Big\{\hat\sigma_s^2 + \big(\hat\theta_s - \tilde\theta\big)^2\Big\} \]

In the 401(k) data, lasso's median within-split standard error of $1,298 becomes $1,558 after adjustment, a 20% increase; for random forests, $1,295 becomes $1,328.

The binary-treatment version is AIPW

For a binary treatment with heterogeneous effects, the "interactive regression model" \(Y = g_0(D, X) + U\) targets the average treatment effect with the orthogonal score

\[ \psi = g(1,X) - g(0,X) + \frac{D\big(Y - g(1,X)\big)}{m(X)} - \frac{(1-D)\big(Y - g(0,X)\big)}{1 - m(X)} - \theta \]

where \(m(X) = P(D = 1 \mid X)\). Term for term this is AIPW: an outcome-model contrast plus a propensity-weighted correction for its residual. Its remainder is the product \(\|\hat m - m_0\|_2 \cdot \|\hat g - g_0\|_2\), so one nuisance may be mediocre if the other is good. DML's addition to the 1994 estimator is the theory: rate conditions for arbitrary learners, and cross-fitting in place of Donsker conditions. Kennedy (2022) reviews DML, AIPW and TMLE as one theory.

Seeing It in Motion

The naive score passes one learner's error straight through; the orthogonal score passes through only a product.

flowchart TB
    subgraph NAIVE["Naive plug-in score"]
        N1["Error in g-hat"] --> N2["Multiplied by D, which tracks m(X)"]
        N2 --> N3["Bias of order sqrt(n) times n^-phi"]
        N3 --> N4["Grows with n, interval misses"]
    end
    subgraph ORTH["Orthogonal score"]
        O1["Error in m-hat"] --> O3["Errors multiply"]
        O2["Error in l-hat"] --> O3
        O3 --> O4["Bias of order sqrt(n) times product"]
        O4 --> O5["Vanishes when product beats n^-1/2"]
    end

    classDef purple fill:#6d28d9,stroke:#a78bfa,stroke-width:1px,color:#fff
    classDef amber fill:#b45309,stroke:#fbbf24,stroke-width:1px,color:#fff
    classDef rose fill:#be123c,stroke:#fb7185,stroke-width:1px,color:#fff
    classDef emerald fill:#047857,stroke:#34d399,stroke-width:1px,color:#fff
    classDef teal fill:#0e7490,stroke:#22d3ee,stroke-width:1px,color:#fff

    class N1,O1,O2 purple
    class N2,N3 amber
    class N4 rose
    class O3,O4 teal
    class O5 emerald

One cross-fitting pass, with repetition over splits summarised in a note:

sequenceDiagram
    participant DS as Dataset
    participant LR as Nuisance learners
    participant FK as Held-out fold k
    participant SC as Orthogonal score
    participant AG as Aggregator

    loop For each fold k from 1 to K
        DS->>LR: Train l-hat and m-hat on other K-1 folds
        LR->>FK: Predict l-hat(X) and m-hat(X) on fold k
        FK->>SC: Residuals U-hat and V-hat for fold k
    end
    SC->>AG: Pooled out-of-fold residuals for all N units
    AG->>AG: Solve the pooled moment equation for theta
    AG->>AG: Sandwich variance from the same residuals
    Note over DS,AG: Repeat for S random partitions, report the median
    AG->>DS: theta and SE inflated for split variation

[IMAGE: 20 squares in 5 fold rows; across five frames one row is teal ("scored") while four are purple ("training"). Caption: "Every unit is scored once, by models that never saw it."]

By the Numbers

401(k) rows are from Table 2 of the paper; the last row is the randomised Pennsylvania bonus experiment (Table 1), where every method should agree, and does.

Setting Nuisance learner Estimate Median SE across splits Split-adjusted SE
401(k), no controls none $19,559 $1,413 n/a
401(k), partially linear, 5-fold Lasso $8,187 $1,298 $1,558
401(k), partially linear, 5-fold Random forest $9,247 $1,295 $1,328
401(k), partially linear, 5-fold Boosting $9,110 $1,314 $1,328
401(k), interactive (AIPW), 5-fold Lasso $7,170 $1,201 $1,398
401(k), interactive (AIPW), 5-fold Random forest $8,105 $1,242 $1,299
401(k), interactive (AIPW), 2-fold Lasso $6,830 $1,282 $1,530
401(k) participation, linear IV, Poterba et al. controls none $13,102 $1,922 n/a
Pennsylvania bonus, log duration, partially linear, 5-fold all seven learners -0.073 to -0.084 0.035 to 0.036 0.035 to 0.037

Sources: Tables 1 and 2 and section 6.2 of Chernozhukov et al. (2018), arXiv:1608.00060. All DML rows use 100 splits and the median method; interactive-model propensities were trimmed at 0.01 and 0.99. 5-fold standard errors are below 2-fold ones "for all methods across all cases", which the authors attribute to more data for the nuisance fits.

[IMAGE: Dot-and-whisker chart of the seven 401(k) partially linear estimates with 95% intervals, beside a reference line at the no-controls 19,559 dollars. Caption: "Different learners, one conclusion."]

A Concrete Example

A linear-Gaussian world, computable by hand, with "machine learning" reduced to the property that matters: shrinkage toward zero.

Setup. \(X \sim N(0,1)\), \(D = X + V\), \(Y = \theta D + 2X + U\) with \(\theta = 1\) and \(U, V \sim N(0,1)\). So \(g_0(x) = 2x\), \(m_0(x) = x\), \(\ell_0(x) = 3x\), and \(Y = 3X + V + U\). Each learner gets the shape right but shrinks its coefficient by a fraction \(s\): \(\hat g(x) = 2(1-s)x\), \(\hat m(x) = (1-s)x\), \(\hat\ell(x) = 3(1-s)x\).

Step 1: moments. \(\mathbb{E}[X^2] = 1\), \(\mathbb{E}[D^2] = 2\), \(\mathbb{E}[DX] = 1\).

Step 2: naive estimate. The learner error is \(2sx\), so

\[ \theta^{\text{naive}} = 1 + \frac{2s \cdot \mathbb{E}[DX]}{\mathbb{E}[D^2]} = 1 + s \]

With \(s = 0.1\) that is 1.10: a 10% bias, linear in the learner's error.

Step 3: orthogonal residuals. \(\hat V = D - (1-s)X = sX + V\) and \(\hat U = Y - 3(1-s)X = 3sX + V + U\).

Step 4: orthogonal estimate. \(\mathbb{E}[\hat V \hat U] = 3s^2 + 1\) and \(\mathbb{E}[\hat V^2] = s^2 + 1\), so

\[ \theta^{\text{orth}} = \frac{1 + 3s^2}{1 + s^2} = 1 + \frac{2s^2}{1 + s^2} \]

With \(s = 0.1\), \(1.03/1.01 = 1.0198\): a bias five times smaller, and quadratic in \(s\). It matches the derivation: with \(h_\ell = -3sx\) and \(h_m = -sx\), the \(r^2\) coefficient is \(3s^2 - s^2 = 2s^2\), divided by \(\mathbb{E}[\hat V^2]\).

Step 5: standard errors. Orthogonal: \(\sqrt{\mathbb{E}[V^2 U^2]}/(\mathbb{E}[V^2]\sqrt n) = 1/\sqrt n\). Naive: \(\sqrt{\mathbb{E}[D^2 U^2]}/(\mathbb{E}[D^2]\sqrt n) = 1/\sqrt{2n}\). The naive estimator is more precise, which makes its bias more dangerous.

Step 6: let the learner improve with \(n\). Take \(s = n^{-1/3}\), a typical nonparametric rate. A nominal 95% interval whose centre is off by \(b\) standard errors covers with probability \(\Phi(1.96 - b) - \Phi(-1.96 - b)\).

\(n\) \(s\) Naive bias Naive bias in SEs Naive coverage Orthogonal bias Orthogonal bias in SEs Orthogonal coverage
1,000 0.100 0.100 4.47 0.6% 0.0198 0.63 90.4%
10,000 0.046 0.046 6.56 under 0.1% 0.0043 0.43 92.9%
100,000 0.022 0.022 9.63 under 0.1% 0.00093 0.29 94.0%
1,000,000 0.010 0.010 14.14 under 0.1% 0.00020 0.20 94.5%

The naive bias shrinks in absolute terms, but slower than its standard error, so in SE units it grows like \(n^{1/6}\) and coverage collapses. The orthogonal bias shrinks like \(n^{-2/3}\), faster than \(n^{-1/2}\), so coverage climbs toward 95%.

Step 7: the boundary case. Slow the learner to \(s = n^{-1/4}\). The orthogonal bias becomes about \(2n^{-1/2}\), between 1.94 and 2.0 standard errors at every sample size, and coverage sticks near 48 to 51%. That is why the product must be \(o(n^{-1/2})\), not merely \(O(n^{-1/2})\).

[IMAGE: Log-x line chart of bias in standard errors against n from 1,000 to 1,000,000: naive (amber) rising 4.5 to 14.1, orthogonal with s = n^-⅓ (teal) falling 0.63 to 0.20, s = n^-¼ (dashed rose) flat near 2.0. Caption: "More data makes the naive estimator more confidently wrong."]

Where It Breaks

Identification is untouched

Omit a confounder and DML estimates the wrong quantity with a correct-looking interval; include a mediator or collider and flexible learners exploit it. Hünermund, Louw and Caspi (2023, Journal of Causal Inference 11(1)) show that "DML is very sensitive to the inclusion of only a few 'bad controls' in the covariate space". The adjustment set must come from a causal argument (see bad controls, colliders and mediators).

For residual confounding, Chernozhukov, Cinelli, Newey, Sharma and Syrgkanis (2021) bound omitted-variable bias for a broad class of causal parameters from "plausibility judgments on the maximum explanatory power of omitted variables".

Overlap problems return through the residuals

If confounders predict treatment almost perfectly, \(\theta\) is identified from a sliver of residual variation in \(\hat V\). The standard error widens, but so does the remainder's relative size, because \(\mathbb{E}[\hat V^2]\) is small while \(\|\hat m - m_0\|^2\) is not. In the interactive model the problem appears as \(1/\hat m(X)\): a propensity of 0.005 gives one unit a weight of 200. Trimming at 0.01 and 0.99, as the 401(k) analysis did, caps weights at 100 but changes the population described. Trimming redefines the estimand.

The partially linear estimand is not the ATE

With heterogeneous effects, \(Y = \tau(X) D + g(X) + U\), the residual regression converges to

\[ \theta^{\text{PLR}} = \frac{\mathbb{E}\big[\operatorname{Var}(D \mid X)\,\tau(X)\big]}{\mathbb{E}\big[\operatorname{Var}(D \mid X)\big]} \]

because \(Y - \ell_0(X) = \tau(X) V + U\). For a binary treatment the weight is \(m(X)(1 - m(X))\): units with propensities near 0.5 dominate, and units almost certain to be treated barely count. In the 401(k) tables the 5-fold partially linear estimates sit $880 to $1,460 above the interactive ones, learner by learner.

Rate conditions are about limits, not your dataset

"Faster than \(n^{-1/4}\)" describes a sequence of sample sizes with an unknown constant. Naimi, Mishler and Kennedy (2023, American Journal of Epidemiology 192(9)) simulated samples of 200, 1,200 and 5,000: under complex confounding, doubly robust estimators "were less biased, but coverage was well below nominal", and only the combination of sample splitting, confounder interactions, a rich learner and a doubly robust estimator "yielded negligible bias and nominal coverage". Zivich and Breskin (2021, Epidemiology 32(3)) favour cross-fit doubly robust estimators but warn they "may require larger sample sizes to avoid finite-sample issues".

Tuning for prediction is not tuning for theta

A learner chosen by prediction error can predict \(Y\) well overall while fitting \(\ell_0\) badly where treatment varies, the only region feeding \(\hat\theta\). The ddml authors recommend stacking several learners, with Monte Carlo support (Ahrens, Hansen, Schaffer and Wiemann, 2023).

Here the literature genuinely disagrees. Commenting on Naimi et al., Balzer and Westling (2023, American Journal of Epidemiology 192(9)) "largely agree" with its recommendations but show that, in both its simulation settings, TMLE without sample splitting improved bias and coverage when the Super Learner library excluded tree-based methods and included regression splines. Whether splitting is needed, they argue, "depends on the specific problem". No finite-sample theorem settles it.

Dependent data

The Chebyshev argument assumes independent units. A customer or store split across folds leaks information, so folds must be formed by cluster; for crossed clusters such as buyers and sellers, Chiang, Kato, Ma and Sasaki (2022, JBES 40(3)) give a multiway cross-fitting scheme. Time series need forward-chained folds.

Alternative Designs

Design How it works Key advantage Key limitation Best when
OLS with hand-picked controls Regress \(Y\) on \(D\) and a few terms of \(X\) Transparent, stable Misspecified form biases \(\theta\) Few confounders, known structure
Post-double-selection lasso Lasso-select controls for \(Y\) and \(D\), OLS on the union Fast, valid under approximate sparsity Linear in the dictionary Many candidate controls, plausibly sparse
IPW with ML propensity Weight outcomes by \(1/\hat m\) One model First-order sensitive to propensity error Simple propensity, intractable outcome
DML Orthogonal score, cross-fitted learners Any fast-enough learner Asymptotic; overlap-sensitive Moderate-to-large \(n\), rich confounders
TMLE Update initial fit along a least-favourable direction Respects outcome bounds; same efficiency Library choice matters Bounded or binary outcomes
R-learner / causal forest Orthogonalised loss for \(\tau(X)\) Estimates heterogeneity Pointwise inference harder The question is "for whom"

The first row is not a strawman: for 401(k) participation, a linear IV model with Poterba et al.'s simple controls gives $13,102, and the paper calls its flexible estimates "broadly consistent with, though somewhat attenuated relative to" it.

The last row is DML's direct descendant. Nie and Wager's R-learner minimises \(\sum_i\big(\hat U_i - \tau(X_i)\hat V_i\big)^2\) over a function class, and with cross-fitted nuisances its error bound matches an oracle that knew \(\ell_0\) and \(m_0\) (Nie and Wager, 2021, Biometrika 108(2)).

How It Is Used in Practice

Software. DoubleML implements the paper's partially linear, interactive and IV models on scikit-learn (Bach et al., 2022, JMLR 23(53)); the Python package stood at 0.11.4 as of September 2026. EconML, now under PyWhy, emphasises heterogeneous effects with LinearDML, CausalForestDML and a panel DynamicDML (Battocchi et al., 2019); v0.17.0 shipped on 31 July 2026.

Compute. Five folds, two learners and 100 splits is 1,000 model fits before tuning; fewer splits is the budget lever, at the cost of a noisier median.

Track record. Baiardi and Naghi (2024, Econometrics Journal 27(2)) re-ran influential economics studies with causal ML. For average effects the gain was credibility; for heterogeneity, discovery: in a teacher-training experiment with an average effect near zero, generic-ML group analysis found a negative effect for the quintile with the lowest predicted benefit and a positive one for the top quintile, both significant at 10%.

A workflow that survives review. Defend the adjustment set with a causal graph; fix trimming before seeing the estimate; fit two structurally different learners plus a stack; use 5 folds, clustered where units repeat, with repeated splits; report every learner; finish with a sensitivity bound. If learners disagree by more than a standard error, that disagreement is the finding.

[IMAGE: A one-page "DML analysis report": adjustment DAG, propensity overlap histogram with trimming lines, forest plot across four learners, split-to-split spread, and an omitted-variable sensitivity contour. Caption: "The estimate is one panel out of five."]

Insights Worth Remembering

  1. Bias in a causal estimate is a property of the estimating equation, not the learner. The same forest is biased in one moment condition and honest in another.

  2. Regularisation bias gets worse, relative to the interval, as data grows. In the worked example the naive estimate drifts from 4.5 to 14.1 standard errors off between a thousand and a million rows.

  3. Confounding is the conduit. The naive bias is weighted by \(m_0(X)\), how much covariates predict treatment. Where they predict nothing, as in a centred randomised experiment, the leak closes.

  4. Orthogonality and cross-fitting fix different biases. One handles systematic learner error, the other the correlation between a unit's noise and its own prediction.

  5. Agreement across learners is the diagnostic, not the goal. When lasso, forests and boosting land within a standard error, the score is working. When they do not, the nuisance rates are probably too slow, whichever learner won cross-validation.

Open Questions

How should nuisance learners be tuned when the target is \(\theta\)? Prediction error is the default, and good predictive fit does not guarantee a good \(\hat\theta\). Open: whether score-targeted criteria reliably do better.

When is cross-fitting worth its cost? Shown: it removes Donsker-type conditions and matters with overfit-prone learners. Contested: Balzer and Westling report good coverage for TMLE without splitting given a well-chosen library.

Can inference stay honest when the rate condition fails? Naimi et al. document well-below-nominal coverage under complex confounding. Higher-order influence-function estimators exist in theory; whether they are practical at analytics scale is open.

What happens when the covariates are learned representations? Adjusting for text or image embeddings is tempting. Whether they carry the confounding information DML needs, or smuggle in post-treatment signal, is largely unstudied; the concern is extrapolated from the bad-controls results.

Sources and Further Reading

  1. Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). "Double/debiased machine learning for treatment and structural parameters." The Econometrics Journal, 21(1), C1-C68. doi:10.1111/ectj.12097; arXiv:1608.00060
  2. Robinson, P. M. (1988). "Root-N-Consistent Semiparametric Regression." Econometrica, 56(4), 931-954. doi:10.2307/1912705
  3. Neyman, J. (1959). "Optimal asymptotic tests of composite statistical hypotheses." In U. Grenander (Ed.), Probability and Statistics, 416-444. Wiley.
  4. Robins, J. M., Rotnitzky, A., & Zhao, L. P. (1994). "Estimation of Regression Coefficients When Some Regressors Are Not Always Observed." JASA, 89(427), 846-866. doi:10.1080/01621459.1994.10476818
  5. van der Laan, M. J., & Rubin, D. (2006). "Targeted Maximum Likelihood Learning." International Journal of Biostatistics, 2(1). doi:10.2202/1557-4679.1043
  6. Belloni, A., Chernozhukov, V., & Hansen, C. (2014). "Inference on Treatment Effects after Selection among High-Dimensional Controls." Review of Economic Studies, 81(2), 608-650. doi:10.1093/restud/rdt044
  7. Chernozhukov, V., Newey, W. K., & Singh, R. (2022). "Automatic Debiased Machine Learning of Causal and Structural Effects." Econometrica, 90(3), 967-1027. doi:10.3982/ECTA18515
  8. Kennedy, E. H. (2022). "Semiparametric doubly robust targeted double machine learning: a review." arXiv:2203.06469
  9. Nie, X., & Wager, S. (2021). "Quasi-oracle estimation of heterogeneous treatment effects." Biometrika, 108(2), 299-319. doi:10.1093/biomet/asaa076
  10. Chiang, H. D., Kato, K., Ma, Y., & Sasaki, Y. (2022). "Multiway Cluster Robust Double/Debiased Machine Learning." Journal of Business & Economic Statistics, 40(3), 1046-1056. doi:10.1080/07350015.2021.1895815
  11. Chernozhukov, V., Cinelli, C., Newey, W., Sharma, A., & Syrgkanis, V. (2021). "Long Story Short: Omitted Variable Bias in Causal Machine Learning." arXiv:2112.13398
  12. Hünermund, P., Louw, B., & Caspi, I. (2023). "Double machine learning and automated confounder selection: A cautionary tale." Journal of Causal Inference, 11(1), 20220078. doi:10.1515/jci-2022-0078
  13. Naimi, A. I., Mishler, A. E., & Kennedy, E. H. (2023). "Challenges in Obtaining Valid Causal Effect Estimates With Machine Learning Algorithms." American Journal of Epidemiology, 192(9), 1536-1544. doi:10.1093/aje/kwab201
  14. Balzer, L. B., & Westling, T. (2023). "Invited Commentary: Demystifying Statistical Inference When Using Machine Learning in Causal Research." American Journal of Epidemiology, 192(9), 1545-1549. doi:10.1093/aje/kwab200
  15. Zivich, P. N., & Breskin, A. (2021). "Machine Learning for Causal Inference: On the Use of Cross-fit Estimators." Epidemiology, 32(3), 393-401. doi:10.1097/EDE.0000000000001332
  16. Ahrens, A., Hansen, C. B., Schaffer, M. E., & Wiemann, T. (2023). "ddml: Double/debiased machine learning in Stata." arXiv:2301.09397
  17. Bach, P., Chernozhukov, V., Kurz, M. S., & Spindler, M. (2022). "DoubleML: An Object-Oriented Implementation of Double Machine Learning in Python." JMLR, 23(53), 1-6. jmlr.org
  18. Baiardi, A., & Naghi, A. A. (2024). "The value added of machine learning to causal inference: evidence from revisited studies." The Econometrics Journal, 27(2), 213-234. doi:10.1093/ectj/utae004; working paper arXiv:2101.00878

Free to read, no ads, no sign-up. If it was useful you can buy me a coffee.