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.
173,513 characters · 28 sections · 163 citation commands
Double/Debiased Machine Learning for Treatment and Structural Parameters
We develop a series of simple results for obtaining root-$N$ consistent estimation, where $N$ is the sample size, and valid inferential statements about a low-dimensional parameter of interest, $\theta_0$, in the presence of a high-dimensional or “highly complex” nuisance parameter, $\eta_0$. The parameter of interest will typically be a causal parameter or treatment effect parameter, and we consider settings in which the nuisance parameter will be estimated using machine learning (ML) methods such as random forests, lasso or post-lasso, neural nets, boosted regression trees, and various hybrids and ensembles of these methods. These ML methods are able to handle many covariates and provide natural estimators of nuisance parameters when these parameters are highly complex. Here, highly complex formally means that the entropy of the parameter space for the nuisance parameter is increasing with the sample size in a way that moves us outside of the traditional framework considered in the classical semi-parametric literature where the complexity of the nuisance parameter space is taken to be sufficiently small. Offering a general and simple procedure for estimating and doing inference on $\theta_0$ that is formally valid in these highly complex settings is the main contribution of this paper.
Regularization Bias. A naive approach to estimation of $\theta_0$ using ML methods would be, for example, to construct a sophisticated ML estimator $D\widehat \theta_0+ \widehat g_0(X)$ for learning the regression function $D\theta_0+ g_0(X)$.\footnote{For instance, we could use lasso if we believe $g_0$ is well-approximated by a sparse linear combination of prespecified functions of $X$. In other settings, we could, for example, use iterative methods that alternate between random forests, for estimating $g_0$, and least squares, for estimating $\theta_0$.} Suppose, for the sake of clarity, that we randomly split the sample into two parts: a main part of size $n$, with observation numbers indexed by $i \in I$, and an auxiliary part of size $N - n$, with observations indexed by $i \in I^c$. For simplicity, we take $n = N/2$ for the moment and turn to more general cases which cover unequal split-sizes, using more than one split, and achieving the same efficiency as if the full sample were used for estimating $\theta_0$ in the formal development in Section (ref). Suppose $\widehat g_0$ is obtained using the auxiliary sample and that, given this $\widehat g_0$, the final estimate of $\theta_0$ is obtained using the main sample:
The estimator $\widehat \theta_0$ will generally have a slower than $1/\sqrt{n}$ rate of convergence, namely,
As detailed below, the driving force behind this “inferior" behavior is the bias in learning $g_0$. Figure 1 provides a numerical illustration of this phenomenon for a naive ML estimator based on a random forest in a simple computational experiment.
To heuristically illustrate the impact of the bias in learning $g_0$, we can decompose the scaled estimation error in $\widehat{\theta}_0$ as $$ \sqrt{n}(\widehat \theta_0 - \theta_0) = \underbrace{\Big (\frac{1}{n}\sum_{i \in I} D_i^2 \Big )^{-1} \frac{1}{\sqrt{n}} \sum_{i \in I}D_i U_i}_{:=a} + \underbrace{\Big (\frac{1}{n}\sum_{i \in I} D_i^2 \Big )^{-1} \frac{1}{\sqrt{n}} \sum_{i \in I}D_i ( g_0(X_i) - \widehat g_0(X_i) )}_{:=b}. $$ The first term is well-behaved under mild conditions, obeying $a \leadsto N(0, \bar\Sigma)$ for some $\bar\Sigma$. Term $b$ is the regularization bias term, which is not centered and diverges in general. Indeed, we have $$ b = ({\mathrm{E}} {D_i}^2 )^{-1} \frac{1}{\sqrt{n}} \sum_{i \in I} m_0(X_i) ( g_0(X_i) - \widehat g_0(X_i) ) +o_P(1) $$ to the first order. Heuristically, $b$ is the sum of $n$ terms that do not have mean zero, $m_0(X_i) ( g_0(X_i) - \widehat g_0(X_i) )$, divided by $\sqrt{n}$. These terms have non-zero mean because, in high dimensional or otherwise highly complex settings, we must employ regularized estimators - such as lasso, ridge, boosting, or penalized neural nets - for informative learning to be feasible. The regularization in these estimators keeps the variance of the estimator from exploding but also necessarily induces substantive biases in the estimator $\widehat g_0$ of $g_0$. Specifically, the rate of convergence of (the bias of) $\widehat g_0$ to $g_0$ in the root mean squared error sense will typically be $n^{-\varphi_g}$ with $\varphi_g < 1/2$. Hence, we expect $b$ to be of stochastic order $\sqrt{n} n^{-\varphi_g} \to \infty$ since $D_i$ is centered at $m_0(X_i) \neq 0$, which then implies ((ref)).
Overcoming Regularization Biases using Orthogonalization. Now consider a second construction that employs an “orthogonalized” formulation obtained by directly partialling out the effect of $X$ from $D$ to obtain the orthogonalized regressor $V= D - m_0(X)$. Specifically, we obtain $\widehat V = D- \widehat m_0(X)$, where $\widehat m_0$ is an ML estimator of $m_0$ obtained using the auxiliary sample of observations. We are now solving an auxiliary prediction problem to estimate the conditional mean of $D$ given $X$, so we are doing “double prediction" or “double machine learning".
After partialling the effect of $X$ out from $D$ and obtaining a preliminary estimate of $g_0$ from the auxiliary sample as before, we may formulate the following “debiased" machine learning estimator for $\theta_0$ using the main sample of observations:
By approximately orthogonalizing $D$ with respect to $X$ and approximately removing the direct effect of confounding by subtracting an estimate of $g_0$, $\check \theta_0$ removes the effect of regularization bias that contaminates ((ref)). The formulation of $\check \theta_0$ also provides direct links to both the classical econometric literature, as the estimator can clearly be interpreted as a linear instrumental variable (IV) estimator, and to the more recent literature on debiased lasso in the context where $g_0$ is taken to be well-approximated by a sparse linear combination of prespecified functions of $X$; see, e.g., BCH2011:InferenceGauss; c.h.zhang:s.zhang; javanmard2014confidence; vandeGeerBuhlmannRitov2013; BelloniChernozhukovHansen2011; and belloni2014pivotal.\footnote{Each of these works differs in terms of detail but can be viewed through the lens of either “debiasing" or “orthogonalization" to alleviate the impact of regularization bias on subsequent estimation and inference.}
To illustrate the benefits of the auxiliary prediction step and estimating $\theta_0$ with $\check\theta_0$, we sketch the properties of $\check \theta_0$ here. We can decompose the scaled estimation error of $\check\theta_0$ into three components: $$ \sqrt{n}(\check \theta_0 - \theta_0) = a^* + b^* +c^*. $$ The leading term, $a^*$, will satisfy $$ a^* = ({\mathrm{E}} V^2)^{-1} \frac{1}{\sqrt{n}} \sum_{i \in I}V_i U_i \leadsto N(0, \Sigma) $$ under mild conditions. The second term, $b^*$, captures the impact of regularization bias in estimating $g_0$ and $m_0$. Specifically, we will have $$ b^* = ({\mathrm{E}} V^2)^{-1} \frac{1}{\sqrt{n}} \sum_{i \in I} (\widehat m_0(X_i) - m_0(X_i))(\widehat g_0(X_i) - g_0(X_i)), $$ which now depends on the product of the estimation errors in $\widehat m_0$ and $\widehat g_0$. Because this term depends only on the product of the estimation errors, it can vanish under a broad range of data-generating processes. Indeed, this term is upper-bounded by $\sqrt{n} n^{-(\varphi_m + \varphi_g )}$, where $n^{-\varphi_m}$ and $n^{-\varphi_g}$ are respectively the rates of convergence of $\widehat m_0$ to $m_0$ and $\widehat g_0$ to $g_0$; and this upper bound can clearly vanish even though both $m_0$ and $g_0$ are estimated at relatively slow rates. Verifying that $\check\theta_0$ has good properties then requires that the remainder term, $c^*$, is sufficiently well-behaved. Sample-splitting will play a key role in allowing us to guarantee that $c^* = o_P(1)$ under weak conditions as outlined below and discussed in detail in Section (ref).
The Role of Sample Splitting in Removing Bias Induced by Overfitting. Our analysis makes use of sample-splitting which plays a key role in establishing that remainder terms, like $c^*$, vanish in probability. In the partially linear model, we have that the remainder $c^*$ contains terms like
that involve $1/\sqrt{n}$ normalized sums of products of structural unobservables from model (ref)-(ref) with estimation errors in learning the nuisance functions $g_0$ and $m_0$ and need to be shown to vanish in probability. The use of sample splitting allows simple and tight control of such terms. To see this, assume that observations are independent and recall that $\widehat g_0$ is estimated using only observations in the auxiliary sample. Then, conditioning on the auxiliary sample and recalling that ${\mathrm{E}}[V_i|X_i] = 0$, it is easy to verify that term ((ref)) has mean zero and variance of order $$ \frac{1}{n} \sum_{i \in I} (\widehat g_0(X_i) - g_0(X_i))^2 \to_P 0. $$ Thus, the term ((ref)) vanishes in probability by Chebyshev's inequality.
While sample splitting allows us to deal with remainder terms such as $c^*$, its direct application does have the drawback that the estimator of the parameter of interest only makes use of the main sample which may result in a substantial loss of efficiency as we are only making use of a subset of the available data. However, we can flip the role of the main and auxiliary samples to obtain a second version of the estimator of the parameter of interest. By averaging the two resulting estimators, we may regain full efficiency. Indeed, the two estimators will be approximately independent, so simply averaging them offers an efficient procedure. We call this sample splitting procedure where we swap the roles of main and auxiliary samples to obtain multiple estimates and then average the results cross-fitting. We formally define this procedure and discuss a $K$-fold version of cross-fitting in Section (ref).
Without sample splitting, terms such as ((ref)) may not vanish and can lead to poor performance of estimators of $\theta_0$. The difficulty arises because model errors, such as $V_i$, and estimation errors, such as $\widehat g_0(X_i) - g_0(X_i)$, are generally related because the data for observation $i$ is used in forming the estimator $\widehat g_0$. The association may then lead to poor performance of an estimator of $\theta_0$ that makes use of $\widehat g_0$ as a plug-in estimator for $g_0$ even when this estimator converges at a very favorable rate, say $N^{-1/2+\epsilon}$.
As an artificial but illustrative example of the problems that may result from overfitting, let $\widehat g_0(X_i) = g_0(X_i) + (Y_i-g_0(X_i))/N^{1/2-\epsilon}$ for any $i$ in the sample used to form estimator $\widehat g_0$, and note that the second term provides a simple model that captures overfitting of the outcome variable within the estimation sample. This estimator is excellent in terms of rates: If the $U_i$'s and $D_i$'s are bounded, $\widehat g_0$ converges uniformly to $g_0$ at the nearly parametric rate $N^{-1/2+\epsilon}$. Despite this fast rate of convergence, term $c^*$ now explodes if we do not use sample splitting. For example, suppose that the full sample is used to estimate both $\widehat g_0$ and $\check\theta_0$. A simple calculation then reveals that term $c^*$ becomes $$ \frac{1}{\sqrt{N}} \sum_{i = 1}^{N} V_i (\widehat g_0(X_i) - g_0(X_i)) \propto N^{\epsilon} \to \infty. $$
This bias due to overfitting is illustrated in the left panel of Figure 2. The histogram in the figure gives a simulated distribution for the studentized $\check\theta$ resulting from using the full sample and the contrived estimator $\widehat{g}(X_i)$ given above. We can see that the histogram is shifted markedly to the left demonstrating substantial bias resulting from overfitting. The right panel of Figure 2 also illustrates that this bias is completely removed by sample splitting. The results the right panel of Figure 2 make use of the two-fold cross-fitting procedure discussed above using the estimator $\check\theta$ and the contrived estimator $\widehat{g}(X_i)$ exactly as in the left panel. The difference is that $\widehat{g}(X_i)$ is formed in one half of the sample and then $\check\theta$ is estimated using the other half of the sample. This procedure is then repeated swapping the roles of the two samples and the results are averaged. We can see that the substantial bias from the full sample estimator has been removed and that the spread of the histogram corresponding to the cross-fit estimator is roughly the same as that of the full sample estimator clearly illustrating the bias-reduction property and efficiency of the cross-fitting procedure.
A less contrived example that highlights the improvements brought by sample-splitting is the sparse high-dimensional instrumental variable (IV) model analyzed in BellChenChernHans:nonGauss. Specifically, they consider the IV model $$ Y = D \theta_0 + \epsilon $$ where ${\mathrm{E}}[\epsilon|D] \ne 0$ but instruments $Z$ exist such that ${\mathrm{E}}[D |Z]$ is not a constant and ${\mathrm{E}}[\epsilon|Z] = 0$. Within this model, BellChenChernHans:nonGauss focus on the problem of estimating the optimal instrument, $\eta_0(Z) = {\mathrm{E}}[D|Z]$ using lasso-type methods. If $\eta_0(Z)$ is approximately sparse in the sense that only $s$ terms of the dictionary of series transformations $B(Z) = (B_1(Z), \ldots, B_p(Z))$ are needed to approximate the function accurately, BellChenChernHans:nonGauss require that $s^2 \ll n$ to establish their asymptotic results when sample splitting is not used but show that these results continue to hold under the much weaker requirement that $s \ll n$ if one employs sample splitting. We note that this example provides a prototypical example where Neyman orthogonality holds and ML methods can usefully be adopted to aid in learning structural parameters of interest. We also note that the weaker conditions required when using sample sample-splitting would also carry over to sparsity-based estimators in the partially linear model cited above. We discuss this in more detail in Section (ref).
While we find substantial appeal in using sample-splitting, one may also use empirical process methods to verify that biases introduced due to overfitting are negligible. For example, consider the problematic term in the partially linear model described previously, $\frac{1}{\sqrt{n}} \sum_{i \in I}V_i (\widehat g_0(X_i) - g_0(X_i))$. This term is clearly bounded by
where $\mathcal{G}_N$ is the smallest class of functions that contains estimators of $g_0$, $\widehat g$, with high probability. In conventional semiparametric statistical and econometric analysis, the complexity of $\mathcal{G}_N$ is controlled by invoking Donsker conditions which allow verification that terms such as ((ref)) vanish asymptotically. Importantly, Donsker conditions require that $\mathcal{G}_N$ has bounded complexity, specifically a bounded entropy integral. Because of the latter property, Donsker conditions are inappropriate in settings using ML methods where the dimension of $X$ is modeled as increasing with the sample size and estimators necessarily live in highly complex spaces. For example, Donsker conditions rule out even the simplest linear parametric model with high-dimensional regressors with parameter space given by the Euclidean ball with the unit radius: $$ \mathcal{G}_N=\{ x\mapsto g(x) = x'\theta; \ \ \theta \in \mathbb{R}^{p_N}: \| \theta\| \leqslant 1\}. $$ The entropy of this model, as measured by the logarithm of the covering number, grows at the rate $p_N$. Without invoking Donsker conditions, one may still show that terms such as ((ref)) vanish as long as $\mathcal{G}_N$'s entropy does not increase with $N$ too rapidly. A fairly general treatment is given in BCFH:Policy who provide a set of conditions under which terms like $c^*$ can vanish making use of the full sample. However, these conditions on the growth of entropy could result in unnecessarily strong restrictions on model complexity, such as very strict requirements on sparsity in the context of lasso estimation as demonstrated in IV example mentioned above. Sample splitting allows one to obtain good results under very weak conditions.
Neyman Orthogonality and Moment Conditions. Now we turn to a generalization of the orthogonalization principle above. The first “conventional” estimator $\widehat\theta_0$ given in (ref) can be viewed as a solution to estimating equations $$ \frac{1}{n} \sum_{i \in I} \varphi(W; \widehat \theta_0, \widehat g_0) = 0, $$ where $\varphi$ is a known “score" function and $\widehat g_0$ is the estimator of the nuisance parameter $g_0$. For example, in the partially linear model above, the score function is $\varphi(W; \theta, g) = (Y- \theta D- g(X)) D.$ It is easy to see that this score function $\varphi$ is sensitive to biased estimation of $g$. Specifically, the Gateaux derivative operator with respect to $g$ does not vanish: $$ \partial_g {\mathrm{E}} \varphi(W; \theta_0, g_0) [g - g_0] \neq 0.\footnote{See Section \ref{sec: orthogonal score construction} for the definition of the Gateaux derivative operator.} $$ The proofs of the general results in Section (ref) show that this term's vanishing is a key to establishing good behavior of an estimator for $\theta_0$.
By contrast the orthogonalized or double/debiased ML estimator $\check\theta_0$ given in (ref) solves $$ \frac{1}{n} \sum_{i \in I} \psi(W; \check \theta_0, \widehat \eta_0) =0, $$ where $\widehat \eta_0$ is the estimator of the nuisance parameter $\eta_0$ and $\psi$ is an orthogonalized or debiased “score” function that satisfies the property that the Gateaux derivative operator with respect to $\eta$ vanishes when evaluated at the true parameter values:
We refer to property ((ref)) as “Neyman orthogonality” and to $\psi$ as the Neyman orthogonal score function due to the fundamental contributions in N59 and Neyman1979, where this notion was introduced. Intuitively, the Neyman orthogonality condition means that the moment conditions used to identify $\theta_0$ are locally insensitive to the value of the nuisance parameter which allows one to plug-in noisy estimates of these parameters without strongly violating the moment condition. In the partially linear model (ref)-(ref), the estimator $\check\theta_0$ uses the score function $\psi(W; \theta, \eta) = (Y- D \alpha - g(X) ) (D-m(X)),$ with the nuisance parameter being $\eta = (m, g)$. It is easy to see that these score functions $\psi$ are not sensitive to biased estimation of $\eta_0$ in the sense that ((ref)) holds. The proofs of the general results in Section (ref) show that this property and sample splitting are two generic keys that allow establishing good behavior of an estimator for $\theta_0$.
Our paper builds upon two important bodies of research within the semiparametric literature. The first is the literature on obtaining $\sqrt{N}$-consistent and asymptotically normal estimates of low-dimensional objects in the presence of high-dimensional or nonparametric nuisance functions. The second is the literature on the use of sample-splitting to relax entropy conditions. We provide links to each of these literatures in turn.
The problem we study is obviously related to the classical semiparametric estimation framework which focuses on obtaining $\sqrt{N}$-consistent and asymptotically normal estimates for low-dimensional components with nuisance parameters estimated by conventional nonparametric estimators such as kernels or series. See, for example, the work by levit:75, IH:book, bickel:1982, robinson, newey90, vaart:1991, andrews94, newey94, newey1998undersmoothing, robins:dr, linton96, bickel:semibook, CLK:EfficientSP, newey2004twicing, vanderlaan:book, and AC2012. Neyman orthogonality ((ref)), introduced by N59, plays a key role in optimal testing theory and adaptive estimation, semiparametric learning theory and econometrics, and, more recently, targeted learning theory. For example, andrews94, newey94 and vdV provide a general set of results on estimation of a low-dimensional parameter $\theta_{0}$ in the presence of nuisance parameters $\eta_0$. andrews94 uses Neyman orthogonality (1.8) and Donsker conditions to demonstrate the key equicontinuity condition $$ \frac{1}{\sqrt{n}} \sum_{i\in I}\Big( \psi(W_i; \theta_0, \widehat \eta) - \int \psi(w; \theta_0, \widehat \eta) d P(w) - \psi(W_i; \theta_0, \eta_0)\Big) \to_P 0, $$ which reduces to ((ref)) in the partially linear regression model. newey94 gives conditions on estimating equations and nuisance function estimators so that nuisance function estimators do not affect the limiting distribution of parameters of interest, providing a semiparametric version of Neyman orthogonality. vdV discusses use of semiparametrically efficient scores to define estimators that solve estimating equations setting averages of efficient scores to zero. He also uses efficient scores to define k-step estimators, where a preliminary estimator is used to estimate the efficient score and then updating is done to further improve estimation; see also comments below on the use of sample-splitting.
There is also a related targeted maximum likelihood learning approach, introduced in SRR1999-rejoinder in the context of treatments effects analysis and substantially generalized by van2006targeted. van2006targeted use maximum likelihood in a least favorable direction and then perform “one-step" or “k-step" updates using the estimated scores in an effort to better estimate the target parameter.\footnote{Targeted minimum loss estimation, which shares similar properties, is also discussed in e.g. vanderlaan:book and vanderlaan:loss.} This procedure is like the least favorable direction approach in semiparametrics; see, for example, severini1992profile. The introduction of the likelihood introduces major benefits such as allowing simple and natural imposition of constraints inherent in the data, such as support restrictions when the outcome is binary or censored, and permitting the use of likelihood cross-validation to choose the nuisance parameter estimator. This data adaptive choice of the nuisance parameter has been dubbed the “super learner” by vanderlaan:super. In subsequent work, vanderlaan:book emphasize the use of ML methods to estimate the nuisance parameters for use with the super learner. Much of this work, including recent work such as luedtke2016optimal, TvdL:Struct, and ZLvdL:Struct, focuses on formal results under a Donsker condition, though the use of sample splitting to relax these conditions has also been advocated in the targeted maximum likelihood setting as discussed below.
The Donsker condition is a powerful classical condition that allows rich structures for fixed function classes $\mathcal{G}$, but it is unfortunately unsuitable for high-dimensional settings. Examples of function classes where a Donsker condition holds include functions of a single variable that have total variation bounded by 1 and functions $x \mapsto f(x)$ that have $r>\dim(x)/2$ uniformly bounded derivatives. As a further example, functions composed from function classes with VC dimensions bounded by $p$ through a fixed number of algebraic and monotone transforms are Donsker. However, this property will no longer hold if we let $\dim(x)$ grow to infinity with the sample size as this increase in dimension would require that the VC dimension also increases with $n$. More generally, Donsker conditions are easily violated once dimensions get large. A major point of departure of the present work from the classical literature on semiparametric estimation is its explicit focus on high-complexity/entropy cases. One way to analyze the problem of estimation in high-entropy cases is to see to what degree equicontinuity results continue to hold while allowing moderate growth of the complexity/entropy of $\mathcal{G}_N$. Examples of papers taking this approach in an approximately sparse settings are BCFH:Policy, belloni2014pivotal, BCY-honest, CHS:AnnRev, JM:ConfidenceIntervals, vandeGeerBuhlmannRitov2013, and c.h.zhang:s.zhang. In all of these examples, entropy growth must be limited in what may be very restrictive ways. The entropy conditions rule out the contrived overfitting example mentioned above, which does approximate realistic examples, and may otherwise place severe restrictions on the model. For example, in BellChernHans:Gauss and BellChenChernHans:nonGauss, the optimal instrument needs to be sparse of order $s \ll \sqrt{n}$.
A key device that we use to avoid strong entropy conditions is cross-fitting via sample splitting. Cross-fitting is a practical, efficient form of data splitting. Importantly, its use here is not simply as a device to make proofs elementary (which it does), but as a practical method to allow us to overcome the overfitting/high-complexity phenomena that commonly arise in data analysis based on highly adaptive ML methods. Our treatment builds upon the sample-splitting ideas employed in BellChernHans:Gauss and BellChenChernHans:nonGauss who considered sample-splitting in a high-dimensional sparse optimal IV model to weaken the sparsity condition mentioned in the previous paragraph to $s \ll n$. This work in turn was inspired by ssiv. We also build on A10 and RZALTLLV13, where ML methods and sample splitting were used in the estimation of a partially linear model of the effects of pollution while controlling for several covariates. We use the term “cross-fitting" to characterize our recommended procedure, partly borrowing the jargon from FGH:refittedCV which employed a slightly different form of sample-splitting to estimate the scale parameter in a high-dimensional sparse regression. Of course, the use of sample-splitting to relax entropy conditions has a long history in semiparametric estimation problems. For example, bickel:1982 considered estimating nuisance functions using a vanishing fraction of the sample, and these results were extended to sample splitting into two equal halves and discretization of the parameter space by schick:1986. Similarly, vdV uses 2-way sample splitting and discretization of the parameter space to give weak conditions for k-step estimators using the efficient scores where sample splitting is used to estimate the “updates"; see also vdl:AdaptiveTarget. robins2008higher and RLMTV17 use sample splitting in the construction of higher-order influence function corrections in semiparametric estimation. Some recent work in the targeted maximum likelihood literature, for example zheng2011cross, also notes the utility of sample splitting in the context of k-step updating, though this sample splitting approach is different from the cross-fitting approach we pursue.
Plan of the Paper. We organize the rest of the paper as follows. In Section (ref), we formally define Neyman orthogonality and provide a brief discussion that synthesizes various models and frameworks that may be used to produce estimating equations satisfying this key condition. In Section (ref), we carefully define DML estimators and develop their general theory. We then illustrate this general theory by applying it to provide theoretical results for using DML to estimate and do inference for key parameters in the partially linear regression model and for using DML to estimate and do inference for coefficients on endogenous variables in a partially linear instrumental variables model in Section (ref). In Section (ref), we provide a further illustration of the general theory by applying it to develop theoretical results for DML estimation and inference for average treatment effects and average treatment effects on the treated under unconfoundedness and for DML estimation of local average treatment effects in an IV context within the potential outcomes framework; see imbens:rubin:book. Finally, we apply DML in three empirical illustrations in Section (ref). In an appendix, we define additional notation and present proofs.
Notation. The symbols ${\mathrm{P}}$ and ${\mathrm{E}}$ denote probability and expectation operators with respect to a generic probability measure that describes the law of the data. If we need to signify the dependence on a probability measure $P$, we use $P$ as a subscript in ${\mathrm{P}}_P$ and ${\mathrm{E}}_P$. We use capital letters, such as $W$, to denote random elements and use the corresponding lower case letters, such as $w$, to denote fixed values that these random elements can take. In what follows, we use $\|\cdot\|_{P,q}$ to denote the $L^q(P)$ norm; for example, we denote $ \|f\|_{P,q} := \|f(W)\|_{P,q} := \left(\int |f(w)|^q d P(w)\right)^{1/q}, $ where $\|f\|_{P,\infty}$ stands for the essential supremum. We use $x'$ to denote the transpose of a column vector $x$. For a differentiable map $x \mapsto f(x)$, mapping $\mathbb{R}^d$ to $\mathbb{R}^k$, we use $\partial_{x'} f$ to abbreviate the partial derivatives $(\partial/\partial x') f$, and we correspondingly use the expression $\partial_{x'} f(x_0)$ to mean $\partial_{x'} f (x) \mid_{x = x_0}$, etc.
\setcounter{equation}{0}
Here we formally introduce the model and discuss several methods for generating orthogonal scores in a wide variety of settings, including the classical Neyman's construction. We also use this as an opportunity to synthesize some recent developments in the literature.
We are interested in the true value $\theta_0$ of the low-dimensional target parameter $\theta \in \Theta$, where $\Theta$ is a non-empty measurable subset of $\mathbb{R}^{d_\theta}$. We assume that $\theta_0$ satisfies the moment conditions
where $\psi = (\psi_1,\dots,\psi_{d_\theta})'$ is a vector of known score functions, $W$ is a random element taking values in a measurable space $(\mathcal{W}, \mathcal{A}_\mathcal{W})$ with law determined by a probability measure $P \in \mathcal{P}_N$, and $\eta_{0}$ is the true value of the nuisance parameter $\eta \in T$, where $T$ is a convex subset of some normed vector space with the norm denoted by $\|\cdot\|_T$. We assume that the score functions $ \psi_{j}\colon \mathcal{W} \times \Theta \times T \to \mathbb{R} $ are measurable once we equip $\Theta$ and $T$ with their Borel $\sigma$-fields, and we assume that a random sample $(W_i)_{i=1}^N$ from the distribution of $W$ is available for estimation and inference.
As discussed in the Introduction, we require the Neyman orthogonality condition for the score $\psi$. To introduce the condition, for $\widetilde T = \{\eta - \eta_0\colon \eta \in T\}$ we define the pathwise (or the Gateaux) derivative map $ \mathrm{D}_{r}\colon \widetilde T \to \mathbb{R}^{d_\theta}$, $$ \mathrm{D}_{r}[\eta - \eta_0]:= \partial_r \bigg\{{\mathrm{E}}_P \Big [ \psi (W; \theta_{ 0 }, \eta_0+ r(\eta - \eta_0) \Big ]\bigg\},\quad \eta\in T, $$ for all $r \in [0,1)$, which we assume to exist. For convenience, we also denote
Note that $\psi(W;\theta_0,\eta_0 + r (\eta - \eta_0))$ here is well-defined because for all $r\in[0,1)$ and $\eta\in T$, $$ \eta_0 + r(\eta - \eta_0) = (1 - r)\eta_0 + r \eta\in T $$ since $T$ is a convex set. In addition, let $\mathcal T_N \subset T$ be a {\em nuisance realization set} such that the estimators $\widehat\eta_0$ of $\eta_0$ specified below take values in this set with high probability. In practice, we typically assume that $\mathcal{T}_N$ is a properly shrinking neighborhood of $\eta_0$. Note that $\mathcal{T}_N - \eta_0$ is the nuisance deviation set, which contains deviations of $\widehat\eta_0$ from $\eta_0$, $\widehat\eta_0 - \eta_0$, with high probability. The Neyman orthogonality condition requires that the derivative in (ref) vanishes for all $\eta\in\mathcal{T}_N$.
We remark here that condition (ref) holds with $\mathcal T_N = T$ when $\eta$ is a finite-dimensional vector as long as $\partial_{\eta}{\mathrm{E}}_P[\psi_j(W;\theta_0,\eta_0)] = 0$ for all $j = 1,\dots,d_{\theta}$, where $\partial_{\eta}{\mathrm{E}}_P[\psi_j(W;\theta_0,\eta_0)]$ denotes the vector of partial derivatives of the function $\eta \mapsto {\mathrm{E}}_P[\psi_j(W;\theta_0,\eta)]$ for $\eta = \eta_0$.
Sometimes it will also be helpful to use an approximate Neyman orthogonality condition as opposed to the exact one given in Definition (ref):
If we start with a score $\varphi$ that does not satisfy the orthogonality condition above, we first transform it into a score $\psi$ that does. Here we outline several methods for doing so.
\break
First, we describe the construction used by N59 to derive his celebrated orthogonal score and $C(\alpha)$-statistic in a maximum likelihood setting.\footnote{The $C(\alpha)$-statistic, or the orthogonal score statistic, has been explicitly used for testing and estimation in high-dimensional sparse models in BCK-LAD.} Such construction also underlies the concept of local unbiasedness in construction of optimal tests in e.g. F67 and was extended to non-likelihood settings by W91. The discussion of Neyman's construction here draws on CHS15.
To describe the construction, let $\theta \in \Theta \subset \mathbb{R}^{d_{\theta}}$ and $\beta \in \mathcal B \subset \mathbb{R}^{d_{\beta}}$, where $\mathcal B$ is a convex set, be the target and the nuisance parameters, respectively. Further, suppose that the true parameter values $\theta_0$ and $\beta_0$ solve the optimization problem
where $\ell(W;\theta,\beta)$ is a known criterion function. For example, $\ell(W;\theta,\beta)$ can be the log-likelihood function associated to observation $W$. More generally, we refer to $\ell(W;\theta,\beta)$ as the quasi-log-likelihood function. Then, under mild regularity conditions, $\theta_0$ and $\beta_0$ satisfy
Note that the original score function $\varphi (W; \theta, \beta) = \partial_\theta \ell (W; \theta, \beta) $ for estimating $\theta_0$ will not generally satisfy the orthogonality condition. Now consider the new score function, which we refer to as the Neyman orthogonal score,
where the nuisance parameter is $$ \eta= (\beta', \textrm{vec}(\mu)')' \in T = \mathcal B\times\mathbb R^{d_{\theta}d_{\beta}}\subset \mathbb{R}^{p}, \quad p=d_{\beta} + d_\theta d_\beta, $$ and $\mu$ is the $d_{\theta} \times d_{\beta}$ orthogonalization parameter matrix whose true value $\mu_0$ solves the equation
for $$ J = \left(
\right) = \partial_{(\theta',\beta')}{\mathrm{E}}_P\Big[\partial_{(\theta',\beta')'}\ell(W;\theta,\beta)\Big]\Big\vert_{\theta = \theta_0; \ \beta = \beta_0}. $$ The true value of the nuisance parameter $\eta$ is
and when $J_{\beta\beta}$ is invertible, (ref) has the unique solution,
The following lemma shows that the score $\psi$ in (ref) satisfies the Neyman orthogonality condition.
If $J_{\beta \beta}$ is not invertible, equation (ref) typically has multiple solutions. In this case, it is convenient to focus on a minimal norm solution, $$ \mu_0 = \arg\min \| \mu \| \text{ such that } \|J_{\theta\beta} - \mu J_{\beta \beta} \|_{q} =0 $$ for a suitably chosen norm $\|\cdot\|_q$ on the space of $d_{\theta}\times d_{\beta}$ matrices. With an eye on solving the empirical version of this problem, we may also consider the relaxed version of this problem,
for some $r_N>0$ such that $r_N\to 0$ as $N\to \infty$. This relaxation is also helpful when $J_{\beta\beta}$ is invertible but ill-conditioned. The following lemma shows that using $\mu_0$ in (ref) leads to Neyman near-orthogonal scores. The proof of this lemma can be found in the Appendix.
\addtocounter{example}{-1}
\break
The construction in the previous section gives a Neyman orthogonal score whenever the moment conditions (ref) hold, and, as discussed in Remark (ref), the resulting score is efficient as long as $\ell(W;\theta,\beta)$ is the log-likelihood function. The question, however, remains about constructing the efficient score when $\ell(W; \theta,\beta)$ is not necessarily a log-likelihood function. In this section, we answer this question and describe a GMM-based method of constructing an efficient and Neyman orthogonal score in this more general case. The discussion here is related to L05, BMS10, and CHS:AnnRev.
Since GMM does not require that the equations (ref) are obtained from the first-order conditions of the optimization problem (ref), we use a different notation for the moment conditions. Specifically, we consider parameters $\theta\in \Theta\subset \mathbb R^{d_{\theta}}$ and $\beta\in\mathcal B\subset \mathbb R^{d_{\beta}}$, where $\mathcal B$ is a convex set, whose true values, $\theta_0$ and $\beta_0$, solve the moment conditions
where $m\colon \mathcal{W} \times \Theta \times \mathcal{B} \to \mathbb{R}^{d_m}$ is a known vector-valued function, and $d_m \geqslant d_{\theta}+ d_{\beta}$ is the number of moment conditions. In this case, a Neyman orthogonal score function is
where the nuisance parameter is $$ \eta= (\beta', \textrm{vec}(\mu)')' \in T = \mathcal B \times\mathbb R^{d_{\theta}d_m}\subset \mathbb{R}^{p}, \quad p=d_{\beta} + d_{\theta} d_m, $$ and $\mu$ is the $d_{\theta} \times d_m$ orthogonalization parameter matrix whose true value is $$ \mu_0= \Big(A' \Omega^{-1} - A' \Omega^{-1} G_\beta (G_\beta' \Omega^{-1} G_\beta)^{-1} G_\beta '\Omega^{-1} \Big), $$ where
for $\gamma = (\theta', \beta')'$ and $\gamma_0 = (\theta_0', \beta_0')'$, $A$ is a $d_m\times d_{\theta}$ moment selection matrix, $\Omega$ is a $d_m \times d_m$ positive definite weighting matrix, and both $A$ and $\Omega$ can be chosen arbitrarily. Note that setting $$ A = G_{\theta} \text{ and } \Omega = \operatorname{Var}_P(m(W; \theta_0, \beta_0) ] ) = {\mathrm{E}}_P\Big[m(W;\theta_0,\beta_0)m(W;\theta_0,\beta_0)'\Big] $$ leads to the efficient score in the sense of yielding an estimator of $\theta_0$ having the smallest variance in the class of GMM estimators (H82), and, in fact, to the semi-parametrically efficient score; see levit:75, N77, and C87. Let $\eta_0 = (\beta_0',\textrm{vec}(\mu_0)')'$ be the true value of the nuisance parameter $\eta = (\beta',\textrm{vec}(\mu)')'$. The following lemma shows that the score $\psi$ in (ref) satisfies the Neyman orthogonality condition.
As in the quasi-likelihood case, we can also consider near-orthogonal scores. Specifically, note that one of the orthogonality conditions that the score $\psi$ in (ref) has to satisfy is that $\mu_0 G_{\beta} = 0$, which can be rewritten as $$ A' \Omega^{-1/2} (I - L (L 'L)^{-1} L' ) L = 0 ,\quad\text{where } L = \Omega^{-1/2} G_\beta $$ Here, the part $A' \Omega^{-1/2} L (L' L)^{-1} L' $ can be expressed as $ \gamma_0 L'$, where $\gamma_0=A' \Omega^{-1/2} L (L' L)^{-1} $ solves the optimization problem $$ \min \| \gamma\|_o \text{ such that } \| A' \Omega^{-1/2}L - \gamma L' L\|_\infty = 0, $$ for a suitably chosen norm $\| \cdot \|_o$. When $L' L$ is close to being singular, this problem can be relaxed:
This relaxation leads to Neyman near-orthogonal scores:
\break
Here we show that the concentrating-out approach described in Remark (ref) for the case of finite-dimensional nuisance parameters can be extended to the case of infinite-dimensional nuisance parameters. Let $\ell (W;\theta ,\beta )$ be a known criterion function, where $\theta$ and $\beta$ are the target and the nuisance parameters taking values in $\Theta$ and $\mathcal B$, respectively and assume that the true values of these parameters, $\theta_0$ and $\beta_0$, solve the optimization problem (ref). The function $\ell(W; \theta,\beta)$ is analogous to that discussed above but now, instead of assuming that $\mathcal B$ is a (convex) subset of a finite-dimensional space, we assume that $\mathcal B$ is some (convex) set of functions, so that $\beta$ is the functional nuisance parameter. For example, $\ell (W;\theta ,\beta )$ could be a semiparametric log-likelihood where $\beta $ is the nonparametric part of the model. More generally, $\ell (W;\theta ,\beta )$ could be some other criterion function such as the negative of a squared residual. Also let
be the “concentrated-out" nonparametric part of the model. Note that $\beta_{\theta}$ is a function-valued function. Now consider the score function
where the nuisance parameter is $\eta\colon \Theta\to\mathcal B$, and its true value $\eta_0$ is given by $$ \eta_0(\theta) = \beta_{\theta},\quad\text{for all }\theta\in\Theta. $$ Here, the symbol $d/d \theta$ denotes the full derivative with respect to $\theta$, so that we differentiate with respect to both $\theta$ arguments in $\ell(W;\theta,\eta(\theta))$. The following lemma shows that the score $\psi$ in (ref) satisfies the Neyman orthogonality condition.
As an example, consider the partially linear model from the Introduction. Let
and let $\mathcal B$ be the set of functions of $X$ with finite mean square. Then $$ (\theta_0,\beta_0) = \arg\max_{\theta\in \Theta,\beta\in \mathcal B}{\mathrm{E}}_P[\ell(W; \theta,\beta)] $$ and
Hence, (ref) gives the following Neyman orthogonal score:
which corresponds to the estimator $\theta_0$ described in the Introduction in (ref).
It is important to note that the concentrating-out approach described here gives a Neyman orthogonal score without requiring that $\ell(W; \theta,\beta)$ is the log-likelihood function. Except for the technical conditions needed to ensure the existence of derivatives and their interchangeability, the only condition that is required is that $\theta_0$ and $\beta_0$ solve the optimization problem (ref). If $\ell(W; \theta,\beta)$ is the log-likelihood function, however, it follows from newey94, p. 1359, that the concentrating-out approach actually yields the efficient score. An alternative, but closely related, approach to derive the efficient score in the likelihood setting would be to apply Neyman's construction described above for a one-dimensional least favorable parametric sub-model; see severini1992profile and Chap. 25 of vdV.
\break
Next we consider the conditional moment restrictions framework studied in C92. To define the framework, let $W$, $R$, and $Z$ be random vectors taking values in $\mathcal W \subset \mathbb R^{d_w}$, $\mathcal R \subset \mathbb R^{d_r}$, and $\mathcal Z \subset \mathbb R^{d_z}$, respectively. Assume that $Z$ is a sub-vector of $R$ and $R$ is a sub-vector of $W$, so that $d_z\leqslant d_r\leqslant d_w$. Also, let $\theta\in\Theta\subset\mathbb R^{d_\theta}$ be a finite-dimensional parameter whose true value $\theta_0$ is of interest, and let $h$ be a vector-valued functional nuisance parameter taking values in a convex set of functions $\mathcal H$ mapping $\mathcal Z$ to $\mathbb R^{d_h}$, with the true value of $h$ being $h_0$. The conditional moment restrictions framework assumes that $\theta_0$ and $h_0$ satisfy the moment conditions
where $m\colon \mathcal W \times \Theta \times \mathbb R^{d_h}\to \mathbb R^{d_m}$ is a known vector-valued function. This framework is of interest because it covers a rich variety of models without having to explicitly rely on the likelihood formulation.
To build a Neyman orthogonal score $\psi(W;\theta,\eta)$ for estimating $\theta_0$, consider the matrix-valued functional parameter $\mu\colon \mathcal R\to \mathbb R^{d_{\theta}\times d_{m}}$ whose true value is given by
where the {\em moment selection} matrix-valued function $A\colon \mathcal R\to \mathbb R^{d_{m}\times d_{\theta}}$ and the {\em weighting} positive definite matrix-valued function $\Omega\colon \mathcal R\to\mathbb R^{d_m\times d_m}$ can be chosen arbitrarily, and the matrix-valued functions $\Gamma\colon \mathcal R\to \mathbb R^{d_{m}\times d_{\theta}}$ and $G\colon \mathcal Z\to \mathbb R^{d_{\theta}\times d_m}$ are given by
Note that $\mu_0$ in (ref) is well-defined even though the right-hand side of (ref) contains both $R$ and $Z$ since $Z$ is a sub-vector of $R$. Then a Neyman orthogonal score is
where the nuisance parameter is $$ \eta = (\mu,h)\in T = \mathcal L^1(\mathcal R;\ \mathbb R^{d_{\theta}\times d_m})\times \mathcal H. $$ Here, $\mathcal L^1(\mathcal R;\ \mathbb R^{d_{\theta}\times d_m})$ is the vector space of matrix-valued functions $f\colon \mathcal R\to\mathbb R^{d_{\theta}\times d_m}$ satisfying ${\mathrm{E}}_P[\|f(R)\|]<\infty$. Also, note that even though the matrix-valued functions $A$ and $\Omega$ can be chosen arbitrarily, setting
leads to an asymptotic variance equal to the semiparametric bound of C92. Let $\eta_0 = (\mu_0,h_0)$ be the true value of the nuisance parameter $\eta = (\mu, h)$. The following lemma shows that the score $\psi$ in (ref) satisfies the Neyman orthogonality condition.
As an application of the conditional moment restrictions framework, let us derive Neyman orthogonal scores in the partially linear regression example using this framework. The partially linear regression model (ref) is equivalent to $$ {\mathrm{E}}_P[Y - D\theta_0 - g_0(X)\mid X,D] = 0, $$ which can be written in the form of the conditional moment restrictions framework (ref) with $W = (Y,D,X')'$, $R = (D,X')'$, $Z = X$, $h(Z) = g(X)$, and $m(W;\theta,v) = Y - D\theta - v$. Hence, using (ref) and (ref) and denoting $\sigma(D,X)^2 = {\mathrm{E}}_P[U^2\mid D,X]$ for $U = Y - D\theta_0 - g_0(X)$, we can take $$ A(R) = -D,\quad \Omega(R) = {\mathrm{E}}_P[U^2\mid D,X] = \sigma(D,X)^2. $$ With this choice of $A(R)$ and $\Omega(R)$, we have $$ \Gamma(R) = - 1,\quad G(Z) = \Big({\mathrm{E}}_P\Big[\frac{D}{\sigma(D,X)^2} \mid X \Big]\Big)\times \Big( {\mathrm{E}}_P\Big[\frac{1}{\sigma(D,X)^2}\mid X\Big] \Big)^{-1}, $$ and so (ref) and (ref) give
By construction, the score $\psi$ above is efficient and Neyman orthogonal. Note, however, that using this score would require estimating the heteroscedasticity function $\sigma(D,X)^2$ which would requires the imposition of some additional smoothness assumptions over this conditional variance function. Instead, if are willing to give up on efficiency to gain some robustness, we can take $$ A(R) = -D,\quad \Omega(R) = 1; $$ in which case we have $$ \Gamma(R) = -1,\quad G(Z) = {\mathrm{E}}_P[D\mid X]. $$ (ref) and (ref) then give
This score $\psi$ is Neyman orthogonal and corresponds to the estimator of $\theta_0$ described in the Introduction in (ref). Note, however, that this score $\psi$ is efficient only if $\sigma(X,D)$ is a constant.
\break
Neyman orthogonality is a joint property of the score $\psi(W;\theta,\eta)$, the true parameter value $\eta_0$, the parameter set $T$, and the distribution of $W$. It is not determined by any particular model for the parameter $\theta$. Nevertheless, it is possible to use semiparametric efficiency calculations to construct the orthogonal score from the original score as in CEINR16. Specifically, an orthogonal score can be constructed by adding to the original score the influence function adjustment for estimation of the nuisance functions that is analyzed in newey94. The resulting orthogonal score will be the influence function of the limit of the average of the original score.
To explain, consider the original score $\varphi(W;\theta,\beta)$, where $\beta$ is some function, and let $\widehat\beta_0$ be a nonparametric estimator of $\beta_0$, the true value of $\beta$. Here, $\beta$ is implicitly allowed to depend on $\theta$, though we suppress that dependence for notational convenience. The corresponding orthogonal score can be formed when there is $\phi(W; \theta, \eta)$ such that
where $\eta$ is a vector of nuisance functions that includes $\beta$. $\phi(W;\theta,\eta)$ is an adjustment for the presence of the estimated function $\widehat\beta_0$ in the original score $\varphi(W; \theta,\beta)$. The decomposition (ref) typically holds when $\widehat\beta$ is either a kernel or a series estimator with a suitably chosen tuning parameter. The Neyman orthogonal score is given by
Here $\psi(W; \theta_0,\eta_0)$ is the influence function of the limit of $n^{-1}\sum_{i=1}^n \varphi(W_i;\theta_0,\widehat\beta_0)$, as analyzed in newey94, with the restriction ${\mathrm{E}}_P[\psi(W; \theta_0,\eta_0)] = 0$ identifying $\theta_0$.
The form of the adjustment term $\phi (W;\theta ,\eta )$ depends on the estimator $\widehat{\beta}_0$ and, of course, on the form of $\varphi (W;\theta ,\beta ).$ Such adjustment terms have been derived for various $\widehat{\beta}_0$ by newey94. Also IN15 show how the adjustment term can be computed from the limit of a certain derivative. Any of these results can be applied to a particular starting score $\varphi (W;\theta ,\beta )$ and estimator $\widehat{\beta}_0$ to obtain an orthogonal score.
For example, consider again the partially linear model with the original score
Here $\widehat{\beta}_0 = \widehat g_0$ is a nonparametric regression estimator. From newey94, we know that we obtain the influence function adjustment by taking the conditional expectation of the derivative of the score with respect to $ g_0(x)$ (obtaining $-m_0(X) = -{\mathrm{E}}_P[D|X]$) and multiplying the result by the nonparametric residual to obtain
The corresponding orthogonal score is then simply
illustrating that an orthogonal score for the partially linear model can be derived from an influence function adjustment.
Influence functions have been used to estimate functionals of nonparametric estimators by HI78 and BR88. newey1998undersmoothing, newey2004twicing showed that $n^{-1/2}\sum_{i=1}^n \psi(W_i; \theta_0, \widehat\eta_0)$ from equation (ref) will have a second order remainder in $\widehat\eta_0$, which is the key asymptotic property of orthogonal scores. Orthogonality of influence functions in semiparametric models follows from vaart:1991, as shown for higher order counterparts in robins2008higher, RLMTV17. CEINR16 point out that in general an orthogonal score can be constructed from an original score and nonparametric estimator $\widehat\beta_0$ by adding to the original score the adjustment term for estimation of $\beta_0$ as described above. This construction provides a way of obtaining an orthogonal score from any initial score $\varphi(W; \theta,\beta)$ and nonparametric estimator $\widehat\beta_0$.
\setcounter{equation}{0}
We assume that we have a sample $(W_i)_{i=1}^N$, modeled as i.i.d. copies of $W$, whose law is determined by the probability measure $P$ on $\mathcal W$. Estimation will be carried out using the finite-sample analog of the estimating equations (ref).
We assume that the true value $\eta_0$ of the nuisance parameter $\eta$ can be estimated by $\widehat\eta_0$ using a part of the data $(W_i)_{i=1}^N$. Different structured assumptions on $\eta_0$ allow us to use different machine-learning tools for estimating $\eta_0$. For instance,
There are performance guarantees for most of these ML methods that make it possible to satisfy the conditions stated below. Ensemble and aggregation methods ensure that the performance guarantee is approximately no worse than the performance of the best method.
We assume that $N$ is divisible by $K$ in order to simplify the notation. The following algorithm defines the simple cross-fitted DML as outlined in the Introduction.
This approach generalizes the 50-50 cross-fitting method mentioned in the Introduction. We now define a variation of this basic cross-fitting approach that may behave better in small samples.
We first consider the case of linear scores, where
Let $c_0 > 0$, $c_1 > 0$, $s> 0$, $q>2$ be some finite constants such that $c_0 \leqslant c_1$; and let $\{\delta_N\}_{N\geqslant 1}$ and $\{\Delta_N\}_{N\geqslant 1}$ be some sequences of positive constants converging to zero such that $\delta_N \geqslant N^{-1/2}$. Also, let $K \geqslant 2$ be some fixed integer, and let $\{\mathcal P_N\}_{N\geqslant 1}$ be some sequence of sets of probability distributions $P$ of $W$ on $\mathcal W$.
Assumption (ref) requires scores to be Neyman orthogonal or near-orthogonal and imposes mild smoothness requirements as well as the canonical identification condition.
Assumptions (ref)(a)-(c) state that the estimator of the nuisance parameter belongs to the realization set $\mathcal{T}_N \subset T$, which is a shrinking neighborhood of $\eta_0$, which contracts around $\eta_0$ with the rate determined by the “statistical” rates $r_N$, $r_N'$, and $\lambda'_N$. These rates are not given in terms of the norm $\| \cdot\|_T $ on $T$, but rather are the intrinsic rates that are most connected to the statistical problem at hand. However, in smooth problems, as discussed below this translates, in the worst cases, to the crude requirement that the nuisance parameters are estimated at the rate $o(N^{-1/4})$.
The conditions in Assumption 3.2 embody refined requirements on the quality of nuisance parameter estimators. In many applications, where $( \theta, \eta) \mapsto \psi(W; \theta, \eta)$ is smooth, we can bound
where $\varepsilon_N$ is the upper bound on the rate of convergence of $\widehat \eta_0$ to $\eta_0$ with respect to the norm $\| \cdot \|_T = \| \cdot\|_{P,2}$: $$ \| \widehat \eta_0 - \eta\|_T \lesssim \varepsilon_N. $$ Note that $\mathcal{T}_N$ can be chosen as the set of $\eta$ that is within a neighborhood of size $\varepsilon_N$ of $\eta_0$, possibly with other restrictions, in this case. If only ((ref)) holds, Assumption (ref), particularly $\lambda'_N = o(N^{-1/2})$, imposes the (crude) rate requirement
This rate is achievable for many ML methods under structured assumptions on the nuisance parameters. See, among many others, BickelRitovTsybakov2009, BvdG:book, BCW-SqLASSO, BC-SparseQR, BellChenChernHans:nonGauss, and BC-PostLasso for $\ell_1$-penalized and related methods in a variety of sparse models; Kozbur:FS for forward selection in sparse models; LuoSpindler:BoostingRates for $L_2$-boosting in sparse linear models; WW:TreeConcentration for concentration results for a class of regression trees and random forests; and CW:NNRates for a class of neural nets.
However, the presented conditions allow for more refined statements than ((ref)). We note that many important structured problems -- such as estimation of parameters in partially linear regression models, estimation of parameters in partially linear structural equation models, and estimation of average treatment effects under unconfoundedness -- are such that some cross-derivatives vanish, allowing more refined requirements than ((ref)). This feature allows us to require much finer conditions on the quality of the nuisance parameter estimators than the crude bound ((ref)). For example, in many problems
because the second derivatives vanish, $$ \partial_{r}^2 {\mathrm{E}}_P[\psi(W; \theta_0,\eta_0 +r(\eta - \eta_0))] =0. $$ This occurs in the following important examples:
If both ((ref)) and ((ref)) hold, Assumption (ref), particularly $r_N = o(1)$ and $r'_N = o(1)$, imposes the weakest possible rate requirement: $$ \varepsilon_N = o(1). $$ We note that similar refined rates have appeared in the context of estimation of treatment effects in high-dimensional settings under sparsity; see Farrell:JMP and AIW:ResidualBalancing and related discussion in Remark (ref). Our refined rate results complement this work by applying to a broad class of estimation contexts, including estimation of average treatment effects, and to a broad set of ML estimators.
The result establishes that the estimator based on the orthogonal scores achieves the root-$N$ rate of convergence and is approximately normally distributed. It is noteworthy that this convergence result, both the rate of concentration and the distributional approximation, holds uniformly with respect to $P$ varying over an expanding class of probability measures $\mathcal{P}_N$. This means that the convergence holds under any sequence of probability distributions $(P_N)_{N\geqslant 1}$ with $P_N \in \mathcal{P}_N$ for each $N$, which in turn implies that the results are robust with respect to perturbations of a given $P$ along such sequences. The same property can be shown to fail for methods not based on orthogonal scores.
Theorems (ref) and (ref) can be used for standard construction of confidence regions which are uniformly valid over a large, interesting class of models:
Indeed, the above theorem implies that $\mathrm{CI}$ obeys $ {\mathrm{P}}_{P_N} ( \ell'\theta_0 \in \mathrm{CI}) \to (1-\alpha)$ under any sequence $\{P_N\} \in \mathcal{P}_N$, which implies that these claims hold uniformly in $P \in \mathcal{P}_N$. For example, one may choose $\{P_N\}$ such that, for some $\epsilon_N \to 0$ $$ \sup_{P \in \mathcal{P}_N} |{\mathrm{P}}_P ( \ell' \theta_0 \in \mathrm{CI}) - (1-\alpha) | \leqslant |{\mathrm{P}}_{P_N} ( \ell' \theta_0 \in \mathrm{CI}) - (1-\alpha) | + \epsilon_N \to 0. $$
Next we note that the estimators need not be semi-parametrically efficient, but under some conditions they can be.
Let $c_0 > 0$, $c_1 > 0$, $a>1$, $v>0$, $s> 0$, and $q>2$ be some finite constants, and let $\{\delta_N\}_{N\geqslant 1}$, $\{\Delta_N\}_{N\geqslant 1}$, and $\{\tau_{N}\}_{N\geqslant 1}$ be some sequences of positive constants converging to zero. To derive the properties of the DML estimator, we will use the following assumptions.
Assumption (ref) is mild and rather standard in moment condition problems. Assumption (ref)(a) requires $\theta_0$ to be sufficiently separated from the boundary of $\Theta$. Assumption (ref)(b) only requires differentiability of the function $(\theta,\eta)\mapsto {\mathrm{E}}_P[\psi(W;\theta,\eta)]$ and does not require differentiability of the function $(\theta,\eta)\mapsto \psi(W;\theta,\eta)$. Assumption (ref)(c) implies sufficient identifiability of $\theta_0$. Assumption (ref)(d) is the orthogonality condition that has already been extensively discussed.
Assumptions (ref)(a)-(c) state that the estimator of the nuisance parameter belongs to the realization set $\mathcal{T}_N \subset T$, which is a shrinking neighborhood of $\eta_0$ that contracts at the “statistical" rates $r_N$ and $r_N'$ and $\lambda'_N$. These rates are not given in terms of the norm $\| \cdot\|_T $ on $T$, but rather are intrinsic rates that are connected to the statistical problem at hand. In smooth problems, these conditions translate to the crude requirement that nuisance parameters are estimated at the $o(N^{-1/4})$ rate as discussed previously in the case with linear scores. However, these conditions can be refined as, for example, when $\lambda'_N=0$ or when some cross-derivatives vanish in $\lambda'_N$; see the linear case in the previous subsection for further discussion. Suitable measurability and pointwise entropy conditions, required in Assumption (ref)(b), are mild regularity conditions that are satisfied in all practical cases. The assumption of a bounded parameter space $\Theta$ in Assumption (ref)(b) is embedded in the entropy condition, but we state it separately for clarity. This assumption was not needed in the linear case, and it can be removed in the nonlinear case with the imposition of more complicated Huber-like regularity conditions. Assumption (ref)(c) is a set of mild growth conditions.
The estimation technique developed in this paper relies on subsamples obtained by randomly partitioning the sample: an auxiliary sample for estimating the nuisance functions and a main sample for estimating the parameter of interest. Although the specific sample partition has no impact on estimation results asymptotically, the effect of the particular random split on the estimate can be important in finite samples. To make the results more robust with respect to the partitioning, we propose to repeat the DML estimator $S$ times, obtaining the estimates $$ \tilde \theta_0^s, \quad s=1, \ldots, S. $$ Features of these estimates may then provide insight into the sensitivity of results to the sample splitting, and we can report results that incorporate features of this set of estimates that should be less driven by any particular sample-splitting realization.
We recommend using medians, reporting $\tilde \theta_0^{\text{median}} $ and $ \widehat \sigma^{2 \text{Median}}$, as these quantities are more robust to outliers.
It would be interesting to investigate the behavior under the regime where $S \to \infty$ as $N \to \infty$.
\setcounter{equation}{0}
Here we revisit the partially linear regression model
The parameter of interest is the regression coefficient $\theta_0$. If $D$ is conditionally exogenous (as good as randomly assigned conditional on covariates), then $\theta_0$ measures the average causal/treatment effect of $D$ on potential outcomes.
The first approach to inference on $\theta_0$, which we described in the Introduction, is to employ the DML method using the score function
where $W = (Y,D,X)$ and $g$ and $m$ are $P$-square-integrable functions mapping the support of $X$ to $\mathbb R$. It is easy to see that $\theta_0$ satisfies the moment condition $ {\mathrm{E}}_P \psi (W; \theta_0, \eta_0) = 0,$ and also the orthogonality condition $ \partial_\eta {\mathrm{E}}_P \psi (W; \theta_0, \eta_0)[\eta - \eta_0] = 0$ where $\eta_0 = (g_0,m_0). $
A second approach employs the Robinson-style “partialling-out" score function
where $W = (Y,D,X)$ and $\ell$ and $m$ are $P$-square-integrable functions mapping the support of $X$ to $\mathbb R$. This gives an alternative parameterization of the score function above, and using this score is first-order equivalent to using the previous score. It is easy to see that $\theta_0$ satisfies the moment condition $ {\mathrm{E}}_P \psi (W; \theta_0, \eta_0) = 0,$ and also the orthogonality condition $\partial_\eta {\mathrm{E}}_P \psi (W; \theta_0, \eta_0)[\eta - \eta_0] = 0,$ for $\eta_0 = (\ell_0,m_0),$ where $\ell_0(X) = {\mathrm{E}}_P[Y|X]$.
In the partially linear model, the DML approach complements BCH2011:InferenceGauss, c.h.zhang:s.zhang, vandeGeerBuhlmannRitov2013, javanmard2014confidence, and BelloniChernozhukovHansen2011, belloni2014pivotal, and BCK-LAD, all of which consider estimation and inference for parameters within the partially linear model using lasso-type methods without cross-fitting. By relying upon cross-fitting, the DML approach allows for the use of a much broader collection of ML methods for estimating the nuisance functions and also allows relaxation of sparsity conditions in the case where lasso or other sparsity-based estimators are used. Both the DML approach and the approaches taken in the aforementioned papers can be seen as heuristically “debiasing" the score function $(Y- D \theta - g(X) ) D$, which does not possess the orthogonality property unless $m_0(X) = 0$.
Let $(\delta_N)_{n=1}^\infty$ and $(\Delta_N)_{n=1}^\infty$ be sequences of positive constants approaching 0 as before. Also, let $c$, $C$, and $q$ be fixed strictly positive constants such that $q > 4$, and let $K \geqslant 2$ be a fixed integer. Moreover, for any $\eta = (\ell_1,\ell_2)$, where $\ell_1$ and $\ell_2$ are functions mapping the support of $X$ to $\mathbb R$, denote $\|\eta\|_{P,q} = \|\ell_1\|_{P,q} \vee \|\ell_2\|_{P,q}$. For simplicity, assume that $N/K$ is an integer.
The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions (ref) and (ref) and will be proven as a special case of Theorem (ref) below.
Here we extend the partially linear regression model studied in Section (ref) to allow for instrumental variable (IV) identification. Specifically, we consider the model
where $Z$ is the instrumental variable. As before, the parameter of interest is $\theta$ and its true value is $\theta_0$. If $Z = D$, the model (ref)-(ref) coincides with (ref)-(ref) but is otherwise different.
To estimate $\theta_0$ and to perform inference on it, we will use the score
where $W = (Y,D,X,Z)$ and $g$ and $m$ are $P$-square-integrable functions mapping the support of $X$ to $\mathbb R$. Alternatively, we can use the Robinson-style score
where $W = (Y,D,X,Z)$ and $\ell$, $m$, and $r$ are $P$-square-integrable functions mapping the support of $X$ to $\mathbb R$. It is straightforward to verify that both scores satisfy the moment condition ${\mathrm{E}}_P\psi(W; \theta_0,\eta_0) = 0$ and also the orthogonality condition $\partial_{\eta}{\mathrm{E}}_P\psi(W; \theta_0, \eta_0)[\eta - \eta_0] = 0$, for $\eta_0 = (g_0, m_0)$ in the former case and $\eta_0 = (\ell_0, m_0, r_0)$ for $\ell_0$ and $r_0$ defined by $\ell_0(X) = {\mathrm{E}}_P[Y \mid X]$ and $r_0(X) = {\mathrm{E}}_P[D\mid X]$, respectively, in the latter case.\footnote{It is interesting to note that the methods for constructing Neyman orthogonal scores described in Section (ref) may give scores that are different from those in (ref) and (ref). For example, applying the method for conditional moment restriction problems in Section (ref) with $\Omega(R) = 1$ gives the score $\psi(W; \theta, \eta) = (Y - D\theta - g(X))(r(Z,X) - f(X))$, where the true values of $r(Z,X)$ and $f(X)$ are $r_0(Z,X) = {\mathrm{E}}_P[D\mid Z,X]$ and $f_0(X) = {\mathrm{E}}_P[D\mid X]$, respectively. It may be interesting to compare properties of the DML estimators $\tilde\theta_0$ based on this score with those based on (ref) and (ref) in future work.}
Note that the score in (ref) has a minor advantage over the score in (ref) because all of its nuisance parameters are conditional mean functions, which can be directly estimated by the ML methods. If one prefers to use the score in (ref), one has to construct an estimator of $g_0$ first. To do so, one can first obtain a DML estimator of $\theta_0$ based on the score in (ref), say $\tilde\theta_0$. Then, using the fact that $g_0(X) = {\mathrm{E}}_P[Y - D \theta_0 \mid X]$, one can construct an estimator $\widehat g_0$ by applying an ML method to regress $Y - D\tilde\theta_0$ on $X$. Alternatively, one can use assumption-specific methods to directly estimate $g_0$, without using the score (ref) first. For example, if $g_0$ can be approximated by a sparse linear combination of a large set of transformations of $X$, one can use the methods of GT11 to obtain an estimator of $g_0$.
Let $(\delta_N)_{n=1}^\infty$ and $(\Delta_N)_{n=1}^\infty$ be sequences of positive constants approaching 0 as before. Also, let $c$, $C$, and $q$ be fixed strictly positive constants such that $q > 4$, and let $K \geqslant 2$ be a fixed integer. Moreover, for any $\eta = (\ell_j)_{j=1}^l$ mapping the support of $X$ to $\mathbb R^l$, denote $\|\eta\|_{P,q} = \max_{1\leqslant j \leqslant l} \|\ell_j\|_{P,q}$. For simplicity, assume that $N/K$ is an integer.
The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions (ref) and (ref).
\setcounter{equation}{0}
In this section, we specialize the results of Section (ref) to estimating treatment effects under the unconfoundedness assumption of RR:prop. Within this setting, there is a large classical literature focused on low-dimensional settings that provides methods for adjusting for confounding variables including regression methods, propensity score adjustment methods, matching methods, and “doubly-robust” combinations of these methods; see, for example, robins:dr, hahn, HIR:PropWeighting, and AI:LargeSampleMatching as well as the textbook overview provided in imbens:rubin:book. In this section, we present results that complement this important classic work as well as the rapidly expanding body of work on estimation under unconfoundedness using ML methods; see, among others, AIW:ResidualBalancing, BCFH:Policy, BelloniChernozhukovHansen2011, Farrell:JMP, and IR:TEHet.
We specifically consider estimation of average treatment effects when treatment effects are fully heterogeneous and the treatment variable is binary, $D \in \{0,1\}$. We consider vectors $(Y,D,X)$ such that
Since $D$ is not additively separable, this model is more general than the partially linear model for the case of binary $D$. A common target parameter of interest in this model is the average treatment effect (ATE), $$ \theta_0 = {\mathrm{E}}_P[ g_0(1,X) - g_0(0,X)].\footnote{Without unconfoundedness/conditional exogeneity, these quantities measure association, and could be referred to as average predictive effect (APE) and average predictive effect for the exposed (APEX). Inferential results for these objects would follow immediately from Theorem 5.1.} $$ Another common target parameter is the average treatment effect for the treated (ATTE), $$ \theta_0 = {\mathrm{E}}_P[ g_0(1,X) - g_0(0,X)|D=1]. $$
The confounding factors $X$ affect the policy variable via the propensity score $m_0(X)$ and the outcome variable via the function $g_0(D,X)$. Both of these functions are unknown and potentially complicated, and we can employ ML methods to learn them.
We proceed to set up moment conditions with scores obeying orthogonality conditions. For estimation of the ATE, we employ
where the nuisance parameter $\eta = (g,m)$ consists of $P$-square-integrable functions $g$ and $m$ mapping the support of $(D,X)$ to $\mathbb R$ and the support of $X$ to $(\varepsilon,1-\varepsilon)$, respectively, for some $\varepsilon\in(0,1/2)$. The true value of $\eta$ is $\eta_0 = (g_0,m_0)$. This orthogonal moment condition is based on the influence function for the mean for missing data from robins:dr.
For estimation of the ATTE, we use the score
where the nuisance parameter $\eta = (\overline g,m,p)$ consists of $P$-square-integrable functions $\overline g$ and $m$ mapping the support of $X$ to $\mathbb R$ and to $(\varepsilon,1-\varepsilon)$, respectively, and a constant $p\in(\varepsilon, 1 - \varepsilon)$, for some $\varepsilon\in(0,1/2)$. The true value of $\eta$ is $\eta_0 = (\overline g_0,m_0,p_0)$, where $\overline g_0(X) = g_0(0,X)$ and $p_0 = {\mathrm{E}}_P[D]$. Note that estimating ATTE does not require estimating $g_0(1,X)$. Note also that since $p$ is a constant, it does not affect the DML estimators $\tilde\theta_0$ based on the score $\psi$ in (ref) but having $p$ simplifies the formula for the variance of $\tilde\theta_0$.
Using their respective scores, it can be easily seen that true parameter values $\theta_0$ for ATE and ATTE obey the moment condition $ {\mathrm{E}}_P \psi(W; \theta_0, \eta_0) = 0,$ and also that the orthogonality condition $\partial_\eta {\mathrm{E}}_P \psi (W; \theta_0, \eta_0)[\eta - \eta_0] = 0$ holds.
Let $(\delta_N)_{n=1}^\infty$ and $(\Delta_N)_{n=1}^\infty$ be sequences of positive constants approaching 0. Also, let $c, \varepsilon, C$ and $q$ be fixed strictly positive constants such that $q > 2$, and let $K \geqslant 2$ be a fixed integer. Moreover, for any $\eta = (\ell_1,\dots,\ell_l)$, denote $\|\eta\|_{P,q} =\max_{1\leqslant j\leqslant l} \|\ell_j\|_{P,q}$. For simplicity, assume that $N/K$ is an integer.
The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions (ref) and (ref).
In this section, we consider estimation of local average treatment effects (LATE) with a binary treatment variable, $D\in\{0,1\}$, and a binary instrument, $Z\in\{0,1\}$.\footnote{Similar results can be provided for the local average treatment effect on the treated (LATT) by adapting the following arguments to use the orthogonal scores for the LATT. See, for example, BCFH:Policy.} As before, $Y$ denotes the outcome variable, and $X$ is the vector of covariates.
Consider the functions $\mu_0$, $m_0$, and $p_0$, where $\mu_0$ maps the support of $(Z,X)$ to $\mathbb R$ and $m_0$ and $p_0$ respectively map the support of $(Z,X)$ and $X$ to $(\varepsilon, 1 - \varepsilon)$ for some $\varepsilon\in (0,1/2)$, such that
We are interested in estimating $$ \theta_0 = \frac{{\mathrm{E}}_P[\mu(1,X)] - {\mathrm{E}}_P[\mu(0,X)]}{{\mathrm{E}}_P[m(1,X)] - {\mathrm{E}}_P[m(0,X)]}. $$ Under the assumptions of IA94 and F07, $\theta_0$ is the LATE - the average treatment effect for compliers which are observations that would have $D = 1$ if $Z$ were $1$ and would have $D = 0$ if $Z$ were $0$. To estimate $\theta_0$, we will use the score
where $W = (Y,D,X,Z)$ and the nuisance parameter $\eta = (\mu, m, p)$ consists of $P$-square-integrable functions $\mu$, $m$, and $p$, with $\mu$ mapping the support of $(Z,X)$ to $\mathbb R$ and $m$ and $p$ respectively mapping the support of $(Z,X)$ and $X$ to $(\varepsilon, 1 - \varepsilon)$ for some $\varepsilon\in (0,1/2)$. It is easy to verify that this score satisfies the moment condition ${\mathrm{E}}_P\psi(W; \theta_0,\eta_0) = 0$ and also the orthogonality condition $\partial_{\eta}{\mathrm{E}}_P\psi(W; \theta_0,\eta_0)[\eta - \eta_0] = 0$ for $\eta_0 = (\mu_0, m_0, p_0)$.
Let $(\delta_N)_{n=1}^\infty$ and $(\Delta_N)_{n=1}^\infty$ be sequences of positive constants approaching 0. Also, let $c$, $C$, and $q$ be fixed strictly positive constants such that $q > 4$, and let $K \geqslant 2$ be a fixed integer. Moreover, for any $\eta = (\ell_1,\ell_2,\ell_3)$, where $\ell_1$ is a function mapping the support of $(Z,X)$ to $\mathbb R$ and $\ell_2$ and $\ell_3$ are functions respectively mapping the support of $(Z,X)$ and $X$ to $(\varepsilon, 1 - \varepsilon)$ for some $\varepsilon\in (0,1/2)$, denote $\|\eta\|_{P,q} = \|\ell_1\|_{P,q}\vee \|\ell_2\|_{P,2} \vee \|\ell_3\|_{P,q}$. For simplicity, assume that $N/K$ is an integer.
The following theorem follows as a corollary to the results in Section 3 by verifying Assumptions (ref) and (ref).
\setcounter{equation}{0}
To illustrate the methods developed in the preceding sections, we consider three empirical examples. The first example reexamines the Pennsylvania Reemployment Bonus experiment which used a randomized control trial to investigate the incentive effect of unemployment insurance. In the second, we use the DML method to estimate the effect of 401(k) eligibility, the treatment variable, and 401(k) participation, a self-selected decision to receive the treatment that we instrument for with assignment to the treatment state, on accumulated assets. In this example, the treatment variable is not randomly assigned and we aim to eliminate the potential biases due to the lack of random assignment by flexibly controlling for a rich set of variables. In the third, we revisit AJR-2001 IV estimation of the effects of institutions on economic growth by estimating a partially linear IV model.
In this example, we re-analyze the Pennsylvania Reemployment Bonus experiment which was conducted by the US Department of Labor in the 1980s to test the incentive effects of alternative compensation schemes for unemployment insurance (UI). This experiment has been previously studied by B:PennSeq and BK:PennQuan. In these experiments, UI claimants were randomly assigned either to a control group or one of five treatment groups.\footnote{There are six treatment groups in the experiments. Following B:PennSeq. we merge the groups 4 and 6.} In the control group, the standard rules of the UI system applied. Individuals in the treatment groups were offered a cash bonus if they found a job within some pre-specified period of time (qualification period), provided that the job was retained for a specified duration. The treatments differed in the level of the bonus, the length of the qualification period, and whether the bonus was declining over time in the qualification period; see BK:PennQuan for further details.
In our empirical example, we focus only on the most generous compensation scheme, treatment 4, and drop all individuals who received other treatments. In this treatment, the bonus amount is high and the qualification period is long compared to other treatments, and claimants are eligible to enroll in a workshop. Our treatment variable, D, is an indicator variable for being assigned treatment 4, and the outcome variable, Y, is the log of duration of unemployment for the UI claimants. The vector of covariates, X, consists of age group dummies, gender, race, the number of dependents, quarter of the experiment, location within the state, existence of recall expectations, and type of occupation.
We report results based on five simple methods for estimating the nuisance functions used in forming the orthogonal estimating equations. We consider three tree-based methods, labeled “Random Forest”, “Reg. Tree”, and “Boosting”, one $\ell_1$-penalization based method, labeled “Lasso”, and a neural network method, labeled “Neural Net". For “Reg. Tree,” we fit a single CART tree to estimate each nuisance function with penalty parameter chosen by 10-fold cross-validation. The results in the “Random Forest” column are obtained by estimating each nuisance function with a random forest which averages over 1000 trees. The results in “Boosting” are obtained using boosted regression trees with regularization parameters chosen by 10-fold cross-validation. To estimate the nuisance functions using the neural networks, we use 2 neurons and a decay parameter of 0.02, and we set activation function as logistic for classification problems and as linear for regression problems.\footnote{We also experimented with “Deep Learning” methods from which we obtained similar results for some tuning parameters. However, we ran into stability and computational issues and chose not to report these results in the empirical section.} “Lasso" estimates an $\ell_1$-penalized linear regression model using the data-driven penalty parameter selection rule developed in BellChenChernHans:nonGauss. For “Lasso”, we use a set of 96 potential control variables formed from the raw set of covariates and all second order terms, i.e. all squares and first-order interactions. For the remaining methods, we use the raw set of covariates as features.
We also consider two hybrid methods labeled “Ensemble" and “Best". “Ensemble" optimally combines four of the ML methods listed above by estimating the nuisance functions as weighted averages of estimates from “Lasso,” “Boosting,” “Random Forest,” and “Neural Net”. The weights are restricted to sum to one and are chosen so that the weighted average of these methods gives the lowest average mean squared out-of-sample prediction error estimated using 5-fold cross-validation. The final column in Table (ref) (“Best”) reports results that combine the methods in a different way. After obtaining estimates from the five simple methods and “Ensemble”, we select the best methods for estimating each nuisance functions based on the average out-of-sample prediction performance for the target variable associated with each nuisance function obtained from each of the previously described approaches. As a result, the reported estimate in the last column uses different ML methods to estimate different nuisance functions. Note that if a single method outperformed all the others in terms of prediction accuracy for all nuisance functions, the estimate in the “Best" column would be identical to the estimate reported under that method.
Table (ref) presents DML2 estimates of the ATE on unemployment duration using the median method described in Section (ref). We report results for heterogeneous effect model in Panel A and for the partially linear model in Panel B. Because the treatment is randomly assigned, we use the fraction of treated as the estimator of the propensity score in forming the orthogonal estimating equations.\footnote{We also estimated the effects using nonparametric estimates of the conditional propensity score obtained from the ML procedures given in the column labels. As expected due to randomization, the results are similar to those provided in Table (ref) and are not reported for brevity.} For both the partially linear model and the interactive model, we report estimates obtained using 2-fold cross-fitting and 5-fold cross-fitting. All results are based on taking 100 different sample splits. We summarize results across the sample splits using the median method. For comparison, we report two different standard errors. In brackets, we report the median standard error from across the 100 splits; and we report standard errors adjusted for variability across the sample splits using the median method in parentheses.
The estimation results are consistent with the findings of previous studies which have analyzed the Pennsylvania Bonus Experiment. The ATE on unemployment duration is negative and significant across all estimation methods at the 5% level regardless of the standard error estimator used. Interestingly, we see that there is no practical difference across the two different standard errors in this example.
The key problem in determining the effect of 401(k) eligibility is that working for a firm that offers access to a 401(k) plan is not randomly assigned. To overcome the lack of random assignment, we follow the strategy developed in pvw:94 and pvw:95. In these papers, the authors use data from the 1991 Survey of Income and Program Participation and argue that eligibility for enrolling in a 401(k) plan in this data can be taken as exogenous after conditioning on a few observables of which the most important for their argument is income. The basic idea of their argument is that, at least around the time 401(k) initially became available, people were unlikely to be basing their employment decisions on whether an employer offered a 401(k) but would instead focus on income and other aspects of the job. Following this argument, whether one is eligible for a 401(k) may then be taken as exogenous after appropriately conditioning on income and other control variables related to job choice.
A key component of the argument underlying the exogeneity of 401(k) eligibility is that eligibility may only be taken as exogenous after conditioning on income and other variables related to job choice that may correlate with whether a firm offers a 401(k). pvw:94 and pvw:95 and many subsequent papers adopt this argument but control only linearly for a small number of terms. One might wonder whether such specifications are able to adequately control for income and other related confounds. At the same time, the power to learn about treatment effects decreases as one allows more flexible models. The principled use of flexible ML tools offers one resolution to this tension. The results presented below thus complement previous results which rely on the assumption that confounding effects can adequately be controlled for by a small number of variables chosen ex ante by the researcher.
In the example in this paper, we use the same data as in CH401k. We use net financial assets - defined as the sum of IRA balances, 401(k) balances, checking accounts, U.S. saving bonds, other interest-earning accounts in banks and other financial institutions, other interest-earning assets (such as bonds held personally), stocks, and mutual funds less non-mortgage debt - as the outcome variable, $Y$, in our analysis. Our treatment variable, $D$, is an indicator for being eligible to enroll in a 401(k) plan. The vector of raw covariates, $X$, consists of age, income, family size, years of education, a married indicator, a two-earner status indicator, a defined benefit pension status indicator, an IRA participation indicator, and a home ownership indicator.
In Table (ref), we report DML2 estimates of ATE of 401(k) eligibility on net financial assets both in the partially linear model as in ((ref)) and allowing for heterogeneous treatment effects using the interactive model outlined in Section (ref). To reduce the disproportionate impact of extreme propensity score weights in the interactive model, we trim the propensity scores at 0.01 and 0.99. We present two sets of results based on sample-splitting as discussed in Section (ref) using 2-fold cross-fitting and 5-fold cross-fitting. As in the previous section, we consider 100 different sample partitions and summarize the results across different sample splits using the median method. For comparison, we report two different standard errors. In brackets, we report the median standard error from across the 100 splits; and we report standard errors adjusted for variability across the sample splits using the median method in parentheses. We consider the same methods with the same tuning choices for estimating the nuisance functions as in the previous example, with one exception, and so do not repeat details for brevity. The one exception is that we implement neural networks with 8 neurons and a decay parameter of 0.01 in this example.
Turning to the results, it is first worth noting that the estimated ATE of 401(k) eligibility on net financial assets is \$19,559 with an estimated standard error of 1413 when no control variables are used. Of course, this number is not a valid estimate of the causal effect of 401(k) eligibility on financial assets if there are neglected confounding variables as suggested by pvw:94 and pvw:95. When we turn to the estimates that flexibly account for confounding reported in Table (ref), we see that they are substantially attenuated relative to this baseline that does not account for confounding, suggesting much smaller causal effects of 401(k) eligibility on financial asset holdings. It is interesting and reassuring that the results obtained from the different flexible methods are broadly consistent with each other. This similarity is consistent with the theory that suggests that results obtained through the use of orthogonal estimating equations and any sensible method of estimating the necessary nuisance functions should be similar. Finally, it is interesting that these results are also broadly consistent with those reported in the original work of pvw:94 and pvw:95 which used a simple intuitively motivated functional form, suggesting that this intuitive choice was sufficiently flexible to capture much of the confounding variation in this example.
As a further illustration, we also report the LATE in this example where we take the endogenous treatment variable to be participating in a 401(k) plan. Even after controlling for features related to job choice, it seems likely that the actual choice of whether to participate in an offered plan would be endogenous. Of course, we can use eligibility for a 401(k) plan as an instrument for participation in a 401(k) plan under the conditions that were used to justify the exogeneity of eligibility for a 401(k) plan provided above in the discussion of estimation of the ATE of 401(k) eligibility.
We report DML2 results of estimating the LATE of 401(k) participation using 401(k) eligibility as an instrument in Table (ref). We employ the procedure outlined in Section (ref) using the same ML estimators to estimate the quantities used to form the orthogonal estimating equation as we employed to estimate the ATE of 401(k) eligibility outlined previously, so we omit the details for brevity. Looking at the results, we see that the estimated causal effect of 401(k) participation on net financial assets is uniformly positive and statistically significant across all of the considered methods. As when looking at the ATE of 401(k) eligibility, it is reassuring that the results obtained from the different flexible methods are broadly consistent with each other. It is also interesting that the results based on flexible ML methods are broadly consistent with, though somewhat attenuated relative to, those obtained by applying the same specification for controls as used in pvw:94 and pvw:95 and using a linear IV model which returns an estimated effect of participation of \$13,102 with estimated standard error of (1922). The mild attenuation may suggest that the simple intuitive control specification used in the original baseline specification is somewhat too simplistic.
Looking at Tables (ref) and (ref), there are other interesting observations that can provide useful insights into understanding the finite sample properties of the DML estimation method. First, the standard errors of the estimates obtained using 5-fold cross-fitting are lower than those obtained from 2-fold cross-fitting for all methods across all cases. This fact suggests that having more observations in the auxiliary sample may be desirable. Specifically, the 5-fold cross-fitting estimates use more observations to learn the nuisance functions than 2-fold cross-fitting and thus likely learn them more precisely. This increase in precision in learning the nuisance functions may then translate into more precisely estimated parameters of interest. While intuitive, we note that this statement does not seem to be generalizable in that there does not appear to be a general relationship between the number of folds in cross-fitting and the precision of the estimate of the parameter of interest; see the next example. Second, we also see that the standard errors of the Lasso estimates after adjusting for variation due to sample splitting are noticeably larger than the standard errors coming from the other ML methods. We believe that this is due to the fact that the out-of-sample prediction errors from a linear model tend to be larger when there is a need to extrapolate. In our framework, if the main sample includes observations that are outside of the range of the observations in the auxiliary sample, the model has to extrapolate to those observations. The fact that the standard errors are lower in 5-fold cross-fitting than in 2-fold cross-fitting for the “Lasso" estimations also supports this hypothesis, because the higher number of observations in the auxiliary sample reduces the degree of extrapolation. We also see that there is a noticeable increase in the standard errors that account for variability due to sample splitting relative to the simple unadjusted standard errors in this case, though these differences do not qualitatively change the results.
To demonstrate DML estimation of partially linear structural equation models with instrumental variables, we consider estimation of the effect of institutions on aggregate output following the work of AJR-2001 (AJR). Estimating the effect of institutions on output is complicated by the clear potential for simultaneity between institutions and output: Specifically, better institutions may lead to higher incomes, but higher incomes may also lead to the development of better institutions. To help overcome this simultaneity, AJR use mortality rates for early European settlers as an instrument for institution quality. The validity of this instrument hinges on the argument that settlers set up better institutions in places where they are more likely to establish long-term settlements; that where they are likely to settle for the long term is related to settler mortality at the time of initial colonization; and that institutions are highly persistent. The exclusion restriction for the instrumental variable is then motivated by the argument that GDP, while persistent, is unlikely to be strongly influenced by mortality in the previous century, or earlier, except through institutions.
In their paper, AJR note that their instrumental variable strategy will be invalidated if other factors are also highly persistent and related to the development of institutions within a country and to the country's GDP. A leading candidate for such a factor, as they discuss, is geography. AJR address this by assuming that the confounding effect of geography is adequately captured by a linear term in distance from the equator and a set of continent dummy variables. Using DML allows us to relax this assumption and replace it by a weaker assumption that geography can be sufficiently controlled by an unknown function of distance from the equator and continent dummies which can be learned by ML methods.
We use the same set of 64 country-level observations as AJR. The data set contains measurements of GDP, settler morality, an index which measures protection against expropriation risk and geographic information. The outcome variable, Y, is the logarithm of GDP per capita and the endogenous explanatory variable, D, is a measure of the strength of individual property rights that is used as a proxy for the strength of institutions. To deal with endogeneity, we use an instrumental variable Z, which is mortality rates for early European settlers. Our raw set of control variables, X, include distance from the equator and dummy variables for Africa, Asia, North America, and South America.
We report results from applying DML2 following the procedure outlined in Section 4.2 in Table (ref). The considered ML methods and tuning parameters are the same as the previous examples except for the Ensemble method, from which we exclude Neural Network since the small sample size causes stability problems in training the Neural Network. We use the raw set of covariates and all second order terms when doing lasso estimation, and we simply use the raw set of covariates in the remaining methods. As in the previous examples, we consider 100 different sample splits and report the “Median" estimates of the coefficient and two different standard error estimates. In brackets, we report the median standard error from across the 100 splits; and we report standard errors adjusted for variability across the sample splits using the median method in parentheses. Finally, we report results from both 2-fold cross-fitting and 5-fold cross-fitting as in the other examples.
In this example, we see uniformly large and positive point estimates across all procedures considered, and estimated effects are statistically significant at the 5% level. As in the second example, we see that adjusting for variability across sample splits leads to noticeable increases in estimated standard errors but does not result in qualitatively different conclusions. Interestingly, we see that the estimated standard errors based on 5-fold cross-fitting are larger than those on twofold cross-fitting in all procedures except lasso, which differs from the finding in the 401(k) example. Further understanding these differences and the impact of the number of folds on inference for objects of interest seems like an interesting question for future research. Finally, although the estimates are somewhat smaller than the baseline estimates reported in AJR - an estimated coefficient of 1.10 with estimated standard error of 0.46 (AJR-2001, Table 4, Panel A, column 7) - the results are qualitatively similar, indicating a strong and positive effect of institutions on output.
Before closing this section we want to emphasize some important conclusions that can be drawn from these empirical examples. First, the choice of the ML method used in estimating nuisance functions does not substantively change the conclusion in any of the examples, and we obtained broadly consistent results regardless of which method we employ. The robustness of the results to the different methods is implied by the theory assuming that all of the employed methods are able to deliver sufficiently high-quality approximations to the underlying nuisance functions. Second, the incorporation of uncertainty due to sample-splitting using the median method increases the standard errors relative to a baseline that does not account for this uncertainty, though these differences do not alter the main results in any of the examples. This lack of variation suggests that the parameter estimates are robust to the particular sample split used in the estimation in these examples.
We would like to acknowledge research support from the National Science Foundation. We also thank participants of the MIT Stochastics and Statistics seminar, the Kansas Econometrics conference, the Royal Economic Society Annual Conference, The Hannan Lecture at the Australasian Econometric Society meeting, The Econometric Theory lecture at the $EC^2$ meetings 2016 in Toulouse, The CORE 50th Anniversary Conference, The Becker-Friedman Institute Conference on Machine Learning and Economics, The INET conferences at USC on Big Data, the World Congress of Probability and Statistics 2016, the Joint Statistical Meetings 2016, the New England Day of Statistics Conference, CEMMAP's Masterclass on Causal Machine Learning, and St. Gallen's summer school on “Big Data", for many useful comments and questions. We would like to thank Susan Athey, Peter Aronow, Jin Hahn, Guido Imbens, Mark van der Laan, Matt Taddy, and Rui Wang for constructive comments. We thank Peter Aronow for pointing us to the literature on targeted learning on which, along with prior works of Neyman, Bickel, and the many other contributions to semiparametric learning theory, we build.