Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
134,865 characters · 16 sections · 167 citation commands
An Introduction to Double/Debiased Machine Learning
\emergencystretch 3em
This paper provides an introduction to Double/Debiased Machine Learning (DML). DML is a general approach to performing inference about a target parameter in the presence of nuisance functions: objects that are needed to identify the target parameter but are not of primary interest. Nuisance functions arise naturally in many settings, such as when controlling for confounding variables or leveraging instruments. The paper describes two biases that arise from nuisance function estimation and explains how DML alleviates these biases. Consequently, DML allows the use of flexible methods, including machine learning tools, for estimating nuisance functions, reducing the dependence on auxiliary functional form assumptions and enabling the use of complex non-tabular data, such as text or images. We illustrate the application of DML through simulations and empirical examples. We conclude with a discussion of recommended practices. A companion website includes additional examples with code and references to other resources.
A large share of empirical research in economics aims to provide insights into the statistical relationships among two or more variables. For example, a common research goal is to understand the causal impact of a policy on economic outcomes. Target parameters summarizing these relationships, including average treatment effects or regression coefficients, frequently depend on nuisances---auxiliary objects that must be accounted for to identify the parameter of interest but are not themselves of primary interest, such as regression coefficients on control variables.
As a concrete example that we will revisit in our empirical illustrations, consider dube2020, who study monopsony power on the online platform MTurk using partially linear regression:
where $Y$ is the logarithm of the time it takes for a posted job to be filled, $D$ is the logarithm of the reward of the job, $X$ denotes observed features of tasks including their type and complexity, and $\varepsilon$ is assumed to be uncorrelated with $D$ and mean independent of $X$, i.e., $\mathrm{E}[D \varepsilon] = 0$ and $\mathrm{E}[\varepsilon|X] = 0$. Their target parameter is the regression coefficient $\theta_0$, which they interpret as a measure of the negative labor supply elasticity.
In this example, $g_0(\cdot)$ is a nuisance function. We are not primarily interested in how task features, outside of reward, relate to the outcome. However, dube2020 emphasize that task heterogeneity needs to be accounted for to meaningfully interpret $\theta_0$.
If $g_0(\cdot)$ were known, estimation of $\theta_0$ could proceed by regressing $Y - g_0(X)$ onto $D$, which is equivalent to estimation based on the moment condition
where $W = (Y,D,X)$ denote observed random variables. Of course, $g_0(\cdot)$ will typically not be known. A common strategy to simplify the problem is to assume that $g_0(X) = X'\beta_0$ with unknown coefficient $\beta_0$, in which case the model reduces to the familiar multiple regression model. When the dimension of $X$ is much smaller than the sample size, estimation would then proceed by ordinary least squares (OLS) of $Y$ on $D$ and $X$.
Even in the simple multiple regression model, we have a nuisance parameter, $\beta_0$. Conventional regression estimates it jointly with the target parameter $\theta_0$. Alternatively, one can partial $X$ out from both $D$ and $Y$ to obtain residuals that isolate the variation identifying $\theta_0$. By the Frisch-Waugh-Lovell Theorem, regressing these residuals on each other yields an estimate of $\theta_0$ that is numerically equivalent to that obtained from regressing $Y$ on $D$ and $X$. This notion of “partialling out” is related to a broader principle for handling nuisance parameters that we emphasize throughout the review.
In the actual dube2020 example, some of the task characteristics are captured as text data, which makes it difficult to justify ad hoc parametric assumptions like those in the linear model. Instead, dube2020 allow $g_0(\cdot)$ to be a flexible, high-dimensional function rather than committing to a low-dimensional functional form. In contrast to the low-dimensional linear case, where $g_0(X)=X'\beta_0$ reduces to estimating a small set of coefficients, estimation of $g_0(\cdot)$ must then also accommodate rich nonlinear variation.
This regression setup illustrates a broader template common in empirical research, in which the target parameter $\theta_0$ is defined as the solution to a moment condition:
Here, $m$ is a score (or moment) function and $W$ again denotes observed random variables. The parameter $\eta_0$ denotes a nuisance object, which is not of direct interest but is used to define $\theta_0$. $\eta_0$ is often high-dimensional. For example, $\eta_0$ will represent a vector of conditional expectation functions in many interesting cases. This semi-parametric structure encompasses the regression example above, where the target parameter is the coefficient on $D$. It more generally applies to the estimation of many other canonical parameters, including average treatment effects, parameters in linear instrumental variables models, local average treatment effects, dynamic treatment effects in staggered adoption designs, and parameters in nonlinear structural models.
High-dimensional nuisance parameters can arise in several ways: (i) when the nuisance function depends on only a few covariates, controls, or instruments, but no parametric model is specified; (ii) when there are many such variables, even under parametric assumptions; or (iii) when numerous variables enter through unknown functions. High-dimensionality is increasingly common in applications using text or image data gentzkow2019b, but it can also arise in simpler settings. Even a single continuous covariate may create a high-dimensional problem---for example, when identification relies on rainfall instruments that can be nonlinearly related to an endogenous regressor hidalgo2010a,gilchrist2016something,dustmann2017.
As rich data become more common in applied research, there is growing appreciation that traditional functional form assumptions are often difficult to justify, motivating the use of more flexible tools for estimating nuisance parameters. As a result, there is increasing interest in using machine learning (ML) methods, which provide flexible tools for estimating high-dimensional nuisance parameters. A natural use of ML would then be to obtain nuisance parameter estimates, $\widehat\eta$, and use these in place of $\eta_0$ in (ref). Specifically, we can define an estimator $\widehat\theta$ of $\theta_0$ as a solution to the sample analog of (ref): $$ \widehat\theta:\: \frac{1}{n} \sum_{i=1}^{n} m(W_i;\widehat\theta,\widehat\eta) = 0 $$ where $W_i$ denote observed variables for observations $i = 1,...,n$. However, the resulting “plug-in” estimator $\widehat\theta$---so-called because it “plugs” $\widehat\eta$ in for $\eta_0$---can behave poorly and lead to misleading conclusions due to errors in estimating $\eta_0$ propagating into $\widehat\theta$.
In settings with high-dimensional nuisance parameters, standard asymptotic approximations may fail due to two distinct forms of sensitivity to nuisance parameter estimation, termed regularization bias and overfitting bias. Both terms describe channels through which using an estimated nuisance parameter, $\widehat{\eta}$, instead of the true but unknown $\eta_0$, can distort the behavior of the plug-in estimator $\widehat\theta$. A leading manifestation of these distortions is bias, as the terminology suggests, but they more broadly invalidate conventional inference methods that fail to account for nuisance estimation.
Roughly speaking, regularization bias refers to the direct impact of estimation error in $\widehat\eta$ on the plug-in estimator $\widehat\theta$ that results from the difference between $m(W;\theta_0,\widehat\eta)$ and $m(W;\theta_0,\eta_0)$. Overfitting bias refers to a more subtle issue. Because $\widehat\eta$ is an estimator, it is itself a random function of the data. $\widehat\eta$ is thus generally correlated with the observations $\{W_i\}_{i=1}^{n}$ also used in the estimating equation $\frac{1}{n}\sum_{i=1}^{n} m(W_i;\theta,\widehat\eta)$. When this dependence is strong, for example due to “overfitting”, it may generate large differences between $\frac{1}{n}\sum_{i=1}^{n} m(W_i;\theta,\widehat\eta)$ and $\frac{1}{n}\sum_{i=1}^{n} m(W_i;\theta,\eta_0)$, which results in poor performance of $\widehat\theta$.
Both regularization and overfitting biases are major concerns in high-dimensional contexts. Accurately estimating high-dimensional $\eta_0$ is inherently difficult and often results in non-negligible errors. Further, estimation in high-dimensional settings typically relies on highly data-adaptive procedures---such as modern ML methods---which amplify the risk of overfitting bias. As such, mitigating these biases is the focus of a large and rapidly growing literature in statistics and econometrics that builds from classic ideas in semiparametric and nonparametric estimation. One method that provides a solution in a wide variety of empirical settings, and that is the topic of this review, is double/debiased machine learning chernozhukov2018.
DML provides a blueprint for alleviating both regularization and overfitting bias. At its core, DML combines two classical ideas from the rich literature on semiparametric inference---using Neyman orthogonal scores\footnote{Such scores are also referred to as orthogonal scores, orthogonal moments, locally robust moments, debiased moments, influence functions, and pathwise derivatives. We follow chernozhukov2018 and use the term “Neyman orthogonal scores” in homage to Neyman's early contributions, e.g., NeymanCalpha,Neyman1979.} to alleviate regularization bias and using cross-fitting to alleviate overfitting bias---in a common methodological framework. Neyman orthogonality ensures that plugging in estimates that are close to, but not exactly equal to, $\eta_0$ does not lead to large changes in the moment condition (ref).\footnote{ Neyman orthogonality is not guaranteed for all scores that serve to identify a parameter of interest. We show how to construct a Neyman orthogonal score from a given score in Appendix (ref).} Cross-fitting, a form of sample splitting, alleviates potential dependence between nuisance estimates $\widehat{\eta}$ and parts of the data used for estimating the target parameter. Used in conjunction, DML's two core components significantly reduce the impact of nuisance estimation on estimates of the target parameter. However, high-quality estimation of the target parameter $\theta_0$ still requires nuisance parameters to be estimated sufficiently well. Consequently, theoretical results for DML assume specific convergence rate conditions on nuisance estimators. Many estimators, including ML methods, can satisfy these conditions.
From a practical perspective, DML enables researchers to leverage a wide range of ML tools, making it particularly valuable in complex data scenarios involving numerous variables, images, or text data. Importantly, the benefits of ML also extend to traditional research settings with fewer covariates and conventional tabular data as they remove the need for researchers to commit beforehand to specific parametric (often linear) models. Furthermore, DML is easy to implement, applicable to a wide range of econometric settings, and readily available in existing software packages, including Stata, R, and Python DoubleML2021R,DoubleML2022Python,Ahrens2023_ddml,ddml_R. DML thus has the potential to enhance the credibility of research findings in a broad spectrum of settings, either when used as a complementary robustness check or when the application necessitates the use of flexible estimation methods for nuisance objects.
While DML offers a framework for combining flexible nuisance estimation with valid asymptotic inference, its implementation raises important challenges. Available theoretical results assume the nuisance functions are estimated with sufficiently high accuracy. Achieving these convergence rates for modern ML methods often demands strong assumptions and special tuning, and they may not hold for off-the-shelf algorithms. These theoretical qualifications manifest in practical problems where empirical results depend on implementation choices, such as selecting and tuning an ML method for nuisance estimation. Assessing the quality of nuisance estimators is often difficult and, in some applications, different seemingly reasonable choices can lead to substantively different conclusions. This paper therefore aims not only to motivate and explain DML, but also to guide its application in empirical research, emphasizing the need for careful diagnostic analysis and robustness checks. To this end, we divide our review into two parts.
First, in Sections (ref) and (ref), we introduce the DML blueprint at a high level. Section (ref) discusses the practical implications of nuisance estimation and the role of DML's two key components in their remedy. Section (ref) summarizes the asymptotic properties of DML and contains algorithmic details on implementation of generic DML estimators.
In the second part of the paper, we turn to simulations and empirical applications, found in Section (ref) through Section (ref). These examples illustrate DML and provide discussion of key issues that arise in its implementation.
In Section (ref), we present results from two simulation examples. The first is a simple linear IV example that demonstrates the importance of cross-fitting. The second compares the benefits of the DML average treatment effect estimator with inverse propensity weighted and regression-based estimators.
In Section (ref), we illustrate how DML can be leveraged to reduce the dependence on functional forms in staggered adoption designs with covariates. We revisit the analysis of dobkin2018economic, who study the economic consequences of hospital admission. We estimate group-time average treatment effects on the treated under a conditional parallel trends assumption, and show how DML inference applies to dynamic average treatment effects. We note that cross-fitting introduces an additional source of randomness induced by sample-splitting. Part of our aim in this example is to illustrate a simple approach to aid in gauging the impact of this source of randomness.
In Section (ref), we apply DML to estimation of regression coefficients in the presence of complex covariates. We specifically revisit dube2020 who apply DML to estimate the labor supply elasticity in online labor markets using textual controls. We have two main goals in this section: to illustrate the use of complex non-tabular data and, more importantly, to illustrate that DML estimates can vary substantially across otherwise reasonable choices of machine learners (i.e., algorithms used to estimate nuisance functions). This sensitivity can lead to qualitatively different conclusions about economically meaningful parameters. Because it seems difficult to know ex ante exactly which learner one should choose in these situations, we use this example to discuss robustness checks and suggest strategies for selecting ML algorithms.
Section (ref) concludes by summarizing takeaways, raising some caveats, and pointing to potential directions for further research.
To complement the present article, we provide additional resources on our regularly updated website \href{https://dmlguide.github.io}{dmlguide.github.io}. The materials include replication files, additional examples with code, references to DML software packages, and links to other resources.
In terms of scope, we emphasize that DML---or any other estimation framework---cannot replace careful reasoning about economic parameters and identifying assumptions. Rather, with a well-defined target parameter and corresponding identifying assumptions, DML can aid in obtaining estimates of the target parameter in the presence of complex data structures and without relying on pre-specified functional form assumptions. In other words, DML is useful only after a target parameter is defined and the assumptions linking observed data to that parameter are well understood.\footnote{This treatment parallels heckman2007econometric, who stress that estimation plays a limited role relative to defining a target parameter and articulating the assumptions that connect it to the data. } With this understanding, we discuss identification assumptions only with the aim of illustrating the economic content of the applications. Further, we will not review specific ML methods. varian2014big, Sendhil:MLAppliedMetrics, athey2019_ml, and dell2024deep provide reviews of ML methods targeted at economists. Hastie2009 and islp are classic textbook treatments of popular ML methods.
\noindentLiterature. Inference about low-dimensional target parameters in the presence of high-dimensional or nonparametric nuisance components has a long history in econometrics and statistics. Classic reviews such as newey:handbook, yatchew1998nonparametric, li:racine:book, chen:handbook76, and ichtodd:handbook74 synthesize early work on semiparametric and nonparametric methods, emphasizing how nuisance parameters can be accommodated without fully specifying the data-generating process. More recent surveys shift the focus toward the use of modern machine learning tools for nuisance estimation and their implications for inference; see, e.g., chernozhukov2018, diaz2020machine, hines2022demystifying, and kennedy:review. Our review complements these contributions by emphasizing the practical consequences of nuisance estimation choices for applied empirical work.
DML combines two ideas with deep roots in the semiparametric inference literature: Neyman orthogonal scores and sample splitting. NeymanCalpha introduced orthogonal scores in the context of efficient parametric hypothesis testing. Orthogonal scores later played a central role in the development of modern semiparametric estimation, especially in settings with high-dimensional or nonparametric nuisance parameters. Key contributions developing these ideas include vaart:1991, andrews94, and newey94, with a comprehensive textbook treatment provided by vdV.
Sample splitting has also played a long-standing role in semiparametric inference. It appears in several early contributions to semiparametric estimation; see, for example, hasminskii1978, bickel:1982, pfanzagl82book, schick1986asymptotically, and bickelritov1988. In economics, sample splitting has long been used in instrumental variable estimation to mitigate bias from many instruments. See, for instance, ssiv and AIK:JIVE for foundational work and NeweyEtAl-JIVE, hansen:kobzur, and chyn2024examiner for current developments.
More recently, sample splitting and variations such as cross-fitting have gained renewed attention in high-dimensional contexts. A growing literature shows how these techniques can mitigate problems introduced by overfitting and improve inference when modern ML methods are used for nuisance estimation. See, for example, robins2008higher, belloni2010lasso, Belloni2012, fan2012variance, robins2013new, vdl:AdaptiveTarget, robins2017minimax, wager2018, athey2019e, and athey2021a.
DML is also related to targeted maximum likelihood (or minimum loss) estimation, which was introduced in SRR1999-rejoinder for treatment effect estimation and generalized by van2006targeted; see also vanderlaan:book. zheng2011cross discuss benefits of sample splitting for targeted maximum likelihood learning. diaz2020machine expands on the difference between DML and targeted maximum likelihood estimation.
Our focus in this review is on the practical implications of nuisance estimation and the core ideas that motivate DML. Accordingly, we do not attempt a comprehensive survey of the rapidly expanding literature that extends DML in a wide range of directions. We nevertheless highlight several representative strands of this work.
A number of papers adapt DML to canonical empirical settings, including panel data and difference-in-differences designs chang2020double,chiang2022multiway, klosin:2023cont, abadie:matrix, HHZ:DMLcontDID, clarkeDoubleMachineLearning2024,chiang2025double. Related contributions study the use of DML in instrumental variables and proxy control settings DMLDensity, deaner2021proxy, singh2024double. A growing body of research also examines treatment effect and policy parameters, including incremental and dynamic treatment effects, nonparametric policy learning, and localized estimands that depend on complex nuisance components BMBK:incremental, Syrgkanis:DTE, klosin:admlcont, nie2021quasi, semenova2021debiased, DL:DMLnp, FS:orthogonallearning, kennedy2023, sasaki:DMLpolicy, kallus2024localized. Other extensions address specific econometric complications. These include settings with generated regressors escanciano:LRgenerated, partial or set identification semenova2023debiased, and sample selection Bia:SampleSelectDML. Complementing these application-driven contributions, several papers develop general frameworks for the automatic construction of Neyman orthogonal moments for broad classes of target parameters chernozhukov2021automatic, FLM:het, chernozhukov2022locally, chernozhukov2022automaticECMA, CNS:globalandlocal.
More broadly, DML allows researchers to avoid auxiliary parametric assumptions. While these assumptions simplify estimation, they are seldom motivated by economics and can be detrimental for applications that aim to estimate causal parameters. This perspective connects DML to the recent literature highlighting that statistically convenient estimands, often based on linear models, may fail to even approximate causal effects. Such failures have been documented in difference-in-differences settings de2020two, goodman2021difference, sun2021estimating, callaway:santanna, baker2022much, Chaisemartin:review, roth2023s, borusyak2024revisiting. Related concerns arise in linear regression with multiple treatments goldsmith-pinkham2022 and in instrumental variables settings blandhol2022tsls. By allowing for flexible estimation of nuisance parameters, DML provides a framework for inference that obviates the need for convenient but potentially detrimental parametric assumptions.
This section describes the two essential components that define DML: Neyman orthogonality and cross-fitting. Together, they help control the sensitivity of the target estimator to nuisance estimation, which can substantially improve both the reliability of point estimates and the quality of conventional asymptotic approximations. By alleviating this sensitivity, DML further opens the door for researchers to use a wide range of flexible estimators, including many modern machine learners, for estimating nuisance parameters.
We frame our discussion of DML within a relatively general semiparametric framework. There are two key elements of the framework. First, we have a target parameter of interest, $\theta_0$, which is low-dimensional; e.g., $\theta_0$ may be an average treatment effect or a fixed vector of regression coefficients. Second, we have a nuisance parameter $\eta_0$ which may be high-dimensional and potentially complex. In many examples, $\eta_0$ is a vector of conditional expectation functions, such as outcome regressions and propensity scores, though it may take other forms.
Throughout this review, we focus on the case where we observe an i.i.d.\ sample $\{W_i : i=1,\ldots,n\}$ from a random vector $W$. Each $W_i$ collects the variables relevant for individual $i$. For example, $W_i$ might include an outcome $Y_i$, a treatment variable $D_i$, a vector of controls $X_i$, and excluded instruments $Z_i$. This structure also extends to cross-sectional and panel settings with arbitrary temporal dependence and fixed $T$.\footnote{To account for cluster dependence, one may simply redefine $W_i$ to include the data of the $i$th individual over multiple time periods---e.g., $W_i = (Y_{i,t}, D_{i,t}, X_{i,t})_{t=1}^T$.}
We assume the target parameter is identified by moment conditions
where $m(\cdot; \theta, \eta)$ is a known score function indexed by $\theta$ and nuisance parameter $\eta$ with true values $\theta_0$ and $\eta_0$. We focus exclusively on the case where the score function $m(\cdot; \theta, \eta)$ defines as many constraints as we have parameters of interest, but note that the framework extends to other settings such as GMM as discussed, e.g., in chernozhukov2018. Throughout, we assume that the target parameter is strongly identified in the sense that (ref) has a unique solution and satisfies regularity conditions such that $\sqrt{n}$-consistent and asymptotically normal inference for $\theta_0$ would be achievable if $\eta_0$ were known.\footnote{Extension to weakly identified settings is possible as in, e.g., chernozhukov2015 and ma:hdwiv.}
This semiparametric framework captures a large range of common parameters of interest in empirical research. We discuss four illustrative examples below.
Example 1. Linear Regression Coefficient. Consider linear regression with a single variable of interest $D$ and a $p \times 1$ vector of controls $X$ that may include a constant:
where $0_p$ denotes a $p \times 1$ vector of zeros. The coefficient $\theta_0$ on $D$ is the target parameter. The vector of coefficients $\beta_0$ on controls $X$ is the nuisance parameter.
The traditional textbook approach is to estimate both $\theta_0$ and $\beta_0$ by applying OLS to equation (ref). This problem can be framed as a semiparametric estimation task by explicitly targeting $\theta_0$ separately from the nuisance parameters. In the linear regression example, this corresponds to another textbook approach: partialling out.
For any random variable $A$, let $\eta_{A,0} = \arg\min_{\eta} \mathrm{E}[(A - X'\eta)^2]$ denote the best linear predictor coefficient of $A$ given $X$. This definition implies the orthogonality condition
A valid score for $\theta_0$ is then
where the nuisance parameter is $\eta = (\eta_Y',\eta_D')'$ with true value $\eta_0 = (\eta_{Y,0}',\eta_{D,0}')'$. By (ref) and the orthogonality condition, (ref), $\mathrm{E}[m_{LM}(W;\theta_0,\eta_0)] = 0$.
Equation (ref) is the population moment condition underlying the partialling out interpretation of linear least squares regression. It corresponds to projecting $Y$ and $D$ onto $X$ and then estimating $\theta_0$ from a regression using the resulting residuals.
By Frisch-Waugh-Lovell, this score yields the same estimator as OLS of $Y$ on $(D,X)$ in the low-dimensional linear setting. The value of writing the problem in this way is therefore primarily conceptual. It makes explicit the construction of a score for the target parameter $\theta_0$ by projecting onto covariates $X$ and appropriate partialling out. This familiar approach generalizes and lies at the core of the key “Neyman orthogonality” property---to be discussed in Section (ref)---that is fundamental for DML. We verify Neyman orthogonality of the score (ref) in Appendix (ref), and discuss a more general construction of Neyman orthogonal scores via “partialling out” in Appendix (ref).
$\qed$
\noindentExample 2. Partially Linear Regression Coefficient. As discussed in the Introduction, partially linear regression (PLR),
is a natural, flexible generalization of multiple linear regression.\footnote{Analogous to linear regression, PLR can be motivated in similar manner as the best partially linear approximation to the conditional expectation function. Arguments for economic interest in the best “fully” linear approximation to the conditional expectation function as outlined, e.g., in angrist2009mostly, also make PLR an attractive baseline choice in many economic analyses.}
There are several moment conditions for identifying the PLR coefficient $\theta_0$. Two leading examples are based on the score functions
where the nuisance parameters are $\eta(X)=g(X)$ in (ref) with true value $g_0(X)$, and $\eta(X)=(\ell(X),r(X))$ in (ref) with true values $\ell_0(X)=\mathrm{E}[Y|X]$ and $r_0(X)=\mathrm{E}[D|X]$.
The first score is equivalent to regressing $Y-g(\cdot)$ against $D$. The second score corresponds to a “partialling out” approach where both $Y$ and $D$ are residualized with respect to $X$ before regressing the residuals on each other. The latter mirrors the Frisch–Waugh–Lovell logic from linear regression but now allows $X$ to enter flexibly. Note that (ref) corresponds to the treatment of the partially linear model in robinson. While both scores identify the target parameter $\theta_0$, only $m_{PLM}$ satisfies the key “Neyman orthogonality” property that is fundamental for DML. We verify this in Section (ref). $\qed$
\noindentExample 3. Linear IV Coefficient. In linear instrumental variable (IV) models, it is often unclear how best to use instruments. For instance, when rainfall or weather variables are employed as instruments, researchers face choices such as whether to use rainfall in levels, logarithms, squared terms, or deviations from historical averages hidalgo2010a,gilchrist2016something,dustmann2017.
Formally, consider the linear structural equation $Y = \theta_0 D + \varepsilon$ where $D$ is an endogenous variable, $Z$ is a vector of excluded instruments, $\mathrm{E}[\varepsilon|Z]=0$ holds, and we abstract from other covariates. A natural score function is then
where the true value of the nuisance function is $\eta_0(Z) = \mathrm{E}[D|Z]$, which corresponds to the optimal instrument under homoskedasticity. $\mathrm{E}[m_{IV}(W;\theta_0,\eta_0)] = 0$ then follows immediately from the exclusion restriction $\mathrm{E}[\varepsilon|Z] = 0$.\footnote{Under mean independence, any function $g(Z)$ serves as a valid instrument in the sense of satisfying the moment condition $\mathrm{E}[(Y-\theta_0 D)g(Z)] = 0$. However, the instrument relevance condition requires that $\mathrm{E}[g(Z) D] \neq 0$. $\mathrm{E}[g(Z) D]$ is also tightly tied to the efficiency of the IV estimator. Under homoskedasticity, the choice of $g(Z) = \eta_0(Z) = \mathrm{E}[D|Z]$ produces an asymptotically efficient estimator of $\theta_0$.}
The IV score satisfies the key “Neyman orthogonality” property that is fundamental for DML. We verify this in Appendix (ref).$ $$\qed$
\noindentExample 4. Average Treatment Effect. A central policy-relevant parameter is the average treatment effect (ATE) of a binary treatment $D$ on an outcome $Y$ defined as
where $Y(d)$ is the potential outcome under treatment status $d\in\{0,1\}$.
In non-experimental settings, identification of the ATE relies on two standard conditions: overlap and unconfoundedness ImbensRubin2015. Overlap requires that the probability of treatment is bounded away from 0 and 1 across all covariate values: $0 < \text{Pr}(D = 1|X=x) < 1$ for all $x$. That is, we should see treatment and control observations at all values of $X$. Unconfoundedness requires that the treatment status is independent of potential outcomes after conditioning on the covariates: $(Y(1),Y(0)) \perp D|X.$ That is, treatment is as good as randomly assigned after conditioning on $X$.
Under these assumptions, the ATE can be identified using moment conditions. We consider two commonly used scores, the inverse propensity weighted (IPW) score and the augmented IPW (AIPW) score newey94,rrz:
where $W = (Y, D, X)$. The true value of the nuisance parameters are $\alpha_0(D, X) = \frac{D}{r_0(X)} - \frac{(1-D)}{1 - r_0(X)} $, $r_0(X) = \mathrm{E}[D\vert X]$, and $\ell_0(D,X) = \mathrm{E}[Y\vert D, X]$. Under overlap and unconfoundedness, it can be shown that both $\mathrm{E}[m_{IPW}(W;\theta_0,\eta_0)] = 0$ and $\mathrm{E}[m_{AIPW}(W;\theta_0,\eta_0)] = 0$.
Importantly, only the AIPW score---also referred to as the “doubly robust” score---satisfies the key “Neyman orthogonality” property that is fundamental for DML. The IPW score is not Neyman orthogonal and should not be used together with generic machine learners. We verify Neyman orthogonality of the AIPW score in Appendix (ref), and illustrate empirical consequences of (non-)orthogonality in Section (ref).$ $$\qed$
Identification of the target parameter in each of these examples depends on nuisance parameters. Outside of special cases such as randomized controlled trials where the propensity score $\mathrm{E}[D|X] = r_0(X)$ is known by design, these nuisance parameters are generally unknown and thus need to be estimated.
In general, many moment conditions will exist for any target parameter. Examples 2 and 4 illustrate this by presenting two different moments that each identify the parameter of interest and could, in principle, be used for estimation. However, in both cases, only one of the proposed moment functions satisfies a key condition---Neyman orthogonality---that is crucial for obtaining reliable estimates in the presence of nuisance parameters.
In the next subsection, we discuss statistical issues stemming from the estimation of these nuisance parameters. Then, in Section (ref), we outline how the combination of DML's two essential components alleviates the impact of nuisance estimation on the main inferential target. Fundamentally, it is this reduction of impact that allows DML to accommodate complex estimators, including ML methods, for nuisance estimation.
Suppose that we have at our disposal a first-step estimator, $\widehat\eta$, for the nuisance parameter $\eta_0$. This might be a parametric estimator such linear regression or a flexible, nonparametric learner from the ML toolbox. A plug-in estimator for the parameter of interest can be constructed as the solution to the sample average of the scores:
Inference about $\theta_0$ follows from a standard asymptotic Taylor expansion of (ref) around the true parameters $(\theta_0, \eta_0)$:\footnote{We use a finite-dimensional expansion here to convey intuition. A formal treatment would require additional technicality to deal with cases where $\eta$ is a high- or infinite-dimensional object such as a function.}
where “higher order terms” capture the impact of squared estimation errors $(\widehat{\eta} - \eta_0)^2$ and other remainders from the linearization. Under standard regularity conditions, the term labeled CLT in (ref) will be approximately normal by a central limit theorem. The focus of DML is addressing the term ($\star$).
The term ($\star$) captures the first-order impact of estimating the nuisance parameter $\eta_0$. Its presence suggests the asymptotic distribution of $\widehat{\theta}$ will generally depend on the asymptotic behavior of the nuisance estimator $\widehat{\eta}$. That is, estimation error in $\widehat{\eta}$ propagates directly into inference about $\theta_0$, making the resulting distribution of $\widehat\theta$ differ from the idealized case in which $\eta_0$ is known. In situations where the nuisance parameter is low-dimensional (e.g., when assuming a linear model with few parameters), this additional uncertainty can be adequately characterized and managed through adjustments to the asymptotic variance; see, e.g., Section 6 of newey:handbook.
However, the first-order dependence of $\widehat\theta$ on $\widehat\eta$ poses substantial complications in settings where the nuisance parameter is high-dimensional and $\widehat{\eta}$ corresponds to a flexible estimator. The complication results because flexible estimators are often associated with non-negligible bias and variance. As a consequence, (ref) may be dominated by the first-order term ($\star$) that involves $\widehat\eta$. In general, ($\star$) diverges due to two issues, referred to as regularization bias and overfitting bias.
Regularization bias refers to the fact that neither term inside the sum in $(\star)$---$\frac{\partial}{\partial \eta} m(W_i; \theta_0, \eta_0)$ and $(\widehat{\eta} - \eta_0)$---is mean zero in general. As a result, $(\star)$ is $\sqrt{n}$ times a sample average of a non-mean zero quantity, which does not converge in general. The term “regularization bias” reflects the fact that high-dimensional or nonparametric methods used to estimate $\eta_0$ often rely on regularization. That is, they introduce bias in order to control variance, implying that $(\widehat{\eta} - \eta_0)$ will not generally be mean zero in finite samples. Importantly, researchers have some control over score functions. The first key ingredient of DML---Neyman orthogonality, discussed in more detail in Section (ref)---is exactly the requirement that estimation be based on scores for which the high-dimensional analog of $\frac{\partial}{\partial \eta} m(W_i; \theta_0, \eta_0)$ is mean zero.
Overfitting bias, also referred to as own-observation bias, arises more subtly. Because $\widehat{\eta}$ is a function of the data used in its estimation, $\widehat{\eta}-\eta_0$ generally depends on the observations $W_i$ that are also used to construct the sample moment condition. This dependence typically occurs when the same dataset is used both to obtain $\widehat{\eta}$ and to evaluate the sample analog of equation (ref), though it can arise more generally in dependent data settings. As a result, the product $\frac{\partial}{\partial \eta} m(W_i; \theta_0, \eta_0)(\widehat{\eta} - \eta_0)$ in ($\star$) will generally not be mean zero, even if either term in the product is mean zero when considered in isolation.
We use the term “overfitting bias” to describe failures in conventional inference about the target parameter that arise from statistical dependence between $(\widehat{\eta} - \eta_0)$ and the observations used in the sample moment condition. We recognize that “overfitting” is typically associated with inflated variance in prediction contexts. Our usage emphasizes that overfitting in nuisance estimation can lead to bias in target parameter estimation, because it inflates the dependence between nuisance estimation error and the data.
The second key ingredient of DML---cross-fitting, also discussed in Section (ref)---addresses overfitting bias by using sample splitting to ensure (approximate) independence between the estimation error in the nuisance function and the observations used in the sample moment condition. This independence implies that the product in the numerator of ($\star$) is mean zero whenever either term is mean zero.
The two main ingredients of DML, discussed in Section (ref), are meant to alleviate the first-order impact of nuisance estimation. Estimation of nuisance parameters also generally affects the target parameter through the “higher order terms” in (ref). A sufficient condition for these terms to be asymptotically ignorable is that the nuisance parameters are estimated accurately enough such that $\sqrt{n} \|\widehat{\eta} - \eta_{0}\|^2 \rightarrow_p 0$ under a suitable norm. A mean-square convergence rate faster than $n^{-1/4}$ is a commonly cited sufficient benchmark.\footnote{Because the treatment of higher-order terms in DML is similar to that in other semiparametric approaches, we do not discuss them further. For more detailed discussion, see, e.g., chernozhukov2018, chernozhukov2022locally, and kennedy:review.}
Establishing such estimation quality guarantees for modern machine learning methods is an active area of research. Available results involve a combination of assumptions on the structure of the underlying data generating process and a choice of estimator that successfully leverages that structure. A canonical example is the lasso, which achieves suitable rates when the true regression function is sparse bickel2009simultaneous,Belloni2012. Related results for other classes of learners similarly require strong restrictions, such as smoothness, low effective dimension, or specific forms of regularization; see, for example, results for neural networks under compositional or smoothness assumptions farrell2021,schmidt2020nonparametric and for random forests under honesty and regularity conditions wager2018,Athey2019_grf.
At the same time, recent work highlights important limitations of these results. In particular, many theoretical guarantees do not apply to off-the-shelf implementations with default tuning choices. For tree-based methods, available analyses show that, absent specific tree depth and subsampling choices, pointwise polynomial convergence rates may fail or that estimators may even be pointwise inconsistent, with only slow $L^2$ consistency established in high-dimensional settings CVFL:hdrandomforests,cattaneo:trees. More broadly, if the nuisance function $\eta_0$ is fully nonparametric and high-dimensional without additional simplifying structure, no known method is able to attain the $n^{-1/4}$ rate required for standard DML asymptotics.
In applications, it is often unclear which assumptions are credible and which estimator is appropriate, making it difficult to judge whether the required convergence conditions plausibly hold. We thus return to the topic of the choice of ML estimator in Section (ref).
To reduce the dependence of estimators for $\theta_0$ on the estimation of high-dimensional nuisances $\eta_0$, DML estimators rely on two essential components---estimation based on Neyman orthogonal scores and cross-fitting---that mitigate regularization and overfitting bias. We now discuss each in turn.
The chief difficulty in using the plug-in estimator $\widehat\theta$ is its dependence on the nuisance parameter estimator $\widehat\eta$. Neyman orthogonality is a local robustness property of the score function in (ref) that decreases sensitivity of $\widehat\theta$ to errors in estimating $\eta_0$.
Formally, in addition to identifying the target parameter by satisfying (ref), a Neyman orthogonal score also satisfies
where $\eta$ denotes a candidate value for the nuisance parameter and $\lambda$ indexes the size of a local deviation away from the true $\eta_0$. Intuitively, (ref) requires that, when we make a small move away from $\eta_0$ in any direction $\eta-\eta_0$, the moment condition remains unchanged.\footnote{This formulation accommodates cases where $\eta_0$ is not finite-dimensional but instead a function or other complex object belonging to an abstract space $\mathcal{T}$.} In other words, small perturbations of the nuisance parameter away from the true value do not create first-order changes in the moment function. As a result, the plug-in estimator of $\theta_0$ is less sensitive to estimation error in $\widehat\eta$, thereby reducing regularization bias.
Returning to the expansion in (ref), consider the case where the estimator $\widehat{\theta}$ is based on a Neyman orthogonal score. If the first-step estimation error $(\widehat{\eta} - \eta_0)$ were independent of the sample $\{W_i\}_{i=1}^{n}$ used to construct the score, Neyman orthogonality would imply that the second term in (ref) vanishes,
under the convergence requirements on $\widehat\eta$ cited above and additional mild regularity conditions. Loosely, condition (ref) means that $\frac{1}{\sqrt{n}}\sum_{i=1}^n\frac{\partial}{\partial \eta} m(W_i; \theta_0, \eta_0)$ behaves like a mean 0 normalized sum, which converges by a central limit theorem. As long as $(\widehat{\eta} - \eta_0)$ converges to 0, the whole term disappears. Maintaining the assumption that higher order terms vanish, it follows that the asymptotic distribution of the plug-in estimator $\widehat\theta$ does not depend on the nuisance estimator $\widehat\eta$. Thus, inference about $\theta_0$ using $\widehat{\theta}$ can proceed as if $\eta_0$ were known. This heuristic derivation can be formalized and underlies the crucial importance of Neyman orthogonality in settings with complex nuisance parameters.
An important practical implication is that estimators based on Neyman orthogonal scores yield inference about $\theta_0$ that does not depend on the detailed statistical properties of the nuisance estimator $\widehat\eta$. This robustness is useful even in classical low-dimensional settings, where it avoids cumbersome variance adjustments or computationally costly resampling approaches to account for estimation of $\widehat\eta$. It becomes particularly important when modern flexible methods are used, as the statistical properties of these methods are still under active development. For example, only coarse rates of convergence are currently available for many promising machine learning methods. By alleviating the first-order impact of nuisance estimation, Neyman orthogonality makes it possible to combine such methods with standard asymptotic approximations, enabling formally valid inference for $\theta_0$ while exploiting modern flexible estimators for nuisance parameters.
We illustrate the importance of using Neyman orthogonal scores for obtaining reliable finite-sample inference in Section (ref). There, we present simulation results demonstrating that estimators of the average treatment effect (ATE) based on non-orthogonal scores may exhibit substantial bias, and the associated confidence intervals are often severely distorted. In contrast, DML estimators, which incorporate Neyman orthogonal scores as a core component, are approximately unbiased, and their confidence intervals achieve coverage rates close to the nominal level.
\paragraph{Neyman orthogonal scores for common parameters.} Neyman orthogonal scores are well-known and readily available for common target parameters. We present Neyman orthogonal scores for six illustrative targets in Table (ref). In Panels (i)-(iv), we present Neyman orthogonal scores corresponding to the target parameters in Examples 1-4 from Section (ref). In Panel (v), we present an orthogonal score for the coefficient $\theta_0$ in the partially linear IV model
where $Z$ is an excluded scalar instrumental variable. This model differs from the simple IV model considered in Example 2 in two key respects: We allow for the presence of controls, and we do not impose full mean independence of the instrument from structural unobservables. Finally, in Panel (vi), we consider a more exotic target parameter---an average derivative corresponding to a continuous variable of interest. We present this case both because continuous variables are practically relevant and, more importantly, to highlight that nuisance parameters are not always conditional expectation functions (or projection coefficients as in the linear regression example).
\paragraph{Deriving and verifying Neyman orthogonal scores.} Given the importance of Neyman orthogonal scores and their relevance for DML, we provide an outline of a general structure for obtaining Neyman orthogonal scores in Appendix (ref).
For a given score, one can generally verify Neyman orthogonality, or a lack thereof, by direct application of its definition in (ref). We illustrate such a derivation for partially linear regression (Example 2) below. Similar derivations for the scores presented in Examples 1, 3, and 4 are provided in Appendix (ref).
\noindentExample 2 (continued). Neyman Orthogonality of Partially Linear Regression Scores. We introduced two score functions for identifying the partially linear regression coefficient $\theta_0$, (ref) and (ref). The first does not satisfy Neyman orthogonality, while the second---which corresponds to flexible partialling out---does.
Consider (ref), which has nuisance parameter $g(X)$. Let $\Delta g(X) = g(X)-g_0(X)$; then
which is generally non-zero when $D$ and $X$ are related. Hence, the score is not orthogonal.
Turning to (ref), we have nuisance parameters $\eta(X) = (\ell(X),r(X))$ with true value $\eta_0(X) = (\ell_0(X) = \mathrm{E}[Y|X],r_0(X) = \mathrm{E}[D|X])$. Letting $\Delta\eta(X) = \eta(X) - \eta_0(X) = (\Delta \ell(X) = \ell(X) - \ell_0(X),\Delta r(X) = r(X)-r_0(X))$, we have
where the last equality follows from $\ell_0(X) = \mathrm{E}[Y|X]$ and $r_0(X) = \mathrm{E}[D|X].$ $ $$\qed$
The intuition for Neyman orthogonality in Example 2 is instructive and corresponds to common intuition provided for partialling out. The orthogonal score uses only variation in $D$ and $Y$ that is (mean) independent of $X$, thereby isolating the identifying variation. As a result, small errors in one nuisance function can be offset by the other, making estimation more robust. In contrast, the non-orthogonal score uses all the variation in $D$ but adjusts only for the effect of controls on $Y$. Any mistakes in estimating $g(X)$ that are correlated with $D$ then directly bias the estimate of $\theta_0$, much like in a classic omitted variable scenario. Appendix (ref) shows that the same partialling out intuition generalizes to average treatment effects and many other canonical target parameters.
Overfitting bias arises due to statistical dependence between the error in the nuisance parameter estimator and the data used in constructing the plug-in estimator. As discussed above, ignoring the first-step estimation of $\widehat\eta$ is justified only if the term ($\star$) in (ref) is asymptotically negligible. In the previous section, we outlined an argument for this term vanishing that depends on Neyman orthogonality and independence between the estimation error in the nuisance function, $\widehat{\eta} - \eta_0$, and the data used to construct the sample moment condition $\{W_i\}_{i=1}^{n}$. If $\widehat{\eta}$ is constructed using these same observations, the independence assumption will be violated. More generally, Neyman orthogonality alone is not sufficient to guarantee that first-step estimation of nuisance parameters can be ignored in inference about low-dimensional target parameters. To address this issue, DML relies on a second key ingredient: cross-fitting.
Cross-fitting is a form of repeated sample splitting intended to reduce the dependence between first-step estimation error and the data used to estimate the target parameter. Intuitively, if two independent datasets were available, one could be used to estimate the nuisance function $\widehat{\eta}$ and the other to estimate $\theta_0$ by plugging in $\widehat{\eta}$. In that case, independence between $\widehat\eta - \eta_0$ and $\frac{\partial}{\partial \eta} m(W_i; \theta_0, \eta_0)$ would follow by construction. With a single sample of independent observations, we can mimic this logic by randomly splitting the sample into two partitions, or “folds,” using one fold to estimate $\widehat{\eta}$ and the other to evaluate the score function and estimate $\theta_0$. Of course, such an approach inefficiently uses the available data because both $\eta_0$ and $\theta_0$ are estimated using only subsets of the data. Cross-fitting restores asymptotic efficiency by rotating which folds are used for each task, so all observations contribute to both steps.
The cross-fit version of the plug-in estimator with a generic score defined in (ref) is
Here, $\{I_k\}_{k=1}^K$ is a random partition of the sample of units $\{1, \ldots, n\}$ into $K$ subsamples of approximately equal size, and $\widehat{\eta}_{-k}$ denotes a first-step nuisance parameter estimate constructed using only observations excluding those in the subsample $I_k$. Because $\widehat{\eta}_{-k}$ uses only observations not in subsample $k$, the estimation error in $\widehat{\eta}_{-k}$ is independent of observations in subsample $k$, which alleviates overfitting bias. By rotating (“crossing”) samples for the estimation of nuisance parameters and the estimation of target parameters, cross-fitting asymptotically avoids losing efficiency relative to the hypothetical full-sample estimator that makes use of the true values of the nuisance parameters.
As a byproduct that comes at no additional computational cost, cross-fitting produces out-of-sample prediction errors associated with the nuisance functions. These enable calculation of diagnostics for evaluating the choice and specification of the nuisance function estimator. For example, in partially linear regression (Example 2), we have $\eta_0=(\ell_0,r_0)$ with $\ell_0(X)=E[Y\vert X]$ and $r_0(X)=E[D\vert X]$. The cross-fitted errors $(Y_i - \widehat{\ell}_{-k}(X_i))$ and $(D_i - \widehat{r}_{-k}(X_i))$---where $\widehat{\ell}_{-k}$ and $\widehat{r}_{-k}$ are estimated using data not containing observation $i$---are equivalent to cross-validated errors from $K$-fold cross-validation. Their availability allows calculating metrics for learner performance such as mean-squared prediction error (MSPE) and out-of-sample $R^2$. We discuss these and other methods for evaluating the nuisance function estimators in Section (ref).
The simulation results in Section (ref) provide finite-sample evidence that cross-fitting plays an important role in delivering reliable inference. The results illustrate that cross-fitting while using non-orthogonal scores delivers relatively little benefit in the considered examples. Similarly, the results show that inference based on Neyman orthogonal scores without cross-fitting often fails to deliver reliable inferential results. In contrast, DML estimators, which incorporate both Neyman orthogonal scores and cross-fitting, deliver the most robust performance across the considered examples.
In this section, we define the DML estimator, present an implementation algorithm, and summarize its asymptotic properties. We remain in the general semiparametric setting of Section (ref). That is, our focus is estimation and inference for a low-dimensional target parameter $\theta_0$, which is defined by moment conditions (ref) and depends on an unknown (potentially high-dimensional) nuisance parameter $\eta_0$.
Within this general setting, a DML estimator is a plug-in estimator that combines both Neyman orthogonal scores and cross-fitting. The DML estimator of $\theta_0$ is then
where $m$ is a Neyman orthogonal score (i.e., a score that satisfies (ref) and (ref)), $\{I_k\}_{k=1}^K$ is a random partition of the sample of individuals $\{1, \ldots, n\}$ into $K$ subsamples of approximately equal size, and $\{\widehat{\eta}_{-k}\}_{k=1}^K$ are cross-fitted nuisance parameter estimators.
chernozhukov2018 provide conditions that allow researchers to make valid inferential statements about $\theta_0$ using the DML estimator $\widehat{\theta}_{DML}$. These conditions involve conventional sampling and regularity conditions along with the assumption that the nuisance parameter estimator, $\widehat\eta$, converges sufficiently quickly, as discussed in Section (ref). Under these assumptions, the DML estimator is asymptotically normal: $$\sqrt{n}\widehat{\Sigma}^{-1/2}(\widehat{\theta}_{DML} - \theta_0) \overset{d}{\to} \mathcal{N}(0, \mathrm{I}),$$ where
Standard errors for $\widehat{\theta}$ are thus given by square-roots of the diagonal values of $\widehat{\Sigma}/n$.
A key practical implication of using both Neyman orthogonal scores and cross-fitting is that standard errors for $\widehat{\theta}_{DML}$ can be computed as if the nuisance functions were known. This approximation result holds in a variety of settings, including standard cross-sectional data, clustered cross-sectional data, and panel data with fixed $T$. In cross-sectional applications, the standard errors coincide with conventional heteroskedasticity-robust formulas. In clustered or panel settings, they are analogous to clustered standard errors that allow for arbitrary dependence across time within units $i$.
Algorithm (ref) illustrates computation of $\widehat{\theta}_{DML}$ in an i.i.d.\ setting. Implementation proceeds in three parts. First, data are randomly split into $K$ subsamples (Step (ref)). Second, the cross-fitted nuisance estimates are computed in each subsample (Steps (ref)--(ref)). Third, estimation and inference about the target parameter takes place (Steps (ref)--(ref)).
Importantly, the algorithm is highly general: it is a blueprint for estimating and performing inference on a broad range of target parameters $\theta_0$. The procedure accommodates a wide range of nuisance parameter estimators $\widehat{\eta}$, including classical approaches as well as modern and emerging ML methods.
Implementing Algorithm (ref) requires researchers to make several design choices. These include selecting an appropriate scheme for generating the cross-fitting folds, choosing the number of cross-fitting folds, and specifying the nuisance function estimator. We discuss these and other implementation choices in the following sections. In Section (ref), we provide a discussion of current best practices for the implementation of DML estimators.
This section employs two simulation exercises to illustrate the consequences of relying on non-Neyman orthogonal scores and of failing to perform cross-fitting when using flexible nuisance estimators such as ML methods. We discuss two simulations to underscore that bias dominates in some contexts, while regularization bias prevails in others
The first simulation focuses on the role of cross-fitting in the linear IV setting with many instruments. This problem has been thoroughly studied since the mid-1990s bekker. One of the recommendations to emerge from this literature is to use sample-splitting to reduce overfitting bias (often termed “many instruments bias” in this context). Our simulation revisits this classical many instrument setup by comparing linear IV estimation with ML-based IV estimation. As each of these approaches relies on Neyman orthogonal scores (see Example 3), the IV setting isolates the impact of overfitting. We show that cross-fitting alleviates overfitting bias and that using first stage estimators other than OLS can lead to efficiency gains.
Our second example turns to the estimation of average treatment effects, where we emphasize the role of Neyman orthogonality in delivering valid inference. Using a calibrated simulation, we compare estimators based on orthogonal and non-orthogonal scores, with and without cross-fitting. The results highlight that only the combination of orthogonality and cross-fitting yields estimators with reliable sampling behavior, thereby illustrating why these two ingredients are central to the DML framework.
We present simulation results from the canonical many-instrument linear IV model (see Example 3). We generate data as i.i.d.\ realizations from
Here $Z$ are $p=200$ simulated instruments, $\theta_0=0$ is the parameter of interest, and $(\varepsilon,\nu)'$ are correlated error terms drawn from the normal distribution. We define the nuisance function as $g_0(Z) = Z'\pi_0$, where the first six elements of $\pi_0$ are set to 0.1 and the remaining entries are set to 0. That is, only the first six instruments carry any signal. In our setting, OLS suffers from an upward bias of around 0.5. All results are based on a sample size of $n=1000$ and 1000 simulation replications.\footnote{Specifically, we generate $Z \sim N(0_p, .4 \mathrm{I}_p + .6 \iota_p \iota_p')$ and $(\varepsilon,\nu)'=N(0_2,\Sigma)$ with $\Sigma_{11}=\Sigma_{22}=1$ and $\Sigma_{21}=\Sigma_{21}=.6$. We use $0_p$ and $\iota_p$ to denote a $p \times 1$ vector of zeros and ones, respectively; $\mathrm{I}_p$ denotes a $p \times p$ identity matrix.}
We compare five estimators, each based on a Neyman orthogonal score. As a benchmark, we report results from the oracle IV estimator that knows the true first-stage function and uses $Z'\pi_0$ as the instrument (Oracle IV). This estimator is infeasible in practice but serves as an efficient baseline. We next consider two-stage least squares using all 200 instruments (2SLS). As is well-known from the many-instruments literature, 2SLS performs poorly when the number of instruments is large. We implement DML with a linear first stage regression using all 200 instruments. Note that this DML estimator is essentially split-sample IV AngristKruegerSplitSample1995,AIK:JIVE.\footnote{More precisely, AngristKruegerSplitSample1995 suggest splitting the sample randomly into two folds, using the first fold to estimate the IV first stage via linear regression, and the second to estimate the target parameter. That is, they do not cross-fit, though note the possibility of doing so in their conclusion. Furthermore, DML with a linear (first-stage) regression and $K=n$ is exactly jackknife IV AIK:JIVE.} One can thus view DML in this context as an extension of SSIV that also accommodates the use of flexible first stage estimators, for example, because the researcher does not want assume linearity and is unsure about which instruments matter. To this end, we also estimate the first stage using gradient boosted trees, both without cross-fitting (Boosted Trees) and with cross-fitting (DML with Boosted Trees). While one should consider alternative learner choices in practice, we focus on gradient boosted trees with a maximum tree depth of four, a learning rate of 0.1, and early stopping based on a 20% validation set. We use 5-fold cross-fitting for SSIV and DML with gradient boosted trees.
We report results in Figure (ref) and Table (ref). As expected and predicted by the theory, the methods that do not use cross-fitting---2SLS and Boosted Trees---perform poorly. Although the distributions of estimates are tightly concentrated, they are centered far from the truth, yielding average biases an order of magnitude larger than estimators relying on cross-fitting. As a result, their 95% confidence intervals exhibit zero coverage across the 1000 simulation replications. In contrast, the estimators with cross-fitting, SSIV and DML with Boosted Trees, perform substantially better. They have relatively small bias and produce coverage near the nominal 95% level. Compared to the infeasible oracle IV, they exhibit a visually larger spread, but are also approximately centered around zero.
Importantly, we observe a meaningful difference between the two cross-fitted estimators. DML with Boosted Trees outperforms SSIV, reflecting the advantage of using a more flexible and regularized first-stage learner in this context. While this is specific to the design considered here, it highlights a broader point: the choice of learner can matter substantially in practice. We return to this theme in the empirical example in Section (ref).
Finally, we emphasize that all five estimators rely on Neyman orthogonal scores. As such, these simulations are designed to isolate the impact of cross-fitting in reducing the impact of overfitting during nuisance parameter estimation. In the next subsection, we shift focus to the role of orthogonality itself by comparing estimators with and without Neyman orthogonal scores in the context of average treatment effect estimation.
This section uses a calibrated simulation to illustrate the importance of Neyman orthogonality for valid inference. We focus on estimation of the average treatment effect (ATE) under unconfoundedness---that is, assuming treatment is as good as randomly assigned after conditioning on covariates. This setting is instructive as it allows comparison of DML to the IPW estimator, which is commonly used despite not being based on a Neyman orthogonal score.
The simulation is based on poterba1995, who study the effect of 401(k) eligibility ($D$) on household net financial assets ($Y$), treating eligibility as random given observed covariates. The observed covariates ($X$), include age, income, education, family size, and indicators for two-earner households, home ownership, and alternative pension coverage. We use this application because it is a standard example in work on treatment effect estimation with machine learning; see, for instance, belloni2017program, chernozhukov2018, wuthrich2021, Ahrens2024_applied.
We base the simulation on the same 9915 observations used in the studies cited above. We calibrate the data-generating process for the simulation by flexibly estimating the propensity score and the conditional outcome model. We then use the covariate values for each individual in the original data to simulate new treatment assignments and outcomes. Under this design, the true value of the ATE is approximately 6,889.\footnote{Specifically, we estimate the probability of 401(k) eligibility given covariates, $r_0(X)={\mathrm{P}}(D=1\mid X)$, and the conditional mean of financial assets given treatment and covariates, $\ell_0(D,X)=\mathrm{E}[Y\mid D,X]$, using random forests trained on the full sample (1,000 trees and minimum node size 10). We estimate treatment-specific residual variances, $\sigma_d^2$, as the sample variance of $y_i-\widehat\ell(d_i,x_i)$ within treatment groups. In each simulation replication, we draw $d_{i,s}\sim\text{Bernoulli}(\widehat r(x_i))$ and generate outcomes as $y_{i,s}=\widehat\ell(d_{i,s},x_i)+\varepsilon_{i,s}$, where $\varepsilon_{i,s}\sim N(0,\widehat\sigma^2_{d_{i,s}})$.}
We compare three scores for estimating the ATE: the IPW score, the Neyman orthogonal AIPW score, and the regression adjustment (RA) score. The IPW and AIPW scores are defined in Example 4. The RA score is given by
where the nuisance parameter is $\eta(D,X) = \ell(D,X)$ with true value $\ell_0(D,X) = \mathrm{E}[Y\vert D, X].$ Like the IPW score, the RA score is not Neyman orthogonal. For each score, we report estimates both with and without cross-fitting. The DML estimator of the ATE corresponds to using the AIPW score with cross-fitting. All estimators use the same tuning parameter choices for nuisance function estimation, so differences in performance reflect only the choice of score and use of cross-fitting. We use random forests with 1000 trees and a maximum depth of 8 and 4 for the outcome and treatment propensity functions, respectively. We set the number of folds to $K=10$.
Figure (ref) is designed to highlight the role of Neyman orthogonality by comparing the DML estimator to two alternatives based on non-orthogonal scores: the cross-fit RA estimator (left panel) and the cross-fit IPW estimator (right panel). The most striking feature in both comparisons is that DML appears approximately unbiased, while the RA and IPW estimators exhibit substantial bias. The distribution of DML is also more concentrated than that of IPW and similar in spread to RA. Taken together, the plots show that DML dominates in this example due to its much smaller bias and comparable or better precision.
Table (ref) reports simulation results across a range of performance metrics. The DML and AIPW estimators, which are each based on Neyman orthogonal scores, dominate the remaining procedures. They exhibit substantially lower bias and mean absolute deviation than the IPW and RA estimators, whether or not cross-fitting is used. While the RA estimators have a slightly smaller standard deviation than DML and AIPW, this is offset by their much larger bias, implying inferior performance for most reasonable loss functions.
DML and AIPW also exhibit superior coverage performance. For these two estimators, we report confidence intervals based on the standard error estimator defined in (ref). This estimator is theoretically justified for DML under relatively weak conditions and would be valid for AIPW if overfitting were sufficiently controlled. In contrast, for the non-orthogonal IPW and RA estimators, we report coverage using the simulation standard deviation, since their first-order behavior depends directly on the nuisance function estimator, which makes valid standard error estimation challenging.
Even using the infeasible simulation standard deviation, the IPW and RA estimators substantially undercover the true ATE. While substantially better than IPW or RA, we see that the AIPW estimator also undercovers. In contrast, DML achieves near-nominal coverage using estimated standard errors, and outperforms AIPW in both bias and coverage.\footnote{If we used simulation standard deviations for DML and AIPW, coverage would be 0.948 and 0.945.} This highlights the central insight of the simulation: valid inference relies on combining Neyman orthogonality with cross-fitting.
To demonstrate the flexibility of DML, we apply it in a staggered adoption panel setting to estimate group-time average treatment effects on the treated and dynamic average treatment effects as discussed in callaway:santanna. Outside of illustrating the application of DML in a canonical panel data setting, we further use this example to discuss the additional randomness introduced in DML due to the use of sample splitting. We show that simply repeating the DML estimation can help us gauge the robustness of conclusions to particular sample splits. We then outline the median aggregation approach presented in chernozhukov2018 as a way to summarize results across repetitions.
Our example builds on dobkin2018economic and sun2021estimating. dobkin2018economic analyze the causal effect of hospital admission on several economic outcomes, including out-of-pocket medical spending, using conventional two-way fixed effects applied to a panel of U.S. households from the Health and Retirement Study (HRS). sun2021estimating extend this analysis by estimating dynamic effects using more flexible methods that allow for treatment effect heterogeneity. We use the same data as sun2021estimating, which consist of a balanced panel of 656 households observed over waves 7 through 11 of the HRS.\footnote{See Table 2 of sun2021estimating for summary statistics of the HRS sample we use.}
We proceed in two steps. In Section (ref), we estimate the change in out-of-pocket medical spending at time period $t$ caused by hospitalization at time $g$ for all potential pairs $(t, g)$ in the data. In Section (ref), we then aggregate these estimates for inference about the dynamic effects of hospitalization.
The group-time average treatment effect on the treated (GT-ATT) measures the effect of hospitalization at time $g$ on medical spending at time $t\leq T$, for individuals first hospitalized in period $g$. While the GT-ATTs are often not of primary interest, they can be aggregated to population-weighted ATTs or dynamic effects. We follow the potential outcome setup of callaway:santanna to formally define the GT-ATT. Let $G_i$ denote the time of first hospital admission for individual $i$. Let $Y_{i,t}(0)$ denote the potential outcome of out-of-pocket medical spending of individual $i$ at time $t$ if they remain untreated throughout, and let $(Y_{i,t}(g))_{g=2}^T$ be the set of potential outcomes at time $t$ for every potential time $g$ of the first hospital admission. Observed outcomes are given by $Y_{i,t}=Y_{i,t}(0) + \sum_{g=2}^T (Y_{i,t}(g) - Y_{i,t}(0))\mathbbm{1}\{G_i=g\}$. Then, the GT-ATT for group $g$ and time $t$ is defined as
Following dobkin2018economic and sun2021estimating, we leverage a parallel trends assumption on not-yet-hospitalized individuals to identify the GT-ATTs. Motivated by the robustness checks of dobkin2018economic, we consider a more flexible version of this assumption: parallel trends are assumed to hold only after conditioning on observed pre-treatment characteristics, such as age, gender, race, and education. Under this conditional parallel trends assumption and additional standard conditions callaway:santanna, the GT-ATT is equivalent to
where $\Delta_g Y_{i,t} = Y_{i,t} - Y_{i,g-1}$ is the difference of outcomes in period $t$ and the baseline period $g-1$. A standard no anticipation assumption implies $\theta_0^{ (g,t)} = 0$ for $t < g$, which is often used to test for pre-trends.
A key insight about the estimation of the GT-ATT is that for each group-time pair, the target parameter is equivalent to a conventional ATT identified under conditional unconfoundedness where the outcome is replaced by an appropriate difference of outcomes. Building on the efficient ATT score of hahn-pp, we employ an augmented inverse probability weighted (AIPW) score for the group–time average treatment effect (GT-ATT) as the basis for DML:\footnote{Our treatment mirrors callaway:santanna. The AIPW score directly builds on chang2020double, who proposes a DML estimator in the canonical two-group, two-period difference-in-differences design. See also santanna2020.}
where $\mathbbm{1}\{G_i=g\}$ denotes the treated group, $\mathbbm{1}\{G_i \ne g\}\mathbbm{1}\{G_i>t\}$ denotes units not yet treated by time $t$ that serve as controls, and the nuisance parameter $\eta = \big(q^{(g,t)}, \ell^{(g,t)}, \pi^{g}\big)$ takes true values $q_0^{(g,t)}(X_i) = {\mathrm{P}}\big(G_i=g \mid X_i,\{G_i = g\} \cup \{G_i > t\}\big),$ $\ell_0^{(g,t)}(X_i)= \mathrm{E}[\Delta_g Y_{i,t}\mid G_i \ne g,\ G_i>t,\ X_i],$ and $\pi^{g}_0 = {\mathrm{P}}(G_i=g)$. One can confirm that the score (ref) satisfies both the moment condition (ref) and Neyman orthogonality (using steps similar to Example 4, Appendix (ref)).
The application of DML to estimation of the GT-ATTs provides a useful robustness check because it seems unlikely that a researcher has prior knowledge of the parametric form of the nuisance functions. Even in this application with only few control variables, it is practically impossible to saturate the model, which would eliminate any additional need for flexible estimation. A fully interacted set of the eight control variables results in 3,072 indicators, which exceeds the number of households in the sample.
We set the number of cross-fitting folds to $K = 15$. Since the number of observations per group-time sample is relatively small, we opt for more folds than the rule of thumb of 5-10 folds often cited in the literature on cross-validation. For ease of exposition, we report only results on a single machine learner: a random forest estimator with 1000 trees and a minimum node size of 10 to estimate the nuisance parameters.
Table (ref) presents DML estimates for all identified group-time average treatment effects. Columns (1) through (5) show results from five different random partitions of the data used for cross-fitting. Although these DML estimators are asymptotically equivalent, they can yield different results in finite samples due to variation across sample splits. This variability may be particularly pronounced in settings with relatively few observations, as is often the case in staggered adoption designs.
To illustrate, consider the first-period treatment effect for group 9. This estimate is based on only 176 treated and 228 untreated observations. The corresponding DML estimates obtained from the five different considered sample splits range from 3204 to 3486---a difference of roughly 30% of the corresponding standard errors. While this variation is not large enough to change the sign or significance of the estimate, it shows that DML can produce meaningful differences in practice depending on the sample split.
We believe that replicating the DML procedure several times is important for understanding the variability induced by sample splitting. While under ideal conditions the DML estimator from any one split is asymptotically equivalent to that from any other, finite-sample differences may arise. More importantly, large variation across splits may signal deeper issues. For example, substantial across-split variation may signal that asymptotic approximations are unreliable in the sample at hand (e.g., due to heavy tails, poor overlap, or unstable nuisance estimates), or that the underlying machine learners are unstable within the available data.
However, reporting all results may lead to information overload or be practically infeasible. We therefore recommend that researchers examine the full set of estimates and standard errors across splits and ensure they are available for transparency. Following chernozhukov2018, we suggest median aggregation as a simple and practical strategy for summarizing results.\footnote{chernozhukov2018 also consider mean aggregation. We are not aware of formal arguments favoring one approach over another in finite samples. Our recommendation is heuristic.} Letting $\widehat{\theta}_s$ and $\widehat{s.e.}_s$ denote the $s^{\textnormal{th}}$ DML estimator and its standard error across $S$ replications, the median-aggregated point estimate and standard error are
where the median is applied element-wise for vector-valued parameters. The standard error accounts for both conventional sampling uncertainty and variability across splits.
Column (6) of Table (ref) presents the median-aggregated estimators for the identified GT-ATTs. The results strengthen the economic conclusions of dobkin2018economic and sun2021estimating. When we use DML to estimate the same effects allowing for parallel trends to hold only conditional on controls and using a flexible method to include those controls, we continue to estimate that hospitalization causes a substantial and significant increase in out-of-pocket medical spending in the period immediately following the event. This conclusion holds for all three cohorts considered in our analysis.
We also compare the DML estimates to those obtained using the parametric AIPW estimator proposed by santanna2020. This estimator relies on the same AIPW score as DML but differs in how the nuisance functions are estimated. Specifically, it assumes that the propensity score ${\mathrm{P}}\big(G_i=g \mid X_i,\{G_i = g\} \cup \{G_i > t\}\big)$ follows a logit model with a linear index and that the outcome regression $\mathrm{E}[\Delta_g Y_{i,t}\mid G_i \ne g,\ G_i>t,\ X_i]$ is linear in $X_i$. If both parametric models are correctly specified, the estimator is first-order equivalent to DML. If either is misspecified, inference remains valid but efficiency may be lost. If both are misspecified, the estimator is inconsistent. Because these parametric forms are typically motivated by convenience rather than theory, DML offers an appealing alternative. In practice, the choice between the two depends on the researcher’s confidence in the parametric models and access to a flexible nuisance estimator of reasonable quality. Regardless of the preferred approach, using the other as a robustness check is a sensible practice.
Columns (7) and (8) of Table (ref) present the santanna2020 estimator without and with controls, respectively. The comparison confirms that adjusting for covariates materially affects the results. Relative to the size of the standard errors, the median-aggregated DML estimates (in Column 6) and the AIPW DiD estimates (Column 8) are qualitatively similar. The most notable exception is group 9's treatment effect in Wave 10. The parametric estimate is 2534 ($s.e.=422$) and thus substantially larger than the DML estimate of 937 ($s.e.=434$). Further, group 10's pre-treatment effect in Wave 7 is negative under DML ($-1492$, $s.e.=2223$) but positive under the parametric estimator ($1222$, $s.e.=1127$), though both are imprecisely estimated due to the small cohort size. While the true effects are unknown, the discrepancy highlights the potential for results to strongly differ based on nuisance function specification. We return to this issue in Section (ref).
In many applications, researchers are primarily interested in summaries of group-time effects. A leading example is the dynamic treatment effect, which underpins many event study designs and provides a way to assess the plausibility of parallel trends.
Formally, dynamic effects are group-weighted averages of the GT-ATTs:
where $\mathcal{G}$ is the set of all treatment initiation periods, and $\theta_0^{(g,t)}$ is the group-time ATT defined earlier. In our context, this corresponds to averaging the effect of hospitalization across all groups observed $e$ periods after their initial admission, weighted by group size. That is, the dynamic effect $\tau_0^{(e)}$ measures the average effect of hospitalization $e$ periods after the initial admission.
Estimating the probabilities ${\mathrm{P}}\left(G_i=g \vert G_i + e \leq T\right)$ in (ref) is straightforward. For example, a natural estimator is given by the binning estimator $\sum_{i=1} \mathbbm{1}\{G_i=g, G_i + e \leq T\}/\sum_{i=1} \mathbbm{1}\{G_i + e \leq T\}$. Given an estimator for these probabilities and the DML estimators for the GT-ATTs, inference for $\tau_0^{(e)}$ then follows from standard results about inference on linear combinations.
Figure (ref) displays the estimated dynamic effects of hospitalization obtained from the different GT-ATT estimators considered previously. Across methods, we find no statistically significant effects in the pre-treatment period and a large, significant increase in out-of-pocket medical spending immediately following hospitalization. DML estimates and the parametric AIPW alternative are quantitatively similar overall, and the variation in DML across different cross-fitting splits is quite small.
One notable difference is at event time $-3$: the DML estimate is substantially less precise than the parametric alternative. Because only group 10 has outcomes observed three periods before hospitalization, this dynamic effect corresponds exactly to group 10's pre-treatment effect in Wave 7, which we discussed in the preceding section. Since this estimate is based solely on 163 treated and 65 never-treated observations, it is not surprising that it is relatively imprecise. That it is markedly less precise than the corresponding parametric estimate highlights how, with limited data, functional form restrictions can meaningfully shape inference. The discrepancy illustrates how DML can serve as a robustness check on parametric inference when data are limited. At a minimum, it cautions against taking the parametric results at face value without careful economic justification for the functional form assumptions.
We report results for other learners and tuning parameter choices in the companion website for this paper, and find similar results. For example, using random forests with higher regularization (minimum node size of 100) yields a first-period treatment effect estimate of 2903 ($s.e.=497$), while lower regularization (minimum node size of 1) yields 2666 ($s.e.=571$). Pre-treatment effects remain statistically insignificant at the 5% level across all learner and tuning parameter choices considered.
The robustness to learner choice is quite different in the example we consider in next section. There we revisit the study on monopsony power of dube2020 and illustrate that the selection and tuning of nuisance estimators can be empirically consequential. We use this application to demonstrate diagnostic analyses that can assist in evaluating and selecting nuisance estimators.
We revisit the dube2020 paper on monopsony power on Amazon MTurk, an online platform for hiring workers to perform tasks. The research falls into the wider literature estimating the labor supply elasticity Sokolova2021,langella2021. This example is useful for at least two reasons. First, it illustrates how DML can incorporate complex, non-tabular data, such as text or images, which are increasingly common in empirical economics ash2023a. Second, it provides an interesting setting for discussing the choice of ML method for nuisance estimation. In particular, it highlights the importance of selecting learners and tuning parameters to support credible inference. We use this example to discuss methods for evaluating such choices.
dube2020 consider multiple identification strategies and datasets for studying monopsony power. We focus here on their analysis of a cross-sectional sample compiled by ipeirotis2010analyzing. The parameter of interest is the partially linear regression coefficient $\theta_0$ introduced in Example 2. The outcome variable, $Y$, is the logarithm of the time required for a posted task to be filled. The treatment variable, $D$, is the logarithm of the payment offered. The covariates, $X$, are a high-dimensional vector of task characteristics.
The authors provide arguments under which $\theta_0$ is the negative labor supply elasticity to the firm and thus a measure of monopsony power. Since tasks vary in complexity, difficulty, and time commitment, it is important to adjust for characteristics that are likely correlated with duration and reward. Some of these task characteristics---such as title, task description, and keywords---come in the form of text. To account for these, the authors include a combination of hand-engineered features and vectorized representations of the task's title, description, and keywords. The hand-engineered features include, among others, the allotted time, common patterns identified using regular expressions, the number of keywords, and the lengths of the title and description. The vectorized text representations used by dube2020 consist of Doc2Vec embeddings, topic distributions, and n-grams.
In our reanalysis, we retain the original hand-engineered covariates but replace the text representations used by dube2020 with fine-tuned embeddings from the large language model DeBERTa v3. Large language models map text into numerical representations, referred to as embeddings, that capture semantic and contextual information. We use fine-tuning to adapt the pre-trained parameters to our setting, thereby improving predictive performance.\footnote{We implement fine-tuning using Low-Rank Adaptation hu2021loralowrankadaptationlarge. Rather than fully retraining all parameters, LoRA optimizes a low-rank decomposition of the parameters updates, effectively modifying the model within a regularized subspace.} Fine-tuning is performed separately for the outcome and treatment models in each cross-fitting iteration, alleviating concerns about inducing over-fitting bias.
We consider 12 candidate learners that are trained on hand-coded controls and fine-tuned embeddings: OLS as a simple unregularized baseline, cross-validated lasso and ridge as regularized linear methods, and, to allow for nonlinearities, three implementations of random forests, three implementations of gradient boosting (XGBoost), and three feed-forward neural network architectures. The tuning parameters are reported in Table (ref). We set the number of folds to $K=3$, repeat cross-fitting $S=5$ times, and report the median aggregate estimates. Finally, to account for possible dependence in the data, we implement cross-fitting by recruiter, assigning all observations from a given recruiter to the same fold, and compute standard errors clustered at the recruiter level.\footnote{The Online Appendix provides additional implementation details, results for $K = 5$, and results obtained when ignoring dependence in the construction of cross-fitting folds.}
Table (ref) shows that DML estimates can be highly sensitive to the choice of nuisance estimator, with point estimates ranging from $-3.9$ to $2.1$. These differences highlight a practical challenge with DML: ML offers a rich set of methods, but different learners may yield qualitatively distinct results. This challenge has been noted in the literature. For example, through calibrated simulations, Ahrens2024_applied illustrate stark differences in inferential results for target parameters estimated with DML with different nuisance estimators; see also wuthrich2021, giannone2021economic, AF:MachineLabor, and bach2024 for related discussions.
A convenient byproduct of cross-fitting is evidence about the out-of-sample performance of nuisance estimators. To provide guidance on the choice of nuisance estimator, we first consider the estimated out-of-sample $R^2$ values calculated using the cross-fitted predicted values and provided in the same table. These values indicate substantial variation across learner specifications.XGBoost (columns 7-9) achieves $R^2$ scores of around $85$% and $74$% for predicting outcome and treatment, respectively, whereas other learners perform markedly worse. The XGBoost-based estimates imply a labor supply elasticity between $0.040$ and $0.061$ ($s.e.\approx 0.022$). By contrast, the neural networks achieve outcome $R^2$ values below $4$% and treatment $R^2$ values below $22$%, and their corresponding point estimates display considerable instability.
The out-of-sample $R^2$ values provide a useful diagnostic but do not directly indicate whether differences across learners are statistically meaningful. To supplement this, Table (ref) also reports $p$-values from the cross-validation with confidence test lei2020. The null hypothesis in the CVC test is that a given learner achieves the lowest predictive risk among the full set of candidates. The alternative is that at least one other learner has a lower risk. lei2020 recommends forming a confidence set of the best learners by retaining those for which the $p$‑value exceeds a pre‑specified threshold, typically 0.1 or 0.2. Using a threshold of 0.2, the CVC tests identify learner specifications 8 and 9 (XGBoost 2 and 3) as the only candidate learners for the outcome equation and treatment equation for which the null cannot be rejected.
Rather than selecting the best-performing nuisance estimators based on the $R^2$ score or the CVC test, we can also combine them through stacking or model averaging approaches wolpert1996,Breiman1996a,vanderlaan:super. Stacking constructs a “super learner” as a weighted average of candidate learners, where the weights are chosen to minimize out-of-sample prediction error. Within the DML framework, the stacking weights can be re-estimated at each step of the cross-fitting procedure, or estimated once using the full sample by regressing the outcome on the learners' cross-fitted predicted values (see Ahrens2024_applied for a discussion of various options).
In this application, we combine all 12 learners using a method commonly employed in stacking applications, namely constrained least squares: We regress outcome and treatment against cross-fitted predicted values while imposing that the coefficients on the learner predicted values are non-negative and sum to one, and use these estimated coefficients as weights to construct “super learners” for the outcome and treatment. We opt to estimate the stacking regression only once for outcome and treatment on the full sample. The weights are reported in Table (ref). The model averaging procedure assigns large weights only to the three XGBoost specifications. The procedure produces a DML estimate of $-0.054$ with $s.e.=0.020$; see Table (ref). We obtain similar point estimates when selecting the single best learner for each nuisance function or when assigning equal weights to all learners selected by the CVC test. Our preferred estimates are thus consistent with dube2020: using several samples, their DML estimates suggest a labor supply elasticity in the 0.0299---0.198 range.
The results in this application highlight that DML estimates can be sensitive to the choice of nuisance estimator. This sensitivity is consistent with conditions in existing DML theory, which indicate that poorly tuned or ill-suited learners can yield misleading results. We emphasize that the choice of which learner to use for a particular nuisance function and data set is rarely, if ever, known ex ante. Further, there is no reason to assume that the same learner is best for every nuisance function, although this is implicitly imposed by many common estimation strategies. In practice, we recommend using a diverse set of nuisance function estimators to increase the credibility of DML estimates. Appropriate tools for evaluating these choices include cross‑fitted performance metrics such as $R^2$, the CVC test, and model averaging approaches.
DML provides a flexible framework for inference in the presence of high-dimensional nuisance parameters. By combining Neyman orthogonal scores with cross-fitting, it mitigates the impact of estimating nuisance parameters, enabling asymptotically valid inference under relatively weak conditions. These conditions are compatible with a wide range of machine learning methods, making DML particularly useful in applications involving complex non-tabular data such as text or images. More broadly, the ability to leverage flexible estimators makes DML attractive as a complementary robustness check or when researchers wish to avoid parametric assumptions imposed for convenience rather than grounded in economic reasoning.
While DML provides robustness to the estimation of nuisance parameters, it is not a panacea. It does not tell the researcher what interesting target parameters are or how to identify them. Rather, it complements careful reasoning about objects of interest and their identification.
\paragraph{Implementation Guidance.} When implementing DML, researchers face several design choices. One key decision is how to partition the data. With independent observations, forming folds at random provides a natural default. More generally, researchers should aim to form approximately independent folds by using partitions that respect any underlying dependence structure in the data, as noted in Remark (ref).
Researchers also need to select the number of folds $K$. Standard asymptotic results apply for any fixed $K$; see, e.g., chernozhukov2018. These results suggest the use of relatively small values for $K$, but do not provide more specific guidance. velez:DML provides a higher-order analysis within a restricted class of nuisance estimators and finds that performance improves with more folds, peaking at $K = n$ but with diminishing returns. Simulation evidence, including that in Section (ref), suggests that conventional choices like $K = 5$ or $10$ work well in many settings. Given this evidence and that computational complexity is often a concern in applications of DML, we recommend choosing a simple round number informed by available computational resources. When the sample size is very small, larger values of $K$ may be desirable, and if resources permit, sensitivity analysis to $K$ is worthwhile.
To address the algorithmic randomness introduced by sample-splitting, we recommend simply repeating the cross-fitting procedure multiple times and reporting summaries of the resulting estimates, as illustrated in Section (ref). This simple procedure reduces unappealing dependence on a particular random split and also serves as a useful diagnostic. Large variation across repetitions may raise concerns about finite-sample behavior or the plausibility of underlying assumptions, while stability across repetitions provides reassurance that results are not sensitive to a particular random split.
Choice of nuisance function estimator is a first-order concern. In Section (ref), we showed that results can vary substantially across learners. Because it is rarely clear ex ante which learner will perform best in a given application, we recommend considering a diverse set of candidate methods, including parametric benchmarks (e.g., linear or logistic regression), regularized regression, tree-based models, and neural networks. From a practical perspective, greater confidence in DML results is warranted when different learners deliver similar estimates, while substantial divergence across learners should be viewed as a warning sign that calls for caution and further investigation. We illustrated how predictive diagnostics, which are a natural byproduct of cross-fitting, can help identify potentially problematic learners. We encourage researchers to report their candidate learners, tuning procedures, and diagnostic metrics to promote transparency and replicability.
\paragraph{Challenges and Caveats} While DML is conceptually straightforward, its practical implementation can be challenging, and results can be highly sensitive to implementation choices. Decisions such as how nuisance functions are estimated, how many folds are used, how often cross-fitting is repeated, and how tuning parameters are selected can materially affect empirical results. Developing a deeper understanding of the finite-sample behavior of DML---including tradeoffs in learner complexity, computational complexity, fold count, and cross-fitting repetition---would be valuable not just for refining our theoretical knowledge but for guiding applied practice.
From a theoretical standpoint, part of DML's appeal is that it delivers asymptotically normal inference even when flexible methods are used to estimate nuisance functions. Crucially, provided the nuisance estimator converges sufficiently quickly, the limiting distribution does not depend on the specific choice of nuisance estimator. Of course, it is unrealistic to expect all nuisance estimators to yield sufficiently accurate estimates in all settings. Many modern algorithms do not have theory establishing sufficiently fast rates without strong assumptions. For example, the results in CVFL:hdrandomforests suggest that random forests, when applied to high-dimensional data with tuning similar to common defaults, may fail to converge quickly enough for DML, and results in cattaneo:trees establish that deep regression trees may fail to be pointwise consistent. Further work establishing combinations of data-generating processes, learners, and tuning strategies that provide rates of convergence compatible with the sufficient conditions for DML is important for clarifying when, where, and how different learners should be used.
Convergence rate conditions may also fail in settings where nuisance functions are highly complex or non-smooth. Extending DML to better handle such settings is an important direction in current research; see, for example, robins2017minimax, SPwsparsity, and zheng2025perturbedDML. Relatedly, the literature on sensitivity analysis in DML settings longstory develops tools that, while motivated by concerns about unobserved confounding, are applicable more broadly. These tools can readily be adapted to assess robustness to other forms of misspecification, including insufficient learner flexibility or excessive nuisance function complexity.
Finally, while accommodating dependent data is conceptually straightforward within the DML framework, its practical implementation raises additional challenges. Forming appropriate partitions, assessing effective sample sizes, and evaluating finite-sample performance become more subtle in dependent settings, further complicating the use of asymptotic approximations as a guide to inference. Given the prevalence of dependent data in economic applications, further development of practical diagnostics and implementation guidance tailored to this setting would be a welcome addition to the literature.
\paragraph{Functional Restrictions.} A recent proposal for improving performance involves empirical calibration of nuisance function estimates, as developed by laanDoublyRobustInference2025. Their approach adapts ideas from predictive modeling, requiring that estimated nuisance functions satisfy known properties. For example, if the nuisance function is $\eta_0(X) = \mathrm{E}[Y|X]$, we should not be able to improve mean squared prediction error by transforming $\eta_0(X)$. laanDoublyRobustInference2025 propose methods to enforce such calibration and show that, for a broad class of target parameters, DML estimators remain asymptotically normal even if only one nuisance estimator satisfies the usual convergence rate. More generally, using known functional restrictions can reduce the complexity of nuisance estimation, potentially improving performance even in complex settings. More work on incorporating economically motivated constraints seems promising, both as a practical aid and as a direction for methodological research.
\paragraph{Conclusion.} In sum, while recognizing that important challenges remain, we believe DML is a valuable addition to the empirical researcher's toolkit. It offers a framework for incorporating machine learning methods into empirical analysis, which is increasingly important as economic data grow in richness and complexity. At the same time, part of our goal in this review is to underscore that DML is not a foolproof, mechanical procedure. Empirical results can depend sensitively on implementation choices, and careful diagnostics, robustness checks, and transparent reporting are essential. Viewed in this light, DML should be understood as a structured approach that, when used thoughtfully, can improve empirical practice, either as a mechanism for obtaining primary conclusions from complex data or as a systematic robustness check alongside more traditional methods.