Skip to main content

A Guide on Data Analysis

45 Machine Learning for Causal Inference

Causal inference from observational data requires us to adjust for confounding variables that jointly influence the treatment and the outcome. The classical toolkit, built around linear regression, propensity score weighting, and parametric matching, assumes that the analyst can write down a correctly specified model for either the outcome surface or the treatment assignment mechanism. When the number of covariates is large, when their functional relationship to the outcome is unknown, or when treatment effects vary in complex ways across the population, these parametric assumptions become fragile. Machine learning offers flexible, data-adaptive estimators for exactly these high-dimensional nuisance functions. The central intellectual contribution of the modern literature is showing how to combine flexible prediction with the structure of semiparametric theory so that causal estimates retain valid confidence intervals despite being built on top of black-box learners.

This chapter develops that synthesis. We begin with the reason naive plug-in machine learning fails for causal questions, namely the regularization bias problem. We then build up the machinery that solves it: Neyman-orthogonal scores, cross-fitting, and the double or debiased machine learning framework of Chernozhukov, Chetverikov, et al. (2018). We connect this to the older doubly robust tradition of augmented inverse probability weighting and targeted maximum likelihood estimation, and to ensemble learning through the SuperLearner. We then take the orthogonality machinery to platform scale with the debiased deep learning framework for large-scale combinatorial experiments, where the treatment is itself a combination of many concurrent A/B tests and a structured deep network with an influence-function correction recovers the effect of combinations that were never run. We then turn to heterogeneity, where causal forests and generalized random forests recover how effects vary across individuals, and where meta-learners offer a complementary modular approach. Finally we describe a time-aware extension, temporal causal forests, and close with practical guidance and pitfalls.

45.1 Why Machine Learning for Causal Inference

45.1.1 High-Dimensional Confounding

Selection-on-observables identification rests on the assumption of unconfoundedness, that conditional on a covariate vector \(X\) the treatment \(W\) is as good as randomly assigned. The credibility of this assumption grows with the richness of \(X\). A study that conditions on only a handful of covariates invites the criticism that some omitted variable still confounds the comparison. Including hundreds or thousands of covariates, their interactions, and flexible transformations strengthens the identifying assumption, but it overwhelms classical parametric estimators. Ordinary least squares with more covariates than observations is not even defined, and even when \(p < n\) the variance explodes as \(p\) approaches \(n\).

Machine learning methods such as the Lasso, random forests, gradient boosting, and neural networks were designed precisely to predict well in high dimensions by trading a small amount of bias for a large reduction in variance. The temptation is therefore to estimate the outcome regression or the propensity score with one of these learners and to plug the result into a standard causal estimator. As the next subsection shows, doing this naively destroys the statistical guarantees we rely on for inference.

45.1.2 The Regularization Bias Problem

Consider the partially linear model \[ Y = \theta_0 W + g_0(X) + \varepsilon, \qquad \mathbb{E}[\varepsilon \mid X, W] = 0, \] where \(\theta_0\) is the causal parameter of interest and \(g_0\) is an unknown, potentially complicated function of the confounders. A direct strategy is to estimate \(g_0\) with a machine learning method \(\hat{g}\) and then regress the residual \(Y - \hat{g}(X)\) on \(W\). Suppose for concreteness that we estimate \(\theta_0\) by \[ \hat{\theta} = \left(\frac{1}{n}\sum_i W_i^2\right)^{-1} \frac{1}{n}\sum_i W_i \bigl(Y_i - \hat{g}(X_i)\bigr). \] Decomposing the scaled error reveals two terms, \[ \sqrt{n}(\hat{\theta} - \theta_0) = \underbrace{a}_{\text{well behaved}} + \underbrace{b}_{\text{regularization bias}}, \] where the second term involves \(\frac{1}{\sqrt{n}}\sum_i W_i \bigl(g_0(X_i) - \hat{g}(X_i)\bigr)\). Because every machine learning estimator regularizes, shrinking coefficients or limiting tree depth to control variance, the error \(g_0 - \hat{g}\) converges to zero more slowly than the parametric rate \(n^{-1/2}\). The bias term \(b\) therefore diverges, and \(\hat{\theta}\) is not even consistent at the rate needed for valid confidence intervals. This is the regularization bias, sometimes called the plug-in or first-order bias. It is not a small-sample nuisance; it is a structural consequence of using a biased, slowly converging nuisance estimator inside a causal functional that is sensitive to errors in that nuisance.

Two ideas rescue us. The first is to construct the estimating equation so that it is insensitive, to first order, to errors in the nuisance functions. This is Neyman orthogonality. The second is to estimate the nuisance functions on data separate from the data used to form the final estimate, breaking the dependence that would otherwise reintroduce bias. This is cross-fitting. Together they are the heart of double machine learning.

45.2 Double and Debiased Machine Learning

The double or debiased machine learning (DML) framework of Chernozhukov, Chetverikov, et al. (2018), published in The Econometrics Journal, provides a general recipe for estimating a low-dimensional causal parameter in the presence of high-dimensional nuisance functions estimated by arbitrary machine learning methods. It rests on two pillars, Neyman-orthogonal scores and cross-fitting, and it delivers estimators that are root-\(n\) consistent and asymptotically normal even though the first-stage learners converge slowly and are biased by regularization. This section develops the framework in full: the orthogonality condition and a complete derivation of why the canonical scores satisfy it, the two cross-fitting algorithms and their asymptotic theory, the four workhorse models of the paper, and then two replications on the paper’s own empirical applications, the effect of 401(k) eligibility on household financial assets and the effect of institutions on economic growth. Throughout we use the DoubleML package (Bach et al. 2024), the object-oriented R implementation maintained by the paper’s authors and collaborators, with learners drawn from the mlr3 ecosystem.

45.2.1 Neyman-Orthogonal Scores

Suppose the target parameter \(\theta_0\) solves a moment condition \(\mathbb{E}[\psi(Z; \theta_0, \eta_0)] = 0\), where \(Z\) collects the observed data, \(\psi\) is a known score function, and \(\eta_0\) denotes the nuisance functions. The score is Neyman orthogonal if its expectation has a vanishing derivative with respect to the nuisance at the truth, \[ \partial_\eta \, \mathbb{E}\bigl[\psi(Z; \theta_0, \eta)\bigr]\Big|_{\eta = \eta_0} = 0. \] Orthogonality means that small errors in \(\hat{\eta}\) have only a second-order effect on the moment condition. Because machine learning errors are slow but not pathological, second-order terms of the form \(\|\hat{\eta} - \eta_0\|^2\) vanish faster than \(n^{-1/2}\) provided each nuisance is estimated at a rate faster than \(n^{-1/4}\), a rate that flexible learners can attain. The first-order regularization bias is thus eliminated by construction.

For the partially linear model, the naive score \(\psi^{\text{naive}}(Z; \theta, g) = W\,(Y - \theta W - g(X))\) is not orthogonal because it is sensitive to errors in \(g\) through its correlation with \(W\). The orthogonal score partials out the confounders from the treatment as well. Define the two conditional means \[ \ell_0(X) = \mathbb{E}[Y \mid X] = \theta_0\, m_0(X) + g_0(X), \qquad m_0(X) = \mathbb{E}[W \mid X], \] and the Robinson-style partialling-out score \[ \psi(Z; \theta, \eta) = \bigl(Y - \ell(X) - \theta\,(W - m(X))\bigr)\,\bigl(W - m(X)\bigr), \] with nuisance \(\eta = (\ell, m)\). This is the estimator that underlies the classic semiparametric partially linear regression of Robinson (1988), now with machine learning for both nuisance functions. Intuitively, we residualize \(Y\) on \(X\) and \(W\) on \(X\) using flexible learners, then regress the \(Y\)-residual on the \(W\)-residual. The mutual residualization is what makes the score robust to errors in either nuisance.

45.2.1.1 Verifying Orthogonality: A Derivation

The orthogonality claims can be checked directly, and doing so once makes the whole framework transparent. Write \(V = W - m_0(X)\) for the treatment residual and \(\varepsilon = Y - \theta_0 W - g_0(X)\) for the structural error, so that \(\mathbb{E}[\varepsilon \mid X, W] = 0\) and \(\mathbb{E}[V \mid X] = 0\) by construction. Note first that at the true nuisance values the partialling-out residual reduces to the structural error, \[ Y - \ell_0(X) - \theta_0\,(W - m_0(X)) = \theta_0 W + g_0(X) + \varepsilon - \theta_0 m_0(X) - g_0(X) - \theta_0 W + \theta_0 m_0(X) = \varepsilon, \] so the moment condition \(\mathbb{E}[\psi(Z; \theta_0, \eta_0)] = \mathbb{E}[\varepsilon V] = 0\) holds.

To check orthogonality, perturb each nuisance in an arbitrary direction and differentiate. For the outcome regression, replace \(\ell_0\) by \(\ell_0 + r \Delta_\ell\) for a scalar \(r\) and any function \(\Delta_\ell(X)\), and evaluate the derivative of the moment at \(r = 0\): \[ \frac{\partial}{\partial r}\, \mathbb{E}\Bigl[\bigl(Y - \ell_0(X) - r \Delta_\ell(X) - \theta_0 V\bigr) V\Bigr]\Big|_{r=0} = -\,\mathbb{E}\bigl[\Delta_\ell(X)\, V\bigr] = -\,\mathbb{E}\bigl[\Delta_\ell(X)\, \underbrace{\mathbb{E}[V \mid X]}_{=\,0}\bigr] = 0. \] For the propensity, replace \(m_0\) by \(m_0 + r \Delta_m\), which perturbs both appearances of \(m\) in the score: \[ \frac{\partial}{\partial r}\, \mathbb{E}\Bigl[\bigl(\varepsilon + r\,\theta_0 \Delta_m(X)\bigr)\bigl(V - r \Delta_m(X)\bigr)\Bigr]\Big|_{r=0} = \theta_0\, \mathbb{E}\bigl[\Delta_m(X)\, V\bigr] - \mathbb{E}\bigl[\varepsilon\, \Delta_m(X)\bigr] = 0, \] where the first term vanishes because \(\mathbb{E}[V \mid X] = 0\) and the second because \(\mathbb{E}[\varepsilon \mid X] = 0\). Both derivatives are exactly zero, so first-order errors in either learner do not move the moment condition. Contrast the naive score: perturbing \(g\) gives \[ \frac{\partial}{\partial r}\, \mathbb{E}\Bigl[W\bigl(Y - \theta_0 W - g_0(X) - r \Delta_g(X)\bigr)\Bigr]\Big|_{r=0} = -\,\mathbb{E}\bigl[W \Delta_g(X)\bigr] = -\,\mathbb{E}\bigl[m_0(X)\, \Delta_g(X)\bigr], \] which is nonzero whenever the treatment is correlated with the confounders, precisely the situation that makes causal inference hard. The regularization bias of the naive plug-in is this nonzero derivative multiplied by the slowly vanishing learner error, inflated by the \(\sqrt{n}\) scaling. Orthogonalization removes the term entirely rather than hoping the learner error is small.

The same calculation explains the second-order remainder. A Taylor expansion of the moment in the nuisances has no linear term, so the leading remainder is quadratic, of order \(\|\hat{\ell} - \ell_0\|\,\|\hat{m} - m_0\|\). This product structure is what the formal theory exploits: it suffices that the product of the two nuisance errors vanish faster than \(n^{-1/2}\), for instance because each is \(o(n^{-1/4})\), a rate attainable by lasso, random forests, boosting, and neural networks under structured sparsity or smoothness conditions. One learner may even converge more slowly if the other compensates, a trade-off familiar from the doubly robust literature discussed below.

45.2.2 Cross-Fitting and the Two DML Algorithms

Even with an orthogonal score, reusing the same observations to estimate the nuisance functions and to evaluate the moment condition introduces an overfitting bias and requires strong empirical-process (Donsker) conditions on the learner, conditions that highly flexible modern learners violate. Cross-fitting removes this by sample splitting. Partition the observation indices \(\{1, \dots, N\}\) into \(K\) folds \(I_1, \dots, I_K\) of roughly equal size, and for each fold \(k\) let \(\hat{\eta}_k\) denote the nuisance estimates computed on the complement \(I_k^c\), that is, on all data outside fold \(k\). Each observation’s score is then evaluated with nuisance functions that never saw that observation, which breaks the dependence driving the overfitting bias, while rotating through all \(K\) folds ensures every observation contributes to the final estimate, so no efficiency is lost. Cross-fitting allows the use of essentially any learner without complexity restrictions, which is what makes the framework practical.

Chernozhukov, Chetverikov, et al. (2018) define two ways to combine the folds. The DML1 estimator solves the moment condition separately within each fold and averages the solutions, \[ \check{\theta}_k \; \text{solves} \;\; \frac{1}{|I_k|} \sum_{i \in I_k} \psi\bigl(Z_i; \check{\theta}_k, \hat{\eta}_k\bigr) = 0, \qquad \hat{\theta}_{\text{DML1}} = \frac{1}{K} \sum_{k=1}^{K} \check{\theta}_k. \] The DML2 estimator pools the folds into a single estimating equation, \[ \frac{1}{N} \sum_{k=1}^{K} \sum_{i \in I_k} \psi\bigl(Z_i; \hat{\theta}_{\text{DML2}}, \hat{\eta}_k\bigr) = 0. \] The two are first-order equivalent, but DML2 is generally preferred and is the default in software: pooling stabilizes the estimating equation in small samples and avoids folds where the within-fold equation is ill-conditioned. For the linear scores used throughout this chapter, which take the form \(\psi(Z; \theta, \eta) = \psi^a(Z; \eta)\,\theta + \psi^b(Z; \eta)\), DML2 has the closed form \[ \hat{\theta} = -\Biggl(\frac{1}{N} \sum_{k}\sum_{i \in I_k} \psi^a\bigl(Z_i; \hat{\eta}_k\bigr)\Biggr)^{-1} \frac{1}{N} \sum_{k}\sum_{i \in I_k} \psi^b\bigl(Z_i; \hat{\eta}_k\bigr). \] The number of folds \(K\) is not a tuning parameter in the usual sense; theory holds for any fixed \(K \ge 2\), and \(K = 5\) is the practical default, with \(K = 10\) sometimes preferred in moderate samples because each nuisance fit then uses more data. The replication below varies \(K\) as a robustness check.

45.2.3 Asymptotic Normality and Valid Confidence Intervals

The main inferential result of Chernozhukov, Chetverikov, et al. (2018) states that under Neyman orthogonality, the product-rate condition on the nuisance errors described above, and mild regularity conditions, the cross-fitted estimator satisfies \[ \sqrt{N}\,\bigl(\hat{\theta} - \theta_0\bigr) \;\xrightarrow{d}\; \mathcal{N}\bigl(0, \sigma^2\bigr), \qquad \sigma^2 = J_0^{-2}\; \mathbb{E}\bigl[\psi\bigl(Z; \theta_0, \eta_0\bigr)^2\bigr], \qquad J_0 = \mathbb{E}\bigl[\psi^a\bigl(Z; \eta_0\bigr)\bigr], \] and the variance is consistently estimated by the plug-in \[ \hat{\sigma}^2 = \hat{J}^{-2}\, \frac{1}{N} \sum_{k}\sum_{i \in I_k} \psi\bigl(Z_i; \hat{\theta}, \hat{\eta}_k\bigr)^2, \qquad \hat{J} = \frac{1}{N} \sum_{k}\sum_{i \in I_k} \psi^a\bigl(Z_i; \hat{\eta}_k\bigr). \] Ordinary Wald confidence intervals \(\hat{\theta} \pm z_{1-\alpha/2}\, \hat{\sigma} / \sqrt{N}\) are therefore valid despite the black-box first stage. Importantly, the result is uniform over large classes of data-generating processes: it does not rely on the perfect model selection that invalidates naive post-selection inference, a failure documented in the post-double-selection literature (Belloni et al. 2014b) that this framework generalizes.

One practical wrinkle remains. The estimate depends on the random fold partition, and two analysts using different seeds will report slightly different numbers. Chernozhukov, Chetverikov, et al. (2018) recommend repeating the entire cross-fitting procedure \(S\) times over independent partitions and reporting the median estimate \(\hat{\theta} = \operatorname{median}_s\, \hat{\theta}^{(s)}\) together with the split-adjusted variance \[ \hat{\sigma}^2 = \operatorname{median}_s \Bigl\{ \hat{\sigma}^2_{(s)} + \bigl(\hat{\theta}^{(s)} - \hat{\theta}\bigr)^2 \Bigr\}, \] which charges the reported uncertainty for the variability induced by the sample splitting itself. This is the n_rep argument in DoubleML, exercised in the robustness analysis below.

45.2.4 The Four Canonical Models

Chernozhukov, Chetverikov, et al. (2018) work out the framework for four workhorse specifications that between them cover most applied selection-on-observables and instrumental-variables problems. Each is defined by an estimand, a set of nuisance functions, and a Neyman-orthogonal score.

