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.
129,273 characters · 53 sections · 115 citation commands
DoubleML: An Object-Oriented Implementation of Double Machine Learning in R
Structural equation models provide a quintessential framework for conducting causal inference in statistics, econometrics, machine learning (ML), and other data sciences. The package \pkg{DoubleML} DoubleMLR for \proglang{R} R implements partially linear and interactive structural equation and treatment effect models with high-dimensional confounding variables as considered in dml2018. Estimation and tuning of the machine learning models is based on the powerful functionalities provided by the \pkg{mlr3} package and the \pkg{mlr3} ecosystem mlr3. A key ingredient of double machine learning (DML) models are score functions identifying the estimates for the target parameter. These functions play an essential role for valid inference with machine learning methods because they have to satisfy a property called Neyman orthogonality. With the score functions as key elements, \pkg{DoubleML} implements double machine learning in a very general way using object orientation based on the \pkg{R6} package R6. Currently, \pkg{DoubleML} implements the double / debiased machine learning framework as established in dml2018 for
The object-oriented implementation of \pkg{DoubleML} is very flexible. The model classes `\code{DoubleMLPLR}', `\code{DoubleMLPLIV}', `\code{DoubleMLIRM}' and `\code{DoubleIIVM}' implement the estimation of the nuisance functions via machine learning methods and the computation of the Neyman-orthogonal score function. All other functionalities are implemented in the abstract base class `\code{DoubleML}', including estimation of causal parameters, standard errors, $t$ tests, confidence intervals, as well as valid simultaneous inference through adjustments of $p$ values and estimation of joint confidence regions based on a multiplier bootstrap procedure. In combination with the estimation and tuning functionalities of \pkg{mlr3} and its ecosystem, this object-oriented implementation enables a high flexibility for the model specification in terms of
It further can be readily extended regarding
Several other \proglang{R} packages for estimation of causal effects based on machine learning methods exist for \proglang{R}. The packages \pkg{grf} grf and \pkg{hdi} hdi implement alternative approaches to causal machine learning. \pkg{grf} implements generalized random forests grfpaper and can be used for forest-based inference methods in different causal models including least-squares regression and estimation of treatment effects with and without instrumental variables. \pkg{hdi} can be used for inference in high-dimensional models with a focus on lasso-based estimation and methods for simultaneous inference.
An alternative approach, which was developed before the double machine learning framework, is the so-called targeted learning framework, and its software implementations and ecosystem (\pkg{tlverse}). For an overview and introduction to this approach and its implementations, we refer to the extensive \pkg{tlverse} handbook tlmanual. Relevant \proglang{R} packages include \pkg{SuperLearner} SuperLearner for flexible estimation using machine learning and \pkg{tmle} tmle which implements estimation of causal parameters using targeted maximum likelihood estimation (TMLE). \pkg{sl3} sl3 and \pkg{tmle3} tmle3 are recent extensions of the \pkg{tlverse} for object-oriented implementation of machine learning algorithms and a unified interface for TMLE.
Previous implementations that are more closely related to the double machine learning framework of dml2018 have been provided by the \proglang{R} packages \pkg{hdm} hdm, \pkg{dmlmt} knaus2018, \pkg{causalDML} knaus2020, \pkg{causalweight} causalweight and \pkg{AIPW} AIPW. \pkg{hdm} offers lasso-based inference methods in a variety of high-dimensional causal models, including estimation of (local) average treatment effects and linear (instrumental variable) regression. The underlying theoretical framework for valid post-selection and post-regularization inference can be considered as a special case of the more generic DML framework of dml2018. Similarly, \pkg{dmlt} knaus2018 provides methods of lasso-based inference on treatment effects of multi-valued treatment variables. knaus2020 focuses on a nonparametric treatment effect model, which is called the interactive regression model (IRM) in dml2018 and also referred to as augmented inverse probability weighting. We adapted the term IRM from dml2018 to denote this causal model and will refer to it accordingly in the following. The model is introduced in Section (ref). Allowing for additional learners that include generalized random forests, ridge and random forests, \pkg{causalDML} knaus2020 provides an implementation of the DML appraoch in an IRM in combination with recent methods for the analysis of heterogeneous treatment effects. Similarly, \pkg{causalweight} causalweight focuses on the IRM and offers estimation methods for various causal quantities in this model as well as extensions thereof, including instrumental variable estimation, mediation analysis and sample selection approaches. In line with dml2018, we denote the instrumental variable extension of the IRM in the following as IIVM. The model is introduced in Section (ref). The \proglang{R} package \pkg{AIPW} AIPW implements estimation of average treatment effects of a binary treatment variable by augmented inverse probability weighting based on machine learning algorithms and integrates well with the \pkg{tlverse} ecosystem discussed later. In \proglang{Python} python, \pkg{EconML} econml offers an implementation of several causal machine learning approaches. The package does not exclusively build on the double machine learning framework by dml2018 and has a focus on heterogeneous effects.
In contrast to existing software packages, the \proglang{R} package \pkg{DoubleML} is intended to be a general implementation of the double machine learning approach of dml2018. An introduction to the three key ingredients of the DML framework is provided in Section (ref). The package can be used to perform inference in basically any causal model that can be characterized in terms of the formal framework of dml2018. For example, it would be straightforward to extend \pkg{DoubleML} to mediation analysis farbmacher2022, sample selection models bia2020 or difference-in-differences chang2020. As we will point out later, a key requirement for new model classes is a Neyman-orthogonal score. The object-oriented implementation makes the package easily extendable in terms of the supported causal models and other features of DML. By building on the \pkg{mlr3} ecosystem estimation can be based on a rich collection of powerful ML methods available in \pkg{mlr3} mlr3, \pkg{mlr3learners} mlr3learners and \pkg{mlr3extralearners} mlr3extralearners. The package \pkg{DoubleML} is available from the Comprehensive \proglang{R} Archive Network (CRAN) at \url{https://CRAN.R-project.org/package=DoubleML}.
We would like to mention that the \proglang{R} package \pkg{DoubleML} was developed together with a \proglang{Python} twin DoubleMLpython that is based on \pkg{scikit-learn} pedregosa2011. The \proglang{Python} package is also available via {GitHub} (\url{https://github.com/DoubleML/doubleml-for-py}), the Python Package Index ({PyPI}, \url{https://pypi.org/project/DoubleML}), and {conda-forge} (\url{https://anaconda.org/conda-forge/doubleml}). Moreover, doubleml_serverless provides a serverless implementation of the \proglang{Python} module \pkg{DoubleML}.
The rest of the paper is structured as follows: In Section (ref), we briefly demonstrate how to install the \pkg{DoubleML} package and give a short motivating example to illustrate the major idea behind the double machine learning approach. Section (ref) introduces the main causal model classes implemented in \pkg{DoubleML}. Section (ref) shortly summarizes the main ideas behind the double machine learning approach and reviews the key ingredients required for valid inference based on machine learning methods. Section (ref) presents the main steps and algorithms of the double machine learning procedure for inference on one or multiple target parameters. Section (ref) provides more detailed insights on the implemented classes and methods of \pkg{DoubleML}. Section (ref) contains real-data and simulation examples for estimation of causal parameters using the \pkg{DoubleML} package. Additionally, this section provides a brief simulation study that illustrates the validity of the implemented methods in finite samples. Section (ref) concludes the paper. The code output that has been suppressed in the main text and further information regarding the simulations are presented in the appendix. To make the code examples fully reproducible, the entire code is available in a supplementary zip file for this paper, as well as at \url{https://github.com/DoubleML/DoubleMLReplicationCode}. We would like to note that minor numerical differences might occur on other platforms when replicating code examples that involve random forest learners (see Appendix (ref) for more information on the infrastructure used).
The latest CRAN release of \pkg{DoubleML} can be installed using the command
Alternatively, the development version can be downloaded and installed from the {GitHub} (\url{https://github.com/DoubleML/doubleml-for-r}) repository using the command remotes
Among others, \pkg{DoubleML} depends on the \proglang{R} package \pkg{R6} for object oriented implementation, \pkg{data.table} datatable for the underlying data structure, as well as the packages \pkg{mlr3} mlr3, \pkg{mlr3learners} mlr3learners and \pkg{mlr3tuning} mlr3tuning for estimation of machine learning methods, model tuning and parameter handling. Moreover, the underlying packages of the machine learning methods that are called in \pkg{mlr3} or \pkg{mlr3learners} must be installed, for example the packages \pkg{glmnet} for lasso estimation glmnet or \pkg{ranger} ranger for random forests.
Load the package after completed installation.
In the following, we provide a brief summary of and motivation to double machine learning methods and show how the corresponding methods provided by the \pkg{DoubleML} package can be applied. The data generating process (DGP) is based on the introductory example in dml2018. We consider a partially linear model: Our major interest is to estimate the causal parameter \(\theta\) in the following regression equation
with covariates \(x_i \sim \mathcal{N}(0, \Sigma)\), where \(\Sigma\) is a matrix with entries \(\Sigma_{kj} = 0.7^{\lvert j-k\rvert}\). In the following, the regression relationship between the treatment variable \(d_i\) and the covariates \(x_i\) will play an important role
The nuisance functions \(m_0\) and \(g_0\) are given by
We construct a setting with \(n=500\) observations and \(p=20\) explanatory variables to demonstrate the use of the estimators provided in \pkg{DoubleML}. Moreover, we set the true value of the parameter \(\theta\) to \(\theta=0.5\). The corresponding data generating process is implemented in the function \code{make_plr_CCDHNR2018()}. We start by generating a realization of a data set as a `\code{data.table}' object, which is subsequently used to create an instance of the data backend of class `\code{DoubleMLData}'.
The data backend implements the causal model: We specify that we perform inference on the effect of the treatment variable \(d_i\) on the dependent variable \(y_i\).
In the next step, we choose the machine learning method as an object of class `\code{Learner}' from \pkg{mlr3}, \pkg{mlr3learners} mlr3learners or \pkg{mlr3extralearners} mlr3extralearners. As we will point out later, we have to estimate two nuisance functions in order to perform valid inference in the partially linear regression model. Hence, we have to specify two learners. Moreover, we split the sample into two folds used for cross-fitting (\code{n_folds\ =\ 2}) in our illustrating examples for simplicity. Two-fold cross-fitting makes it necessary to estimate the ML models only twice, i.e., once per fold, which reduces the computational costs. In practice, it is generally recommended to choose a larger number of folds, cf. Remark 3. The default for the number of folds is \code{n_folds\ =\ 5}.
Load \pkg{mlr3} and \pkg{mlr3learners} packages and suppress output during estimation.
Initialize a random forests learner with specified parameters.
Initialize a causal model object, here a PLR.
To estimate the causal effect of variable \(d_i\) on \(y_i\), we call the \code{fit()} method. To summarize the estimation output, we call the \code{summary()} method.
The output shows that the estimated coefficient is close to the true parameter \(\theta=0.5\). Moreover, there is evidence to reject the null hypotheses \(H_0: \theta=0\) at all common significance levels.
\pkg{DoubleML} provides estimation of causal effects in four different models: Partially linear regression models (PLR), partially linear instrumental variable regression models (PLIV), interactive regression models (IRM) and interactive instrumental variable regression models (IIVM). We will shortly introduce these models.
Partially linear regression models (PLR), which encompass the standard linear regression model, play an important role in data analysis robinson1988. Partially linear regression models take the form
where \(Y\) is the outcome variable and \(D\) is the policy variable of interest. The high-dimensional vector \(X=(X_1, \dots ,X_p)\) consists of other confounding covariates, and \(\zeta\) and \(V\) are stochastic errors. Equation (ref) is the equation of interest, and \(\theta_0\) is the main regression coefficient that we would like to infer. If \(D\) is conditionally exogenous (randomly assigned conditional on \(X\)), \(\theta_0\) has the interpretation of a structural or causal parameter. The causal diagram supporting such interpretation is shown in Figure (ref). The second equation keeps track of confounding, namely the dependence of \(D\) on covariates/controls. The characteristics \(X\) affect the policy variable \(D\) via the function \(m_0(X)\) and the outcome variable via the function \(g_0(X)\). The partially linear model generalizes both linear regression models, where functions \(g_0\) and \(m_0\) are linear with respect to a collection of basis functions with respect to \(X\), and approximately linear models.
An applied example from the economics literature is the analysis of the causal effect of 401(k) pension plans on employees' net financial assets by 401k1 and 401k2. In these studies, which are based on observational data, it is argued that eligibility for 401(k) pension plans can be assumed to be conditionally exogenous, once it is controlled for a set of confounders \(X\), for example income. Following this argumentation and modelling approach, the estimate on \(\theta_0\) as obtained by a PLR can be interpreted as the average treatment effect of 401(k) eligibility on net financial assets. A reassessment and summary of the 401(k) example is available in dml2018 as well as on the \pkg{DoubleML} website (\url{https://docs.doubleml.org/stable/examples/R_double_ml_pension.html}).
We next consider the partially linear instrumental variable regression model (PLIV)
Note that this model is not a regression model unless \(Z=D\). Model (ref)--(ref) is a canonical model in causal inference, going back to wright1928, with the modern difference being that \(g_0\) and \(m_0\) are nonlinear, potentially complicated functions of high-dimensional \(X\). The idea of this model is that there is a structural or causal relation between \(Y\) and \(D\), captured by \(\theta_0\), and \(g_0(X) + \zeta\) is the stochastic error, partly explained by covariates \(X\). \(V\) and \(\zeta\) are stochastic errors that are not explained by \(X\). Since \(Y\) and \(D\) are jointly determined, we need an external factor, commonly referred to as an instrument, \(Z\), to create exogenous variation in \(D\). Note that \(Z\) should affect \(D\). The \(X\) here serve again as confounding factors, so we can think of variation in \(Z\) as being exogenous only conditional on \(X\).
A simple contextual example is from biostatistics permutt1989, where \(Y\) is a health outcome and \(D\) is an indicator of smoking. Thus, \(\theta_0\) captures the effect of smoking on health. Health outcome \(Y\) and smoking behavior \(D\) are treated as being jointly determined. \(X\) represents patient characteristics, and \(Z\) could be a doctor's advice not to smoke (or another behavioral treatment) that may affect the outcome \(Y\) only through shifting the behavior \(D\), conditional on characteristics \(X\).
We consider estimation of average treatment effects when treatment effects are fully heterogeneous, i.e., the response curves under control and treatment can be different nonparametric functions, 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).
Without unconfoundedness/conditional exogeneity, these quantities measure association, and could be referred to as average predictive effects (APE) and average predictive effect for the exposed (APEX). Inferential results for these objects would follow immediately from Theorem (ref).
Another common target parameter is the average treatment effect for the treated (ATTE)
In business applications, the ATTE is often the main interest, as it captures the treatment effect for those who have been affected by the treatment. A difference of the ATTE from the ATE might arise if the characteristics of the treated individuals differ from those of the general population.
The confounding factors \(X\) affect the policy variable via the propensity score \(m_0(X)\) and the outcome variable via the function \(g_0(X)\). Both of these functions are unknown and potentially complex, and we can employ ML methods to learn them.
Taking up the 401(k) example from Section (ref), the general idea for identification of \(\theta_0\) using the IRM is similar. Once we are able to account for all confounding variables \(X\) in our analysis, we can consistently estimate the causal parameter \(\theta_0\). A difference to the PLR refers to assumptions on the functional form of the main regression equation in (ref) and (ref), respectively. Whereas it is assumed that the effect of \(D\) on \(Y\) in the PLR model is additively separable, the IRM model comes with less restrictive assumptions. For example, it is possible that treatment effects are heterogeneous, i.e., vary across the population.
We consider estimation of the local average treatment effect (LATE) with a binary treatment variable \(D\in\{0,1\}\), and a binary instrument, \(Z\in\{0,1\}\). As before, \(Y\) denotes the outcome variable, and \(X\) is the vector of covariates. In a setting where unobserved factors drive the take-up of the treatment \(D\), the average treatment effect is no longer identified. However, if a valid instrumental variable is available that changes individuals' decision to take up the treatment, it is possible to identify the LATE. The LATE measures the average causal effect for the subgroup of compliers, i.e., those individuals who receive the treatment only if the instrument takes value \(Z = 1\). Hence, the LATE is of interest in many studies, where the treatment assignment cannot be assumed to be conditionally independent. For a more detailed treatment of the LATE and the key assumptions required for its identification, we would like to refer to LATE, mixtape and mostlyharmless.
The structural equation for the IIVM is
Consider the functions \(g_0\), \(r_0\), and \(m_0\), where \(g_0\) maps the support of \((Z,X)\) to \(\mathbb{R}\) and \(r_0\) and \(m_0\) map the support of \((Z,X)\) and \(X\) to \((\epsilon, 1-\epsilon)\) for some \(\epsilon \in (0, 1/2)\), such that
We are interested in estimating
Under the well-known assumptions of LATE, \(\theta_0\) is the LATE -- the average causal effect for compliers, in other words, those observations that would have \(D=1\) if \(Z\) were \(1\) and would have \(D=0\) if \(Z\) were \(0\).
In the smoking example from Section (ref), the setting is similar to the section before, but now the binary treatment variable (“smoking”) is endogenous and is instrumented by a binary instrument variable \(Z\) (“doctor's advice”). In this example, the group of compliers would comprise those individuals who quit smoking once their doctor advises them to do so and would otherwise continue to smoke. Similar to the comparison of the IRM model and the PLR model, the IIVM model does not impose the assumptions of linearity and additive separability that are maintained in the PLIV.
Here we provide an intuitive discussion of how double machine learning works in the first model, the partially linear regression model. Naive application of machine learning methods directly to Equations (ref)--(ref) may have a very high bias. Indeed, it can be shown that small biases in estimation of \(g_0\), which are unavoidable in high-dimensional estimation, create a bias in the naive estimate of the main effect, \(\hat{\theta}_{0}^{\text{naive}}\), which is sufficiently large to cause failure of conventional inference. The left panel in Figure (ref) illustrates this phenomenon. The histogram presents the empirical distribution of the studentized estimator, \(\hat{\theta}_{0}^{\text{naive}}\), as obtained in \(1000\) independent repetitions of the data generating process presented in Section (ref). The functions \(g_0\) and \(m_0\) in the PLR model are estimated with random forest learners and corresponding predictions are then plugged into a non-orthogonal score function. The regularization performed by the random forest learner leads to a bias in estimation of \(g_0\) and \(m_0\). Due to non-orthogonality of the score, this translates into a considerable bias of the main estimator \(\hat{\theta}_{0}^{\text{naive}}\): The distribution of the studentized estimator \(\hat{\theta}_{0}^{\text{naive}}\) is shifted to the right of the origin and differs substantially from a normal distribution that would be obtained if the regularization bias was negligible as shown by the red curve.
The PLR model above can be rewritten in the following residualized form
The variables \(W\) and \(V\) represent original variables after taking out or partialling out the effect of \(X\). Note that \(\theta_0\) is identified from this equation if \(V\) has a non-zero variance.
Given identification, double machine learning for a PLR proceeds as follows
The reason we work with this residualized form is that it eliminates the bias arising from solving the prediction problems in stage (1). The estimates \(\hat{\ell}_0\) and \(\hat{m}_0\) carry a regularization bias due to having to solve prediction problems well in high-dimensions. However, the nature of the estimating equation for \(\theta_0\) are such that these biases are eliminated to the first order, as explained below. This results in a high-quality low-bias estimator \(\tilde{\theta}_0\) of \(\theta_0\), as illustrated in the right panel of Figure (ref). The estimator is adaptive in the sense that the first stage estimation errors do not affect the second stage errors.
Our goal is to construct high-quality point and interval estimators for \(\theta_0\) when \(X\) is high-dimensional and we employ machine learning methods to estimate the nuisance functions such as \(g_0\) and \(m_0\). Example ML methods include lasso, random forests, boosted trees, deep neural networks, and ensembles or aggregated versions of these methods.
We shall use a method-of-moments estimator for \(\theta_0\) based upon the empirical analog of the moment condition
where we call \(\psi\) the score function, \(W = (Y,D,X,Z)\), \(\theta_0\) is the parameter of interest, and \(\eta\) denotes nuisance functions with population value \(\eta_0\).
The first key input of the inference procedure is using a score function \(\psi(W; \theta; \eta)\) that satisfies (ref), with \(\theta_0\) being the unique solution, and that obeys the Neyman orthogonality condition
Neyman orthogonality (ref) ensures that the moment condition (ref) used to identify and estimate \(\theta_0\) is insensitive to small pertubations of the nuisance function \(\eta\) around \(\eta_0\). The derivative \(\partial_{\eta}\) denotes the pathwise (Gateaux) derivative operator.
In general, it is important to distinguish whether machine learning methods are used for prediction or in the context of statistical inference. An accurate prediction rule for the nuisance parameters \(\eta_0\) does not necessarily lead to a consistent estimator for the causal parameter \(\theta_0\). Replacing the true value of \(\eta_0\) by an ML estimator \(\hat{\eta}_0\) likely introduces a bias, for example, due to heavy regularization in high-dimensional settings. If this bias is not taken into account, the estimator \(\hat{\theta}_0\) will generally be inconsistent and not have an asymptotically normal distribution. Using a Neyman-orthogonal score makes estimation of the causal parameter \(\theta_0\) robust against first order biases that arise from regularization. The Neyman orthogonality property is responsible for the adaptivity of the DML estimator -- namely, the approximate distribution of \(\hat{\theta}_0\) will not depend on the fact that the estimate \(\hat{\eta}_0\) contains error, if the latter is mild. Other approaches, as targeted maximum likelihood and semiparametric sieves estimation recognize this as well. For a more detailed treatment of Neyman orthogonality we refer to dml2018.
The right panel of Figure (ref) presents the empirical distribution of the studentized DML estimator \(\tilde{\theta}_0\) that is based on an orthogonal score. Note that estimation is performed on the identical simulated data sets and with the same machine learning method as for the naive learner, which is displayed in the left panel. The histogram of the studentized estimator \(\tilde{\theta}_0\) illustrates the favorable performance of the double machine learning estimator, which is based on an orthogonal score: The DML estimator is robust to the bias that is generated by regularization. The estimator is approximately unbiased, is concentrated around \(0\) and the distribution is well-approximated by the normal distribution.
PLR score: In the PLR model, we can employ two alternative score functions. We will shortly indicate the option for initialization of a model object in \pkg{DoubleML} to clarify how each score can be implemented. Using the option \code{score\ =\ "partialling\ out"} leads to estimation of the score function
where \(W=(Y,D,X)\) and \(\ell\) and \(m\) are \(P\)-square-integrable functions mapping the support of \(X\) to \(\mathbb{R}\), whose true values are given by
Alternatively, it is possible to use the following score function for the PLR via the option \code{score\ =\ "IV-type"}
with \(g\) and \(m\) being \(P\)-square-integrable functions mapping the support of \(X\) to \(\mathbb{R}\) with values given by
The scores above are Neyman-orthogonal by elementary calculations. Now, it is possible to see the connections to the residualized system of equations presented in Section (ref).
PLIV score: In the PLIV model, we can employ two alternative score functions. Using the option \code{score\ =\ "partialling\ out"} leads to estimation of the score function
where \(W=(Y,D,X,Z)\) and \(\ell\), \(m\), and \(r\) are \(P\)-square integrable functions mapping the support of \(X\) to \(\mathbb{R}\), whose true values are given by
Alternatively, it is possible to use the following score function for the PLIV via the option \code{score\ =\ "IV-type"}
with \(g\) and \(m\) being \(P\)-square-integrable functions mapping the support of \(X\) to \(\mathbb{R}\) with values given by
IRM score: For estimation of the ATE parameter of the IRM model, we employ the score (\code{score\ =\ "ATE"})
where \(W=(Y,D,X)\) and \(g\) and \(m\) map the support of \((D,X)\) to \(\mathbb{R}\) and the support of \(X\) to \((\epsilon, 1-\epsilon)\), respectively, for some \(\epsilon \in (0, 1/2)\), whose true values are given by
This orthogonal score is based on the influence function for the mean for missing data from robins1995. For estimation of the ATTE parameter in the IRM, we use the score (\code{score\ =\ "ATTE"})
where \(p_0=\Prob(D=1)\). Note that this score does not require estimating \(g_0(1,X)\).
IIVM score: To estimate the LATE parameter in the IIVM, we will use the score (\code{score\ =\ "LATE"})
where \(W=(Y,D,X,Z)\) and the nuisance parameter \(\eta=(g, m, r)\) consists of \(P\)-square integrable functions \(g\), \(m\), and \(r\), with \(g\) mapping the support of \((Z,X)\) to \(\mathbb{R}\) and \(m\) and \(r\), respectively, mapping the support of \((Z,X)\) and \(X\) to \((\epsilon, 1-\epsilon)\) for some \(\epsilon \in (0,1/2)\).
The second key input is the use of high-quality machine learning estimators for the nuisance parameters.
For instance, in the PLR model with \code{score = "IV-type"}, we need to have access to consistent estimators of \(g_0\) and \(m_0\) with respect to the \(L^2(P)\) norm \(\lVert \cdot \lVert_{P,2}\), such that
In the PLIV model, the sufficient condition is
These conditions are plausible for many ML methods. Different structured assumptions on \(\eta_0\) lead to the use of different machine-learning tools for estimating \(\eta_0\) as listed in dml2018:
For most of these ML methods, performance guarantees are available that make it possible to satisfy the theoretical requirements. For deep learning results can be found in farrell2021deep, for lasso in buhlmann2011statistics. Moreover, if \(\eta_0\) can be well approximated by at least one model mentioned in the list above, ensemble or aggregated methods Wolpert1992, Breiman1996 can be used. Ensemble and aggregation methods ensure that the performance guarantee is approximately no worse than the performance of the best method van2007super, dudoit2005asymptotics.
The third key input is to use a form of sample splitting at the stage of producing the estimator of the main parameter \(\theta_0\), which allows to avoid biases arising from overfitting.
Biases arising from overfitting could result from using highly complex fitting methods such as boosting, random forests, ensemble, and hybrid machine learning methods. We specifically use cross-fitted forms of the empirical moments, as detailed below in Algorithms 1 and 2, in estimation of \(\theta_0\). If the same samples would be used to estimate \(\eta_0\) and the causal parameter \(\theta_0\), we may end up with very large bias, which we refer to as an overfitting bias. While sample splitting is key for the DML approach, other approaches, like target maximum likelihood, allow for the use of arbitrary machine learning methods for the estimation of the nuisance parameters without sample splitting. The overfitting bias is illustrated in Figure (ref). The left panel shows the histogram of a studentized estimator \(\hat{\theta}^{\text{nosplit}}_0\) with \(\hat{\theta}^{\text{nosplit}}_0\) being obtained from solving the orthogonal score of Equation (ref) without sample splitting. All observations are used to learn functions \(g_0\) and \(m_0\) in the PLR model and to solve the score \(\frac{1}{N}\sum_i^{N} \psi(W_i; \hat\theta^{\text{nosplit}}_0, \hat{\eta}_0)\). Consequently, this overfitting bias leads to a considerable shift of the empirical distribution to the left. The double machine learning estimator underlying the histogram in the right panel is obtained with cross-fitting according to Algorithm 2. The sample-splitting procedure makes it possible to completely eliminate the bias induced by overfitting.
We assume that we have a sample \((W_i)^N_{i_1}\), modeled as i.i.d. copies of \(W=(Y,D,Z,X)\), whose law is determined by the probability measure \(P\). We assume that \(N\) is divisible by \(K\) in order to simplify the notation. Let \(\E_N\) denote the empirical expectation
Both Algorithm 1 and 2 use out-of-sample predictions generated by ML learners in order to solve an orthogonal moment condition and, hence, share the same steps (1) and (2). However, the algorithms differ in the way the nuisance predictions are plugged into the score function and in the subsequent solution for \(\theta_0\). In Algorithm 1, the score is solved on each of the \(K\) folds and the estimate \(\tilde{\theta}_0\) is obtained by averaging the \(K\) preliminary estimators, \(\check{\theta}_{0,k}\) with \(k=1, \ldots, K\). According to Algorithm 2, the out-of-sample predictions \(\hat{\eta}_0\) are all plugged into one score function, which is then solved to obtain the estimate \(\tilde{\theta}_0\).
The score for the models PLR, PLIV, IRM and IIVM are linear in \(\theta\), having the form
hence the estimator \(\tilde{\theta}_{0,k}\) for DML2 (\(\check{\theta}_{0,k}\) for DML1) takes the form
The linear score function representations of the PLR, PLIV, IRM and IIVM are
PLR with \code{score\ =\ "partialling\ out"}
PLR with \code{score\ =\ "IV-type"}
PLIV with \code{score\ =\ "partialling\ out"}
PLIV with \code{score\ =\ "IV-type"}
IRM with \code{score\ =\ "ATE"}
IRM with \code{score\ =\ "ATTE"}
IIVM with \code{score\ =\ "LATE"}
In Step (2) of the Algorithm DML1 and DML2, the estimator \(\hat{\eta}_{0,k}\) can generally be an ensemble or aggregation of several estimators as long as we only use the data \((W_i)_{i\not\in I_k}\) outside the \(k\)-th fold to construct the estimators.
We have found that \(K=4\) or \(K=5\) to work better than \(K=2\) in a variety of empirical examples and in simulations. The default for the option \code{n_folds} that implements the value of \(K\) is \code{n_folds=5}. Moreover, we generally recommend to repeat the estimation procedure multiple times and use the estimates and standard errors as aggregated over multiple repetitions as described in dml2018. This aggregation will be automatically executed if the number of repetitions \code{n_rep} is set to a value larger than 1.
The properties of the estimator are as follows.
The presented double machine learning method was developed in dml2018. The idea of using property (ref) to construct estimators and inference procedures that are robust to small mistakes in nuisance parameters can be traced back to neyman and has been used explicitly or implicitly in the literature on debiased sparsity-based inference belloni2011, belloni2014pivotal, javanmard2014hypothesis, van2014asymptotically, zhang2014, chernozhukov2015valid as well as (implicitly) in the classical semi-parametric learning theory with low-dimensional \(X\) levit1975, hasminskii1978, bickel1993efficient, newey1994asymptotic, van2000asymptotic, van2011targeted. These references also explain that if we use scores \(\psi\) that are not Neyman-orthogonal in high dimensional settings, then the resulting estimators of \(\theta_0\) are not \(1/\sqrt{N}\) consistent and are generally heavily biased.
Sample splitting has been used in the traditional semiparametric estimation literature to establish good properties of semiparametric estimators under weak conditions klaassen1987, schick1986, van2000asymptotic, zheng2011cross. In sparse learning problems with high-dimensional \(X\), sample splitting was employed in belloni2012sparse. There and here, the use of sample splitting results in weak conditions on the estimators of nuisance parameters, translating into weak assumptions on sparsity in the case of sparsity-based learning.
The presented approach builds upon and generalizes the approach of belloni2011, zhang2014, javanmard2014hypothesis, javanmard2014hypothesis, javanmard2018debiasing, restud, bck2014, buhlmann2015high, which considered estimation of the special case ((ref))--((ref)) using lasso without cross-fitting. This generalization, by relying upon cross-fitting, opens up the use of a much broader collection of machine learning methods and, in the case the lasso is used to estimate the nuisance functions, allows relaxation of sparsity conditions. All of these approaches can be seen as “debiasing” the estimation of the main parameter by constructing, implicitly or explicitly, score functions that satisfy the exact or approximate Neyman orthogonality.
In addition to estimation of target causal parameters, standard errors, and confidence intervals, the package \pkg{DoubleML} provides methods to perform valid simultaneous inference based on a multiplier bootstrap procedure introduced in cck2013 and cck2014 and suggested in high-dimensional linear regression models in bck2014. Accordingly, it is possible to (i) construct simultaneous confidence bands for a potentially large number of causal parameters and (ii) adjust $p$ values in a test of multiple hypotheses based on the inferential procedure introduced above.
We consider a causal PLR with \(p_1\) causal parameters of interest \(\theta_{0,1}, \ldots, \theta_{0,p_1}\) associated with the treatment variables \(D_1, \ldots, D_{p_1}\). The parameter of interest \(\theta_{0,j}\) with \(j=1, \ldots, p_1\) solves a corresponding moment condition
as for example considered in zestim. To perform inference in a setting with multiple target coefficients \(\theta_{0,j}\), the double machine learning procedure implemented in \pkg{DoubleML} iterates over the target variables of interest. During estimation of the effect of treatment \(D_j\) on \(Y\) as measured by the coefficient \(\theta_{0,j}\), the remaining treatment variables enter the nuisance terms by default, i.e., they are added to the set of control variables \(X\).
The multiplier bootstrap procedure of cck2013 and cck2014 is computationally efficient because it does not require resampling and reestimation of the causal parameters. Instead, it is sufficient to introduce a random pertubation of the score \(\psi\) and solve for \(\theta_0\), accordingly.
To construct simultaneous \((1-\alpha)\)-confidence bands, the multiplier bootstrap presented in Algorithm 4 can be used to obtain a constant \(c_{1-\alpha}\) that will guarantee asymptotic \((1-\alpha\)) coverage
The constant \(c_{1-\alpha}\) is obtained in two steps.
Moreover, it is possible to derive an adjustment method for $p$ values obtained from a test of multiple hypotheses, including classical adjustments such as the Bonferroni correction as well as the Romano-Wolf stepdown procedure rw2005, rw22005. The latter is implemented according to the algorithm for adjustment of $p$ values as provided in rw2016 and adapted to high-dimensional linear regression based on the lasso in bach2018valid.
In this section, we briefly provide information on the implementation details such as the class structure, the data backend and the use of machine learning methods. Section (ref) provides a demonstration of \pkg{DoubleML} in real-data and simulation examples. More information on the implementation can be found in the DoubleML User Guide, that is available online at \url{https://docs.doubleml.org/stable/}. All class methods are documented in the documentation of the corresponding class, which can be browsed online at \url{https://docs.doubleml.org/r/stable/} or, for example, by using the commands \code{help(DoubleML)}, \code{help(DoubleMLPLR)}, or \code{help(DoubleMLData)} in \proglang{R}. For an introduction to \pkg{R6} we refer to the introduction of the online book for \pkg{mlr3}, available at \url{https://mlr3book.mlr-org.com/intro.html}.
As pointed out in the previous sections, the double machine learning framework provides a general inferential framework in that it covers a plethora of causal models that can be characterized in terms of a Neyman-orthogonal score function $\psi$. In order to design an implementation that is similarly general, the implementation of \pkg{DoubleML} for \proglang{R} is based on object orientation as enabled by the the \pkg{R6} package R6. The choice of the object orientation provided by \pkg{R6} as compared to alternative approaches (e.g., \proglang{S}3 or \proglang{S}4 classes) has been motivated by mainly three reasons: First, we would like to obtain an optimal compatibility with the \pkg{mlr3} ecosystem that is built with \pkg{R6} classes as well. Second, \pkg{R6} makes it possible to use encapsulation, inheritance, active bindings and to distinguish between private and public methods which are important features required in our implementation. Third, the object-oriented implementation of \pkg{DoubleML} makes it possible to achieve a high degree of comparability with its \proglang{Python} twin, which will likely facilitate and accelerate the continuous development of both packages in the future. For an introduction to object orientation in \proglang{R} and the \pkg{R6} package, we refer to the vignettes of the \pkg{R6} package that are available online at \url{https://r6.r-lib.org/articles/}, Chapter 2.1 of mlr3book, and the chapters on object orientation in advr. The structure of the classes are presented in Figure (ref). Moreover, the flowchart in Figure (ref) illustrates the main steps of an analysis in \pkg{DoubleML} and links them to the provided object classes. Figure (ref) provides a short code demonstration, too. The abstract class `\code{DoubleML}' provides all methods for estimation and inference, for example the methods \code{fit()}, \code{bootstrap()}, \code{confint()}. All key components associated with estimation and inference are implemented in `\code{DoubleML}', for example the sample splitting, the implementation of Algorithm 1 (DML1) and Algorithm 2 (DML2), the estimation of the causal parameters, and the computation of the scores \(\psi(W;\theta, \eta)\). Only the model-specific properties and methods are allocated at the classes `\code{DoubleMLPLR}' (implementing the PLR), `\code{DoubleMLPLIV}' (PLIV), `\code{DoubleMLIRM}' (IRM), and `\code{DoubleMLIIVM}' (IIVM). For example, each of the models has one or several Neyman-orthogonal score functions that are implemented for the specific child classes.
The `\code{DoubleMLData}' class serves as the data backend and implements the causal model of interest. The user is required to specify the roles of the variables in a data set at hand. Depending on the causal model considered, it is necessary to declare the dependent variable, the treatment variable(s), confounding variables(s), and, in the case of instrumental variable regression, one or multiple instruments. The data backend can be initialized from a \code{data.table} datatable. \pkg{DoubleML} provides wrappers to initialize from `\code{data.frame}' and `\code{matrix}' objects, as well.
Generally, all learners provided by the packages \pkg{mlr3}, \pkg{mlr3learners} and \pkg{mlr3extralearners} can be used for estimation of the nuisance functions of the structural models presented above. An interactive list of supported learners is available at the \pkg{mlr3extralearners} website (\url{https://mlr3extralearners.mlr-org.com/articles/learners/list_learners.html}). The \pkg{mlr3extralearners} package makes it possible to add new learners, as well. The performance of the double machine learning estimator \(\tilde\theta_0\) will depend on the predictive quality of the used machine learning method. Machine learning methods usually have several (hyper-)parameter that need to be adapted to the specific application. Tuning of model parameters can be either performed externally or internally. The latter is implemented in the method \code{tune()} and is further illustrated in an example in Section (ref). Both cases build on the functionalities provided by the package \pkg{mlr3tuning}.
The flexible architecture of the \pkg{DoubleML} package allows users to modify the estimation procedure in many regards. Among others, users can provide customized sample splitting rules after initialization of the causal model via the method \code{set_sample_splitting()}. Moreover, it is possible to adjust the Neyman-orthogonal score function by externally providing a customized function via the \code{score} option during initialization of the causal model object. Short examples for both of these potential extensions are presented in Section (ref).
In this section, we will first demonstrate the use of \pkg{DoubleML} in a real-data example, which is based on data from the Pennsylvania Reemployment Bonus experiment bilias2000. This empirical example has been used in dml2018, as well. The goal in the empirical example is to estimate the causal parameter in a partially linear and an interactive regression model. We further provide a short example on how valid simultaneous inference can be performed with \pkg{DoubleML}. Finally, we present results from a short simulation study as a brief assessment of the finite-sample performance of the implemented estimators. Here we want to stress that in real world applications modelling choices of the estimation of the nuisance parameters and proper tuning of the parameters are very important. We would like to mention that the presented examples are mainly included for the purpose of illustration. In practice, we recommend to carefully choose and tune the ML learners in terms of their hyperparamaters.
We begin our real-data example by downloading the Pennsylvania Reemployment Bonus data set. To do so, we use the call \code{fetch_bonus()} (a connection to the internet is required).
Load data as \code{data.table}.
The output is suppressed for the sake of brevity.
The data backend `\code{DoubleMLData}' can be initialized from a `\code{data.table}' object by specifying the dependent variable \(Y\) via a character in \code{y_col}, the treatment variable(s) \(D\) in \code{d_cols}, and the confounders \(X\) via \code{x_cols}. Moreover, in IV models, an instrument can be specified via \code{z_cols}. In the next step, we assign the roles to the variables in the data set: \code{y_col\ =\ \textquotesingleinuidur1\textquotesingle} serves as outcome variable \(Y\), the column \code{d_cols\ =\ \textquotesingletg\textquotesingle} serves as treatment variable \(D\) and the columns \code{x_cols} specify the confounders.
The data backend can be printed to obtain a summary of the main attributes of the `\code{DoubleMLData}' object.
Print the data set (output suppressed).
To initialize an instance of the class `\code{DoubleMLData}' from a `\code{data.frame}' or a collection of `\code{matrix}' objects, \pkg{DoubleML} provides the convenient wrappers \code{double_ml_data_from_}- \code{data_frame()} and \code{double_ml_data_from_matrix()}. Although the data backend does not provide a formula interface, `\code{DoubleMLData}' objects can be initialized from a `\code{model.matrix}' object. The following example demonstrates how users may proceed to specify the causal model by using a `\code{formula}'. We load the Pennsylvania Reemployment Bonus data set as a `\code{data.frame}' and replicate a flexible model specification used in the empirical analysis of dml2018. To flexibly model the nuisance function, we generate all two-way interactions of the control variables.
Load the data as a `\code{data.frame}'.
Print the names of the variables.
Specify a `\code{formula}' object to generate all two-way interactions of the control variables.
Create a `\code{model.matrix}' based on the `\code{formula}' object.
Initialize using the wrapper \code{double_ml_data_from_data_frame()}.
Alternatively, initialize via the wrapper \code{double_ml_data_from_matrix()}.
To initialize a PLR model, we have to provide a learner for each nuisance part in the model in Equation (ref)--(ref). In \proglang{R}, this is done by providing learners to the arguments \code{ml_m} for nuisance part \(m\) and \code{ml_l} for nuisance part \(\ell\). We can pass a learner as instantiated in \pkg{mlr3} and \pkg{mlr3learners}, for example a random forest as provided by the \proglang{R} package \pkg{ranger} ranger. Previous installation of \pkg{ranger} is required. Moreover, we can specify the score (allowed choices for PLR are \code{"partialling\ out"} or \code{"IV-type"}) and the algorithm via the option \code{dml_procedure} (allowed choices \code{"dml1"} and \code{"dml2"}) . Optionally, it is possible to change the number of folds used for sample splitting through \code{n_folds} and the number of repetitions via \code{n_rep}, if the sample splitting and estimation procedure should be repeated.
Set a seed for replication of the sample split.
To perform estimation, call the \code{fit()} method. The output can be summarized using the method \code{summary()}.
Hence, there is evidence to reject the null hypothesis that \(\theta_{0, tg}=0\) at the 5% significance level. The estimated coefficient and standard errors can be accessed via the attributes \code{coef} and \code{se} of the object \code{doubleml_bonus}.
After completed estimation, we can access the resulting score \(\psi(W_i; \tilde{\theta}_0, \hat{\eta}_0)\) or the components \(\psi_a(W_i; \hat{\eta}_0)\) and \(\psi_b(W_i; \hat{\eta}_0)\). The estimated score for the first 5 observations can be obtained via the public field \code{psi}. \code{psi} is an \code{array} with \code{dim = c(n_obs, n_rep, n_treat)} with
Similarly, the components of the score \(\psi_a(W_i; \hat{\eta}_0)\) and \(\psi_b(W_i; \hat{\eta}_0)\) are available as fields.
To construct a \((1-\alpha)\) confidence interval, we use the \code{confint()} method.
The treatment variable \(D\) in the Pennsylvania Reemployment Bonus example is binary. Accordingly, it is possible to estimate an IRM model. Since the IRM requires estimation of the propensity score \(\Prob(D\mid X)\), we have to specify a classifier for the nuisance part \(m_0\).
Initialize a classifier for estimation of the propensity score and create a new instance of a causal model, here an IRM.
The output is suppressed for the sake of brevity.
To perform estimation, call the \code{fit()} method. The output can be summarized using the method \code{summary()}.
The estimated coefficient is very similar to the estimate of the PLR model and our conclusions remain unchanged.
We consider a simulated example of a PLR model to illustrate the use of methods for simultaneous inference. First, we will generate a sparse linear model with only three variables having a non-zero effect on the dependent variable.
Generate a \code{data.frame} and use the corresponding wrapper.
We use the wrapper \code{double_ml_data_from_data_frame()} to specify a data backend that assigns the first 10 columns of \(X\) as treatment variables and declares the remaining columns as confounders.
The output is suppressed for the sake of brevity.
A sparse setting suggests the use of the lasso learner. Here, we use the lasso estimator with cross-validated choice of the penalty parameter \(\lambda\) as provided in the \pkg{glmnet} package for \proglang{R} glmnet.
Output messages during fitting are suppressed.
The multiplier bootstrap procedure can be executed using the \code{bootstrap()} method where the option \code{method} specifies the choice of the random pertubations and \code{n_rep_boot} the number of bootstrap repetitions.
The resulting bootstrapped $t$ statistics are available via the field \code{boot_t_stat}. To construct a simultaneous confidence interval, we set the option \code{joint\ =\ TRUE} when calling the \code{confint()} method.
The correction of the $p$ values of a joint hypotheses test on the considered causal parameters is implemented in the method \code{p_adjust()}. By default, the adjustment procedure specified in the option \code{method} is the Romano-Wolf stepdown procedure.
Alternatively, the correction methods provided in the \pkg{stats} function \code{p.adjust} can be applied, for example the Bonferroni, Bonferroni-Holm, or Benjamini-Hochberg correction. For example a Bonferroni correction could be performed by specifying \code{method\ =\ "bonferroni"}.
The performance of the final double machine learning estimator depends on the predictive performance of the underlying ML method. First, we briefly show how externally tuned parameters can be passed to the learners in \pkg{DoubleML}. Second, it is demonstrated how the parameter tuning can be done internally by \pkg{DoubleML}.
Section 3 of the \pkg{mlr3} book mlr3book provides a step-by-step introduction to the powerful tuning functionalities of the \pkg{mlr3tuning} package. Accordingly, it is possible to manually reconstruct the \pkg{mlr3} regression and classification problems, which are internally handled in \pkg{DoubleML}, and to perform parameter tuning accordingly. One advantage of this procedure is that it allows users to fully exploit the powerful benchmarking and tuning tools of \pkg{mlr3} and \pkg{mlr3tuning}.
Consider the sparse regression example from above. We will briefly consider a setting where we explicitly set the parameter \(\lambda\) for a \pkg{glmnet} estimator rather than using the interal cross-validated choice with \code{cv_glmnet}.
Suppose for simplicity, some external tuning procedure resulted in an optimal value of \(\lambda=0.1\) for nuisance part \(m\) and \(\lambda=0.09\) for nuisance part \(\ell\) for the first treatment variable and \(\lambda=0.095\) and \(\lambda=0.085\) for the second variable, respectively. After initialization of the model object, we can set the parameter values using the method \code{set_ml_nuisance_params()}.
To set the values, we have to specify the treatment variable and the nuisance part. If no values are set, the default values are used. Note that variable names are overwritten by the wrapper for the matrix interface.
All externally specified parameters can be retrieved from the field \code{params}. The output is omitted for the sake of brevity.
An alternative to external tuning and parameter provisioning is to perform the tuning internally. The advantage of this approach is that users do not have to specify the underlying prediction problems manually. Instead, \pkg{DoubleML} uses the underlying data backend to ensure that the machine learning methods are tuned for the specific model under consideration and, hence, to possibly avoid mistakes. We initialize our structural model object with the learner. At this stage, we do not specify any parameters.
Load required packages for tuning and set logger to omit messages during tuning and fitting.
To perform parameter tuning, we provide a grid of values used for evaluation for each of the nuisance parameters. To set up a grid of values, we specify a named list with names corresponding to the learner names of the nuisance part (see method \code{learner_names()}). The elements in the list are objects of the class `\code{ParamSet}' of the \pkg{paradox} package paradox.
The hyperparameter tuning is performed according to options passed through a named list \code{tune_settings}. The entries in the list specify options during parameter tuning with \pkg{mlr3tuning}:
In the next code chunk, the value of the parameter \(\lambda\) is tuned via grid search in the range 0.05 to 0.1 at a resolution of 11. The resulting grid has 11 equally spaced values ranging from a minimum value of 0.05 to a maximum value of 0.1. Type \code{generate_design_grid(par_grids\$ml_l,\ resolution\ =\ 11)} to access the grid for nuisance function \code{ml_l}. To evaluate the predictive performance in both nuisance functions, the cross-validated mean squared error is used.
Provide tune settings.
With these parameters we can run the tuning by calling the \code{tune()} method for `\code{DoubleML}' objects.
Execution might take around 50 seconds.
Output omitted for the sake of brevity, available in the appendix. Access tuning results for target variable \code{X1}.
Access tuned parameters (output suppressed).
Estimate model and call the \code{summary()} method.
By default, the parameter tuning is performed on the whole sample, for example in the case of \(K_{\text{tune}}\)-fold cross-validation, the entire sample is split into \(K_{\text{tune}}\) folds for evaluation of the cross-validated error. Alternatively, each of the \(K\) folds used in the cross-fitting procedure could be split up into \(K_{\text{tune}}\) subfolds that are then used for evaluation of the candidate models. As a result, the choice of the tuned parameters will be fold-specific. To perform fold-specific tuning, users can set the option \code{tune_on_folds\ =\ TRUE} when calling the method \code{tune()}.
The flexible architecture of the \pkg{DoubleML} package allows users to modify the estimation procedure in many regards. We will shortly present two examples on how users can adjust the double machine learning framework to their needs in terms of the sample splitting procedure and the score function.
By default, \pkg{DoubleML} performs cross-fitting as presented in Algorithms 1 and 2. Alternatively, all implemented models allow a partition to be provided externally via the method \code{set_sample_splitting()}. Note that by setting \code{draw_sample_splitting\ =\ FALSE} one can prevent that a partition is drawn during initialization of the model object. The following calls are equivalent. In the first sample code, we use the standard interface and draw the sample-splitting with \(K=4\) folds during initialization of the `\code{DoubleMLPLR}' object.
First generate some data and initialize ML learners and a data backend.
In the second sample code, we manually specify a sampling scheme using the `\code{mlr3::Resampling}' class. Alternatively, users can provide a nested list that has the following structure:
Set up a task and cross-validation resampling scheme in \pkg{mlr3}.
Structure of the specified sampling scheme.
Fit the model and summarize.
Setting the option \code{apply_cross_fitting\ =\ FALSE} at the instantiation of the causal model allows double machine learning being performed without cross-fitting. It results in randomly splitting the sample into two parts. The first half of the data is used for the estimation of the nuisance models with the machine learning methods and the second half for estimating the causal parameter, i.e., solution of the score. Note that cross-fitting performs well empirically and is recommended to remove bias induced by overfitting. Moreover, cross-fitting allows to exploit full efficiency: Every fold is used once for training the ML methods and once for estimation of the score dml2018. A short example on the efficiency gains associated with cross-fitting is provided in Figure (ref).
Users may want to adjust the score function $\psi(W;\theta_0, \eta_0)$, for example, to adjust the DML estimators in terms of a re-weighting, e.g., to adjust for missing outcome via inverse probability of censoring weight (IPCW). An alternative to the choices provided in \pkg{DoubleML} is to pass a function via the argument \code{score} during initialization of the model object. The following examples are equivalent. In the first example, we use the score option \code{"partialling out"} for the PLR model whereas in the second case, we explicitly provide a function that implements the same score. The arguments used in the function refer to the internal objects that implement the theoretical quantities in Equation (ref).
Use score \code{"partialling out"}.
We define the function that implements the same score and specify the argument \code{score} accordingly. The function must return a named list with entries \code{psi_a} and \code{psi_b} to pass values for computation of the score.
Required input:
To illustrate the validity of the implemented double machine learning estimators, we perform a brief simulation study.
As mentioned before the use of the cross-fitting Algorithms 1 (DML1) and 2 (DML2) makes it possible to use sample splitting and exploit full efficiency at the same time. To illustrate the superior performance due to cross-fitting, we compare the double machine learning estimator with and without a cross-fitting procedure in the simulation setting that was presented in Section (ref). Figure (ref) illustrates that efficiency gains can be achieved if the role of the random partitions is swapped in the estimation procedure. Using cross-fitting makes it possible to obtain smaller standard errors for the DML estimator: The empirical distribution of the double machine learning estimator that is based on the cross-fitting Algorithm 2 (DML2) exhibits a more pronounced concentration around zero.
We provide simulation results for double machine learning estimators in the presented key causal models in Figure (ref). In a replication of the simulation example in Section (ref), we show that the confidence intervals for the DML estimator in the partially linear regression model achieves an empirical coverage (= \(0.952\)) close to the specified level of \(1-\alpha=0.95\). The estimator is, again, based on a random forest learner. The corresponding results are presented in the top-left panel of Figure (ref).
In a simulated example of a PLIV model, the DML confidence interval that is based on a lasso learner (\code{regr.cv_glmnet} of \pkg{mlr3}) achieves a coverage of 95.6%. The underlying data generating process is based on a setting considered in CHSAERpp with one instrumental variable. Moreover for simulations of the IRM model, we make use of a DGP of belloni2017program. The DGP for the IIVM is inspired by a simulation run in latest. We present the formal DGPs in the appendix. To perform estimation of the nuisance functions in the interactive models, we employ the regression and classification predictors \code{regr.cv_glmnet} and \code{classif.cv_glmnet} as provided by the \pkg{mlr3} package. In all cases, we employ the cross-validated \code{lambda.min} choice of the penalty parameter with five folds, in other words, that \(\lambda\) value that minimizes the cross-validated mean squared error. Figure (ref) shows that the empirical distribution of the centered estimators as obtained in finite sample settings is relatively well-approximated by a normal distribution. In all models the empirical coverage that is achieved by the constructed confidence bands is close to the nominal level.
To verify the finite-sample performance of the implemented methods for simultaneous inference, we perform a small simulation study in a regression setup which is similar as the one used in bach2018valid. We would like to perform valid simultaneous inference on the coefficients \(\theta\) in the regression model
with \(n=1000\) and \(p_1=42\) regressors. The errors \(\varepsilon_i\) are normally distributed with \(\varepsilon_i \sim N(0,\sigma^2)\) and variance \(\sigma^2=3\). The regressors \(d_i\) are generated by a joint normal distribution \(d_i \sim N(\mu, \Sigma)\) with \(\mu = \mathbf{0}\) and \(\Sigma_{j,k} = 0.5^{\lvert j-k\rvert}\). The model is sparse in that only the first \(s=12\) regressors have a non-zero effect on outcome \(y_i\). The \(p_1\) coefficients \(\theta_1, \ldots, \theta_{p_1}\) are generated as
for \(j=1, \ldots, s\) with \(\theta^{\max}=9\), \(\theta^{\min}=0.75\), and \(a=0.99\). All other coefficients have values exactly equal to \(0\). Estimation of the nuisance components has been performed by using the lasso as provided by \code{regr.cv_glmnet} in \pkg{mlr3}.
We report the empirical coverage as achieved by a joint \((1-\alpha)\)-confidence interval for all \(p_1=42\) coefficients and the realized family-wise error rate of the implemented $p$ value adjustments in \(R=500\) repetitions in Table (ref). The finite sample performance of the Romano-Wolf stepdown procedure that is based on the multiplier bootstrap as well as the classical Bonferroni and Bonferroni-Holm correction are evaluated. Table (ref) shows that all methods achieve an empirical FWER close to the specified level of \(\alpha = 0.1\). In all cases, the double machine learning estimators reject all 12 false null hypotheses in every repetition.
In this paper, we provide an overview on the key ingredients and the major structure of the double/debiased machine learning framework as established in dml2018 together with an overview on a collection of structural models. Moreover, we introduce the \proglang{R} package \pkg{DoubleML} that serves as an implementation of the double machine learning approach. A brief simulation study provides insights on the finite sample performance of the double machine learning estimator in the key causal models.
The structure of \pkg{DoubleML} is intended to be flexible with regard to the implemented structural models, the resampling scheme, the machine learning methods and the underlying algorithm, as well as the Neyman-orthogonal scores considered. By providing the \proglang{R} package \pkg{DoubleML} together with its \proglang{Python} twin DoubleMLpython, we hope to make double machine learning more accessible to users in practice. Finally, we would like to encourage users to add new structural models, scores and functionalities to the package.
This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -- Project Number 431701914.