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.
73,854 characters · 22 sections · 22 citation commands
Double Machine Learning based Program Evaluation under Unconfoundedness
\doublespacing
The adaptation of so-called machine learning to causal inference has been a productive area of methodological research in recent years. The resulting new methods complement the existing econometric toolbox for program evaluation along at least two dimensions \cite<see for recent overviews>{Athey2017,Athey2019MachineAbout,Abadie2018EconometricEvaluation}. On the one hand, they provide flexible methods to estimate standard average effects. In particular, they provide a data-driven approach to variable and model selection in studies that rely on an unconfoundedness assumption\footnote{Also known as exogeneity, selection on observables, ignorability, or conditional independence assumption.} for identification. On the other hand, they enable a more comprehensive evaluation by providing new methods for the flexible estimation of heterogeneous effects and of treatment assignment rules.
This paper considers Double Machine Learning (DML) Chernozhukov2018 as a framework for flexible and comprehensive program evaluation. The DML framework seems attractive because (i) it can be combined with a variety of standard supervised machine learning methods, (ii) it covers average effects for binary \cite<e.g.>{Belloni2014InferenceControls,Belloni2017,Chernozhukov2018}, multiple \cite<e.g.>{Farrell2015} as well as continuous treatments \cite<e.g.>{Kennedy2017Non-parametricEffects,Colangelo2019DoubleTreatments,Semenova2021DebiasedFunctions}, (iii) it naturally extends to the estimation of heterogeneous treatment effects of different forms like canonical subgroup effects, the best linear prediction of effect heterogeneity, or nonparametric effect heterogeneity \cite<e.g>{fan2020EstimationData,Zimmert2019NonparametricConfounding,Foster2019OrthogonalLearning,Oprescu2019OrthogonalInference,Semenova2021DebiasedFunctions,Kennedy2020OptimalEffects,Curth2021NonparametricAlgorithms}, and (iv) it can be used to estimate optimal treatment assignment rules \cite<e.g.>{Dudik2011DoublyLearning,Athey2021PolicyData,Zhou2018OfflineOptimization}. All these DML based methods have favourable statistical properties and allow the use of standard tools like t-tests, OLS, kernel regression, series regression, or supervised machine learning for estimating causal parameters of interest after flexibly adjusting for confounding.
This paper starts with a review of DML based methods, then applies these methods in a standard labour economic setting, and comes back to the methods by proposing the normalised DR-learner as a potential fix to a finite sample problem encountered in the application. Thus, it contributes to the steadily growing literature of causal machine learning for program evaluation in three ways. First, the review highlights that methods for different parameters build on the same doubly robust score. The construction of this score might be computationally expensive because it requires the estimation of outcomes and treatment probabilities via machine learning methods. However, once constructed the score can be reused to estimate a variety of interesting parameters. This paper focuses on methods that build on the doubly robust score because they allow to leverage conceptual and computational synergies. The result is a comprehensive pipeline for program evaluation within the same framework as Figure (ref) illustrates. This is currently not possible with the variety of more specialised alternatives that integrate machine learning in the estimation of average treatment effects \cite<e.g.>{vanderLaan2006TargetedLearning,Athey2018ApproximateDimensions,Avagyan2017HonestEstimation,Tan2020Model-assistedData,Ning2020RobustScore}, heterogeneous treatment effects \cite<e.g.>{Tian2014,Athey2016,Chernozhukov2017GenericExperiments,Wager2017,Athey2017a,Kunzel2017,Nie2021} and optimal treatment assignment \cite<e.g.>{Bansak2018ImprovingAssignment,Kallus2018BalancedLearning}.
Second, we use DML based methods to provide a comprehensive and computationally convenient evaluation of four programs of the Swiss Active Labour Market Policy (ALMP) in a standard dataset Lechner2020SwissDataset. The evaluation in this paper illustrates the potential of DML based methods for program evaluations under unconfoundedness and provides a potential blueprint for similar analyses. This adds to a small but steadily growing literature that applies causal machine learning to program evaluation in general \cite<e.g.>{Bertrand2017,Strittmatter2018WhatEvaluation,Gulyas2019UnderstandingApproach,Knittel2019UsingUse,Davis2020RethinkingJobs,Baiardi2021TheStudies,Farbmacher2021HeterogeneousCognition} and to evaluations based on unconfoundedness in particular \cite<e.g>{Kreif2019MachineInference,Cockx2020PriorityBelgium,Knaus2020HeterogeneousApproach,Knaus2021ASkills}.
Third, we contribute to the methodological literature on the flexible estimation of individualised treatment effects \cite<see for a recent overview>{Knaus2021} by proposing the normalised DR-learner (NDR-learner), which builds on the recent DR-learner of \citeA{Kennedy2020OptimalEffects}. The application reveals that the plain DR-learner produces few extreme effect estimates. It turns out that individualised effect estimates can be stabilised by an individualised normalisation of inverse probability weights. Thus, the NDR-learner can be considered as a generalisation of the popular \citeA{Hajek1971CommentOne} normalisation for inverse probability weighting estimators for average effects. The increased stability comes at the price that the NDR-learner limits the class of permissible machine learning methods for effect heterogeneity estimation to methods that form predictions as convex combination of outcomes (e.g. Random Forests).
Overall, we find that DML based methods provide a promising set of methods for program evaluation. The estimated average program effects are in line with the previous literature. We find that computer, vocational and language courses increase employment in the 31 months after programs start, while the effects of job search trainings are mostly negative. The heterogeneity analysis additionally reveals substantial heterogeneities by gender, nationality, previous labour market success and qualification. These are picked up by the estimated optimal assignment rules.
The paper proceeds as follows. Section (ref) defines the estimands of interest and their identification under unconfoundedness. Section (ref) reviews DML based methods for estimation and introduces the NDR-learner. Section (ref) presents the application. Section (ref) describes the implementation of the methods. Section (ref) reports the results. Section (ref) concludes. The Appendix provides additional explanations and results. The R-package \href{https://github.com/MCKnaus/causalDML}{\color{blue}causalDML} implements the applied estimators. An \href{https://mcknaus.github.io/assets/code/Notebook_DML_ALMP_MCK2020.html}{\color{blue}R notebook} replicating the analysis is provided.
We define the estimands of interest in the multiple treatment version of the potential outcomes framework Rubin1974,Imbens2000TheFunctions,Lechner2001. Let $\mathcal{W} = \{0,...,T\}$ denote a set of multiple programs and $D_i(w) = \mathds{1}(W_i = w)$ a binary variable indicating in which program individual $i$ ($i=1, ..., N$) is actually observed.\footnote{For DML based estimation with continuous treatments see, e.g \citeA{Kennedy2017Non-parametricEffects}, \citeA{Colangelo2019DoubleTreatments}, and \citeA{Semenova2021DebiasedFunctions}.} We assume that each individual has a potential outcome $Y_i(w)$ for all $w \in \mathcal{W}$. Without loss of generality, the discussion below assumes that higher outcome values are desirable.
The first estimand of interest is the average potential outcome (APO), $\gamma_w = E[Y_i(w)]$. It answers the question about the average outcome if the whole population was assigned to program $w$. However, the more interesting question is usually to compare different programs $w$ and $w'$. To this end, we take the difference of the according individual potential outcomes, $Y_i(w) - Y_i(w')$,\footnote{This would be $Y_i(1) - Y_i(0)$ in the canonical binary treatment setting.} and aggregate them to different estimands: First, the average treatment effect (ATE), $\delta_{w,w'} = E[Y_i(w) - Y_i(w')]$. Second, the average treatment effect on the treated (ATET), $\theta_{w,w'} = E[Y_i(w) - Y_i(w') \mid W_i = w]$. Third, the conditional average treatment effect (CATE), $\tau_{w,w'}(z) = E[Y_i(w) - Y_i(w') \mid Z_i=z]$, where $Z_i \in \mathcal{Z}$ is a vector of observed pre-treatment variables.\footnote{We focus in this study on expectations of the individual treatment effects. DML based methods for quantile treatment effects can be found, e.g. in \citeA{Belloni2017} and \citeA{Kallus2019LocalizedBeyond}.}
The different aggregations accommodate the notion that treatment effects might be heterogeneous. ATE represents the average effect in the population, while ATET shows it for the subpopulation that is actually observed in program $w$. Thus, the comparison of ATE and ATET can be informative about the quality of the program assignment mechanism. For example, ATET being larger than ATE indicates that the observed program assignment is better than random.
The ATET is defined by the observed program assignment and thus not subject to the choice of the researcher. In contrast, the conditioning variables $Z_i$ of the CATE are specified by the researcher to investigate potentially heterogeneous effects across the groups of individuals that are defined by different values of $Z_i$. Such heterogeneous effects can be indicative for underlying mechanisms. Further, CATEs characterise which groups win and which lose by how much by receiving program $w$ instead of $w'$.
The different average effects above provide a comprehensive evaluation of programs under the current program assignment policy. In many applications, however, we want to conclude the analysis with a recommendation how the assignment policy could be improved. This can either be done using the evidence on the different average effects defined above or by formally defining the objective of an optimal assignment rule. The latter is pursued by the literature on statistical treatment rules \cite<e.g.>[and references therein]{Manski2004StatisticalPopulations,Hirano2009AsymptoticsRules,Stoye2009MinimaxSamples,Stoye2012MinimaxExperiments,Kitagawa2018WhoChoice,Athey2021PolicyData}. Here we focus on the case with multiple treatment options as considered by \citeA{Zhou2018OfflineOptimization}.
Let $\pi(Z_i)$ be a policy that assigns individuals to programs according to their characteristics $Z_i$ or, put more formally, the function $\pi(Z_i)$ maps observable characteristics to a program: $\pi: \mathcal{Z} \rightarrow \mathcal{W}$. In principle, the policy rule can be completely flexible and in the ideal world we would assign each individual to the program with the highest conditional APO, $E[Y_i(w) \mid Z_i =z]$. However, in many cases we want to restrict the set of candidate policy rules denoted by $\Pi$ to be interpretable for the communication with decision makers or to incorporate costs or fairness constraints. Each of these candidate policy rules has a policy value function denoted by $Q(\pi) = E[Y_i(\pi(Z_i))] = E \left[ \sum_w \mathds{1}(\pi(Z_i) = w) Y_i(w) \right]$. $Q(\pi)$ quantifies the average population outcome if policy rule $\pi$ would be used to assign programs. The estimand of interest is then the optimal policy rule $\pi^*$ with the highest value function for the set of candidate policy rules, or formally $\pi^* = \operatorname*{\arg\max}_{\pi \in \Pi} Q(\pi)$.
The previous section defined the estimands of interest in terms of potential outcomes. However, each individual is only observed in one program. Thus, only one potential outcome per individual is observable and the other potential outcomes remain latent. This is the fundamental problem of causal inference Holland1986StatisticsInference and we need further assumptions to identify the estimands of interest. In this paper, we consider the unconfoundedness assumption that assumes access to a vector of pre-treatment variables $X_i \in \mathcal{X}$ containing $Z_i$ such that the following standard assumptions hold \cite<e.g.>{Imbens2015CausalSciences}:
The unconfoundedness assumption requires that $X_i$ contains all confounding variables that jointly affect program assignment and the outcome. Common support states that it must be possible to observe each individual in all programs. SUTVA rules out interference. These assumptions allow the identification of the average potential outcome (APO) conditional on confounders in three common ways:
Equation (ref) shows that the conditional APO is identified as a conditional expectation of the observed outcome. Equation (ref) shows that it is identified by reweighting the observed outcome with the inverse treatment probability. Finally, Equation (ref) adds the reweighted outcome residual to the conditional outcome representation of Equation (ref). This seems redundant because we can check that the reweighted residual has expectation zero under unconfoundedness. However, this identification result is doubly robust in the sense that it still holds if we replace either $\mu(w,x)$ or $e_w(x)$ in Equation (ref) by arbitrary functions of $x$.\footnote{Appendix (ref) reviews identification and identification double robustness of Equation (ref) for completeness.} This doubly robust structure plays a crucial role for the estimation procedures that we discuss in the next section.
From an identification perspective, $\Gamma(w,x)$ defined in Equation (ref) suffices to identify all estimands of interest stated in the previous subsection:
All Double Machine Learning (DML) based estimators for the estimands of interest build on the doubly robust scores of \citeA{Robins1994,Robins1995AnalysisData} and their Augmented Inverse Probability Weighting (AIPW) estimator in particular. In the following, large Greek letters denote the scores corresponding to the small Greek letters used to define the estimands in Section (ref).
The construction of the doubly robust scores requires the input of so-called nuisance parameters that are usually of secondary interest and considered as tool to eventually obtain the parameters of interest. In our case, the two nuisance parameters are ${\mu(w,x) = E[Y_i \mid W_i = w, X_i = x]}$ and $e_w(x) = P[W_i = w \mid X_i = x]$ for all $w$. $\mu(w,x)$ is the conditional outcome mean for the subgroup observed in program $w$. $e_w(x)$ is the conditional probability to be observed in program $w$, also known as the propensity score. Usually these functions are unknown and need to be estimated. Following \citeA{Chernozhukov2018} they are estimated based on $K$-fold cross-fitting: (i) randomly divide the sample in $K$ folds of similar size, (ii) leave out fold $k$ and estimate models for the nuisance parameters in the remaining $K-1$ folds, (iii) use these models to predict $\hat{\mu}^{-k}(w,x)$ and $\hat{e}^{-k}_w(x)$ in the left out fold $k$, and (iv) repeat (i) to (iii) such that each fold is left out once. This procedure avoids overfitting in the sense that no observation is used to predict its own nuisance parameters. To avoid notational clutter, we ignore the dependence on the specific fold in the following notation and refer to the cross-fitted nuisance parameters as $\hat{\mu}(w,x)$ and $\hat{e}_w(x)$.
The main building block of the following estimators is the doubly robust score of the APO, which replaces the true nuisance parameters in Equation (ref) by their cross-fitted predictions:
The ATE score for the comparison of treatment $w$ and $w'$ is then constructed as the difference of the respective APO scores:
The only estimator we consider that uses the same nuisance parameter but plugs them into a different score is the ATET estimator. Although the identification result with the doubly robust APO score in the previous section holds, it is not doubly robust. However, the doubly robust score for the ATET exists and is defined as
where $\hat{e}_w = N_w / N$ is the unconditional treatment probability with $N_w$ counting the number of individuals observed in program $w$ \cite<see also, e.g.>{Farrell2015}.
The estimation of the APOs, ATEs and ATETs boils down to taking the means of the previously defined doubly robust scores. For statistical inference, we can rely on standard one-sample t-tests. Thus, the score's mean and the variance of this mean are the point and the variance estimate of the respective estimand of interest:
Note that the estimated variances require no adjustment for the fact that we have estimated the nuisance parameters in a first step. The resulting estimators are consistent, asymptotically normal and semiparametrically efficient under the main assumption that the estimators of the cross-fitted nuisance parameters are consistent and converge sufficiently fast Belloni2014InferenceControls,Farrell2015,Belloni2017,Chernozhukov2018. In particular, the product of the convergence rates of the outcome and propensity score estimators must be faster than $n^{1/2}$. This allows to apply machine learning to estimate the nuisance parameters.\footnote{Further results, regularity conditions and discussions can be found in section 5.1 of \citeA{Chernozhukov2018}.} Flexible machine learning estimators converge usually slower than the parametric rate $n^{1/2}$ but several are known to be able to achieve $n^{1/4}$ and faster, which would be sufficient if both nuisance parameter estimators achieve it.\footnote{For example, versions of Lasso Belloni2013LeastModels, Boosting Luo2016High-DimensionalConvergence, Random Forests Wager2015AdaptiveForests,Syrgkanis2020EstimationDimensions, Neural Nets Farrell2021DeepInference, forward model selection Kozbur2020AnalysisSelection or ensembles of those can be shown to achieve the required rates under conditions stated in the original papers.}
It is well known that estimators using doubly robust scores and parametric models for the nuisance parameters are doubly robust in the sense that they remain consistent if one of the parametric models is misspecified \cite<see, e.g.>{Glynn2009AnEstimator}. The difference of the DML version is that it exploits what \citeA{Smucler2019AContrasts} call 'rate double robustness'. This robustness allows to estimate the parameters of interest at the parametric rate $n^{1/2}$ even if the nuisance parameters are estimated at slower rates using machine learning methods that do not require the specification of an actual parametric model.
The rate double robustness is the consequence of the so-called Neyman orthogonality of the doubly robust score. Neyman orthogonality is at the heart of the general DML framework of Chernozhukov2018. Scores with this orthogonality are immune against small errors in the estimation of nuisance parameters and thus allow them to be estimated via machine learning. Appendix (ref) revisits what this means in formal terms.
We can reuse the ATE score of Equation (ref) to estimate conditional effects. The so-called DR-learner was introduced for binary treatments but directly translates also to multiple treatments settings. It exploits that the conditional expectation of the score with known nuisance parameters equals CATE: $\tau_{w,w'}(z) = E[\Delta_{i,w,w'} \mid Z_i = z]$.\footnote{Note that this does not work for the ATET score in Equation (ref) and suitable adaptations are beyond the scope of this paper.} Thus, a natural way to estimate CATEs is to use the score with estimated nuisance parameters, $\hat{\Delta}_{i,w,w'}$, as pseudo-outcome in a general regression framework:
From a conceptual and estimation perspective it is instructive to distinguish two special cases of CATEs at this point \cite<see also>{Knaus2021}: (i) Group average treatment effects (GATE) provide the average effects for pre-specified, usually low-dimensional, groups.\footnote{Note that the GATE is different to the Sorted Group Average Treatment Effect (GATES) of \citeA{Chernozhukov2017GenericExperiments}.} This covers standard subgroup analysis comparing, e.g., effects of men and women, or heterogeneity along pre-specified continuous variables like age. (ii) Individualised average treatment effects (IATEs) aim for the most detailed effect heterogeneity considering all confounders as heterogeneity variables, i.e. $Z_i = X_i$ and thus $IATE(x) = \tau_{w,w'}(x) = E[Y_i(w) - Y_i(w') \mid X_i=x]$.
OLS, series or kernel regressions of the pseudo-outcome on low-dimensional heterogeneity variables estimate GATEs. The outputs of such regressions can be interpreted in the standard way. The only difference is that instead of modelling the level of an outcome, they now model the level of a causal effect. Most importantly standard statistical inference applies as is shown for OLS and series regression by \citeA{Semenova2021DebiasedFunctions} as well as for kernel regression by \citeA{fan2020EstimationData} and \citeA{Zimmert2019NonparametricConfounding}. Similar to the discussion in the previous section, the Neyman orthogonality of $\Delta_{i,w,w'}$ allows to ignore that nuisance parameters are estimated with flexible methods potentially converging slower than $n^{1/2}$ when calculating standard errors. The details about the required convergence rates are discussed in the referenced papers.
IATEs may be estimated using the pseudo-outcome in supervised machine learning regressions with the full set of confounders as predictors. As discussed by \citeA{Chernozhukov2017GenericExperiments} statistical inference is not yet well understood for low-dimensional $Z_i$ and even harder for high-dimensional $Z_i$ when machine learning is used to solve Equation (ref). However, \citeA{Kennedy2020OptimalEffects} shows that the doubly robust structure of the ATE score results in favourable bounds on the mean squared error for the estimated IATEs that would not be attainable by outcome regression or IPW based methods alone.\footnote{The Orthogonal Random Forest of \citeA{Oprescu2019OrthogonalInference} is another estimator that is based on the pseudo-outcome idea and can be asymptotically normal under the assumption of parameteric nuisance parameters. We focus in this paper on the more general DR-learner. See also \citeA{Curth2021NonparametricAlgorithms} for a more nuanced analysis of the DR-learner in comparison to other alternatives.}
We consider two variants of the DR-learner for IATEs. First, we reuse the pseudo-outcome in one supervised machine learning regression to estimate IATEs in-sample. This full sample procedure is computationally convenient but prone to overfitting. Thus, the second variant produces out-of-sample IATE predictions for each individual in the sample. Following Algorithm 1 of \citeA{Kennedy2020OptimalEffects}, this requires a four-fold cross-fitting scheme that is detailed in Algorithm 2 of Appendix (ref). The computational downside of this procedure is that we cannot reuse the same nuisance parameter predictions as for the average estimator and need to estimate them for the IATE only. However, the results below suggest that this computational effort is important to avoid severe overfitting.
Note that the point estimates of the plain DR-learner can be expressed as $\hat{\tau}^{dr}(z) = \sum_{i=1}^N \alpha_i \hat{\Delta}_{i,w,w'}$ if the weight $\alpha_i$ that each observation receives can be calculated. For example, the ATE estimator as special case of the DR-learner with $Z_i$ being a constant uses $\alpha_i = 1/N$, the least squares regression uses $\alpha_i = z (\bm{Z'}\bm{Z})^{-1} \bm{Z'}$ with $\bm{Z}$ being the stacked covariate matrix, and the kernel regression uses $\alpha_i = \frac{\mathcal{K}_h \left(Z_i - z \right) }{ \sum_{i=1}^N \mathcal{K}_h \left( Z_i - z \right) }$ with $\mathcal{K}_h$ representing a proper kernel function. The class of estimators with a known weighted representation is called linear smoothers \cite<see e.g.>{Buja1989LinearModels}. Popular machine learners like tree-based methods (regression trees, Random Forests or boosted trees), Ridge or any method that runs OLS after variable selection like Post-Lasso Belloni2013LeastModels have this structure.\footnote{In practice most of these methods are applied with data-driven selection of tuning parameters, which makes them strictly speaking non-linear smoothers Buja1989LinearModels. However, this does not affect our results.} Also for these methods we know the weight $\alpha_i(x)$ that each observation receives in predicting the (pseudo-)outcome at $x$. These weights usually sum up to one, i.e. $\sum_{i=1}^N \alpha_i(x) = 1$. Using such outcome weighting predictors in the final step allows to express the DR-learner estimated IATE as
where $\Tilde{Y_i}(w,X_i) = Y_i - \hat{\mu}(w,X_i)$ denotes the outcome residual.
The DR-learner shares the problem of all estimators that involve reweighting by the inverse of the propensity score. In finite samples, ${\lambda^w_i}(x)$ and ${\lambda^{w'}_i}(x)$ usually do not sum to one, i.e. $\sum_{i=1}^N {\lambda^w_i}(x) \neq 1$ and $\sum_{i=1}^N {\lambda^{w'}_i}(x) \neq 1$. This is especially problematic if it sums to something much greater than one. In this case the weighted residuals receive much more weight than the outcome regressions. This might result in implausibly large effect estimates that even could fall outside of the possible bounds of a given outcome variable Kang2007,Robins2007Comment:Variable.\footnote{For bounded outcomes, the effects must lie in the interval $[Y_{min} - Y_{max},Y_{max} - Y_{min}]$, with $Y_{min}$ and $Y_{max}$ denoting the minimum and maximum values of the outcome, respectively.}
For average effects the \citeA{Hajek1971CommentOne} normalisation is recommended to stabilise estimators using inverse probability weights \cite<e.g.>{Imbens2004NonparametricReview,Lunceford2004StratificationStudy,Robins2007Comment:Variable,busso2014new}. However, Equation (ref) highlights that the inverse probability weights become $x$-specific and such a one time normalisation that targets the average effect does not solve the problem for the individualised effect. This can be problematic as finite sample imbalances are more likely to occur on the individualised level. Thus, we propose the normalised DR-learner (NDR-learner) as a stabilised complement to the DR-learner.
The NDR-learner normalises the weighted residuals by the sum of weights:
This ensures that the weights of the residuals sum up to one under the condition that weights $\alpha_i(x)$ are non-negative. Thus, methods like Ridge or Post-Lasso with potentially negative weights might not be applicable.
The NDR-learner is more demanding from a computational point of view because it requires to calculate the weights $\alpha_i(x)$ and the normalisation for each $x$ of interest (Algorithm 2 in Appendix (ref) provides the details of the implementation). However, the application below shows that the normalisation deals well with the cases where outcome residuals receive high weights leading to implausibly large effect estimates. Thus, the NDR-learner is an interesting alternative to the DR-learner if effect sizes become suspicious.
The APO score of Section (ref) can also be reused to estimate optimal treatment assignment. To this end, note that the value function of any policy rule $\pi(Z_i)$ can be estimated as \[ \hat{Q}(\pi) = N^{-1} \sum_{i=1}^N \sum_{w=0}^{T} \mathds{1}(\pi(Z_i) = w) \hat{\Gamma}_{i,w}. \]
This means each individual contributes the score of the treatment that she is assigned to under this policy rule. However, we are not necessarily interested in the value function of some policy rule, but want to estimate the optimal policy rule that maximises this value function, $\hat{\pi}^* = \operatorname*{\arg\max}_{\pi \in \Pi} \hat{Q}(\pi)$. This requires to search over all candidate policy rules to find the optimum as there exists no closed form solution.
Example: Consider the case where $Z_i$ is a binary covariate and $W_i$ is a binary treatment. We have four different policy rules: treat nobody ($\pi^1$), treat only those with $Z_i=1$ ($\pi^2$), treat only those with $Z_i=0$ ($\pi^3$), or treat everybody ($\pi^4$). We illustrate this using two representative observations, $i=1$ with $Z_1=0$, and $i=2$ with $Z_2=1$ in Table (ref). The columns three to six show the assignments under the four potential assignment rules. For example, the first observation receives no treatment under policy rules $\pi^1$ and $\pi^2$, but is treated under policy rules $\pi^3$ and $\pi^4$. To find the optimal rule, we compare the means of the APO scores in the last four columns and pick the policy rule that corresponds to the largest mean. The number of policy values to compare increases dramatically in settings with multiple treatments and $Z_i$ being a vector of potentially non-binary variables.
We expect that the estimated policy in finite samples and with estimated nuisance parameters does not coincide with the true optimal policy rule. This is conceptualised as the 'regret' defined as the difference between the true and the estimated optimal value function, $R(\hat{\pi}^*) = Q(\pi^*) - Q(\hat{\pi}^*)$.
\citeA{Zhou2018OfflineOptimization} show that the DML based procedure minimises the maximum regret asymptotically under two main conditions: First, the same convergence conditions for the nuisance parameters that are required for ATE estimation (the product of the nuisance parameter convergence rates achieves $n^{1/2}$). Second, the set of candidate policy rules $\Pi$ is not too complex. In particular, \citeA{Zhou2018OfflineOptimization} show that decision trees with fixed depth are a suitable class of policy rules. Again the double robustness of the used scores results in statistical guarantees that are not achievable for methods based on outcome regressions or IPW alone.
We use a standard observational dataset of Swiss Active Labour Market Policy (ALMP) that is already basis of previous studies Huber2017,Lechner2018,Knaus2020HeterogeneousApproach to estimate the effect of different programs on employment.\footnote{\citeA{gerfin2002microeconometric}, \citeA{Lalive2008TheUnemployment} and \citeA{Knaus2020HeterogeneousApproach} among others provide a more detailed description of the surrounding institutional setting.} In particular, we start with the sample of 100,120 unemployed individuals of \citeA{Huber2017} that consists of 24 to 55 year old individuals registered unemployed in 2003.\footnote{The dataset is available as restricted use file via the platform \href{https://forsbase.unil.ch/project/study-public-overview/17035/1/}{FORSbase} Lechner2020SwissDataset.} We consider non-participants and participants of four different program types: job search, vocational training, computer programs and language courses.\footnote{The dataset contains also participants of an employment program and personality training. However, we leave them out to keep the number of obtained results manageable.} As the assignment policies differ substantially across the three language regions, we focus only on individuals living in the German speaking part and remove those in the French and Italian speaking part to avoid common support problems.
We evaluate the first program participation within the first six months after the begin of the unemployment spell. One problem of this definition is that non-participants comprise people that quickly come back into employment before they would be assigned to a training program. This could result in an overly optimistic evaluation of non-participation. We follow \citeA{lechner1999earnings} and \citeA{lechner2007value} and assign pseudo program starting points to the non-participants and keep only those who are still unemployed at this point.\footnote{The assignment of the pseudo starting point is based on estimated probabilities to start a program at a specific time. The probability depends also on covariates and is estimated using the same random forest specification that is discussed later in Section (ref).} This results in a final sample size of 62,497 observations.
The outcome of interest is the cumulated number of months in employment in the 31 months after program start, which is the maximum available time span in the dataset. Row one of Table (ref) provides the number of observations in each group. Roughly 75% participate in no program. By far the largest program is the job search program, which is also called basic program. The more specific programs are much smaller with roughly 1000 observations each. Row two shows that the average outcomes substantially differ by different groups. However, it is not clear whether this is only due to selection effects because the observable characteristics are not comparable across groups, as the remaining rows show. Especially the share of females, the share of foreigners and past income differ quite substantially across programs. The confounders comprise 45 variables and are reported in Table (ref) of Appendix (ref). They consist of socio-economic characteristics of the unemployed individuals, caseworker characteristics, information about the assignment process, information about the previous job, and regional economic indicators.
The nuisance parameters are estimated via Random Forest Breiman2001 using the implementation with honest splitting in the grf R-package Athey2017a and 5-fold cross-fitting. The tuning parameters in each regression are selected by out-of-bag validation. All regressions apply the full set of confounders. We run the outcome regressions for each treatment group separately to obtain $\hat{\mu}(w,x)$. Also the propensity scores are separately estimated for each treatment using a treatment indicator as outcome in the random forest. The propensity scores are then normalised to sum to one within an individual.
We estimate CATEs at different granularity. First, we investigate GATEs for subgroups by gender, foreigners and three categories of employability. These are regularly used in the program evaluation literature and usually investigated by re-estimating everything in the subgroups. However, it can be performed at very low computational costs after DML for average effects using only a standard OLS regression with the pseudo-outcome as described in Section (ref) and using dummy variables for all groups but the reference group as covariates. Second, we estimate kernel regression and spline regression GATEs for the continuous variables age and past income based on the R-packages np Hayfield2008NonparametricPackage and crs Racine2021Crs:Splines, respectively. The kernel regressions apply a second-order Gaussian kernel function and use 0.9 of the cross-validated bandwidth for undersmoothing as suggested by \citeA{Zimmert2019NonparametricConfounding}. The spline regressions use B-splines with cross-validated degree and number of knots. Third, we specify an OLS regression in which all the five previously used variables enter linearly. Finally, we go beyond the handpicked variables and estimate the IATEs using all 45 confounders in the DR-learner and the NDR-learner. Both are implemented with the honest Random Forest because the grf package allows to extract the prediction weights $\alpha_i(x)$ required for the NDR-learner. We apply both variants described in Section (ref). Once we estimate the IATE for each observation using DR- and NDR-learner in the full sample and once we predict them out-of-sample. For the latter, Appendix (ref) provides a detailed description of the underlying DR- and NDR-learner algorithms.
The optimal treatment assignment rule is estimated as decision trees of depth one, two and three. We follow Algorithm 2 for exact tree-search of \citeA{Zhou2018OfflineOptimization} that is implemented in the policytree R-package Sverdrup2020Policytree:Trees. We estimate the trees first with the five handpicked variables. However, these variables include gender and foreigner status that might be too sensitive to include in practice. Thus, we investigate another set of 16 variables that includes only the objective measures of education and labour market history of the unemployed persons that would be available for recommendations from the administrative records.
Table (ref) summarises all required implementation steps. It highlights that a comprehensive DML based program evaluation can be run with few lines of code in any statistical software program that is capable of the operations in the third column. Thus, researchers can build their customised analyses in a modular fashion based on established code. Alternatively, the R-package \href{https://github.com/MCKnaus/causalDML}{causalDML} already implements the required steps as showcased in the \href{https://mcknaus.github.io/assets/code/Notebook_DML_ALMP_MCK2020.html}{replication notebook} accompanying this paper. Most importantly the package provides a fast implementation of the individualised normalisation required for the NDR-learner in C++.
We focus here on the effect estimates and discuss the nuisance parameters in Appendix (ref). Throughout this section, we compare the four programs to non-participation.\footnote{The underlying APOs are shown in Figure (ref) of Appendix (ref).} Recall that the outcome of interest is the cumulated number of months employed in the 31 months after program start. Figure (ref) depicts ATE and ATET estimates and shows substantial differences in the effectiveness of programs. The job search program decreases the months in employment on average by about one month. In contrast, other programs that teach hard skills show substantial improvements with roughly three additional months in employment on average.\footnote{For a better understanding of the underlying dynamics, Figure (ref) of Appendix (ref) reports and discusses the effects of program participation on the employment probabilities over time.}
Comparing ATE and ATET shows no big differences for most programs. This suggests that there is either no effect heterogeneity correlated with observables or that the assignment does not take advantage of this heterogeneity. We would expect to see ATETs being higher than ATEs if program assignment is well targeted. However, we find only evidence for the opposite as the actual participants of a language course show a 1.5 months lower treatment effect compared to the population. This difference suggests that there is substantial effect heterogeneity to uncover and the potential to improve treatment assignment.
This subsection studies effect heterogeneity at different granularity. We start by estimating group average treatment effects (GATEs) for discrete subgroups. Panel A of Table (ref) shows the result of an OLS regression with a female dummy as covariate, $\hat{\Delta}_{i,w,w'} = \beta_0 + \beta_1 female_i + error_i$. The constant ($\beta_0$) provides the GATE for the reference group men and the female coefficient ($\beta_1$) describes how much the GATE differs for women. The results show substantial gender differences in the effectiveness of programs. Women significantly suffer less or profit more from job search and computer program participation. This gender gap in the effectiveness of ALMPs is also well-documented in the literature crepon2016active,Card2018WhatEvaluations. In contrast to this, we find that women profit on average significantly less from language courses than men.
Panel B replaces the female dummy in the regression by a foreigner dummy. Strikingly, Swiss citizens as reference group show a big positive effect for participating in language courses but the effect disappears for foreigners. After adding the coefficient for foreigners to the constant, the foreigners' GATE is only 0.72 ($3.56 - 2.84$, standard error: 0.69). A crucial information to better understand this finding would be to know which languages they learn.\footnote{See \citeA{Heiler2021EffectTreatments} for a discussion about how treatment heterogeneity could drive effect heterogeneity.} However, this information is unfortunately not available in the dataset.
Panel C shows the results of a similar regression but now with two dummies indicating medium and high employability such that low employability becomes the reference group. The F-statistic in the last line tests the joint significance of the two dummies. It is statistically significant at least at the 10%-level for the programs in the first three columns. They all show a common gradient that individuals with low employability benefit substantially more or at least suffer less from program participation.
While subgroup analyses are standard in program evaluations, the estimation of nonparametric GATEs using kernel or series regression is rarely pursued. We estimate such GATEs along two continuous variables past income and age. We find no notable heterogeneity for the latter.\footnote{Figure (ref) shows the according results.} However, effect sizes are clearly associated with past income. Figure (ref) shows on the left the results of the kernel regression and on the right of the series regression. Both estimators agree by and large. They document that effects decrease with higher past income for all but for language programs. The latter have only a small positive effect for individuals with low past income but it increases with higher income. One potential explanation for these findings is that the value of language skills is larger for high-skilled workers in multilingual countries like Switzerland because they reduce information costs across language borders \cite<see, e.g.>{Isphording2014LanguageSuccess}.
The GATEs considered so far were nonparametric but only univariate. Now we model the GATE by specifying a multivariate OLS regression with the previously used covariates entering linearly. It is most likely misspecified and thus estimates the best linear predictor (BLP) of GATEs with respect to these variables. However, it provides a compact and accessible summary of the effect heterogeneities. Additionally, it holds the other included variables constant. Consider for example the coefficients for being female in Table (ref). Compared to Table (ref), the coefficients in the first three columns are smaller and the one for language courses is larger (for example for job search it is 0.25 instead of 0.60). The reason is that it represents a partial effect that holds other variables like past income fixed. The subgroup female coefficient in Table (ref) partly picks up that women have lower past income and that lower income is associated with higher treatment effects for all but language courses. This example illustrates that the same strategies that are usually applied to interpret an outcome OLS regression can now be used to interpret the effect OLS regression.
We focus on the results based on the out-of-sample variant of the DR- and NDR-learner as the full sample variant leads to severe overfitting with predicted IATEs ranging from -265 to 161 that are up to eight times larger than what is possible given that the outcome is bounded between zero and 31.\footnote{See Appendix (ref) for results and discussion of the full sample.} However, Figure (ref) shows that the DR-learner produces impossible effect sizes even out-of-sample, which motivates the proposal of the NDR-learner as stabilised variant. Figure (ref) provides boxplots of the predicted IATEs and shows outliers lying below the smallest possible value of -31. However, the descriptive statistics provided in Table (ref) and the joint and marginal distributions depicted in Figure (ref) document that besides the outliers the distributions are quite similar and correlate with at least 0.90. Not surprisingly, the impact of normalised weights is much larger for the three smaller programs and nearly negligible for job search programs. Still, we base the following discussion for all programs on the more stable results of the NDR-learner.
We conduct a classification analysis as proposed by \citeA{Chernozhukov2018TheAverages} to understand which variables are most predictive of effect sizes. To this end, we split the predicted IATE distributions in quintiles and compare the covariate means of the observations falling into the fifth and first quintile. For comparability, we normalise all covariates to have mean zero and variance one. Table (ref) shows the seven variables that have at least one absolute difference between the highest and lowest quintile that is larger than one standard deviation. For example, we observe that the group with the highest effects (the fifth quintile) of a job search program has a 1.32 standard deviations lower past income compared to the lowest IATE group (the first quintile). Also the other variables confirm the patterns that we document already in previous subsections. The effects of job search, vocational and computer training are higher for unskilled workers with lower previous labour market success and foreigners, while the opposite holds for language programs.\footnote{Table (ref) shows the classification analysis for all variables.}
\singlespacing
\doublespacing
The previous section documented substantial heterogeneities in the program effects. To leverage this heterogeneity for better targeting, we apply the DML based optimal policy algorithm of Section (ref). Figure (ref) shows the simplest decision tree with only one split for the five handpicked covariates. It would allocate men to vocational training and women to computer courses. This split is probably similar to what we would have suggested given the evidence presented in Table (ref). For a tree of depth two, such an eyeballing approach has its limits and the algorithmic approach provides a systematic way to arrive at an estimated optimal decision tree. The depth two tree in Figure (ref) splits first on past income and then recommends to send low earning men to vocational and low earning women to computer training, while high earners (more than CHF 55k) are recommended to participate in language training. In the absence of the possibility to split on gender, the depth one tree in Figure (ref) splits on past income roughly at the same value where the nonparametric GATEs of computer and language training intersect in Figure (ref).\footnote{Appendix (ref) provides the trees of depth three.}
Panel A of Table (ref) summarises the results of the different trees. It shows the percentage of individuals that are placed in the different programs. Not surprisingly, all individuals are recommended to be placed into one of the three positively evaluated hard skill enhancing programs.
One yet unsolved challenge is how to conduct statistical inference about the quality and stability of the decision trees. \citeA{Athey2021PolicyData} propose a form of cross-validation. To this end, we take the same folds that were used in the cross-fitting procedure to estimate the nuisance parameters. We build the decision tree in four folds and evaluate the value in the left out fold. First, we inspect how often the recommendations based on these trees coincide with the full sample policy rules. Figures (ref) and (ref) of Appendix (ref) show that the cross-validated trees are not identical to the full sample ones.
\singlespacing
\doublespacing
\citeA{Zhou2018OfflineOptimization} propose another validation idea and test whether the optimal policy rules perform significantly better than sending all individuals to the same program. This is achieved by taking the difference of the APO score of the cross-validated policy rule and the APO score of the program $w$: $ \hat{\Delta}^{cv}_{i,w}(\pi) = \sum_{t=0}^{T} \mathds{1}(\hat{\pi}^{cv}(Z_i) = t) \hat{\Gamma}_{i,t} - \hat{\Gamma}_{i,w}$, where $\hat{\pi}^{cv}(Z_i)$ is the policy rule that is estimated without individual $i$. A standard t-test on the mean of $\hat{\Delta}^{cv}_{i,w}(\pi)$ tests then whether the cross-validated policy rules are significantly better than sending everybody to the same program. Note that the cross-validated policy rules do not necessarily coincide with the trees in the full sample and the cross-validation estimates not the value function for that specific tree. This requires to hold out a test set, which would be viable for an application with bigger programs.
The results are provided in Panel B of Table (ref). We can interpret the mean of $\hat{\Delta}^{cv}_{i,w}(\pi)$ as average treatment effect comparing a regime under the estimated assignment rule or a regime where everybody is sent to the program $w$. This effect is always positive indicating that the estimated rules can leverage the effect heterogeneities to improve the allocation. However, the cross-validated policy rules perform not significantly better than sending just everybody into vocational or computer programs. This would probably change if we could take costs or capacity constraints into account. However, we do not observe costs in this dataset and the optimal decision tree algorithm is currently not capable of incorporating capacity constraints in a systematic way. We leave both extensions for future research using a more detailed database on both costs and capacity constraints.
This paper considers recent methodological developments based on Double Machine Learning (DML) through the lens of a standard program evaluation under unconfoundedness. DML based methods provide a convenient toolbox for a comprehensive program evaluation as different parameters of interest can be estimated using the same framework and a combination of standard statistical software. The application to an Active Labour Market Policy evaluation shows that the methods also produce plausible results in practice. The only exception is the DR-learner that required a modification, the newly introduced NDR-learner, before producing stable results for all individualised treatment effects. However, several conceptual and implementational issues remain open for investigation and refinement.
In general, we know little about how to choose the estimator for the nuisance parameters. The pool of potential machine learning algorithms and their combinations is large and little is known, e.g., about the trade-off between high prediction performance and computation time in the causal setting. Also clear recommendations for the implementation of cross-fitting are missing. Another open question is how to deal with common support in general and for each estimand specifically. The literature on trimming rules is well developed for propensity score based methods estimating average effects. However, we are not only interested in average effects and the propensity score is not the only nuisance parameter of DML. It remains an open question whether the established trimming methods are also sensible for DML when common support becomes an issue.
The estimators for flexible heterogeneous treatment effects provide interesting new tools. However, it is currently not clear to what extent we can actually explore heterogeneity or to what extent we need to pre-define the heterogeneity of interest. The possibility to summarise pre-defined heterogeneity of interest using OLS, kernel or series regressions provide clearly valuable and easy to use options in applications. The instability of methods that aim for individualised heterogeneous effects shows that they should be used with caution and more research is required to investigate whether adjustments like the proposed NDR-learner are useful beyond the application of this paper.
The estimation of optimal treatment assignment rules is mostly unexplored in practice and many interesting issues in applications regarding inference, the implementation of different constraints, more flexible rules than decision trees, or the choice of variables that could or should enter the set of policy variables, which could be explored in future research.
The investigation of these DML specific questions but also the comparison with other more specialised causal machine learning methods for each estimand provides another interesting direction of future research. Such evidence would help to understand and guide which choices are critical in applications similar to the one in this paper.