The partially linear regression (PLR) model is the specification developed above. It imposes a constant treatment effect \(\theta_0\) but leaves the confounding function unrestricted, and its orthogonal score is the partialling-out score already derived, with nuisances \(\ell_0(X) = \mathbb{E}[Y \mid X]\) and \(m_0(X) = \mathbb{E}[W \mid X]\).

The partially linear instrumental variable (PLIV) model handles the case where the treatment \(D\) is endogenous even conditional on \(X\), but an instrument \(Z\) is available that is relevant and excluded conditional on \(X\): \[ Y = \theta_0 D + g_0(X) + \zeta, \qquad \mathbb{E}[\zeta \mid X, Z] = 0. \] The orthogonal score residualizes all three observables on the confounders, \[ \psi(Z; \theta, \eta) = \bigl(Y - \ell(X) - \theta\,(D - r(X))\bigr)\,\bigl(Z - m(X)\bigr), \] with nuisances \(\ell_0(X) = \mathbb{E}[Y \mid X]\), \(r_0(X) = \mathbb{E}[D \mid X]\), and \(m_0(X) = \mathbb{E}[Z \mid X]\). This is exactly the logic of two-stage least squares after flexibly partialling out the controls, and it reduces to the classical IV estimator when the nuisances are estimated by linear regression. The orthogonality verification mirrors the PLR derivation, using \(\mathbb{E}[Z - m_0(X) \mid X] = 0\) in place of the treatment residual.

The interactive regression model (IRM) drops the constant-effect restriction for a binary treatment and targets the average treatment effect \[ \theta_0 = \mathbb{E}\bigl[\mu_0(1, X) - \mu_0(0, X)\bigr], \] where \(\mu_0(w, x) = \mathbb{E}[Y \mid W = w, X = x]\). Its orthogonal score is the augmented inverse probability weighting (AIPW) score \[ \psi(Z; \theta, \eta) = \mu(1, X) - \mu(0, X) + \frac{W\,\bigl(Y - \mu(1, X)\bigr)}{e(X)} - \frac{(1 - W)\,\bigl(Y - \mu(0, X)\bigr)}{1 - e(X)} - \theta, \] with the outcome regressions and the propensity score \(e_0(X) = \mathbb{P}(W = 1 \mid X)\) as nuisances. This score is the efficient influence function for the ATE, and the next section develops its doubly robust interpretation in detail. Double machine learning for the interactive model is AIPW with cross-fitted, machine-learned nuisances. In practice the propensity denominators require trimming, and the replication below trims estimated propensities at 0.01.

The interactive instrumental variable model (IIVM) combines heterogeneity and endogeneity for a binary treatment \(D\) and binary instrument \(Z\), targeting the local average treatment effect (LATE) of Imbens and Angrist (1994), the effect among compliers. Its orthogonal score is a ratio of two AIPW-type constructions, one for the outcome and one for the treatment, \[ \psi(Z; \theta, \eta) = \mu(1, X) - \mu(0, X) + \frac{Z\,\bigl(Y - \mu(1, X)\bigr)}{\pi(X)} - \frac{(1 - Z)\,\bigl(Y - \mu(0, X)\bigr)}{1 - \pi(X)} - \theta \Biggl[ r(1, X) - r(0, X) + \frac{Z\,\bigl(D - r(1, X)\bigr)}{\pi(X)} - \frac{(1 - Z)\,\bigl(D - r(0, X)\bigr)}{1 - \pi(X)} \Biggr], \] where \(\mu_0(z, x) = \mathbb{E}[Y \mid Z = z, X = x]\), \(r_0(z, x) = \mathbb{E}[D \mid Z = z, X = x]\), and \(\pi_0(x) = \mathbb{P}(Z = 1 \mid X = x)\). The numerator part is the intent-to-treat effect of the instrument on the outcome and the bracketed part is its effect on takeup, so the estimator is a doubly robust, machine-learned Wald ratio.

Table 45.1: The four canonical double machine learning models, their estimands, and the nuisance functions each requires.
Model Estimand Endogeneity Nuisance functions
PLR Constant effect \(\theta_0\) None given \(X\) \(\mathbb{E}[Y \mid X]\), \(\mathbb{E}[W \mid X]\)
PLIV Constant effect \(\theta_0\) Instrument \(Z\) \(\mathbb{E}[Y \mid X]\), \(\mathbb{E}[D \mid X]\), \(\mathbb{E}[Z \mid X]\)
IRM ATE (or ATT) None given \(X\) \(\mathbb{E}[Y \mid W, X]\), \(\mathbb{P}(W = 1 \mid X)\)
IIVM LATE Instrument \(Z\) \(\mathbb{E}[Y \mid Z, X]\), \(\mathbb{E}[D \mid Z, X]\), \(\mathbb{P}(Z = 1 \mid X)\)

The framework has an important precursor in the post-double-selection estimator of Belloni et al. (2014b). There the nuisances are estimated by lasso regressions of \(Y\) on \(X\) and of \(W\) on \(X\), the union of the selected controls is taken, and \(\theta_0\) is estimated by ordinary least squares of \(Y\) on \(W\) and the selected union. Selecting controls relevant to either equation is precisely what makes the procedure immune to the model-selection mistakes that break single-selection inference, and it is the special case of the PLR score with lasso nuisances. Double machine learning generalizes the insight from lasso to arbitrary learners and from the partially linear model to any orthogonal moment problem. The hdm package implements post-double-selection directly, and the 401(k) replication below uses it as a cross-check on the DML estimates.

45.2.5 A Runnable Partially Linear Example

The following simulation illustrates the orthogonalization and cross-fitting logic from scratch, using the built-in random forest only through base R so that it runs in a clean session. We generate data from a partially linear model with a nonlinear confounding surface, then compare the naive plug-in estimator with the cross-fitted orthogonal estimator.

set.seed(123)

n <- 2000
p <- 10
X <- matrix(runif(n * p), n, p)

# Nonlinear confounding for both treatment and outcome
g0 <- function(X) sin(2 * pi * X[, 1]) + X[, 2]^2 + 0.5 * X[, 3]
m0 <- function(X) plogis(1.5 * X[, 1] - X[, 2])

theta_true <- 1.0
W <- m0(X) + rnorm(n, sd = 0.3)          # continuous treatment
Y <- theta_true * W + g0(X) + rnorm(n, sd = 1)

dat <- data.frame(Y = Y, W = W, X)

# Naive plug-in: estimate g on full data, regress residual on W
naive_fit <- lm(Y ~ poly(X1, 3) + poly(X2, 3) + X3, data = dat)
resid_Y <- dat$Y - predict(naive_fit)
theta_naive <- coef(lm(resid_Y ~ dat$W - 1))

# Cross-fitted orthogonal (partialling-out) estimator
K <- 5
folds <- sample(rep(1:K, length.out = n))
vY <- numeric(n); vW <- numeric(n)
for (k in 1:K) {
    tr <- folds != k
    te <- folds == k
    gY <- lm(Y ~ poly(X1, 3) + poly(X2, 3) + X3, data = dat[tr, ])
    gW <- lm(W ~ poly(X1, 3) + poly(X2, 3) + X3, data = dat[tr, ])
    vY[te] <- dat$Y[te] - predict(gY, newdata = dat[te, ])
    vW[te] <- dat$W[te] - predict(gW, newdata = dat[te, ])
}
theta_dml <- sum(vW * vY) / sum(vW * vW)

c(naive = theta_naive, dml = theta_dml, truth = theta_true)
#> naive.dat$W         dml       truth 
#>   0.2112158   1.0079888   1.0000000

The orthogonal, cross-fitted estimate sits close to the true value of one, while the naive plug-in is biased because it fails to residualize the treatment. With genuine machine learning learners in place of the polynomial regressions, the same logic delivers valid inference in far higher dimensions. The DoubleML package (Bach et al. 2024) implements this machinery for any learner in the mlr3 ecosystem. Rerunning the same problem with random forest nuisances reproduces the from-scratch estimate and adds the score-based standard error from the asymptotic theory above.

library(DoubleML)
library(mlr3)
library(mlr3learners)
lgr::get_logger("mlr3")$set_threshold("warn")

set.seed(123)
dml_data <- DoubleMLData$new(
    dat,
    y_col = "Y",
    d_cols = "W",
    x_cols = paste0("X", 1:p)
)

learner <- lrn("regr.ranger", num.trees = 200, min.node.size = 20)
dml_plr <- DoubleMLPLR$new(
    dml_data,
    ml_l = learner$clone(),
    ml_m = learner$clone(),
    n_folds = 5
)
dml_plr$fit()
dml_plr$summary()
#> Estimates and significance testing of the effect of target variables
#>   Estimate. Std. Error t value Pr(>|t|)    
#> W   0.98981    0.07542   13.12   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The point estimate agrees with the hand-rolled cross-fitting loop, and the reported standard error is the plug-in estimator \(\hat{\sigma}/\sqrt{N}\) from the asymptotic normality result, so the printed confidence interval is the honest Wald interval that the theory licenses.

45.2.6 Replication: 401(k) Eligibility and Household Financial Assets

The flagship empirical application of Chernozhukov, Chetverikov, et al. (2018) asks what effect eligibility for a 401(k) pension plan has on a household’s net financial assets, using data on 9,915 households from the 1991 Survey of Income and Program Participation. The question, whether 401(k) saving crowds out other saving or represents new wealth accumulation, is a classic in public economics (Poterba et al. 1995). The identification argument is that eligibility is determined by employers, so once we condition on the job and income characteristics that drive both employer plan offerings and saving behavior, eligibility is plausibly unconfounded, even though actual participation is a household choice and clearly endogenous. The covariates are age, income, education, family size, marital status, two-earner status, defined-benefit pension coverage, IRA participation, and homeownership. Income enters saving behavior nonlinearly, which is exactly why flexible nuisance estimation matters here: the sample includes households with net financial assets ranging from heavy debt to more than half a million dollars, and misspecifying the income-assets relationship contaminates a parametric estimate.

The same data support all three selection-on-observables and instrumental-variables estimands from the model taxonomy above. The PLR and IRM models estimate the effect of eligibility itself, an intention-to-treat quantity. The IIVM model estimates the LATE of actual participation, instrumented by eligibility. The instrument is relevant because eligibility is a precondition for participation, satisfies monotonicity trivially for the same reason (there are no always-takers, since nobody can participate without being eligible), and inherits its conditional exogeneity from the eligibility argument. We tell DoubleML about the absent always-takers through the subgroups argument, which fixes the takeup regression at zero in the uninstrumented arm rather than asking a classifier to fit a degenerate outcome.

The dataset ships with the hdm package (Chernozhukov et al. 2016), so the analysis is fully reproducible offline. Following the paper, we give the lasso a flexible dictionary of transformed controls, quadratics in the continuous variables plus income interactions, while the tree-based learners receive the raw covariates and are left to find the nonlinearities themselves.

library(DoubleML)
library(mlr3)
library(mlr3learners)
library(data.table)
library(hdm)
lgr::get_logger("mlr3")$set_threshold("warn")

data(pension, package = "hdm")
pen <- data.table(pension)

# Raw controls for the tree-based learners
xs <- c("age", "inc", "educ", "fsize", "marr", "twoearn", "db", "pira", "hown")

# Flexible dictionary for the lasso: quadratics in the continuous controls,
# main effects for the indicators, income interactions with the indicators
flex <- model.matrix(
    ~ 0 + poly(age, 2, raw = TRUE) + poly(inc, 2, raw = TRUE) +
        poly(educ, 2, raw = TRUE) + poly(fsize, 2, raw = TRUE) +
        marr + twoearn + db + pira + hown +
        inc:(marr + twoearn + db + pira + hown),
    data = pen
)
colnames(flex) <- paste0("x", seq_len(ncol(flex)))
pen_flex <- data.table(net_tfa = pen$net_tfa, e401 = pen$e401,
                       p401 = pen$p401, flex)

c(households = nrow(pen), eligible = sum(pen$e401),
  participants = sum(pen$p401))
#>   households     eligible participants 
#>         9915         3682         2594

About 37 percent of households are eligible and 26 percent participate, and participation implies eligibility, the no-always-takers structure noted above. We now run the full grid: three estimands (PLR, IRM ATE, IIVM LATE) by three learners (lasso, random forest, gradient boosting), each with five-fold cross-fitting and propensity trimming at 0.01.

fit_401k <- function(kind) {
    if (kind == "lasso") {
        reg <- lrn("regr.cv_glmnet", s = "lambda.min")
        cls <- lrn("classif.cv_glmnet", s = "lambda.min")
        xcols <- colnames(flex); dat_use <- pen_flex
    } else if (kind == "forest") {
        reg <- lrn("regr.ranger", num.trees = 200, min.node.size = 20)
        cls <- lrn("classif.ranger", num.trees = 200, min.node.size = 20)
        xcols <- xs; dat_use <- pen[, c("net_tfa", "e401", "p401", xs), with = FALSE]
    } else {
        reg <- lrn("regr.xgboost", nrounds = 150, eta = 0.1, max_depth = 4,
                   objective = "reg:squarederror", verbose = 0)
        cls <- lrn("classif.xgboost", nrounds = 150, eta = 0.1, max_depth = 4,
                   objective = "binary:logistic", verbose = 0)
        xcols <- xs; dat_use <- pen[, c("net_tfa", "e401", "p401", xs), with = FALSE]
    }

    d_elig <- DoubleMLData$new(dat_use[, c("net_tfa", "e401", xcols), with = FALSE],
                               y_col = "net_tfa", d_cols = "e401", x_cols = xcols)
    plr <- DoubleMLPLR$new(d_elig, ml_l = reg$clone(), ml_m = reg$clone(),
                           n_folds = 5)
    plr$fit()

    irm <- DoubleMLIRM$new(d_elig, ml_g = reg$clone(), ml_m = cls$clone(),
                           n_folds = 5, trimming_threshold = 0.01)
    irm$fit()

    d_late <- DoubleMLData$new(dat_use, y_col = "net_tfa", d_cols = "p401",
                               z_cols = "e401", x_cols = xcols)
    iiv <- DoubleMLIIVM$new(d_late, ml_g = reg$clone(), ml_m = cls$clone(),
                            ml_r = cls$clone(), n_folds = 5,
                            trimming_threshold = 0.01,
                            subgroups = list(always_takers = FALSE,
                                             never_takers = TRUE))
    iiv$fit()

    data.table(
        learner  = kind,
        model    = c("PLR (eligibility)", "IRM ATE (eligibility)",
                     "IIVM LATE (participation)"),
        estimate = c(plr$coef, irm$coef, iiv$coef),
        se       = c(plr$se, irm$se, iiv$se)
    )
}

set.seed(3141)
res_401k <- rbindlist(lapply(c("lasso", "forest", "boost"), fit_401k))
res_401k[, `:=`(lower = estimate - 1.96 * se, upper = estimate + 1.96 * se)]
Table 45.2: Double machine learning estimates of the effect of 401(k) eligibility and participation on net financial assets in dollars, SIPP 1991, across nuisance learners. Replication of the empirical application in Chernozhukov et al. (2018).
Model Learner Estimate SE Lower Upper
IIVM LATE (participation) boost 11,577 2,140 7,384 15,771
IIVM LATE (participation) forest 10,782 1,764 7,324 14,239
IIVM LATE (participation) lasso 11,724 1,834 8,130 15,319
IRM ATE (eligibility) boost 8,729 1,465 5,859 11,600
IRM ATE (eligibility) forest 8,405 1,189 6,074 10,736
IRM ATE (eligibility) lasso 8,053 1,325 5,456 10,650
PLR (eligibility) boost 9,819 1,553 6,776 12,862
PLR (eligibility) forest 8,833 1,309 6,266 11,399
PLR (eligibility) lasso 10,007 1,414 7,236 12,777
Estimates and 95 percent confidence intervals for the three 401(k) estimands across nuisance learners. The estimates are stable across learners within each estimand, the pattern the orthogonality theory predicts.

Figure 45.1: Estimates and 95 percent confidence intervals for the three 401(k) estimands across nuisance learners. The estimates are stable across learners within each estimand, the pattern the orthogonality theory predicts.

Three findings line up with the published results. First, eligibility raises net financial assets by roughly nine thousand dollars on the partially linear specification and about eight thousand dollars on the interactive (ATE) specification, and the gap between the two is informative: the IRM model reweights toward the covariate distribution of the whole population, so the difference reflects effect heterogeneity that the constant-effect PLR model averages differently. Second, the LATE of actual participation is larger, on the order of eleven to twelve thousand dollars, which is mechanical once one recalls that only about seven in ten eligible households participate: diluting the participation effect over all eligibles yields the intention-to-treat magnitude. Third, and most important for the methodological message, the estimates are stable across the three very different learners within each estimand. Orthogonality plus cross-fitting makes the causal estimate insensitive to the details of the first stage, which is visible here as the flat profile within each panel of the figure.

As a cross-check, the post-double-selection estimator of Belloni et al. (2014b), the lasso-specific precursor to double machine learning, gives a same-ballpark answer on the identical flexible dictionary.

set.seed(3141)
pds <- rlassoEffect(x = flex, y = pen$net_tfa, d = pen$e401,
                    method = "double selection")
summary(pds)
#> [1] "Estimates and significance testing of the effect of target variables"
#>    Estimate. Std. Error t value Pr(>|t|)    
#> d1      6979       1415   4.931 8.18e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The point estimate sits somewhat below the tree-based DML estimates, a reminder that with a fixed dictionary the lasso can only capture the nonlinearities the analyst wrote into it, while forests and boosting adapt their functional form to the data.

45.2.7 Replication: Institutions and Growth with Partially Linear IV

The second empirical application in Chernozhukov, Chetverikov, et al. (2018) revisits the institutions-and-growth question of Acemoglu et al. (2001) with the PLIV model. The outcome is log GDP per capita in 1995 across 64 former colonies, the endogenous treatment is a measure of protection against expropriation risk, and the instrument is the log of early European settler mortality: where settlers faced deadly disease environments, colonial powers built extractive institutions rather than settler institutions, and those institutional differences persisted. The exclusion restriction requires that settler mortality affect modern output only through institutions, conditional on geography. That conditionality is where machine learning earns its place: rather than betting on a particular parametric form for latitude and continent effects, the PLIV model partials all three variables, outcome, treatment, and instrument, on a flexible function of the geographic controls.

data(AJR, package = "hdm")
ajr <- data.table(AJR)[, .(GDP, Exprop, logMort, Latitude, Latitude2,
                           Africa, Asia, Namer, Samer)]

d_ajr <- DoubleMLData$new(
    ajr, y_col = "GDP", d_cols = "Exprop", z_cols = "logMort",
    x_cols = c("Latitude", "Latitude2", "Africa", "Asia", "Namer", "Samer")
)

set.seed(42)
pliv_lasso <- DoubleMLPLIV$new(
    d_ajr,
    ml_l = lrn("regr.cv_glmnet", s = "lambda.min"),
    ml_m = lrn("regr.cv_glmnet", s = "lambda.min"),
    ml_r = lrn("regr.cv_glmnet", s = "lambda.min"),
    n_folds = 5
)
pliv_lasso$fit()

set.seed(42)
pliv_forest <- DoubleMLPLIV$new(
    d_ajr,
    ml_l = lrn("regr.ranger", num.trees = 500, min.node.size = 5),
    ml_m = lrn("regr.ranger", num.trees = 500, min.node.size = 5),
    ml_r = lrn("regr.ranger", num.trees = 500, min.node.size = 5),
    n_folds = 5
)
pliv_forest$fit()

rbind(
    data.frame(learner = "lasso", estimate = pliv_lasso$coef, se = pliv_lasso$se),
    data.frame(learner = "forest", estimate = pliv_forest$coef, se = pliv_forest$se)
)
#>         learner  estimate        se
#> Exprop    lasso 0.8872468 0.2396606
#> Exprop1  forest 0.9944699 0.3906057

Both learners put the effect of institutions on log output near one, close to the two-stage least squares estimates in Acemoglu et al. (2001), with confidence intervals that comfortably exclude zero. Two cautions are in order. With 64 observations the machine learning is doing modest work, and the exercise is best read as a demonstration that the framework nests and reproduces a canonical IV study rather than as an improvement on it. And the small sample makes the estimate more sensitive to the particular fold partition than in the 401(k) application, which is precisely the situation where the repeated cross-fitting of the next subsection matters.

45.2.8 Robustness: Folds, Sample Splits, and Learner Choice

The asymptotic theory holds for any fixed number of folds, but in finite samples the choice can matter, so a serious application reports a sensitivity analysis. We rerun the forest-based PLR estimate of the 401(k) eligibility effect across two, five, and ten folds, and then address the remaining seed dependence with repeated cross-fitting.

d_elig <- DoubleMLData$new(pen[, c("net_tfa", "e401", xs), with = FALSE],
                           y_col = "net_tfa", d_cols = "e401", x_cols = xs)
forest_reg <- lrn("regr.ranger", num.trees = 200, min.node.size = 20)

set.seed(7)
fold_sweep <- rbindlist(lapply(c(2, 5, 10), function(K) {
    m <- DoubleMLPLR$new(d_elig, ml_l = forest_reg$clone(),
                         ml_m = forest_reg$clone(), n_folds = K)
    m$fit()
    data.table(folds = K, estimate = m$coef, se = m$se)
}))
fold_sweep
#>    folds estimate       se
#>    <num>    <num>    <num>
#> 1:     2 8984.912 1307.626
#> 2:     5 9046.653 1336.966
#> 3:    10 9004.840 1329.179

The estimate moves by well under one standard error as the fold count varies, which is the expected behavior in a sample of nearly ten thousand observations. The last source of arbitrariness is the random partition itself. Following the recommendation of Chernozhukov, Chetverikov, et al. (2018), we repeat the entire procedure over independent partitions and aggregate by the median, with the split-adjusted variance from the theory section.

set.seed(7)
m_rep <- DoubleMLPLR$new(d_elig, ml_l = forest_reg$clone(),
                         ml_m = forest_reg$clone(),
                         n_folds = 5, n_rep = 5)
m_rep$fit()
data.frame(
    aggregated = m_rep$coef, se = m_rep$se,
    min_split = min(m_rep$all_coef), max_split = max(m_rep$all_coef)
)
#>      aggregated       se min_split max_split
#> e401   9186.753 1317.844  8810.604  9507.073

The spread across partitions is small relative to the sampling uncertainty here, but it need not be in smaller samples or with more variable learners, and reporting the median-aggregated estimate with the adjusted standard error costs nothing. In published work a larger number of repetitions, fifty or one hundred, is standard.

Taken together, the two replications illustrate the practical workflow this section has built up piece by piece: choose the model class by the estimand and the endogeneity structure, pick several dissimilar learners, give dictionary-based learners a rich feature expansion, cross-fit with five folds and repeat over partitions, trim propensities, and report the full grid so readers can see that the causal conclusion does not hinge on any single machine learning choice.

45.3 Doubly Robust Estimation and Targeted Learning

The double machine learning framework for the interactive model is closely tied to a longer tradition in biostatistics built around doubly robust estimators. Understanding that tradition clarifies why the AIPW score is orthogonal and motivates targeted maximum likelihood estimation.

45.3.1 Augmented Inverse Probability Weighting

For the average treatment effect under unconfoundedness, two simple estimators are available. Outcome regression imputes both potential outcomes from \(\hat{\mu}(w, x)\) and averages the difference. Inverse probability weighting reweights observed outcomes by the inverse propensity score \(\hat{e}(x)\). Each is consistent only if its single nuisance model is correct. The augmented inverse probability weighting estimator combines them into the influence-function form \[ \hat{\tau}_{\text{AIPW}} = \frac{1}{n}\sum_{i=1}^{n}\Biggl[\hat{\mu}(1, X_i) - \hat{\mu}(0, X_i) + \frac{W_i\bigl(Y_i - \hat{\mu}(1, X_i)\bigr)}{\hat{e}(X_i)} - \frac{(1 - W_i)\bigl(Y_i - \hat{\mu}(0, X_i)\bigr)}{1 - \hat{e}(X_i)}\Biggr]. \] This estimator is doubly robust: it is consistent if either the outcome model or the propensity model is correctly specified, not necessarily both. The augmentation term is exactly the correction that makes the moment condition Neyman orthogonal, which is why AIPW with cross-fitting is the interactive-model double machine learning estimator. The doubly robust property has its roots in the missing-data and marginal structural model literature of Robins et al. (2000) and the augmentation theory of Scharfstein et al. (1999), where the efficient influence function for the treatment effect was first characterized.

45.3.2 Targeted Maximum Likelihood Estimation

Targeted maximum likelihood estimation (TMLE), developed by Laan and Rubin (2006), is a substitution estimator that achieves the same double robustness and semiparametric efficiency through a different route. Rather than plug machine learning predictions directly into the target parameter, TMLE proceeds in two stages. First it obtains an initial estimate of the outcome regression, typically from a flexible learner. Second it updates that initial fit through a parametric fluctuation submodel whose direction is chosen so that the updated estimate solves the efficient influence function estimating equation. The fluctuation uses a clever covariate built from the estimated propensity score. The final substitution estimator then respects the natural bounds of the parameter, for example remaining a valid probability for binary outcomes, and inherits the efficiency and double robustness of the influence-function approach.

The advantages of TMLE over plain AIPW are largely practical. Because TMLE is a substitution estimator, it tends to be more stable when propensity scores approach zero or one, and it never produces estimates outside the natural parameter space. Conceptually, AIPW and TMLE are two members of the same family of asymptotically efficient, doubly robust estimators. The tmle package implements the method in R.

library(tmle)
library(SuperLearner)

tmle_fit <- tmle(
    Y = Y, A = W, W = X,
    Q.SL.library = c("SL.glm", "SL.ranger", "SL.glmnet"),
    g.SL.library = c("SL.glm", "SL.ranger", "SL.glmnet")
)
tmle_fit$estimates$ATE

45.3.3 SuperLearner Ensembling

A natural worry is how to choose the machine learning method for the nuisance functions. The SuperLearner of Laan et al. (2007) answers this by not choosing at all. It builds a library of candidate learners, ranging from simple parametric models to flexible ensembles, and forms an optimal convex combination of their cross-validated predictions. The weights are chosen to minimize cross-validated risk, and the oracle results show that the resulting ensemble performs asymptotically as well as the best single algorithm in the library, or the best weighted combination. Because the SuperLearner is itself a valid prediction method, it slots directly into AIPW, TMLE, or double machine learning as the nuisance estimator, and it is the recommended default when the analyst has no strong prior about which learner is appropriate.

library(SuperLearner)

sl_lib <- c("SL.mean", "SL.glm", "SL.glmnet", "SL.ranger", "SL.xgboost")
sl_fit <- SuperLearner(
    Y = Y, X = as.data.frame(X),
    family = gaussian(),
    SL.library = sl_lib,
    cvControl = list(V = 10)
)
sl_fit  # prints the learned ensemble weights

45.3.4 Flexible Nuisance Learners

Several learners recur as components of these ensembles and as standalone nuisance estimators. The Lasso of Tibshirani (1996) performs variable selection and shrinkage, and its high-dimensional inference theory underpins the post-double-selection approach of Belloni et al. (2014b) and Belloni et al. (2014a), an early and influential instance of orthogonalized estimation that prefigures double machine learning. Random forests, introduced by Breiman (2001) and based on the recursive partitioning of Breiman (2017), provide flexible, low-tuning regression for both the outcome and the propensity score. Gradient boosting machines build additive ensembles of shallow trees and are common for propensity score estimation. Bayesian additive regression trees offer a flexible, regularized sum-of-trees prior with automatic uncertainty quantification and are widely used for the outcome surface in causal applications, including the endogeneity-aware formulation of Hill et al. (2021). The unifying message is that any of these can serve as a nuisance learner, and the orthogonal-score plus cross-fitting machinery makes the downstream causal estimate robust to the particular choice.

45.4 Deep Learning for Large-Scale Combinatorial Experiments

The double machine learning machinery of this chapter was built for a single treatment. Large online platforms face a harder version of the problem: they run hundreds of randomized experiments concurrently, and any given user is enrolled in many of them at once. Ye et al. (2026), published in Management Science, study this setting, which they call combinatorial experimentation, and develop a debiased deep learning (DeDL) framework that extends the Neyman-orthogonality and cross-fitting ideas of double machine learning to the case where the treatment itself is a high-dimensional combination of binary A/B tests. The paper pairs a complete semiparametric theory (identification, influence functions, asymptotic normality, and best-arm inference) with one of the cleanest large-scale empirical validations of debiased machine learning to date, a field deployment in which the ground-truth effect of every treatment combination is observable. This section develops the framework in full: the structural model and its link functions, the identification argument that turns an exponential design burden into a linear one, the derivation of the influence function by the pathwise-derivative method, the cross-fitting algorithm and its asymptotic theory including inference on the best arm, the platform evidence, and an end-to-end replication on synthetic data built to the authors’ own specification, using a genuine deep network trained with torch, followed by a Monte Carlo study of coverage and best-arm selection under controlled first-stage bias.

45.4.1 The Combinatorial Experimentation Problem

To use user traffic efficiently, a large platform typically randomizes each of its many concurrent experiments independently, an orthogonal traffic assignment design. Orthogonality guarantees a credible causal estimate for each individual experiment, but it means that the joint effect of any particular combination of treatments is essentially never directly observed: with \(m\) concurrent binary experiments there are \(2^m\) possible treatment combinations, and testing all of them exhaustively with a classical full factorial design requires traffic that grows exponentially in \(m\). Practitioners’ default workaround is to assume treatment effects are linearly additive, so that the effect of running two experiments together is simply the sum of their solo effects, and to make expansion decisions for each experiment independently of the others.

Ye et al. (2026) show empirically that this assumption is not innocuous. Collaborating with a large-scale short-video platform that, for one unusual set of three experiments, tested every one of the \(2^3 = 8\) treatment combinations on the same population, they divide the experimented users into 1,254 subgroups by pre-treatment covariates (gender, age, location, and degree of activeness) and document subgroup joint effects that diverge sharply from the linear-addition prediction, with the direction and magnitude of the gap varying across subgroups: some segments show increasing marginal returns to combined treatment and others show decreasing marginal returns, exactly the heterogeneity that a flexible nuisance model can capture but a linear rule cannot.

Formally, let \(T \in \{0,1\}^m\) denote a user’s combination of treatments across \(m\) concurrent experiments, \(X\) her pre-treatment covariates, and \(Y\) her outcome. In potential-outcomes language each user carries \(2^m\) potential outcomes \(Y(t)\), one per combination, which is the same combinatorial explosion discussed for interference in Chapter 30, now arising within a single user rather than across users. The treatment assignment mechanism \(\nu(t \mid x) = \mathbb{P}[T = t \mid X = x]\) is known by design; under the orthogonal assignment used in practice it is a product of independent Bernoulli draws, one per experiment, so every combination receives positive probability but most receive far too little traffic for cell-by-cell estimation. The framework assumes a semiparametric data-generating process \[ \mathbb{E}[Y \mid X = x, T = t] = G(\theta^*(x), t), \] where \(G\) is a known link function and \(\theta^*(\cdot): \mathbb{R}^{d_X} \to \mathbb{R}^{d_\theta}\) is an unknown, potentially high-dimensional nuisance function of the covariates, estimated with a deep neural network. The link function is the interpretable, structural piece of the model; the network supplies the flexibility needed to capture rich covariate heterogeneity in how users respond. Training minimizes the empirical mean squared error, whose population minimizer is \(\theta^*(\cdot)\) because the true regression function solves the least-squares problem by construction.

Two estimands carry all the managerial content. The average treatment effect of combination \(t\) against the all-control condition \(t_0 = \mathbf{0}\) is \[ \mu(t) := \mathbb{E}\bigl[G(\theta^*(X), t)\bigr] - \mathbb{E}\bigl[G(\theta^*(X), t_0)\bigr], \] and the platform wants \(\mu(t)\) for all \(2^m\) combinations, including the never-observed ones. The best arm is \(t^* := \arg\max_{t} \mu(t)\), and the regret of deploying any other combination is \(\tau(t) := \mu(t^*) - \mu(t)\), so identifying the best arm and testing whether its lead over the runners-up is real amounts to inference on \(\tau(\cdot)\). Choosing the combination that maximizes an outcome is a policy learning problem in the sense of Kitagawa and Tetenov (2018) and Athey and Wager (2021), with the policy space being the \(2^m\) treatment combinations themselves.

45.4.2 Where the Framework Sits in the Literature

The statistical problem of learning all \(2^m\) joint effects is old, and the framework sits at the junction of three literatures. The classical answer is the factorial experiment, the workhorse of industrial experimentation. Dasgupta et al. (2015) recast \(2^K\) factorial designs in the potential-outcomes framework, defining factorial effects as contrasts among each unit’s \(2^K\) potential outcomes and deriving randomization-based inference for them, and Li et al. (2020) extend rerandomization to factorial settings so that covariate balance can be enforced jointly across many factors at design time. On a platform running hundreds of concurrent tests, however, the factorial toolkit stops scaling: a full factorial design needs traffic exponential in \(m\), and a fractional factorial design identifies only on the order of \(m\) effects while deliberately aliasing away the rest, an acceptable sacrifice when high-order interactions are ignorable but not when the interactions are the business question. Egami and Imai (2019) confront the same combinatorial explosion in conjoint experiments, where the number of potential causal interactions forces both regularization and a careful choice of estimand. The DeDL move is different in kind: rather than restricting which interactions exist, it restricts how effects combine, through a structural link function, and in exchange recovers the entire \(2^m\) surface from \(m + 2\) observed cells.

The second stream is experimentation practice on online platforms, treated in Chapter 29. The overlapping, orthogonally hashed traffic assignment that makes hundreds of concurrent tests possible was formalized at Google by Tang et al. (2010), and precisely this infrastructure produces the combinatorial blind spot: every user experiences a random combination of many treatments, yet each test is analyzed marginally and joint effects are filled in by linear addition. The platform literature has documented other failures of the one-test-at-a-time lens, including marketplace interference in two-sided platforms (Johari, Li, et al. 2022) and temporal carryover addressed by switchback designs (Bojinov et al. 2023); combinatorial experimentation is the complementary failure, arising even with perfectly isolated users because treatments interact within a user. Because the end goal is to ship the estimated best combination, the selection problem also looms: picking the arm with the largest estimate and then reporting that same estimate overstates the winner, the winner’s curse for which Andrews et al. (2024) develop exact corrected inference. The best-arm theory below is the asymptotic answer to the same problem.

The third stream is the semiparametric and debiased machine learning theory developed earlier in this chapter. The influence-function calculus behind DeDL is the pathwise-derivative method of Newey (1994); the orthogonal score plus cross-fitting recipe is Chernozhukov, Chetverikov, et al. (2018); and the guarantee that a deep network can estimate the nuisance functions at the \(o(n^{-1/4})\) rate the second stage requires comes from the deep learning estimation theory of Farrell et al. (2021), built on approximation results for ReLU networks of the kind established by Schmidt-Hieber (2020). Where deriving the correction term analytically is impractical, Chernozhukov, Newey, et al. (2022) show how to learn it from data with a Lasso-based Riesz representer; DeDL is the opposite design point, a hand-derived correction made possible because both the link function and the assignment mechanism are known. Finally, the empirical stakes are set by Gordon et al. (2023), who show across a large collection of advertising experiments that debiased machine learning applied to observational data can miss the experimental benchmark badly when unconfoundedness fails. The DeDL validation is the mirror image: with experimental data, where unconfoundedness holds by design, the debiased estimator tracks ground truth almost exactly, locating the failure mode of observational applications in the identification assumptions rather than in the estimation machinery.

45.4.4 Identification with m + 2 Observed Combinations

The structural layer has a striking payoff: the entire \(2^m\) response surface is identified from \(m + 2\) observed combinations. The formal requirements are mild. On the data side, observations are i.i.d. and bounded, and each component of \(\theta^*(\cdot)\) is uniformly bounded and \(p\)-smooth (it lies in a Hölder ball), the standard smoothness currency of nonparametric estimation. On the design side, two conditions on the assignment mechanism do all the work. First, writing \(\tilde T = (1, T')'\), the matrix \(\mathbb{E}[\tilde T \tilde T' \mid X = x]\) must be positive definite uniformly in \(x\), which holds as soon as the control condition and each of the \(m\) single-treatment conditions receive positive probability. Second, for the generalized sigmoid form II the ceiling must be pinned down, which requires \(|\theta_{m+1}(x)| > 0\) and one overlapping pair of treatment sets: a triplet \((i, S_1, S_2)\) with \(i \notin S_1 \cup S_2\) and \(S_1 \neq S_2\) such that the four combinations \(t(S_1)\), \(t(S_1 \cup \{i\})\), \(t(S_2)\), \(t(S_2 \cup \{i\})\) all occur with positive probability, plus a non-degeneracy condition ruling out a knife-edge coincidence among their conditional means. Observing the control, the \(m\) singletons, and one pair such as \(t(\{1,2\})\) satisfies everything, taking \(i = 1\), \(S_1 = \{2\}\), \(S_2 = \emptyset\): hence \(m + 2\) cells.

The identification argument is a ratio trick worth seeing. Under form II, \[ \frac{\theta_{m+1}(x)}{G(\theta(x), t)} - 1 = \exp\bigl(-(\theta_0(x) + \theta_1(x)t_1 + \cdots + \theta_m(x)t_m)\bigr), \] so adding treatment \(i\) to any set \(S\) multiplies the left side by the same factor \(e^{-\theta_i(x)}\) regardless of \(S\). Equating that factor across the two sets, \[ \frac{\theta_{m+1}/G(\theta, t(S_1 \cup \{i\})) - 1}{\theta_{m+1}/G(\theta, t(S_1)) - 1} = \frac{\theta_{m+1}/G(\theta, t(S_2 \cup \{i\})) - 1}{\theta_{m+1}/G(\theta, t(S_2)) - 1}, \] gives one equation in the single unknown \(\theta_{m+1}(x)\) whose coefficients are four observable conditional means, and the non-degeneracy condition guarantees a unique solution. With the ceiling known, dividing it out and inverting the logistic function recovers the linear index \(\theta_0(x) + \theta(x)'t\) at every observed cell, and the rank condition on \(\mathbb{E}[\tilde T \tilde T' \mid x]\) identifies \((\theta_0, \dots, \theta_m)\) exactly as in a linear regression of the inverted index on \(\tilde T\). Nothing in this argument uses more than \(m + 2\) cells. For ten concurrent experiments that is 12 combinations instead of \(1{,}024\), about one percent of the full factorial burden, and the traffic requirement grows linearly rather than exponentially in \(m\), which is what makes the framework usable on a real platform.

Convergence comes packaged with identification. Ye et al. (2026) prove that with a structured network of width \(O\bigl(n^{d_X / (2(p + d_X))} \log^2 n\bigr)\) and depth \(O(\log n)\), the empirical-loss minimizer satisfies, with probability approaching one, \[ \bigl\|\hat\theta_k - \theta^*_k\bigr\|_{L^2(X)}^2 \;\le\; C\left(n^{-\frac{p}{p + d_X}} \log^8 n + \frac{\log\log n}{n}\right) \] for each component \(k\). The proof sandwiches the nuisance error between multiples of the excess prediction risk, \(c_1\,\mathbb{E}\|\theta - \theta^*\|_2^2 \le \mathbb{E}[\ell(\theta)] - \mathbb{E}[\ell(\theta^*)] \le c_2\,\mathbb{E}\|\theta - \theta^*\|_2^2\), so the deep learning theory that controls how well the network predicts outcomes (Farrell et al. 2021; Schmidt-Hieber 2020) transfers to how well it recovers the structural parameters. The lower inequality, called sufficient curvature, is exactly where the identification conditions and the assignment design earn their keep. The resulting rate is \(o(n^{-1/4})\), as the second stage requires, whenever the smoothness \(p\) exceeds the covariate dimension \(d_X\) up to logarithmic factors, the familiar semiparametric trade-off between dimension and smoothness.

45.4.5 Deriving the Influence Function

A network trained this way delivers an accurate predictor of \(\mathbb{E}[Y\mid X,T]\), but the plug-in estimator of the average treatment effect, \[ \hat\mu(t) = \frac{1}{n}\sum_{i=1}^n \Bigl(G(\hat\theta(x_i), t) - G(\hat\theta(x_i), t_0)\Bigr), \] inherits the same regularization bias problem that motivated double machine learning earlier in this chapter, now amplified by the fact that \(\theta(\cdot)\) is vector-valued and the link function is nonlinear in \(\theta\). The remedy is the influence function for \(\mu(t)\), which Ye et al. (2026) derive by the pathwise-derivative method of Newey (1994). Define the advantage function \(H(x, \theta(x); t, t_0) := G(\theta(x), t) - G(\theta(x), t_0)\), so that \(\mu(t) = \mathbb{E}[H(X, \theta^*(X); t, t_0)]\), and let \(\ell_\theta = 2\,(G - y)\,G_\theta\) denote the gradient of the squared loss with respect to \(\theta\). Differentiating the loss twice and taking conditional expectations gives the generalized information matrix \[ \Lambda(x) := \mathbb{E}\bigl[\ell_{\theta\theta} \mid X = x\bigr] = 2\,\mathbb{E}\bigl[G_\theta G_\theta' \mid X = x\bigr] + 2\,\mathbb{E}\bigl[G_{\theta\theta}\,(G - Y) \mid X = x\bigr] = 2\,\mathbb{E}\bigl[G_\theta(\theta^*(x), T)\,G_\theta(\theta^*(x), T)' \mid X = x\bigr], \] where the middle term vanishes because \(\mathbb{E}[Y \mid X, T] = G\) at the truth. Because the assignment mechanism \(\nu(\cdot \mid x)\) is known, \(\Lambda(x)\) is a finite weighted sum over the observed support and requires no additional estimation, one of the two places where the experimental setting pays off (the other being that unconfoundedness holds by design).

The derivation runs along a smooth parametric submodel indexed by \(\eta\), with score \(S(z)\) and truth at \(\eta = 0\), under which both the covariate density \(f(x;\eta)\) and the nuisance \(\theta(x; \eta)\) vary, and asks which function \(\psi\) satisfies \(\partial_\eta\, \mu(\eta)\vert_{\eta=0} = \mathbb{E}[\psi(Z) S(Z)]\) for every such submodel; by the Riesz-representation logic of Newey (1994) that \(\psi\), centered at \(\mu(t)\), is the influence function. Differentiating \(\mu(\eta) = \int H(x, \theta(x; \eta); t, t_0)\, f(x; \eta)\, dx\) by the product rule gives a term through the density, which the usual score identities turn into \(\mathbb{E}[H\,S(Z)]\), and a term through the nuisance path, \(\mathbb{E}[H_\theta' \,\partial_\eta \theta(X; \eta)]\). The key step eliminates the unknowable path derivative \(\partial_\eta \theta\): because \(\theta(\cdot; \eta)\) minimizes squared loss at every \(\eta\), the conditional first-order condition \(\mathbb{E}_\eta[\ell_\theta(Z, \theta(X;\eta)) \mid X = x] \equiv 0\) holds identically in \(\eta\), and differentiating this identity in \(\eta\) and solving yields \[ \partial_\eta \theta(x; \eta)\big\vert_{\eta=0} = -\,\Lambda(x)^{-1}\, \mathbb{E}\bigl[\ell_\theta\, S(Y, T \mid X)\ \big\vert\ X = x\bigr]. \] Substituting back and collecting both pieces as a single expectation against the score produces \[ \psi(z,\theta,\Lambda; t, t_0) = \underbrace{H(x,\theta(x); t,t_0)}_{\text{plug-in}} - \underbrace{H_\theta(x,\theta(x);t,t_0)' \Lambda(x)^{-1} \ell_\theta(y,\check t,\theta(x))}_{\text{debiasing term}}, \] where \(\check t\) is the treatment combination the user actually received and \(H_\theta = G_\theta(\theta(x), t) - G_\theta(\theta(x), t_0)\). The debiasing term is a product of three computable gradients, and for the generalized sigmoid form II every piece is available in closed form: writing \(s = 1/(1 + e^{-u})\) with \(u = \theta_0(x) + \theta(x)'t\), the gradient is \(G_{\theta_k} = \theta_{m+1}\, s(1-s)\,\tilde t_k\) for \(k = 0, \dots, m\) (with \(\tilde t = (1, t')'\)) and \(G_{\theta_{m+1}} = s\).

Two properties make \(\psi\) work. First, at the truth the debiasing term has conditional mean zero given \(X\), because it is a linear transformation of \(\ell_\theta\) and the first-order condition sets \(\mathbb{E}[\ell_\theta \mid X] = 0\); so \(\mathbb{E}[\psi] = \mu(t)\), and the correction adds no bias, only the sampling variability that the plug-in estimator wrongly ignores. Second, \(\psi\) is Neyman orthogonal. Perturb the nuisance \(\theta^* \to \theta^* + r\delta\) and differentiate at \(r = 0\): the plug-in piece moves by \(\mathbb{E}[H_\theta'\delta]\), while the correction piece moves by \(\mathbb{E}[H_\theta' \Lambda^{-1}\, \mathbb{E}(\ell_{\theta\theta} \mid X)\, \delta] = \mathbb{E}[H_\theta' \Lambda^{-1} \Lambda\, \delta] = \mathbb{E}[H_\theta'\delta]\), and the two movements cancel exactly. First-order errors in the trained network therefore do not propagate into the estimate, the same cancellation that drove the partially linear model earlier in this chapter, now executed for a vector-valued nuisance inside a nonlinear structural layer. Invertibility of \(\Lambda(x)\) reduces, under the \(m+2\) design, to the same kind of knife-edge exclusion as in the identification argument: the gradient vectors across the observed cells form a matrix of full rank \(m + 2\) except on a measure-zero coincidence. In practice one adds a small ridge, \(\hat\Lambda(x) + 0.0005\, I\), the moral equivalent of propensity-score trimming, which is exactly what the authors’ own implementation does.

45.4.6 Cross-Fitting, Confidence Intervals, and Best-Arm Inference

Combined with cross-fitting, the influence function yields the DeDL algorithm. Split the sample into \(S\) folds. For each fold \(s\), train the structured network on the other folds to obtain \(\hat\theta_s(\cdot)\), compute \(\hat\Lambda_s(\cdot)\) directly from \(\hat\theta_s\) and the known assignment mechanism, and average the influence function over the held-out fold: \[ \hat\mu_{\text{DeDL}}(t) = \frac{1}{S}\sum_{s=1}^S \frac{1}{|S_s|}\sum_{i \in S_s} \psi(z_i, \hat\theta_s(x_i), \hat\Lambda_s(x_i); t, t_0), \qquad \hat\Psi_{\text{DeDL}}(t;\mu) = \frac{1}{S}\sum_{s=1}^S \frac{1}{|S_s|}\sum_{i \in S_s}\bigl(\psi_i - \hat\mu_{\text{DeDL}}(t)\bigr)^2, \] with the Wald confidence interval \(\hat\mu_{\text{DeDL}}(t) \pm z_{1-\alpha/2}\sqrt{\hat\Psi_{\text{DeDL}}(t;\mu)/n}\). Because the treatment distribution is known, two-way splitting suffices; with an unknown assignment mechanism a third split would be needed to estimate \(\Lambda\) separately. Under the conditions above, with the network converging at \(\|\hat\theta - \theta^*\|_{L^2} = o(n^{-1/4})\), Ye et al. (2026) prove that \(\sqrt{n}\,\hat\Psi_{\text{DeDL}}(t;\mu)^{-1/2}\,(\hat\mu_{\text{DeDL}}(t) - \mu(t)) \to_d \mathcal{N}(0,1)\) for every treatment combination \(t\), observed or not, so Wald intervals are valid across the whole lattice of combinations.

Best-arm identification rides on the same machinery with one genuine subtlety. The natural estimator \(\hat t^* := \arg\max_t \hat\mu_{\text{DeDL}}(t)\) and the regret estimator \(\hat\tau_{\text{DeDL}}(t) := \hat\mu_{\text{DeDL}}(\hat t^*) - \hat\mu_{\text{DeDL}}(t)\) are built from an arm that was itself chosen with the data, so fixed-parameter asymptotics do not apply directly. The proof introduces an infeasible companion anchored at the true best arm, \(\hat\zeta(t) := \hat\mu(t^*) - \hat\mu(t)\), which is asymptotically normal by the fixed-arm theory since \(\psi(z, \theta, \Lambda; t^*, t)\) is just the influence function with the control condition replaced by \(t\). If the best arm is unique, consistency of every \(\hat\mu(t)\) plus a union bound over the \(2^m - 1\) competitors gives \(\mathbb{P}[\hat t^* = t^*] \to 1\), so the feasible and infeasible estimators coincide with probability approaching one and, by Slutsky’s theorem, \(\hat\tau_{\text{DeDL}}(t)\) inherits the normal limit. Testing whether the winner’s lead is real then amounts to checking which combinations have \(\hat\tau\) significantly above zero, and the set of arms whose regret interval covers zero is the honest statement of what the experiment ruled in. The caveat is that the selection argument is asymptotic: when the top arms are nearly tied relative to sampling noise, the event \(\hat t^* \neq t^*\) retains nonvanishing probability at any finite \(n\) and the winner’s estimate is biased upward, the winner’s curse. In that regime the conditional and hybrid corrections of Andrews et al. (2024) are the finite-sample complement to the DeDL asymptotics, and reporting the full set of statistically indistinguishable arms rather than a single point winner is the sensible practice, especially as \(m\) grows and near-ties multiply.

45.4.7 Evidence from a Large-Scale Platform

Ye et al. (2026) validate the framework on three concurrent A/B tests run between January 10 and February 1, 2021 on Platform O, a short-video platform serving billions of users. Each test changed the video recommendation algorithm on a different product surface: the discovery page (first coordinate), the live-stream page (second), and the personalized feed (third). Every user was independently assigned to each test’s treatment with probability \(0.6\), an orthogonal design whose product-Bernoulli assignment mechanism is known exactly, and the three tests happened to target the same population of \(4.45\) million users, so all \(2^3=8\) treatment combinations are genuinely observed, the rare architecture that makes ground-truth validation possible. The outcome is each user’s total daily screen time across the three pages, the platform’s primary engagement metric. Because some of the \(26\) pre-treatment covariates (\(16\) discrete and \(10\) continuous, spanning demographics, device type, social-graph size, and pre-period watch behavior) were imperfectly balanced in the raw data, the authors stratify: they discretize the continuous covariates by quartiles, partition users into \(69{,}111\) strata, and sample equal numbers per treatment cell within each stratum, yielding a stratified sample of \(2{,}066{,}606\) users, roughly \(258{,}000\) per combination, that passes randomization checks on demographics and pre-period behavior (all \(p > 0.05\)). For validation they then treat five combinations as observed (the control, the three individual experiments, and the fully combined treatment, the cells a platform would realistically see from independent A/B tests followed by one back-test) and three as artificially unobserved, and ask whether DeDL, fit on the roughly \(1.29\) million observed-cell users with four-fold cross-fitting, recovers the true, independently measured effects of the held-out cells.

Table 45.3: Ground-truth relative average treatment effect on total daily screen time for all eight combinations of three concurrent experiments on Platform O, and which combinations were treated as observed versus artificially held out for validation. The best arm is \(t^*=(1,0,1)\): treat the discovery page and the personalized feed, leave the live page at control. Source: Ye et al. (2026), Table 2.
Combination Relative ATE Significance Observed?
(0,0,0) 0.000% reference Yes
(0,0,1) +1.091% p < 0.01 Yes
(0,1,0) -0.267% not significant Yes
(1,0,0) +0.758% p < 0.05 Yes
(1,1,1) +2.121% p < 0.0001 Yes
(1,1,0) +0.689% not significant No
(1,0,1) +2.299% p < 0.0001 No
(0,1,1) +1.387% p < 0.001 No

The ground truth alone already refutes linear additivity. Adding the individual effects predicts \(1.091 + 0.758 = 1.849\) percent for the combination \((1,0,1)\), understating its true effect of \(2.299\) percent by roughly a fifth: the discovery-page and feed treatments are synergistic. In the other direction, the fully combined deployment \((1,1,1)\), the one configuration platforms do routinely back-test, delivers \(2.121\) percent, less than \((1,0,1)\) alone, so layering the live-page treatment on top of the best pair destroys about \(0.18\) points of engagement even though that treatment is close to neutral on its own. A platform that shipped everything and trusted the back-test would leave roughly eight percent of the attainable engagement gain on the table, and discovering that requires estimating precisely the combinations nobody tested. The linear-addition heuristic does point at the right arm in this particular configuration, because the live-page treatment’s solo estimate is already negative, but it misprices that arm by twenty percent and misjudges the direction or significance of other combinations, and in the synthetic experiments reported below its best-arm accuracy collapses as \(m\) grows.

Table 45.4: Performance of DeDL against linear-addition, linear-regression, pure-DNN, and undebiased structured-DNN benchmarks in recovering the ATE of held-out treatment combinations and of all eight combinations. Source: Ye et al. (2026), Table 3.
Estimator Correct direction (held out) MAPE (held out) Correct direction (all 8) MAPE (all 8)
Linear addition 2/3 30.06% 7/8 12.02%
Linear regression 2/3 4.90% 7/8 17.37%
Pure DNN 2/3 6.86% 6/8 14.76%
Structured DNN (no debias) 2/3 14.71% 6/8 14.03%
DeDL 3/3 1.75% 8/8 3.07%

The comparison in the second table is the paper’s central empirical result: DeDL dominates every benchmark on every metric, correctly signing and identifying the significance of all three held-out combinations while the linear-addition and linear-regression benchmarks each miss one, and cutting the mean absolute percentage error on held-out combinations by roughly two-thirds relative to the next-best method. Comparing the undebiased structured DNN (SDL) against DeDL isolates the value of the influence-function correction specifically: without debiasing, the confidence intervals are too narrow and understate the true sampling variability, which would inflate false discovery rates in a platform that used them to decide which combination to ship, while DeDL’s intervals have the correct coverage. The authors also report a revealing interaction with training quality: the advantage of debiasing over the plug-in SDL estimator only appears once the network is trained to convergence, tracking the same \(o(n^{-1/4})\) rate condition that the theory requires; with an undertrained network, or with a badly misspecified link function, debiasing can inject more noise than it removes, so in practice one should monitor the network’s cross-validation loss, and compare it to an unrestricted DNN that also takes treatment as a direct input, before trusting the debiased estimates.

Table 45.5: Best-arm identification on Platform O: accuracy of each estimator for the regret of every combination relative to the identified best arm. Correct direction counts combinations whose regret is correctly signed with the correct significance. Source: Ye et al. (2026), Table 4.
Estimator Correct direction MAPE
Linear addition 7/8 21.92%
Linear regression 7/8 11.86%
Pure DNN 7/8 12.83%
Structured DNN (no debias) 8/8 17.45%
DeDL 8/8 5.97%

Best-arm identification sharpens the same message. Every estimator points at \((1,0,1)\) as the winner in this three-experiment problem, but they differ sharply in how well they measure each arm’s regret \(\tau(t)\), which is what a manager needs in order to know how much worse any alternative deployment would be. DeDL signs every regret correctly with the right significance and estimates regret magnitudes with a mean absolute percentage error near six percent, half the error of the next-best benchmark; the undebiased structured network also signs everything correctly but misprices regrets by seventeen percent. The regret scale is the platform’s decision currency, and in a closer contest a mispricing of that size is the difference between ranking two adjacent arms correctly and shipping the wrong one.

The fold-level detail behind these averages carries a practical lesson of its own. In the four-fold cross-fitting, individual folds can be far off: for the fully combined arm \((1,1,1)\), three folds estimate between \(2.30\) and \(2.46\) percent while one fold lands at \(1.25\) percent, an absolute percentage error above forty percent, yet the cross-fitted average comes in at \(2.118\) percent against a truth of \(2.121\), an error of about a tenth of a percent. The pattern repeats for the held-out arm \((0,1,1)\), where fold estimates range from \(0.70\) to \(2.13\) percent but the average misses the truth by under three percent. Single-fold estimates from a trained network are noisy in exactly the way the theory predicts, and averaging over folds is not a cosmetic refinement but the mechanism that delivers the accuracy, one more reason never to skip or truncate the cross-fitting step.

45.4.8 Official Replication Package

Ye et al. (2026) have released an official code and data package at github.com/zikunye2/deep_learning_based_causal_inference_for_combinatorial_experiments, which is worth knowing about independently of the replication below. The repository ships two Jupyter notebooks, main.ipynb and figure1.ipynb, that reproduce every table and figure in the paper, including the heterogeneity plot in Figure 1, the ground-truth, benchmark, and fold-level results reproduced above, and the CATE analyses that space does not permit covering here. Because Platform O’s real data are proprietary, the authors instead release a synthetic stand-in, synthetic_data.csv, with 50,000 rows carrying the same schema as the real deployment: 87 pre-treatment covariate columns (device, browsing, geography, follower and fan counts, and activity level, dummy coded exactly as in the deployed model’s 87-node input layer), the three binary treatment indicators is_DP, is_LP, is_FYP, an outcome y, and a pre-assigned kfold column for the four-way cross-fitting behind Table 7. The synthetic outcome does not reproduce the literal proprietary effect sizes in Table 2, so it validates the pipeline rather than the published numbers, but it is a substantially more realistic covariate structure than a hand-rolled simulation, and readers who want to run the authors’ own PyTorch implementation on data with genuine platform-scale dimensionality should start there.

45.4.9 Replicating the Framework End to End

Because the released synthetic data intentionally does not carry the paper’s true effect sizes, we build our own synthetic environment with a known, checkable ground truth to validate the DeDL machinery, generated to the authors’ own specification for their simulation study (their Section 5). This has the advantage of letting us confirm that the estimator, implemented independently in R with a real neural network rather than the authors’ Python code, recovers a ground truth we control exactly, which is a complementary check to running their notebooks: an independent reimplementation in a different language reaching the same qualitative conclusions is stronger evidence that the result is about the method and not about a particular codebase. We generate \(m\) concurrent binary experiments, covariates \(X \sim U(0,1)^{10}\), and nuisance functions \(\theta_j^*(x) = (A_j'x)^3\) for a random coefficient matrix \(A\), combined through the generalized sigmoid form II link. Consistent with the platform setting, only \(m+2\) of the \(2^m\) treatment combinations are ever observed during training, so the estimator never sees most of the combinations whose effects it is asked to recover.

library(torch)

set.seed(2025); torch_manual_seed(2025)

dX <- 10L; m <- 4L

# True nuisance functions: cubic transforms of a random linear index.
A <- matrix(runif((m + 1) * dX, -0.5, 0.5), nrow = m + 1)
theta_m1_true <- runif(1, 10, 20)
theta_star <- function(X) { raw <- X %*% t(A); sign(raw) * abs(raw)^3 }
G_true <- function(X, Tmat) {
  th <- theta_star(X)
  u <- th[, 1] + rowSums(th[, -1, drop = FALSE] * Tmat)
  theta_m1_true / (1 + exp(-u))
}

# Partial-observation support: control, each single treatment, and one pair (m + 2 combinations).
all_combos <- as.matrix(expand.grid(rep(list(0:1), m)))
obs_support <- all_combos[rowSums(all_combos) %in% c(0, 1), , drop = FALSE]
extra <- rep(0, m); extra[1:2] <- 1
obs_support <- rbind(obs_support, extra)
storage.mode(obs_support) <- "double"
n_support <- nrow(obs_support)
t0 <- rep(0, m)

sample_partial <- function(n) {
  X <- matrix(runif(n * dX), n, dX)
  idx <- sample.int(n_support, n, replace = TRUE)
  Tmat <- obs_support[idx, , drop = FALSE]
  y <- G_true(X, Tmat) + runif(n, -0.05, 0.05)
  list(X = X, T = Tmat, y = y)
}

train <- sample_partial(500L * m)

# Structured network: two-layer MLP maps X to theta_0(x),...,theta_m(x);
# theta_{m+1} is a single trainable scale parameter, as in the paper's simulation.
net <- nn_module(
  initialize = function(dX, hidden, out_dim) {
    self$fc1 <- nn_linear(dX, hidden)
    self$fc2 <- nn_linear(hidden, hidden)
    self$fc3 <- nn_linear(hidden, out_dim)
  },
  forward = function(x) x |> self$fc1() |> nnf_relu() |> self$fc2() |> nnf_relu() |> self$fc3()
)
model <- net(dX, hidden = 10L, out_dim = m + 1L)
theta_m1_hat <- torch_tensor(15, requires_grad = TRUE)

Xtr <- torch_tensor(train$X, dtype = torch_float())
Ttr <- torch_tensor(train$T, dtype = torch_float())
ytr <- torch_tensor(train$y, dtype = torch_float())
opt <- optim_adam(c(model$parameters, list(theta_m1_hat)), lr = 0.01)

forward_G <- function(X, Tt) {
  th <- model(X)
  u <- th[, 1] + torch_sum(th[, 2:(m + 1)] * Tt, dim = 2)
  theta_m1_hat / (1 + torch_exp(-u))
}
for (epoch in seq_len(500L)) {
  opt$zero_grad()
  loss <- nnf_mse_loss(forward_G(Xtr, Ttr), ytr)
  loss$backward(); opt$step()
}

theta_hat_fn <- function(X) {
  with_no_grad({ th <- model(torch_tensor(X, dtype = torch_float())) })
  as.matrix(th)
}
theta_m1_val <- as.numeric(theta_m1_hat)

# Closed-form sigmoid gradients (G is a sigmoid, so its theta-gradient has a simple form).
G_hat_row <- function(th_row, t) theta_m1_val / (1 + exp(-(th_row[1] + sum(th_row[-1] * t))))
Gtheta_hat_row <- function(th_row, t) {
  u <- th_row[1] + sum(th_row[-1] * t)
  s <- 1 / (1 + exp(-u)); ds <- s * (1 - s)
  c(theta_m1_val * ds, theta_m1_val * ds * t, s)
}

inf_data <- sample_partial(500L * m)
th_inf <- theta_hat_fn(inf_data$X)
d <- m + 2L

# Per-unit debiasing correction Lambda(x_i)^{-1} ell_theta(y_i, t_i, theta_hat(x_i));
# this does not depend on the target combination t, so it is computed once and reused.
correction <- matrix(0, nrow(inf_data$X), d)
for (i in seq_len(nrow(inf_data$X))) {
  th_i <- th_inf[i, ]
  Lam <- matrix(0, d, d)
  for (k in seq_len(n_support)) {
    g <- Gtheta_hat_row(th_i, obs_support[k, ])
    Lam <- Lam + outer(g, g)
  }
  Lam <- Lam * 2 / n_support + diag(0.0005, d)
  g_ti <- Gtheta_hat_row(th_i, inf_data$T[i, ])
  ell_th <- 2 * g_ti * (G_hat_row(th_i, inf_data$T[i, ]) - inf_data$y[i])
  correction[i, ] <- solve(Lam, ell_th)
}

# Linear-addition and linear-regression benchmarks, fit on the same observed data.
lr_dat <- data.frame(y = train$y, train$X, train$T)
names(lr_dat) <- c("y", paste0("x", seq_len(dX)), paste0("t", seq_len(m)))
lr_fit <- lm(y ~ ., data = lr_dat)

y0 <- mean(train$y[rowSums(train$T) == 0])
la_effects <- sapply(seq_len(m), function(k) {
  sel <- train$T[, k] == 1 & rowSums(train$T) == 1
  mean(train$y[sel]) - y0
})

Xmc <- matrix(runif(10000L * dX), 10000L, dX)
mu_true <- function(t) mean(G_true(Xmc, matrix(t, nrow(Xmc), m, byrow = TRUE))) -
                        mean(G_true(Xmc, matrix(t0, nrow(Xmc), m, byrow = TRUE)))

combo_results <- data.frame()
for (k in seq_len(nrow(all_combos))) {
  t_target <- as.numeric(all_combos[k, ])
  if (all(t_target == t0)) next
  n_inf <- nrow(inf_data$X)

  Hx  <- sapply(seq_len(n_inf), function(i) G_hat_row(th_inf[i, ], t_target) - G_hat_row(th_inf[i, ], t0))
  Hth <- t(sapply(seq_len(n_inf), function(i) Gtheta_hat_row(th_inf[i, ], t_target) - Gtheta_hat_row(th_inf[i, ], t0)))
  psi <- Hx - rowSums(Hth * correction)

  newx  <- data.frame(matrix(runif(2000L * dX), 2000L, dX)); names(newx) <- paste0("x", seq_len(dX))
  newt1 <- as.data.frame(matrix(t_target, 2000L, m, byrow = TRUE)); names(newt1) <- paste0("t", seq_len(m))
  newt0 <- as.data.frame(matrix(t0, 2000L, m, byrow = TRUE)); names(newt0) <- paste0("t", seq_len(m))

  combo_results <- rbind(combo_results, data.frame(
    combo = paste(t_target, collapse = ""),
    truth = mu_true(t_target),
    LA    = sum(la_effects[t_target == 1]),
    LR    = mean(predict(lr_fit, cbind(newx, newt1))) - mean(predict(lr_fit, cbind(newx, newt0))),
    SDL   = mean(Hx),
    DeDL  = mean(psi),
    DeDL_se = sd(psi) / sqrt(n_inf),
    observed = any(apply(obs_support, 1, function(s) all(s == t_target)))
  ))
}
Table 45.6: Estimated average treatment effect of every non-control combination of four synthetic concurrent experiments, against ground truth, for four benchmark estimators trained on only \(m+2=6\) of the \(2^4=16\) possible combinations.
Combination True ATE Linear addition Linear regression Structured DNN (no debias) DeDL DeDL SE Observed?
1000 0.9312 1.0422 1.0106 0.9675 0.9465 0.0239 TRUE
0100 -1.5259 -1.4903 -1.5134 -1.4760 -1.4907 0.0278 TRUE
1100 -0.5892 -0.4481 -0.5028 -0.6289 -0.6305 0.0375 TRUE
0010 1.0352 1.0722 1.0693 1.0464 1.0751 0.0240 TRUE
1010 1.8614 2.1144 2.0799 2.0078 1.9873 0.0322 FALSE
0110 -0.5170 -0.4182 -0.4440 -0.6203 -0.6035 0.0317 FALSE
1110 0.4259 0.6240 0.5665 0.2940 0.3072 0.0409 FALSE
0001 0.2444 0.2821 0.3250 0.2860 0.2650 0.0161 TRUE
1001 1.1581 1.3243 1.3355 1.2596 1.2091 0.0299 FALSE
0101 -1.2933 -1.2082 -1.1884 -1.2268 -1.2595 0.0337 FALSE
1101 -0.3579 -0.1660 -0.1778 -0.3606 -0.3840 0.0431 FALSE
0011 1.2510 1.3543 1.3943 1.3387 1.3436 0.0308 FALSE
1011 2.0532 2.3965 2.4049 2.2935 2.2419 0.0379 FALSE
0111 -0.2839 -0.1360 -0.1191 -0.3432 -0.3471 0.0389 FALSE
1111 0.6450 0.9062 0.8915 0.5809 0.5680 0.0474 FALSE
Table 45.7: Mean absolute error against ground truth across all fifteen non-control treatment combinations.
Estimator MAE
LA LA 0.1474
LR LR 0.1396
SDL SDL 0.0788
DeDL DeDL 0.0677 
library(ggplot2)

combo_results$combo <- factor(combo_results$combo, levels = combo_results$combo[order(combo_results$observed, decreasing = TRUE)])

ggplot(combo_results, aes(x = combo)) +
  geom_point(aes(y = truth, color = "Truth"), size = 2.5, shape = 18) +
  geom_pointrange(aes(y = DeDL, ymin = DeDL - 1.96 * DeDL_se, ymax = DeDL + 1.96 * DeDL_se, color = "DeDL")) +
  labs(x = "Treatment combination", y = "Average treatment effect", color = NULL) +
  causalverse::ama_theme() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
DeDL point estimates and 95 percent confidence intervals against ground truth for every treatment combination in the synthetic replication. Combinations marked as observed were part of the training support; the rest were never seen during training.

Figure 45.3: DeDL point estimates and 95 percent confidence intervals against ground truth for every treatment combination in the synthetic replication. Combinations marked as observed were part of the training support; the rest were never seen during training.

The DeDL confidence intervals cover the true effect for essentially every combination, including those the network never saw a single training example of, which is the payoff of Neyman orthogonality: the debiasing term corrects for the network’s residual approximation error well enough to deliver valid inference on combinations entirely outside the training support.

45.4.10 The Cost of Linear Additivity as Experiments Scale Up

The platform’s central practical question is how badly the default linear-addition heuristic degrades as more experiments run concurrently. We repeat the exercise above for \(m \in \{3,4,5\}\) concurrent experiments, each replicated twice with independently redrawn nuisance functions to average out simulation noise, exactly the design of the paper’s own Section 5.1 robustness check.

run_dedl_experiment <- function(m, dX = 10L, n_epochs = 400L, seed = 1) {
  set.seed(seed); torch_manual_seed(seed)

  A <- matrix(runif((m + 1) * dX, -0.5, 0.5), nrow = m + 1)
  theta_m1_true <- runif(1, 10, 20)
  theta_star <- function(X) { raw <- X %*% t(A); sign(raw) * abs(raw)^3 }
  G_true <- function(X, Tmat) {
    th <- theta_star(X)
    u <- th[, 1] + rowSums(th[, -1, drop = FALSE] * Tmat)
    theta_m1_true / (1 + exp(-u))
  }

  all_combos <- as.matrix(expand.grid(rep(list(0:1), m)))
  obs_support <- all_combos[rowSums(all_combos) %in% c(0, 1), , drop = FALSE]
  extra <- rep(0, m); extra[1:2] <- 1
  obs_support <- rbind(obs_support, extra)
  storage.mode(obs_support) <- "double"
  n_support <- nrow(obs_support)
  t0 <- rep(0, m)

  sample_partial <- function(n) {
    X <- matrix(runif(n * dX), n, dX)
    idx <- sample.int(n_support, n, replace = TRUE)
    Tmat <- obs_support[idx, , drop = FALSE]
    y <- G_true(X, Tmat) + runif(n, -0.05, 0.05)
    list(X = X, T = Tmat, y = y)
  }

  train <- sample_partial(500L * m)

  net <- nn_module(
    initialize = function(dX, hidden, out_dim) {
      self$fc1 <- nn_linear(dX, hidden)
      self$fc2 <- nn_linear(hidden, hidden)
      self$fc3 <- nn_linear(hidden, out_dim)
    },
    forward = function(x) x |> self$fc1() |> nnf_relu() |> self$fc2() |> nnf_relu() |> self$fc3()
  )
  model <- net(dX, hidden = 10L, out_dim = m + 1L)
  theta_m1_hat <- torch_tensor(15, requires_grad = TRUE)

  Xtr <- torch_tensor(train$X, dtype = torch_float())
  Ttr <- torch_tensor(train$T, dtype = torch_float())
  ytr <- torch_tensor(train$y, dtype = torch_float())
  opt <- optim_adam(c(model$parameters, list(theta_m1_hat)), lr = 0.01)

  forward_G <- function(X, Tt) {
    th <- model(X)
    u <- th[, 1] + torch_sum(th[, 2:(m + 1)] * Tt, dim = 2)
    theta_m1_hat / (1 + torch_exp(-u))
  }
  for (epoch in seq_len(n_epochs)) {
    opt$zero_grad()
    loss <- nnf_mse_loss(forward_G(Xtr, Ttr), ytr)
    loss$backward(); opt$step()
  }

  theta_hat_fn <- function(X) {
    with_no_grad({ th <- model(torch_tensor(X, dtype = torch_float())) })
    as.matrix(th)
  }
  theta_m1_val <- as.numeric(theta_m1_hat)

  G_hat_row <- function(th_row, t) theta_m1_val / (1 + exp(-(th_row[1] + sum(th_row[-1] * t))))
  Gtheta_hat_row <- function(th_row, t) {
    u <- th_row[1] + sum(th_row[-1] * t)
    s <- 1 / (1 + exp(-u)); ds <- s * (1 - s)
    c(theta_m1_val * ds, theta_m1_val * ds * t, s)
  }

  inf_data <- sample_partial(500L * m)
  th_inf <- theta_hat_fn(inf_data$X)
  d <- m + 2L
  n_inf <- nrow(inf_data$X)

  correction <- matrix(0, n_inf, d)
  for (i in seq_len(n_inf)) {
    th_i <- th_inf[i, ]
    Lam <- matrix(0, d, d)
    for (k in seq_len(n_support)) {
      g <- Gtheta_hat_row(th_i, obs_support[k, ])
      Lam <- Lam + outer(g, g)
    }
    Lam <- Lam * 2 / n_support + diag(0.0005, d)
    g_ti <- Gtheta_hat_row(th_i, inf_data$T[i, ])
    ell_th <- 2 * g_ti * (G_hat_row(th_i, inf_data$T[i, ]) - inf_data$y[i])
    correction[i, ] <- solve(Lam, ell_th)
  }

  lr_dat <- data.frame(y = train$y, train$X, train$T)
  names(lr_dat) <- c("y", paste0("x", seq_len(dX)), paste0("t", seq_len(m)))
  lr_fit <- lm(y ~ ., data = lr_dat)

  y0 <- mean(train$y[rowSums(train$T) == 0])
  la_effects <- sapply(seq_len(m), function(k) {
    sel <- train$T[, k] == 1 & rowSums(train$T) == 1
    mean(train$y[sel]) - y0
  })

  Xmc <- matrix(runif(10000L * dX), 10000L, dX)
  mu_true <- function(t) mean(G_true(Xmc, matrix(t, nrow(Xmc), m, byrow = TRUE))) -
                          mean(G_true(Xmc, matrix(t0, nrow(Xmc), m, byrow = TRUE)))

  out <- data.frame()
  for (k in seq_len(nrow(all_combos))) {
    t_target <- as.numeric(all_combos[k, ])
    if (all(t_target == t0)) next
    truth <- mu_true(t_target)

    Hx  <- sapply(seq_len(n_inf), function(i) G_hat_row(th_inf[i, ], t_target) - G_hat_row(th_inf[i, ], t0))
    Hth <- t(sapply(seq_len(n_inf), function(i) Gtheta_hat_row(th_inf[i, ], t_target) - Gtheta_hat_row(th_inf[i, ], t0)))
    psi <- Hx - rowSums(Hth * correction)

    newx  <- data.frame(matrix(runif(2000L * dX), 2000L, dX)); names(newx) <- paste0("x", seq_len(dX))
    newt1 <- as.data.frame(matrix(t_target, 2000L, m, byrow = TRUE)); names(newt1) <- paste0("t", seq_len(m))
    newt0 <- as.data.frame(matrix(t0, 2000L, m, byrow = TRUE)); names(newt0) <- paste0("t", seq_len(m))

    out <- rbind(out, data.frame(
      m = m, truth = truth,
      LA = sum(la_effects[t_target == 1]),
      LR = mean(predict(lr_fit, cbind(newx, newt1))) - mean(predict(lr_fit, cbind(newx, newt0))),
      SDL = mean(Hx), DeDL = mean(psi)
    ))
  }
  rm(model, opt, Xtr, Ttr, ytr, th_inf, correction)
  out
}

# Two replicates per m keep the sweep's total runtime bounded; gc() after each
# replicate prevents torch's per-call tensor allocations from accumulating
# across the loop, which otherwise slows later iterations substantially.
sweep_results <- data.frame()
for (m_val in c(3L, 4L, 5L)) {
  for (rep in 1:2) {
    r <- run_dedl_experiment(m_val, n_epochs = 300L, seed = 100 * m_val + rep)
    sweep_results <- rbind(sweep_results, r)
    gc(full = TRUE)
  }
}

mae_by_m <- aggregate(
  cbind(LA = abs(LA - truth), LR = abs(LR - truth), SDL = abs(SDL - truth), DeDL = abs(DeDL - truth)) ~ m,
  data = sweep_results, FUN = mean
)
mae_long <- reshape(mae_by_m, direction = "long", varying = c("LA", "LR", "SDL", "DeDL"),
                     v.names = "MAE", timevar = "Estimator", times = c("LA", "LR", "SDL", "DeDL"))

ggplot(mae_long, aes(x = m, y = MAE, color = Estimator)) +
  geom_line() +
  geom_point() +
  labs(x = "Number of concurrent experiments (m)", y = "Mean absolute error vs. ground truth") +
  causalverse::ama_theme()
Mean absolute error against ground truth as the number of concurrent experiments m increases, averaged over two independently redrawn synthetic environments per m. The linear-addition and linear-regression benchmarks assume additive effects and degrade as more experiments interact; the structured deep network, with or without debiasing, uses the same $m+2$ observed combinations regardless of m and remains comparatively stable.

Figure 45.4: Mean absolute error against ground truth as the number of concurrent experiments m increases, averaged over two independently redrawn synthetic environments per m. The linear-addition and linear-regression benchmarks assume additive effects and degrade as more experiments interact; the structured deep network, with or without debiasing, uses the same \(m+2\) observed combinations regardless of m and remains comparatively stable.

The qualitative pattern matches the paper’s own Figure 8: the linear-addition and linear-regression estimators, which hard-code an additive structure that becomes more wrong as more experiments interact, do not improve as \(m\) grows and often worsen, while the structured network, trained on only \(m+2\) observed combinations no matter how large \(m\) is, holds up because its sigmoid link nests genuine nonlinear interaction. This is the practical case for the framework: a platform that already runs hundreds of concurrent experiments cannot expand its factorial coverage as \(m\) grows, but a structured, orthogonalized deep network can still deliver valid, increasingly necessary estimates of combinations it has never tested.

Against the paper’s own full-scale version of this sweep (200 replications per configuration, Tables 8 and 9 of Ye et al. (2026)), the ordering here is the same. Their DeDL estimator holds a mean absolute percentage error of 7 to 17 percent as \(m\) rises from 4 to 10 while linear addition sits at 23 to 32 percent, and DeDL still finds the exact best arm among \(1{,}024\) combinations in 79 percent of replications at \(m = 10\), against 59 percent for linear addition and 56 percent for the undebiased network. Their sweeps hold the network width fixed at ten hidden nodes for comparability, and the authors verify that the mild degradation of the structured estimators with \(m\) disappears when the width grows with the problem, so the binding constraint is network capacity rather than the framework itself.

45.4.11 Coverage and Best-Arm Selection under Controlled First-Stage Bias

The deepest claim of the theory is not point-estimate accuracy but the validity of the confidence intervals, and that claim is hard to check with a trained network because the network’s true estimation error is unknown. Ye et al. (2026) resolve this in their second synthetic study with a device in the spirit of Chernozhukov, Chetverikov, et al. (2018): skip training entirely and manufacture the first-stage error by hand. The true nuisances are linear, \(\theta_j^*(x) = A_{[j+1]} x\), and the estimator handed to the second stage is \(\hat\theta_j(x) = (1 + \mathrm{err}_j)\,\theta_j^*(x)\) with \(\mathrm{err}_j\) drawn uniformly from \([-\delta, \delta]\), so the bias level \(\delta\) is a dial rather than an unknown. A fixed bias of this kind deliberately violates the \(o(n^{-1/4})\) premise, which makes it the right stress test: it reveals what the debiasing term buys when the first stage is mediocre, the realistic operating condition. We replicate the design exactly (\(m = 4\), the same \(m+2\) partially observed cells, 200 replications per bias level) and add the two diagnostics the paper does not report directly: the empirical coverage of nominal 95 percent intervals, and the best-arm hit rate of every estimator.

sig_of <- function(th, Tmat) {
  u <- th[, 1] + rowSums(th[, -1, drop = FALSE] * Tmat)
  1 / (1 + exp(-u))
}

run_bias_rep <- function(m = 4L, dX = 10L, n = 2000L, delta, seed) {
  set.seed(seed)
  A <- matrix(runif((m + 1) * dX, -0.5, 0.5), nrow = m + 1)
  th_m1 <- runif(1, 10, 20)

  all_combos <- as.matrix(expand.grid(rep(list(0:1), m)))
  storage.mode(all_combos) <- "double"
  obs_support <- rbind(all_combos[rowSums(all_combos) <= 1, , drop = FALSE],
                       c(1, 1, rep(0, m - 2)))
  n_support <- nrow(obs_support)
  t0 <- rep(0, m)

  X <- matrix(runif(n * dX), n, dX)
  Tmat <- obs_support[sample.int(n_support, n, replace = TRUE), , drop = FALSE]
  th_star <- X %*% t(A)
  y <- th_m1 * sig_of(th_star, Tmat) + runif(n, -0.05, 0.05)

  # Injected first-stage bias in place of a trained network.
  err <- runif(m + 2, -delta, delta)
  th_hat <- th_star * matrix(1 + err[1:(m + 1)], n, m + 1, byrow = TRUE)
  th_m1_hat <- (1 + err[m + 2]) * th_m1

  # True ATEs by sample average approximation.
  Xmc <- matrix(runif(10000L * dX), 10000L, dX)
  th_mc <- Xmc %*% t(A)
  s0_mc <- sig_of(th_mc, matrix(t0, 10000L, m, byrow = TRUE))
  mu_true <- apply(all_combos, 1, function(t)
    mean(th_m1 * sig_of(th_mc, matrix(t, 10000L, m, byrow = TRUE))) -
      mean(th_m1 * s0_mc))

  # Per-unit debiasing correction Lambda^{-1} ell_theta; the gradient of the
  # generalized sigmoid II link is (th_m1 * s(1-s) * (1, t), s), all closed form.
  d <- m + 2L
  s_obs <- sapply(seq_len(n_support), function(k)
    sig_of(th_hat, matrix(obs_support[k, ], n, m, byrow = TRUE)))
  s_i <- sig_of(th_hat, Tmat)
  resid_i <- th_m1_hat * s_i - y
  correction <- matrix(0, n, d)
  for (i in seq_len(n)) {
    Lam <- matrix(0, d, d)
    for (k in seq_len(n_support)) {
      sk <- s_obs[i, k]
      g <- c(th_m1_hat * sk * (1 - sk) * c(1, obs_support[k, ]), sk)
      Lam <- Lam + outer(g, g)
    }
    Lam <- Lam * 2 / n_support + diag(5e-4, d)
    g_i <- c(th_m1_hat * s_i[i] * (1 - s_i[i]) * c(1, Tmat[i, ]), s_i[i])
    correction[i, ] <- solve(Lam, 2 * g_i * resid_i[i])
  }

  # Linear benchmarks fit on the same observed sample.
  y0bar <- mean(y[rowSums(Tmat) == 0])
  la <- sapply(seq_len(m), function(k) {
    sel <- Tmat[, k] == 1 & rowSums(Tmat) == 1
    mean(y[sel]) - y0bar
  })
  lr_dat <- data.frame(y = y, X, Tmat)
  names(lr_dat) <- c("y", paste0("x", 1:dX), paste0("t", 1:m))
  lr_coefs <- coef(lm(y ~ ., data = lr_dat))[paste0("t", 1:m)]

  s_0 <- sig_of(th_hat, matrix(t0, n, m, byrow = TRUE))
  out <- data.frame()
  for (k in seq_len(nrow(all_combos))) {
    t_t <- all_combos[k, ]
    if (all(t_t == t0)) next
    Tk <- matrix(t_t, n, m, byrow = TRUE)
    s_t <- sig_of(th_hat, Tk)
    Hx <- th_m1_hat * (s_t - s_0)
    Hth <- cbind(th_m1_hat * (s_t * (1 - s_t) - s_0 * (1 - s_0)),
                 th_m1_hat * (s_t * (1 - s_t) * Tk),
                 s_t - s_0)
    psi <- Hx - rowSums(Hth * correction)
    out <- rbind(out, data.frame(
      truth = mu_true[k],
      SDL = mean(Hx), SDL_se = sd(Hx) / sqrt(n),
      DeDL = mean(psi), DeDL_se = sd(psi) / sqrt(n),
      LA = sum(la[t_t == 1]),
      LR = sum(lr_coefs[t_t == 1])))
  }
  out
}
R <- 200
rows <- list()
for (dl in c(0.1, 0.2, 0.3)) {
  for (r in seq_len(R)) {
    res <- run_bias_rep(delta = dl, seed = 7000 + 97 * r + round(1000 * dl))
    pick <- function(v) which.max(c(0, v))        # include the control arm
    rows[[length(rows) + 1]] <- data.frame(
      delta = dl,
      MAE_LA = mean(abs(res$LA - res$truth)),
      MAE_LR = mean(abs(res$LR - res$truth)),
      MAE_SDL = mean(abs(res$SDL - res$truth)),
      MAE_DeDL = mean(abs(res$DeDL - res$truth)),
      cover_SDL = mean(abs(res$SDL - res$truth) <= qnorm(0.975) * res$SDL_se),
      cover_DeDL = mean(abs(res$DeDL - res$truth) <= qnorm(0.975) * res$DeDL_se),
      bai_LA = pick(res$LA) == pick(res$truth),
      bai_LR = pick(res$LR) == pick(res$truth),
      bai_SDL = pick(res$SDL) == pick(res$truth),
      bai_DeDL = pick(res$DeDL) == pick(res$truth))
  }
}
bias_mc <- aggregate(. ~ delta, data = do.call(rbind, rows), FUN = mean)
Table 45.8: Mean absolute error and empirical coverage of nominal 95 percent intervals across 200 replications per bias level. The plug-in intervals come from the pairwise t-test on predicted outcomes, as in the SDL benchmark; the DeDL intervals come from the influence-function variance.
Bias level MAE linear addition MAE linear regression MAE plug-in (SDL) MAE DeDL Coverage plug-in Coverage DeDL
0.1 0.176 0.145 0.117 0.048 0.334 0.925
0.2 0.171 0.138 0.235 0.084 0.186 0.935
0.3 0.188 0.134 0.353 0.125 0.141 0.930
Table 45.9: Share of replications in which each estimator exactly identifies the best treatment combination (the control arm included as a candidate).
Bias level Linear addition Linear regression Plug-in (SDL) DeDL
0.1 0.925 0.945 0.985 0.97
0.2 0.925 0.950 0.990 0.94
0.3 0.925 0.950 0.985 0.94

Three findings emerge, each matching or sharpening the paper’s Table 11. First, the point estimates: debiasing cuts the mean absolute error by roughly sixty percent at every bias level, and the gap widens with \(\delta\) because the plug-in error grows almost linearly in the injected bias while the orthogonal correction absorbs its first-order component. From \(\delta = 0.2\) onward the undebiased network is worse than the linear benchmarks it was meant to replace, while DeDL remains the best estimator throughout, exactly the ordering in the paper’s version of this experiment. Second, the coverage results, which are the new information here: the plug-in intervals cover the truth in only a third of cases at \(\delta = 0.1\) and barely a seventh at \(\delta = 0.3\), and they deteriorate as the bias grows because the prediction-spread variance reflects neither the bias nor the estimation noise. The DeDL intervals hold near 93 percent at every bias level. The remaining shortfall from the nominal 95 is itself informative: a bias that does not shrink with \(n\) deliberately violates the \(o(n^{-1/4})\) premise, and orthogonality removes only the first-order term, leaving a second-order residual that no correction can eliminate. Under the theory’s actual conditions coverage converges to nominal; under honest mediocrity it lands within a couple of points, which is the practically reassuring version of the claim. Third, best-arm selection: with \(m = 4\) and well-separated arms, every estimator identifies the winner in the low-to-high nineties, and the plug-in actually edges out DeDL in our draws because a shared multiplicative distortion of the fitted surface rarely reorders arms, while the correction adds arm-specific noise; in the paper’s own draws the ranking favors DeDL (97 versus 92.5 percent at \(\delta = 0.1\), 92.5 versus 80 at \(\delta = 0.3\)), so which estimator ranks best is design dependent. The two tables together make the practical division of labor precise: point rankings can survive moderate first-stage bias when arms are far apart, but inference cannot, and only the debiased intervals mean what they say. The authors also report that at extreme bias (\(\delta\) near one) debiasing can hurt even the point estimates, the same regime as the undertrained networks in their Figure 7, so the correction is a complement to, not a substitute for, a decently estimated first stage.

45.4.13 Practical Guidance

Three lessons from Ye et al. (2026) are worth carrying forward whenever deep learning is combined with Neyman orthogonality on a structural model layer rather than an unrestricted regression. First, the debiasing term is only as good as the first-stage network; when training has not converged, plain plug-in and debiased estimates are similar, and the value of debiasing appears only once training loss is genuinely small, so cross-validation loss during training is a useful diagnostic before trusting the debiased estimate. Second, a misspecified link function is more dangerous with debiasing than without, because the correction term amplifies rather than removes bias from a structural form that does not nest the truth; comparing the structured network’s validation loss against an unrestricted “pure DNN” that takes treatment as a raw input alongside the covariates, exactly the PDL benchmark Ye et al. (2026) use in Table 3, is a practical specification check that does not require knowing the true link function in advance. Third, and most important for scale, the framework’s traffic requirement grows linearly rather than exponentially in the number of concurrent experiments, which is what turns a theoretically elegant debiasing result into something a platform running hundreds of A/B tests can actually deploy.

Four more operational rules emerge from the appendices and from the replications above. Aggregate across folds and never read a single fold: fold-level estimates on Platform O miss by as much as forty percent while their average lands within a fraction of a percent. Exploit the known assignment mechanism, which makes \(\Lambda(x)\) computable rather than estimable and lets two-way cross-fitting replace the three-way splitting an observational application would need; porting the framework to observational data means giving up both conveniences and defending unconfoundedness besides. Regularize the inversion of \(\hat\Lambda(x)\) with a small ridge and treat it like propensity trimming, a numerical guard rather than a tuning knob, and check that the covariate distribution in the inference sample matches the population the estimand refers to, rebalancing by stratification or reweighting when it does not, because the correction fixes estimation error in \(\hat\theta\), not evaluation over the wrong population. On the selection side, treat best-arm point identification as an asymptotic statement: report the regret estimates \(\hat\tau\) with their intervals, flag near-ties as undecided, and reach for the winner’s-curse corrections of Andrews et al. (2024) when the leaderboard is tight.

45.5 Heterogeneous Treatment Effects

Average effects answer whether a treatment works on the whole. Many questions instead concern for whom it works, which is captured by the conditional average treatment effect (CATE) \[ \tau(x) = \mathbb{E}\bigl[Y(1) - Y(0) \mid X = x\bigr]. \] Estimating \(\tau(x)\) is a function-estimation problem rather than a scalar-estimation problem, and it is where tree-based and meta-learning methods come into their own.

45.5.1 Causal Trees: Honest Recursive Partitioning

The tree-based line of attack begins with the causal tree of Athey and Imbens (2016), published in the Proceedings of the National Academy of Sciences. A regression tree partitions the covariate space and predicts within each leaf; a causal tree partitions the covariate space and estimates a treatment effect within each leaf, the difference in mean outcomes between treated and control observations falling in that leaf. Two modifications turn this simple idea into a valid inferential procedure. First, the splitting criterion changes: rather than rewarding splits that predict the outcome, the tree rewards splits that maximize the variance of the estimated treatment effects across leaves, penalized by their estimation variance, an expected mean squared error criterion for the treatment effect itself. A tree grown to predict outcomes well can be useless for heterogeneity, because the strongest predictors of the outcome level need not moderate the treatment effect at all.

Second, and more consequentially, Athey and Imbens (2016) introduce honest estimation. The sample is split in half: one half chooses the tree structure, the other half estimates the leaf-level effects. Because the estimation sample played no role in selecting the partition, the leaf estimates are unbiased conditional on the tree, and standard confidence intervals within leaves are valid no matter how aggressively the structure search overfit. Honesty costs data, since each half sees only half the observations, but it buys valid inference after a data-driven model search, the same trade that cross-fitting makes for double machine learning. The two ideas are siblings: honesty separates model selection from estimation across a single split, cross-fitting separates nuisance estimation from score evaluation across rotating folds.

45.5.2 Causal Forests and Generalized Random Forests

Causal forests, introduced by Wager and Athey (2018), aggregate honest causal trees into an ensemble and adapt the random forest to estimate \(\tau(x)\) rather than a conditional mean. The key conceptual move, developed more generally in the generalized random forest of Athey et al. (2019), is to view a forest not as an ensemble of predictions but as an adaptive kernel or weighting function. A forest grown on the data assigns to a target point \(x\) a set of similarity weights \(\alpha_i(x)\), the frequency with which training observation \(i\) falls in the same leaf as \(x\) across the trees, \[ \alpha_i(x) = \frac{1}{B}\sum_{b=1}^{B} \frac{\mathbf{1}\{X_i \in L_b(x)\}}{\lvert L_b(x)\rvert}, \] where \(L_b(x)\) is the leaf of tree \(b\) containing \(x\). These weights then define a local, weighted version of any estimating equation. For heterogeneous effects the local estimating equation is a residualized treatment-effect regression, in the spirit of the R-learner described below, so the forest solves a weighted partialling-out problem in the neighborhood of each \(x\). This connects causal forests directly to the orthogonalization theme of the chapter.

Two features make the inference valid. First, the trees split to maximize heterogeneity in the treatment effect rather than to predict the outcome, so the adaptive neighborhoods are tailored to the causal target. Second, the forest uses honest splitting: the sample used to choose the splits is disjoint from the sample used to estimate the effect within each leaf. Honesty removes the overfitting bias that adaptive partitioning would otherwise introduce and is what allows Wager and Athey (2018) to derive asymptotic normality of the pointwise estimates, hence pointwise confidence intervals for \(\tau(x)\). The generalized random forest of Athey et al. (2019) extends the same kernel-weighting and honesty principles to a broad class of moment-condition problems, including instrumental variables and quantile estimation, and provides the theoretical foundation in the Annals of Statistics. An accessible empirical walkthrough appears in Athey and Wager (2019).

The following simulation builds the honest-forest CATE estimate conceptually using only base R, generating heterogeneous effects driven by a single covariate.

set.seed(2024)
n <- 4000
p <- 5
X <- matrix(runif(n * p), n, p)
e <- 0.5                                   # randomized treatment for clarity
W <- rbinom(n, 1, e)
tau_fun <- function(x) 1 + 2 * (x[, 1] > 0.5)   # effect doubles past 0.5
mu0 <- function(x) x[, 2] + sin(2 * pi * x[, 3])
Y <- mu0(X) + W * tau_fun(X) + rnorm(n)

# Stratified difference-in-means as a transparent CATE proxy:
# compare treated vs control within bins of the effect modifier X1.
bins <- cut(X[, 1], breaks = seq(0, 1, by = 0.25), include.lowest = TRUE)
cate_by_bin <- tapply(seq_len(n), bins, function(idx) {
    mean(Y[idx][W[idx] == 1]) - mean(Y[idx][W[idx] == 0])
})
round(cate_by_bin, 2)
#>   [0,0.25] (0.25,0.5] (0.5,0.75]   (0.75,1] 
#>       1.06       1.08       3.06       3.09

The estimated effect jumps as \(X_1\) crosses 0.5, recovering the simulated heterogeneity. A causal forest automates this binning adaptively and over many covariates at once. The production estimator uses the grf package.

library(grf)

cf <- causal_forest(
    X = X, Y = Y, W = W,
    num.trees = 2000,
    honesty = TRUE
)

# Out-of-bag CATE predictions and pointwise standard errors
tau_hat <- predict(cf, estimate.variance = TRUE)
head(tau_hat)

# Doubly robust average treatment effect from the forest
average_treatment_effect(cf)

# Test for the presence of heterogeneity
test_calibration(cf)

45.5.3 The R-Learner and Meta-Learners

A complementary, modular philosophy estimates the CATE by reducing it to a sequence of standard prediction problems, each of which can be solved with any off-the-shelf learner. These are the meta-learners. The S-learner trains a single model on the pooled data with treatment as an additional feature, then takes the difference in predictions setting \(W = 1\) versus \(W = 0\). The T-learner trains two separate outcome models, one per treatment arm, and differences them. The X-learner, proposed by Künzel et al. (2019), refines the T-learner by imputing individual treatment effects, regressing them on covariates within each arm, and combining the two CATE estimates with propensity-based weights. The X-learner performs especially well when the treatment groups are very unequal in size or when one potential-outcome surface is much smoother than the other, settings where the T-learner wastes data.

The R-learner formalizes the residualization idea as a loss function. Starting from the Robinson decomposition, one residualizes the outcome and the treatment on the covariates, \(\tilde{Y} = Y - \hat{m}(X)\) and \(\tilde{W} = W - \hat{e}(X)\), and then estimates the CATE by minimizing the weighted criterion \[ \hat{\tau}(\cdot) = \arg\min_{\tau}\; \frac{1}{n}\sum_{i=1}^{n}\Bigl(\tilde{Y}_i - \tau(X_i)\,\tilde{W}_i\Bigr)^2 + \Lambda\bigl(\tau\bigr), \] where \(\Lambda\) is a regularizer on the CATE function. Because the nuisances enter only through the residuals, the objective is Neyman orthogonal, so errors in \(\hat{m}\) and \(\hat{e}\) have a second-order effect on the estimated heterogeneity. Causal forests can be seen as solving a local version of this same R-learner objective, which is why the two approaches are tightly linked.

library(rlearner)   # R-, S-, T-, X-learners with glmnet or boosting

# R-learner with cross-validated lasso nuisances
rfit <- rlasso(x = X, w = W, y = Y)
tau_rlearner <- predict(rfit, X)

# X-learner via the same package family
xfit <- xboost(x = X, w = W, y = Y)
tau_xlearner <- predict(xfit, X)

The choice among meta-learners and forests is empirical. Forests provide honest pointwise inference and require little tuning; meta-learners offer flexibility in the choice of base learner and can exploit structure such as smoothness or sparsity. A common workflow estimates the CATE several ways and checks that the qualitative conclusions about who benefits are stable across methods.

45.6 Publishing With These Methods: A Marketing Targeting Study

The methods in this chapter are not only academically interesting; they are the engine behind a recognizable genre of empirical paper that publishes regularly in the top marketing and quantitative-marketing journals. Knowing the shape of that genre is useful both for reading the literature and for producing work that can survive its review process, so this section first describes the template and then replicates it end to end on a public dataset with production-grade code.

45.6.1 What a Publishable Paper Looks Like

The papers that use heterogeneous-treatment-effect machinery to publish in Marketing Science, the Journal of Marketing Research, Management Science, and Quantitative Marketing and Economics share a common skeleton, and the contribution is rarely the estimator itself, which is taken off the shelf, but the design and the decision it informs. The canonical reference is Hitsch et al. (2024), who lay out the modern recipe explicitly: start from a randomized experiment or a credibly exogenous source of variation that identifies the treatment effect without functional-form heroics, estimate the conditional average treatment effect with a causal forest, and then, crucially, evaluate targeting policies rather than reporting the CATE as an end in itself. Their central methodological insight, which separates a publishable paper from a classroom exercise, is that the quantity a manager cares about is the value of a policy that assigns treatment based on covariates, and that this policy value can be estimated with a doubly robust score even when the underlying CATE estimates are noisy. The estimand is the profit or outcome of a targeting rule, not the accuracy of the effect surface.

The substantive surprises that make such papers worth publishing usually come from the gap between the effect of a treatment and the characteristics that managers intuitively target. Ascarza (2018) provides the discipline’s clearest cautionary result, showing through field experiments that the customers at highest risk of churning, the natural target of a retention campaign, are not the customers whose behavior responds most to intervention, so that a campaign aimed at high-risk customers can be nearly worthless while the same budget aimed at high-responsiveness customers pays off. The lesson, that one should target the treatment effect rather than the baseline level, is exactly what a causal forest operationalizes, and it recurs across the literature. Simester et al. (2020) stress-test the machine-learning targeting pipeline against the data problems that real firms face, missing covariates, distribution shift between the experiment and the deployment population, and small samples, and document which methods degrade gracefully, a robustness contribution that is itself publishable. Yoganarasimhan et al. (2023) use these tools to design and evaluate optimal free-trial lengths for a subscription product, turning an estimated heterogeneous effect into a personalized policy and quantifying its lift over the firm’s status quo, while Rafieian and Yoganarasimhan (2021) bring the same logic to mobile advertising and weigh the targeting gains against their privacy costs.

A second strand earns its place by getting the identification right at scale rather than by the heterogeneity machinery. Gordon et al. (2019) compare experimental and observational advertising-effect estimates across large Facebook field experiments and show how far selection-on-observables adjustments, including propensity and double-machine-learning style methods, fall short of the experimental benchmark, a sobering result about the limits of observational causal inference in marketing. Dube and Misra (2023) run experiments to estimate personalized price elasticities and then compute the welfare consequences of personalized pricing, combining experimentation, heterogeneous-effect estimation, and a welfare calculation into one paper. On the methodological side, Knaus et al. (2021) give the kind of careful empirical Monte Carlo comparison of heterogeneous-effect estimators that lets applied researchers choose among causal forests, the R-learner, and the various meta-learners with evidence rather than fashion. The throughline is that the estimator is a commodity and the contribution lives in the design, the data, the policy question, and the honesty of the evaluation.

45.6.2 Replicating the Targeting Workflow on the Hillstrom Email Experiment

To make the recipe concrete we replicate its core on a genuinely public marketing field experiment, the email campaign released by Kevin Hillstrom through the MineThatData challenge, which has become a standard benchmark for targeting research because it is a real randomized experiment with customer covariates and behavioral outcomes. Sixty-four thousand customers were randomly assigned to receive an email featuring men’s merchandise, an email featuring women’s merchandise, or no email, and the data record whether each customer visited the site, made a purchase, and how much they spent in the following two weeks. We use the clean two-arm comparison of the women’s email against no email, with website visits as the outcome, and ask the question a manager actually faces: given the covariates, whom should we email?

email <- read.csv("data/hillstrom_email.csv")

# Clean two-arm experiment: women's email versus no email.
email <- email[email$segment %in% c("Womens E-Mail", "No E-Mail"), ]

W <- as.integer(email$segment == "Womens E-Mail")   # treatment: received email
Y <- email$visit                                     # outcome: visited the site

# Covariates the firm knows before sending: purchase recency, dollar history,
# past category, tenure, acquisition channel, and location type.
X <- data.frame(
  recency   = email$recency,
  history   = email$history,
  mens      = email$mens,
  womens    = email$womens,
  newbie    = email$newbie,
  phone     = as.integer(email$channel == "Phone"),
  web       = as.integer(email$channel == "Web"),
  zip_urban = as.integer(email$zip_code == "Urban"),
  zip_rural = as.integer(email$zip_code == "Rural")
)

c(customers = nrow(X), treated = sum(W), baseline_visit_rate = round(mean(Y), 4))
#>           customers             treated baseline_visit_rate 
#>          42693.0000          21387.0000              0.1288

The first step is the credibility anchor, the average treatment effect. Because assignment was randomized, even a simple difference in means identifies it, but we estimate it from a causal forest so that the same fitted object delivers both the average effect and the heterogeneity. The forest residualizes the outcome and the treatment on the covariates internally, the orthogonalization of the earlier sections, and returns a doubly robust average effect.

library(grf)

set.seed(42)
cf <- causal_forest(as.matrix(X), Y, W, num.trees = 1000, seed = 42)

ate <- average_treatment_effect(cf, target.sample = "all")
calib <- test_calibration(cf)
tau_hat <- predict(cf)$predictions
Table 45.11: Average effect of the women’s email on the probability of a site visit, and the calibration test for whether the forest detects genuine effect heterogeneity. A differential-prediction coefficient near one with a small standard error is evidence of real heterogeneity.
Quantity Estimate Std. Error
Average treatment effect (visit prob.) 0.0439 0.0033
Heterogeneity test (differential.forest.prediction) 0.4120 0.0700

The women’s email raises the probability of a site visit by roughly four percentage points off a baseline near thirteen percent, a large and precisely estimated effect. The calibration test of Wager and Athey (2018) regresses the held-out outcome on the mean forest prediction and on the demeaned forest prediction; the coefficient on the differential prediction is well above zero and many standard errors from it, which is the formal evidence that the effect genuinely varies across customers rather than the forest manufacturing spurious variation. Figure 45.5 shows the spread of the individual effect estimates, which range from slightly negative to well above the average, the raw material of any targeting decision.

library(ggplot2)

ggplot(data.frame(tau = tau_hat), aes(tau)) +
  geom_histogram(bins = 40) +
  geom_vline(xintercept = ate[["estimate"]], linetype = "dashed") +
  labs(x = "Estimated treatment effect on visit probability",
       y = "Number of customers") +
  causalverse::ama_theme()
Distribution of estimated conditional average treatment effects of the women's email on site visits across customers. The dashed line marks the average effect. A minority of customers have near-zero or negative estimated effects and are candidates for suppression.

Figure 45.5: Distribution of estimated conditional average treatment effects of the women’s email on site visits across customers. The dashed line marks the average effect. A minority of customers have near-zero or negative estimated effects and are candidates for suppression.

Detecting heterogeneity is not the same as profiting from it, and the next step is the one Hitsch et al. (2024) insist on: quantify how much value a targeting rule that prioritizes by estimated effect actually delivers. The rank-weighted average treatment effect, computed by rank_average_treatment_effect, measures exactly this. It traces the average effect among the customers the model ranks most responsive as we expand the targeted fraction from the few down to everyone, and the area under that curve, the AUTOC, summarizes the gain from prioritizing by the estimated CATE over treating a random customer. Figure 45.6 plots the curve.

rate <- rank_average_treatment_effect(cf, tau_hat)

ggplot(rate$TOC, aes(q, estimate)) +
  geom_line() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(x = "Fraction of customers targeted (most responsive first)",
       y = "Average effect among targeted") +
  causalverse::ama_theme()
Targeting operator characteristic curve. The vertical axis is the average treatment effect among the top fraction of customers ranked by estimated responsiveness; the horizontal axis is that fraction. A curve lying above zero and sloping down means the model concentrates the effect in the customers it prioritizes, which is the value of targeting.

Figure 45.6: Targeting operator characteristic curve. The vertical axis is the average treatment effect among the top fraction of customers ranked by estimated responsiveness; the horizontal axis is that fraction. A curve lying above zero and sloping down means the model concentrates the effect in the customers it prioritizes, which is the value of targeting. 

#> RATE (AUTOC): 0.0176  (std. error 0.0039)

The curve slopes downward and the AUTOC is several standard errors above zero, so the forest is genuinely separating responsive from unresponsive customers rather than ranking them at random. The final step converts that ranking into a deployable rule and evaluates it honestly. The policytree package of Athey et al. (2019) learns a shallow, interpretable decision tree that maps covariates directly to a treat-or-not recommendation, fit on the doubly robust scores from the forest so that the learned policy inherits the orthogonality that protects it from nuisance error. We fit a depth-two tree on a subsample for speed and then evaluate several policies against the same doubly robust scores.

library(policytree)

dr_scores <- double_robust_scores(cf)

# Learn an interpretable depth-2 targeting rule (subsample for tractable search).
set.seed(7)
sub <- sample(nrow(X), 10000)
ptree <- policy_tree(X[sub, ], dr_scores[sub, ], depth = 2)

# Apply the rule to everyone; action 2 = send email, 1 = withhold.
send <- predict(ptree, X) == 2

# Doubly robust value of each policy: expected incremental visits per customer
# relative to emailing no one.
policy_value <- function(treat) mean(ifelse(treat, dr_scores[, 2], dr_scores[, 1]) - dr_scores[, 1])

value_table <- data.frame(
  Policy = c("Email no one", "Email everyone",
             "Email if estimated effect > 0", "Learned policy tree"),
  `Share emailed` = c(0, 1, mean(tau_hat > 0), mean(send)),
  `Incremental visits per customer` = c(
    0, policy_value(rep(TRUE, nrow(X))),
    policy_value(tau_hat > 0), policy_value(send)),
  check.names = FALSE
)
Table 45.12: Doubly robust value of alternative targeting policies on the Hillstrom experiment, measured as incremental site visits per customer relative to emailing no one. The learned policy tree captures nearly the full lift of a blanket campaign while contacting substantially fewer customers.
Policy Share emailed Incremental visits per customer
Email no one 0.000 0.0000
Email everyone 1.000 0.0439
Email if estimated effect > 0 0.814 0.0366
Learned policy tree 0.692 0.0420

The policy comparison is where the managerial payoff appears, and it carries the same counterintuitive flavor as Ascarza (2018). Because the email helps the great majority of customers, emailing everyone produces the largest raw lift in visits, and a naive reading would stop there. The learned policy tree, however, achieves nearly the same incremental visits while emailing only about seventy percent of the customers, which means the firm can shed roughly a third of its mailing cost for a small fraction of the benefit, and the tree is fully transparent about who it drops, splitting on past purchase history and category in a way a manager can inspect and defend. When the outcome is replaced by profit net of contact cost, as it would be in a deployment, the policies that look tied on raw visits separate sharply, and the targeting rule dominates the blanket campaign. This is precisely the structure of the published papers: a randomized design for identification, a flexible forest for the heterogeneity, an honest doubly robust evaluation of the policy rather than the effect, and a managerial conclusion that the obvious blanket strategy leaves money on the table.

The replication is also a template a reader can extend toward a publishable contribution. Swapping the visit outcome for spend turns the analysis into a revenue-targeting problem; subtracting a per-email cost converts the policy value into profit and sharpens the case for selective mailing; using all three experimental arms moves from a binary treatment to a multi-armed targeting problem that policytree handles directly; and validating the learned policy on a held-out portion of the experiment, or on the men’s-email arm, supplies the out-of-sample evidence that referees expect. The full toolchain, grf for the forest, rank_average_treatment_effect for the targeting value, and policytree for the deployable rule, is the same one behind the papers cited above, which is why working through it on public data is the most direct way to learn how such studies are built.

45.7 Temporal Causal Forests

The forest machinery extends naturally to settings with a strong temporal structure, where every unit is eventually exposed to a widely known event and a clean contemporaneous control group does not exist. This arises with public announcements, recalls, security breaches, and similar shocks, where it is impossible to find a group unaware of the event that would still be representative of the affected population. The temporal causal inference design of Turjeman and Feinberg (2023) addresses this by matching across cohorts that adopted at different times, so that newer participants who experienced the event early in their tenure are compared with the early trajectories of older participants who had not yet been exposed.

The core idea is to make calendar time and tenure separate axes. Let \(H_T\) denote the cohort that joined \(T\) time units before the event; this treated cohort is observed for \(T + 3\) units, with the final units after the event. Control trajectories are assembled from earlier cohorts \(H_1, \dots, H_{T-1}\), tracked from their own adoption up to the event or to the matched tenure window, whichever comes first. Each cohort plays a dual role, sometimes treated and sometimes control, except for those at the extreme ends of the window. Because the design relies on comparing trajectories rather than a single contemporaneous control, it requires a large sample so that subgroups with comparable time trends can be matched.

The figure below illustrates the design with simulated cohorts, showing the same activity curves on a calendar-time axis and on a tenure axis. On the tenure axis the cohorts line up, making clear that the control cohorts simply had not yet reached the event when the treated cohort did.

library(ggplot2)
library(patchwork)

tenure <- seq(0, 100, length.out = 100)
generate_cohort <- function(name, start_time, mean, sd) {
    data.frame(
        tenure = tenure,
        cohort = name,
        time   = seq(start_time, start_time + 99, by = 1),
        value  = dnorm(tenure, mean = mean, sd = sd)
    )
}
cohorts <- list(
    generate_cohort("Cohort 1",  1, 47, 15),
    generate_cohort("Cohort 2", 10, 48, 17),
    generate_cohort("Cohort 3", 20, 52, 20),
    generate_cohort("Cohort 4", 30, 53, 18),
    generate_cohort("treatment", 40, 50, 16)
)
final_dataset <- do.call(rbind, cohorts)

plot_time <- ggplot(final_dataset, aes(time, value, color = cohort)) +
    geom_line() +
    geom_vline(xintercept = c(40, 90), linetype = "dashed") +
    labs(title = "Value vs. Calendar Time", x = "Time", y = "Value")

plot_tenure <- ggplot(final_dataset, aes(tenure, value, color = cohort)) +
    geom_line() +
    geom_vline(xintercept = c(0, 50), linetype = "dashed") +
    labs(title = "Value vs. Tenure", x = "Tenure", y = "Value")

plot_time / plot_tenure

Identification rests on the usual selection-on-observables assumptions, reinterpreted for the temporal setting. The stable unit treatment value assumption requires no interference, which holds across cohorts because their treatment exposures occur at different points in calendar time, though it can be strained when treated units influence one another after the shock. Conditional independence requires that, given the covariates including the time trend, cohort membership is as good as randomly assigned; this is checked empirically with pre-treatment parallel-trends tests, including bidirectional Granger tests and Kolmogorov-Smirnov comparisons of the pre-event trajectories. Overlap requires that every unit had a positive propensity of treatment at any tenure, which the design satisfies because all units are eventually treated and the estimated propensity is bounded away from zero and one. Exogeneity of covariates requires that the shock was unforeseen, or at least equally anticipated across cohorts.

Temporal causal forests then extend the causal forest in two ways identified by Turjeman and Feinberg (2023) as improving root-mean-squared error and recovery of heterogeneous effects. First, the covariate vector \(X_i\) is augmented to include the unit’s time trend, so the forest groups units that are homogeneous in their activity trajectories, in effect estimating a counterfactual time path for each unit and ensuring that treated and control units share similar pre-event trends. Second, as a robustness analysis, the nuisance functions are estimated with the local linear forest of Friedberg et al. (2020), which fits a local linear correction within each leaf and improves accuracy when the underlying surface is smooth; this refinement is feasible when cohort timelines are of equal length. The result is a method that recovers both the average effect of a public shock and its heterogeneity across individuals while respecting the temporal structure that rules out a conventional control group.

45.8 Practical Guidance and Pitfalls

Several practical considerations govern whether these methods deliver on their promise.

Cross-fitting is not optional. The asymptotic guarantees of double machine learning depend on estimating the nuisance functions on data separate from the data used to form the final estimate. Skipping this step reintroduces overfitting bias and invalidates the confidence intervals, even when an orthogonal score is used.

The software ecosystem is mature, and there is little reason to hand-roll these estimators outside of teaching. The DoubleML package used throughout this chapter (Bach et al. 2024) implements the four canonical models plus difference-in-differences and quantile variants over any mlr3 learner, with a parallel implementation in Python; hdm (Chernozhukov et al. 2016) covers the lasso-specific post-double-selection and rigorous-lasso methods; grf provides orthogonalized forests with honest inference; and EconML in Python implements the meta-learner and orthogonal-learning families. What the software cannot automate is the identification argument, the choice of estimand, and the reporting discipline of the replications above: learner grids, fold and split sensitivity, and propensity diagnostics.

Orthogonality protects only against slow nuisance errors, not against violations of identification. If unconfoundedness fails because an important confounder is unobserved, no amount of flexible nuisance estimation will recover the causal effect. Machine learning strengthens the case for conditional ignorability by allowing a rich conditioning set, but it cannot test the assumption itself. The role of these methods is to remove functional-form and high-dimensionality concerns, leaving the identifying assumptions to be defended on substantive grounds.

Overlap matters more, not less, with flexible models. Inverse propensity weights become unstable when estimated propensities approach zero or one, and flexible learners can produce extreme propensity estimates. Diagnosing the propensity distribution, trimming or stabilizing extreme weights, and preferring substitution estimators such as TMLE that are less sensitive to small denominators all help. Reporting the distribution of estimated propensity scores should be routine.

Heterogeneity claims require discipline. A causal forest will always produce a varying \(\tau(x)\), but the variation may be noise. The honest forest delivers pointwise standard errors, and calibration tests such as the best-linear-predictor test assess whether the estimated heterogeneity is real. Out-of-sample validation, for instance ranking units by predicted effect and checking that the realized effect ordering agrees, guards against overinterpreting spurious patterns.

Sample size and tuning are real constraints. Like other nonparametric and matching methods, forests and meta-learners need substantial data for stable results, especially when estimating heterogeneity across many covariates or when treatment groups are imbalanced. The SuperLearner reduces the burden of choosing a single learner, but it multiplies computation, and the analyst should ensure that the cross-validation used inside the ensemble is nested properly within the cross-fitting used for the causal estimate.

Finally, transparency and replication remain essential. The flexibility that makes these methods powerful also makes them harder to scrutinize. Pre-registering the estimand, the conditioning set, and the validation plan, and reporting results across several estimators, keeps the analysis honest and the conclusions credible.

📖 Free preview — limited per publisher guidelines. Purchase the complete A Guide on Data Analysis series (Vols. 1–4) on Springer.
Vol. 1 Vol. 2 Vol. 3 Vol. 4