Double machine learning for sample selection models

\doublespacing \pagestyle{plain}

Abstract: { This paper considers the evaluation of discretely distributed treatments when outcomes are only observed for a subpopulation due to sample selection or outcome attrition. For identification, we combine a selection-on-observables assumption for treatment assignment with either selection-on-observables or instrumental variable assumptions concerning the outcome attrition/sample selection process. We also consider dynamic confounding, meaning that covariates that jointly affect sample selection and the outcome may (at least partly) be influenced by the treatment. To control in a data-driven way for a potentially high dimensional set of pre- and/or post-treatment covariates, we adapt the double machine learning framework for treatment evaluation to sample selection problems. We make use of (a) Neyman-orthogonal, doubly robust, and efficient score functions, which imply the robustness of treatment effect estimation to moderate regularization biases in the machine learning-based estimation of the outcome, treatment, or sample selection models and (b) sample splitting (or cross-fitting) to prevent overfitting bias. We demonstrate that the proposed estimators are asymptotically normal and root-n consistent under specific regularity conditions concerning the machine learners and investigate their finite sample properties in a simulation study. We also apply our proposed methodology to the Job Corps data for evaluating the effect of training on hourly wages which are only observed conditional on employment. The estimator is available in the causalweight package for the statistical software R. }

{ Keywords: sample selection, double machine learning, doubly robust estimation, efficient score.}

{ JEL classification: C21. \quad }

{ {\tiny We have benefited from comments by Alyssa Carlson, David Kaplan, Peter Mueser, and seminar participants at the University of Missouri. Addresses for correspondence: Michela Bia, Luxembourg Institute of Socio-Economic Research, 11 Porte des Sciences, Maison des Sciences Humaines, 4366 Esch-sur-Alzette/Belval, Luxembourg, [email removed], [email removed]; Martin Huber, University of Fribourg, Bd.\ de P\'{e}rolles 90, 1700 Fribourg, Switzerland, [email removed]; Luk\'{a}\v{s} Laff\'{e}rs, Matej Bel University, Tajovskeho 40, 97411 Bansk\'{a} Bystrica, Slovakia, [email removed]. Laff\'{e}rs acknowledges support provided by the Slovak Research and Development Agency under contract no. APVV-17-0329 and VEGA-1/0692/20. }\thispagestyle{empty} }

{ \setcounter{footnote}{0} \setcounter{footnote}{0} \setcounter{page}{1} }

Introduction

In many studies aiming at evaluating the causal effect of a treatment or policy intervention, the empirical analysis is complicated by non-random outcome attrition or sample selection. Examples include the estimation of the returns to education when wages are only observed for the selective subpopulation of working individuals or the effect of educational interventions like vouchers for private schools on college admissions tests when students non-randomly abstain from the test. Furthermore, in observational studies, treatment assignment is typically not random, implying that the researcher faces a double selection problem, namely selection into the treatment and observability of the outcome. A large literature addresses treatment selection by assuming a selection-on-observables assumption, implying that treatment is as good as randomly assigned conditional on observed pre-treatment covariates, see for instance the reviews by Im04 and ImWo08. Furthermore, a growing number of studies addresses the question of how to control for the crucial confounders in a potentially high-dimensional vector of covariates in a data-driven way based on machine learning algorithms, see for instance the double machine learning framework of Chetal2018.

In this paper, we adapt the double machine learning framework to the evaluation of binary or multiply discrete treatments in the presence of sample selection or outcome attrition. In terms of identifying assumptions, we combine a selection-on-observables assumption for the treatment assignment with either selection-on-observables or instrumental variable assumptions concerning the outcome attrition/sample selection process. Such assumptions have previously been considered in Hu11 and Hu11b for the estimation of the average treatment effect (ATE) based on inverse probability weighting, however, for pre-selected (or fixed) covariates. As methodological advancement, we derive doubly robust and efficient score functions for evaluating treatment effects under double selection and demonstrate that they satisfy so-called Neyman1959 orthogonality. The latter property permits controlling for covariates in a data-driven way by machine learning-based estimation of the treatment, outcome, and attrition models under specific conditions. Therefore, the subset of important confounders need not be known a priori (but must be contained in the total set of covariates), which is particularly useful in high dimensional data with a vast number of covariates that could potentially serve as control variables. We also consider dynamic confounding based on a sequential selection-on-observables assumption that is closely related to assumptions found in the dynamic treatment effect literature as e.g.\ in Ro86, Robins1998, and Lech09. This assumption permits that covariates that jointly affect sample selection and the outcome may themselves be a function of the treatment, a scenario widely neglected in sample selection models despite its likely relevance in empirical applications. In particular when there is a substantial time lag between treatment assignment and the sample selection process, exploiting post-treatment covariates to tackle selection-outcome confounding seems more convincing than solely relying on pre-treatment covariates (as in conventional selection-on-observables assumptions) for addressing both treatment endogeneity and sample selection.

Following Chetal2018, we show that treatment effect estimation based on our score functions (that are tailored to the various identifying assumptions) is root-$n$ consistent and asymptotically normal under particular regularity conditions, in particular the $n^{-1/4}$-convergence of the machine learners. A further condition in the double machine learning framework is the prevention of overfitting bias due to correlations between the various estimation steps. This is obtained by estimating the treatment, outcome, and selection models on the one hand and the treatment effect on the other hand in different parts of the data. As in Chetal2018, we subsequently swap the roles of the data parts and average over treatment effects in order to prevent asymptotic efficiency losses, a procedure known as cross-fitting. We also provide a simulation study suggesting that our estimators perform decently in terms of the root mean squared error and coverage (by confidence intervals) in the simulation designs with several thousand observations considered. Finally, we present an empirical illustration considering the female sample of a study on Job Corps, a large training program for disadvantaged youth in the U.S. We apply our DML estimators to assess the effects of academic and vocational training on hourly wage, which is only observed conditional on employment, one and four years after program assignment and find some statistical evidence for positive longer-run impacts.

Our paper is related to a range of studies tackling sample selection and selective outcome attrition. One strand of the literature models the attrition process based on a selection-on-observables assumption also known as missing at random (MAR) condition. The latter imposes conditional independence of sample selection and the outcome given observed information like the covariates and the treatment. Examples include Ru76b, LittleRubin87, Carroll+95, Shah+97, FitzgeraldGottschalkMoffitt98a, Abowd+01, Wooldridge02, and Wo07. Robins+94, RobinsRotnitzkyZhao95, and BaRo05 discuss doubly robust estimators of the outcome that are consistent under MAR when either the conditional outcome or the attrition model are correctly specified. This approach satisfies Neyman orthogonality as required for double machine learning.\footnote{Relatedly, BarnwellChaudhuri2020 consider several outcome periods under a monotonic MAR assumption (i.e.\ outcome attrition being an absorbing state weakly increasing over time) and also discuss the evaluation of randomly assigned treatments in this context based on the efficient influence function. In contrast, our framework considers a single outcome period and permits selection into treatment to be related to observed confounders.} However, their framework does not consider double selection into treatment and the observability of the outcome at the same time as we do in this paper.

Negi2020 suggests an alternative estimator under double selection that falls into the weighted M-estimation framework described in slonwooldridge2018 and also satisfies doubly robustness, i.e.\ remains consistent under parametric misspecification of either the conditional outcome model or the treatment and selection models. This approach based on reweighting outcome models is nevertheless different to ours making use of efficient influence functions and to the best of our knowledge, Neyman1959 orthogonality (as required for double machine learning) has not been shown for weighted M-estimation (while we prove this property for our proposed estimators). A further difference is that Negi2020 focuses on treatment evaluation when controlling for pre-treatment covariates to tackle double selection, while we in addition consider identification based on both pre- and post-treatment covariates (dynamic confounding) or an instrument for sample selection.

In contrast to MAR-based identification, so-called sample selection or nonignorable nonresponse models allow for unobserved confounders of the attrition process and the outcome. Unless strong functional form assumptions as in Heck76, Heckman79, Hausman79, and Little1995 are imposed, identification requires an instrumental variable (IV) for sample selection. We refer to Daneve03, Ne07, Hu11, and Hu11b for nonparametric estimation approaches in this context. To the best of our knowledge, this study is the first one to propose a doubly robust treatment effect estimator under nonignorable outcome attrition and to consider machine learning techniques to control for (possibly high-dimensional) covariates in this context. Our estimators are available in the causalweight package for R by BodoryHuber2018.

This paper proceeds as follows. Using the potential outcome framework, Section (ref) discusses the identification of the average treatment effect when outcomes are assumed to be missing at random (i.e.\ selection is on observables, as for the treatment) conditional on pre-treatment covariates. Section (ref) considers identification when outcome attrition is related to unobservables, known as nonignorable nonresponse, and an instrument is available for tackling this issue. Section (ref) demonstrates identification under a sequential selection-on-observables which allows for dynamic confounding, meaning that outcomes are assumed to be missing at random conditional on pre- and post-treatment covariates. Section (ref) proposes an estimator based on double machine learning and shows root-n consistency and asymptotic normality under specific regularity conditions. Section (ref) provides a simulation study. Section (ref) presents an empirical application to data from the US Job Corps Study. Section (ref) concludes.

Identification under missingness at random

Our target parameter is the average treatment effect (ATE) of a binary or multiply discretely distributed treatment variable $D$ on an outcome variable $Y$. To define the effect of interest, we use the potential outcome framework, see Rubin74. Let $Y(d)$ denote the potential outcome under hypothetical treatment assignment $d$ $\in$ $\{0,1,...,Q\}$, with $0$ indicating non-treatment and $1,...,Q$ the different treatment choices (where $Q$ denotes the number of non-zero treatments). The ATE when comparing two distinct treatment $d\neq d'$ then corresponds to $\Delta=E[Y(d)-Y(d')]$. Furthermore, let $Y$ denote the outcome realized under the treatment (f)actually assigned to a subject, i.e.\ $Y=Y(D)$. Therefore, $Y$ corresponds to the potential outcome under the treatment received, while the potential outcome under any counterfactual treatment assignment remains unknown. A further complication in our evaluation framework is that $Y$ is assumed to be only observed for a subpopulation, i.e. conditional on $S=1$, where $S$ is a binary variable indicating whether $Y$ is observed/selected, or not.

Empirical examples with partially observed outcomes include wage regressions, with $S$ being an employment indicator, see for instance Gr74, or the evaluation of the effects of policy interventions in education on test scores, with $S$ being participation in the test, see AnBeKr04. Throughout our discussion, $S$ is permitted to be a function of $D$ and $X$, i.e.\ $S=S(D,X)$. However, $S$ must neither be affected by nor affect $Y$.\footnote{See for instance Imai09 for alternative assumptions, which imply that selection is associated with the outcome but is independent of the treatment conditional on the outcome and other observable variables.} Therefore, selection per se does not causally influence the outcome. The following nonparametric outcome and selection models satisfy this framework:

where $U,V$ are unobserved characteristics and $\phi, \psi$ are general functions.\footnote{Note that $Y(d)=\phi(d,X,U)$, which means that fixing the treatment yields the potential outcome.} Throughout the paper we assume that the stable unit treatment value assumption (SUTVA, Rubin80) holds such that $\Pr(D = d \implies Y = Y(d))=1$ This rules out interaction or general equilibrium effects and implies that the treatment is uniquely defined.

We subsequently formalize the assumptions that permit identifying the average treatment effect when both selection into the treatment and outcome attrition is related to observed characteristics. \newline Assumption 1 (conditional independence of the treatment):\newline $Y(d) \bot D | X=x$ for all $d \in \{0,1,...,Q\}$ and $x$ in the support of $X$.\newline By Assumption 1, there are no unobservables jointly affecting the treatment and the outcome conditional on covariates $X$. For model (ref), this implies that $U$ is not associated with unobserved terms affecting $D$ given $X$. In observational studies, the plausibility of this assumption crucially hinges on the richness of the data, while in experiments, it is satisfied if the treatment is randomized within strata defined by $X$ or randomized independently of $X$. \newline Assumption 2 (conditional independence of selection):\newline $ Y \bot S | D=d, X=x $ for all $d \in \{0,1,...,Q\}$ and $x$ in the support of $X$. \newline By Assumption 2, there are no unobservables jointly affecting selection and the outcome conditional on $D,X$, such that outcomes are missing at random (MAR) in the denomination of Ru76b. Put differently, selection is assumed to be selective w.r.t.\ observed characteristics only. For model (ref), this implies that $U$ and $V$ are conditionally independent given $D,X$. \newline Assumption 3 (common support):\newline (a) $\Pr(D=d| X=x)>0$ and (b) $\Pr(S=1| D=d, X=x)>0$ for all $d \in \{0,1,...,Q\}$ and $x$ in the support of $X$.\newline Assumption 3(a) is a common support restriction requiring that the conditional probability to receive a specific treatment given $X$, henceforth referred to as treatment propensity score, is larger than zero in either treatment state. Assumption 3(b) requires that for any combination of $D,X$, the conditional probability to be observed, henceforth referred to as selection propensity score, is larger than zero. Otherwise, the outcome is not observed for some specific combinations of these variables implying yet another common support issue.

Figure (ref) provides a graphical illustration of our identification setup using a directed acyclic graph, with arrows representing causal effects. Each of $D$, $S$, and $Y$ might be causally affected by distinct and statistically independent sets of unobservables not displayed in Figure (ref), but none of these unobservables may jointly affect $D$ and $Y$ given $X$ or $S$ and $Y$ given $D$ and $X$.\newline

Our identifying assumptions imply that

where the first equality follows from Assumption 1 and the second equality from Assumption 2. Therefore, the mean potential outcome identified by

or, using the fact that $E[Y| D=d,S=1, X]=E[I\{D=d\} \cdot S \cdot Y | X]/\Pr(D=d, S=1|X)$, by

where the second equality follows from the law of iterated expectations. $I\{\cdot\}$ denotes the indicator function, which is equal to one if its argument is satisfied and zero otherwise. Division by $\Pr(D=d|X) \cdot \Pr(S=1| D=d, X)$ in (ref) also demonstrates the importance of Assumption 3 for nonparametric identification. For the sake of brevity, we henceforth denote by $\mu(D,S,X)=E[Y| D ,S , X]$ the conditional mean outcome and by $p^d(X)=\Pr(D=d|X)$ and $\pi(D,X)=\Pr(S=1| D, X)$ the propensity scores. Expressions (ref) and (ref) suggest that the mean potential outcomes (and thus, the ATE) are identified, either based on conditional mean outcomes or inverse probability weighting using the treatment and selection propensity scores.

Following the literature on doubly robust methods, see e.g.\ RobinsMarkNewey1992, Robins+94, and RobinsRotnitzkyZhao95, we combine both approaches to obtain the following identification result:

The result in (ref) is based on the so-called efficient score function, which is formally derived in Appendix (ref) following the approach outlined in Levy2019. By noting that

it is easy to see that (ref) is equivalent to (ref) and thus, (ref). In contrast to (ref) and (ref), however, expression (ref) is doubly robust in the sense that it identifies $E[Y(d)]$ if either the conditional mean outcome $\mu(d,1,X)$ or the propensity scores $p_d(X)$ and $\pi (d,X)$ are correctly specified. Furthermore, it satisfies the so-called Neyman1959 orthogonality, i.e.\ is first-order insensitive to perturbations in $\mu(D,S,X)$, $p_d(X)$, and $\pi (D,X)$, see Appendix (ref). This entails desirable robustness properties when using machine learning to estimate the outcome, treatment, and selection models in a data-driven way.

Identification under nonignorable nonresponse

When sample selection or outcome attrition is related to unobservables even conditional on observables, identification generally requires an instrument for $S$. We therefore replace Assumptions 2 and 3, but maintain Assumption 1 (i.e.\ selection into treatment is on observables). \newline Assumption 4 (Instrument for selection):\newline (a) There exists an instrument $Z$ that may be a function of $D$, i.e. $Z=Z(D)$, is conditionally correlated with $S$, i.e.\ $E[Z\cdot S | D,X]\neq 0$, and satisfies (i) $Y(d,z)=Y(d)$ and (ii) $Y \bot Z | D=d, X=x$ for all $d \in \{0,1,...,Q\}$ and $x$ in the support of $X$,\\ (b) $S=I\{V \leq \chi(D,X,Z)\}$, where $\chi$ is a general function and $V$ is a scalar (index of) unobservable(s) with a strictly monotonic cumulative distribution function conditional on $X$,\\ (c) $V \bot (D,Z)|X$. \newline Assumption 4 no longer imposes the conditional independence of $Y$ and $S$ given $D,X$. As the unobservable $V$ in the selection equation is allowed to be associated with unobservables affecting the outcome, Assumptions 1 and 2 generally do not hold conditional on $S=1$ due to the endogeneity of the post-treatment variable $S$. In fact, $S=1$ implies that $\chi(D,X,Z)> V$ such that conditional on $X$, the distribution of $V$ generally differs across values of $D$. This entails a violation of the conditional independence of $D$ and $Y(d)$ given $S=1$ and $X$ if the potential outcome distributions differ across values of $V$. We therefore require an instrumental variable denoted by $Z$, which must not affect $Y$ or be associated with unobservables affecting $Y$ conditional on $D$ and $X$, as invoked in 4(a).\footnote{As an alternative set of IV restrictions in the context of selection, Haultfoeuille2010 permits the instrument to be associated with the outcome, but assumes conditional independence of the instrument and selection given the outcome.} We apply a control function approach based on this instrument,\footnote{Control function approaches have been applied in semi- and nonparametric sample selection models, e.g. Ahpo93, Daneve03, Ne07, Hu11, and Hu11b, as well as in nonparametric instrumental variable models, see for example NeweyPowellVella99, BlPo04, and ImNe09.} which requires further assumptions.

By the threshold crossing model postulated in 4(b), $\Pr(S = 1 | D, X , Z )= \Pr(V \leq \chi(D, X , Z ) ) = F_{V}(\chi(D, X , Z ))$, where $F_{V}(v)$ denotes the cumulative distribution function of $V$ evaluated at $v$. We will henceforth use the notation $\Pi=\pi(D,X,Z)=\Pr(S=1|D,X,Z)$ for the sake of brevity. Again by Assumption 4(b), the selection probability $\Pi$ increases strictly monotonically in $\chi$, such that there is a one-to-one correspondence between the distribution function $F_{V}$ and specific values $v$ given $X$. By Assumption 4(c), $V$ is independent of $(D,Z)$ given $X$, implying that the distribution function of $V$ given $X$ is (nonparametrically) identified. By comparing individuals with the same $\Pi$, we control for $F_V$ and thus for the confounding associations of $V$ with $D$ and $Y(d)$ that occur conditional on $S=1,X$. In other words, $\Pi$ serves as control function where the exogenous variation comes from $Z$. Controlling for the distribution of $V$ based on the instrument is thus a feasible alternative to the (infeasible) approach of directly controlling for levels of $V$.

Figure (ref) provides an acyclic graph of a causal model that can satisfy Assumptions 1 and 4. $U$ denotes unobservables affecting the outcome, which may be arbitrarily associated with $V$, the unobservable affecting selection. Note that the dashed lines indicate that $V,U$ are not observed in the data. Identification relies on instrument $Z$, which must not be associated with outcome $Y$ conditional on $D$ and $X$.\newline

Furthermore, identification requires the following common support assumption, which is similar to Assumption 3(a), but in contrast to the latter also includes $\Pi$ as a conditioning variable. \newline Assumption 5 (common support):\newline $\Pr(D=d| X=x, \Pi=\pi )>0$ for all $d \in \{0,1,...,Q\}$ and $x,z$ in the support of $X,Z$.\newline This means that in fully nonparametric contexts, the instrument $Z$ must in general be continuous and strong enough to importantly shift the selection probability $\Pi$ conditional on $D,M,X$ in the selected population. Assumptions 1, 4, and 5 are sufficient for the identification of mean potential outcomes and the ATE in the selected population, denoted as $\Delta_{S=1}=E[Y(1)-Y(0)| S=1]$.

To see this, note that the identifying assumptions imply

The first equality follows from $\Pi=F_V$ under Assumption 4, the second from the fact that when controlling for $F_V$, conditioning on $S=1$ does not result in an association between $Y(d)$ and $D$ given $X$ such that $Y(d) \bot D | X, \Pi, S=1$ holds by Assumptions 1 and 4. Therefore

Denoting by $p_d(X, \Pi)=\Pr(D=d |X, \Pi )$ and $\mu (D,S,X,\Pi)=E[Y |D,S,X, \pi(D,X,Z)]$, an alternative expression for the mean potential outcome among the selected is obtained by

where division by $p_d(X, \Pi)$ makes the reliance on Assumption 5 explicit. By applying the law of iterated expectations to replace $[Y- \mu(d,1,X,\Pi)]$ with $E[Y- \mu(d,1,X,\Pi)|D=d,S=1,X, \Pi]$ and noting that the latter expression is zero, one can see that (ref) is equivalent to (ref). But in contrast to the latter, the identification result in (ref) satisfies Neyman orthogonality and is based on the efficient influence function,\footnote{While the efficient influence function associated with (ref) is technically speaking doubly robust, i.e.\ consistent if either $\mu(d,1,X,\Pi)$ or $p_d(X, \Pi)$ is correctly specified, it is worth noting that this property can generally only hold if $\Pi$ is correctly specified because it enters both $\mu(d,1,X,\Pi)$ and $p_d(X, \Pi)$ as first step estimator. However, our approach does not rely on (global) doubly robustness but on Neyman orthogonality, which implies that DML is robust to local perturbations in $\Pi$ under specific regularity conditions.} see Appendix (ref).

The identification of the ATE in the total (rather than the selected) population is not feasible without further assumptions. The reason is that effects among selected observations cannot be extrapolated to the non-selected population if the effect of $D$ interacts with unobservables affecting the outcome, i.e.\ $U$ in (ref), as the latter are in general distributed differently across $S=1,0$ even conditional on $(X, \Pi)$ or $(D, X, \Pi)$. To see this, note that conditional on $\Pi=\Pr(V \leq \chi(D, X , Z ) )$, the distribution of $V$ differs across the selected (satisfying $V \leq \chi(D, X , Z )$) and the non-selected (satisfying $V > \chi(D, X , Z )$), such that the distribution of $U$ differs, too, if $V$ and $U$ are associated. This generally implies that $E[ Y(1)-Y(0)|S=1, X, \Pi]\neq E[ Y(1)-Y(0)|S=0, X, \Pi]$. While control function $\Pi$ ensures (together with $X$) that the treatment is unconfounded in the selected subpopulation, it does not permit extrapolating effects to the non-selected population with unobserved outcomes, see also HuberMelly for further discussion.

Assumption 6 therefore imposes homogeneity in the average treatment effect across selected and non-selected populations conditional on $X,V$. A sufficient condition for effect homogeneity is the separability of observed and unobserved components in the outcome equation, i.e.\ $Y=\eta(D,X)+\nu(U)$, where $\eta,\nu$ are general functions. Furthermore, common support as postulated in Assumption 5 needs to be strengthened to hold in the entire population. In addition, the selection probability $\Pi$ must be larger than zero for any $d,x,z$ in their support. Otherwise, outcomes are not observed for some values of $D,X$. Assumption 7 formalizes this common support restriction.\newline Assumption 6 (conditional effect homogeneity):\newline $E[Y(d)-Y(d')|S=1,X=x,V=v]=E[Y(d)-Y(d')|X=x,V=v]$ for all $d\neq d' \in \{0,1,...,Q\}$ and $x,v$ in the support of $X,V$.\newline Assumption 7 (common support):\newline $\pi(d,x,z)>0$ for all $d \in \{0,1,...,Q\}$ and $x,z$ in the support of $X,Z$.\newline

Under Assumptions 1,4,5,6, and 7, it follows that

where the first equality follows from Assumptions 1 and 4, see (ref), and the second one from Assumption 6. Therefore, the ATE is identified by

An alternative expression for the ATE that is based on the efficient influence function and respects Neyman orthogonality is given by

where division by $p_d(X,\Pi)\cdot \pi (d,X,Z)$ relies on the satisfaction of Assumptions 5 and 7.

Identification under sequential conditional independence

In many applications, it might appear unrealistic that one can control for all variables jointly affecting the sample selection indicator by conditioning only on baseline covariates measured prior to treatment assignment, in particular when no instrument is at hand. This is particularly the case when there is a substantial time lag between treatment assignment and sample selection/attrition, which raises concerns about dynamic confounding. The latter implies that some confounders influencing both the outcome and sample selection are themselves a function of the treatment. We subsequently reconsider the MAR framework, but not modify the identifying assumptions such that observed post-treatment confounders of $Y$ and $S$ are permitted. We will subsequently refer to observed post-treatment variables by $M$, in order to distinguish them from pre-treatment covariates $X$. Identification is based on a sequential conditional independence, which is based on maintaining Assumption 1 (conditional independence of $D$ given $X$), but replacing Assumption 2 by a modified conditional independence assumption for the selection indicator $S$ that allows for dynamic confounding due to $M=M(D)$, i.e.\ covariates possibly influenced by the treatment.\newline Assumption 8 (conditional independence of selection):\newline $ Y \bot S | D=d, X=x, M=m $ for all $d \in \{0,1,...,Q\}$ and $x,m$ in the support of $X$ and $M$. \newline By Assumption 8, there are no unobservables jointly affecting selection and the outcome conditional on $D,X,M$, such that sample selection is selective w.r.t.\ observed characteristics only. When modifying the nonparametric outcome and selection models in (ref) to $Y=\phi(D,X,M,U)$ and $S=\psi(D,X,M,V)$, Assumption 8 is satisfied if unobservables $U$ and $V$ are independent. \newline Assumption 9 (common support):\newline (a) $\Pr(D=d| X=x)>0$ and (b) $\Pr(S=1| D=d, X=x, M=m)>0$ for all $d \in \{0,1,...,Q\}$ and $x,m$ in the support of $X,M$.\newline Part (b) in Assumption 9 imposes a somewhat stronger common support restriction than part (b) in Assumption 3, as it requires the selection propensity score to be larger than zero for any combination of $D,X,M$ (rather $D,X$ only).

Figure (ref) provides an acyclic graph in which Assumptions 1 and 8 hold. Post-treatment covariates $M$ may be influenced by $D,X$ and might jointly affect $S$ and $Y$. Conditional on $D,X,M$, there are, however, no unobservables jointly influencing $S$ and $Y$.\newline

Our identifying assumptions imply that

where the first and third equalities follow from the law of iterated expectations, the second from Assumption 1, and the fourth from Assumption 8. Alternatively to this regression-based result using nested conditional mean outcomes, an IPW-based expression can be obtained, in which we use $\pi(D,X,M)=\Pr(S=1| D , X , M )$ as shortcut notation for the selection propensity score.

where the first and third equalities follow from basic probability theory and the second and last ones from the law of iterated expectations. Combining regression and IPW yields the following doubly robust identification result based on the efficient influence function, in which $\mu(d,1,X,M)=E[Y| D=d, S=1, X,M]$ and $\nu(d,1,X)=E[E[Y| D=d, S=1, X,M]|D=d,X]$ denote the conditional mean outcome and the nested conditional mean outcome, respectively:

where division by $p_d(X)\cdot \pi (d,X,M)$ relies on Assumption 9. The derivation of the efficient influence function is provided in Appendix (ref).

Estimation of the counterfactual with K-fold Cross-Fitting

We subsequently propose an estimation strategy for the counterfactual $E[Y(d)]$ under MAR as discussed in Section (ref) based on identification result (ref) and show its root-n consistency under specific regularity conditions. Let to this end $\mathcal{W} = \{W_i|1\leq i \leq n\}$ with $W_i = (Y_{i}\cdot S_{i}, D_{i}, S_{i}, X_{i})$ for all $ i $ denote the set of observations in an i.i.d.\ sample of size $n$. $\eta$ denotes the plug-in (or nuisance) parameters, i.e.\ the conditional mean outcome, mediator density and treatment probability. Their respective estimates are referred to by $\hat{\eta} = \{\hat\mu(D,1,X),\hat{p}_d(X), \hat{\pi} (D,X) \}$ and the true parameters by $\eta_0 = \{\mu_0(D,1,X), p_{d0}(X), \pi_0(D,X) \}$. Finally, $\Psi_{d0}=E[Y(d)]$ denotes the true counterfactual.

We estimate $\Psi_{d0}$ by the following algorithm that combines the estimation of Neyman-orthogonal scores with sample splitting or cross-fitting and is root-$n$ consistent under conditions outlined further below. \newline Algorithm 1: Estimation of $E[Y(d)]$ based on equation (ref)

In order to obtain root-n consistency for counterfactual estimation, we make the following assumption about the prediction qualities of machine learning for estimating the nuisance parameters. Following Chetal2018, we introduce some further notation: let $(\delta_n)_{n=1}^{\infty}$ and $(\Delta_n)_{n=1}^{\infty}$ denote sequences of positive constants with $\lim_{N\rightarrow \infty} \delta_n = 0 $ and $\lim_{N\rightarrow \infty} \Delta_n = 0.$ Furthermore, let $c, \epsilon, C$ and $q$ be positive constants such that $q>2,$ and let $K \geq 2$ be a fixed integer. Also, for any random vector $R = (R_1,...,R_l)$, let $\left\| R \right\|_{q} = \max_{1\leq j \leq l}\left\| R_l \right\|_{q},$ where $\left \| R_l \right\|_{q} = \left( E\left[ \left| R_l \right|^q \right] \right)^{\frac{1}{q}}$. In order to ease notation, we assume that $n/K$ is an integer. For the sake of brevity we omit the dependence of probability $\Pr_P,$ expectation $E_P(\cdot),$ and norm $\left\| \cdot \right\|_{P,q}$ on the probability measure $P.$ \newline

Assumption 10 (regularity conditions and quality of plug-in parameter estimates): \newline For all probability laws $P \in \mathcal{P},$ where $\mathcal{P}$ is the set of all possible probability laws the following conditions hold for the random vector $(Y,D,S,X)$ for $d \in \{0,1,...,Q\}$:

The only non-primitive condition is the condition (d), which puts restrictions on the quality of the nuisance parameter estimators. Condition (a) states that the distribution of the outcome does not have unbounded moments. (b) refines the common support condition such that the treatment and selection propensity scores are bounded away from $0$and $1$ and $0$, respectively. (c) states that covariates $X$ do not perfectly predict the conditional mean outcome.

For demonstrating the root-n consistency of our estimator of the mean potential outcome, we show that it satisfies the requirements of the DML framework in Chetal2018 by first verifying linearity and Neyman orthogonality of the score (see Appendix (ref)). As $ \psi_{d}(W, \eta, \Psi_{d0}) $ is smooth in $ (\eta, \Psi_{d0}) $, it then suffices that the plug-in estimators converge with rate $ n^{-1/4} $ for achieving $ n^{-1/2} $-convergence in the estimation of $ \hat{\psi}$, see Theorem 1. A rate of $ n^{-1/4} $ is achievable by many commonly used machine learners under specific conditions, such as lasso, random forests, boosting and neural nets, see for instance Bellonietal2014, LuoSpindler2016, WagerAthey2018, and FarrellLiangMisra2018.\newline Theorem 1\\ Under Assumptions 1-3 and 10, it holds for estimating $\Psi_{d0}=E[Y(d)]$ based on Algorithm 1: \\ $\sqrt{n} \Big(\hat \Psi_{d} - \Psi_{d0} \Big) \rightarrow N(0,\sigma^2_{\psi_d})$, where $\sigma^2_{\psi_{d}}= E[(\psi_d-\Psi_{d0})^2]$. \\ The proof is provided in Appendix (ref).

We subsequently discuss the estimation of $\Psi_{d0}$ based on (ref). We note that in this case, one needs to estimate the nested nuisance parameters $\mu(d,1,X,\Pi)$ and $p_d(X,\Pi)$, because they require the first-step estimation of $\Pi=\pi(D,X,Z)$. To avoid overfitting in the nested estimation procedure, the models for $\Pi$ on the one hand and $\mu(d,1,X,\Pi), p_d(X,\Pi)$ on the other hand are estimated in different subsamples. The plug-in estimates are now denoted by $\hat{\eta} = \{\hat\mu(D,1,X,\Pi),\hat{p}_d(X,\Pi), \hat{\pi} (D,X,Z) \}$ and the true plug-ins by $\eta_0 = \{\mu_0(D,1,X,\Pi), p_{d0}(X,\Pi), \pi_0(D,X,Z) \}$.\newline Algorithm 2: Estimation of $E[Y(d)]$ based on equation (ref)

An estimator of $\Psi_{d0}^{S=1}=E[Y(d)|S=1]$ based on (ref) is obtained by two modifications in Algorithm 2. First, rather than relying on the total sample $n$, one merely uses the subsample with observed outcomes which of size $\sum_{i=1}^n S_i$ to split it into $K$ subsamples. Second, in step 3, $\hat \phi_{d,i}^{k}$ is to be replaced by

to estimate $\Psi_{d0}^{S=1}$ in step 4 by $\hat \Phi_d^{S=1}=\frac{1}{\sum_{i=1}^n S_i} \sum_{k=1}^{K} \sum_{i=1}^{n_k} \hat \phi_{d,S=1,i}^{k}$. As $\sum_{i=1}^n S_i$ is an asymptotically fixed proportion of $n$, also this approach can be shown to be root-$n$ consistent under particular regularity conditions outlined in Assumption 11, that are in analogy to those in Assumption 10, but now adapted our IV-dependent identifying assumptions. \newline Assumption 11 (regularity conditions and quality of plug-in parameter estimates): \newline For all probability laws $P \in \mathcal{P},$ where $\mathcal{P}$ is the set of all possible probability laws the following conditions hold for the random vector $(Y,D,S,X,Z)$ for $d \in \{0,1,...,Q\}$:

Theorems 2 and 3 postulate the root-n consistency and asymptotic normality of the estimators of the mean potential outcomes in the selected and total populations, respectively. \newline Theorem 2\\ Under Assumptions 1, 4, 6, 7, and 11, it holds for estimating $\Psi_{d0}=E[Y(d)]$ based on Algorithm 2: \\ $\sqrt{n} \Big(\hat \Phi_{d} - \Psi_{d0} \Big) \rightarrow N(0,\sigma^2_{\phi_d})$, where $\sigma^2_{\phi_{d}}= E[(\phi_d-\Psi_{d0})^2]$. \\ Theorem 3\\ Under Assumptions 1, 4, 5, and 11, it holds for estimating $\Psi_{d0}^{S=1}=E[Y(d)|S=1]$ based on Algorithm 2: \\ $\sqrt{n} \Big(\hat \Phi_{d}^{S=1} - \Psi_{d0}^{S=1} \Big) \rightarrow N(0,\sigma^2_{\phi_{d,S=1}})$, where $\sigma^2_{\phi_{d,S=1}}= E[(\phi_{d,S=1}-\Psi_{d0}^{S=1})^2]$. \\ The proofs are provided in Appendices (ref) and (ref).

Next, we consider the estimation of $\Psi_{d0}$ based on (ref). Similarly to estimation based on (ref), we are required to estimate a nested nuisance parameter, namely $\nu(d,1,X)=E[\mu(d,1,X,M)|D=d,X]$. To avoid overfitting in the nested estimation procedure, the models for $\mu(d,1,X,M)$ and $\nu(d,1,X)$ estimated in different subsamples.\newline Algorithm 3: Estimation of $E[Y(d)]$ based on equation (ref)

To show root-$n$ consistency for this estimation approach, we impose the following regularity conditions, where we again assume that $n/K$ is an integer and omit the dependence of probability $\Pr_P,$ expectation $E_P(\cdot),$ and norm $\left\| \cdot \right\|_{P,q}$ on the probability measure $P$: \newline Assumption 12 (regularity conditions and quality of plug-in parameter estimates): \newline For all probability laws $P \in \mathcal{P}$ the following conditions hold for the random vector $(Y,D,S,X,M)$ for all $d \in \{0,1,...,Q\}$:

Under these regularity conditions and the sequential conditional independence assumption, estimation based on Algorithm 3 is asymptotically normal, as postulated in Theorem 4. \newline Theorem 4\\ Under Assumptions 1, 8, 9, and 12, it holds for estimating $E[Y(d)]$ based on Algorithm 3: \\ $\sqrt{n} \Big(\hat \Theta_{d} - \Psi_{d0} \Big) \rightarrow N(0,\sigma^2_{\theta_{d}})$, where $\sigma^2_{\theta_{d}}= E[(\theta_{d}-\Psi_{d0})^2]$. \\ The proof of Theorem 4 is provided in Appendix (ref).

Simulation study

This section provides a simulation study to investigate the finite sample behavior of our estimation approaches either relying on a MAR assumption of an instrument for selection based on the following data generating process:

Outcome $Y$ is a linear function of $D$ (whose treatment effect is one), covariates $X$ (for $\beta\neq 0$), and the unobservable $U$ and is only observed if the selection indicator $S$ is equal to one. Selection is a function of $D$, $X$, the unobservable $V$, and of instrument $Z$ if $\gamma\neq 0$. The treatment $D$ is a function of $X$ and the unobservable $W$. Both $Z$ and $W$ are random, standard normally distributed variables that are uncorrelated with $X$ or $(U,V)$. The correlation between the mean zero and normally distributed covariates in $X$ is determined by the covariance matrix $\sigma^2_X$. Similarly, $\sigma^2_{U,V}$ determines the correlation between the mean zero and normally distributed unobservables in the outcome and selection equation. In this setup, MAR is violated if the covariance between $U$ and $V$ is non-zero. We consider the perfromance of our estimators in $1000$ simulations with two sample sizes of $n=2000$ and $8000$.

In our simulations, we set the number of covariates $p$ to 100. $\sigma^2_X$ is defined based on setting the covariance of the $i$th and $j$th covariate in $X$ to $0.5^{|i-j|}$. $\beta$ gauges the impacts of the covariates on $Y$, $S$, and $D$, respectively, and thus, the magnitude of confounding. The $i$th element in the coefficient vector $\beta$ is set to $0.4/i^2$ for $i=1,...,p$, implying a squared decay of covariate importance in terms of confounding. In our first simulation design, we set $\gamma=0$ and $\sigma^2_{U,V}=\left(

\right)$ such that MAR as discussed in Section \ref{marident} holds. We consider the performance of DML based on Theorem 1 (henceforth DML MAR), which does not make use of the instrument $Z$, as well as based on Theorem 2 (DML IV), which exploits the instrument despite the satisfaction of MAR.

The nuisance parameters, i.e.\ the linear and probit specifications of the outcome, selection, and treatment equations, are estimated by lasso regressions using the default options of the SuperLearner package provided by vanderLaanetal2007 for the statistical software R. We use 3-fold cross-fitting for the estimation of the treatment effects. We drop observations whose products of estimated treatment and selection propensity scores are close to zero, namely smaller than a trimming threshold of $0.01$ (or 1%). This avoids an explosion of the propensity score-based weights and, thus, of the variance when estimating the mean potential outcomes or ATE by the sample analogues of (ref) and (ref), where the product of the propensity scores enters the respective denominators for reweighing the outcome. Our estimation procedure is available in the treatselDML command of the causalweight package for R by BodoryHuber2018.

Table (ref) presents the simulation results. The biases (bias) of both DML MAR and DML IV are rather close to zero independent of the sample size. Furthermore, the estimators have virtually the same variance, despite the fact that DML IV unnecessarily relies on the control function approach and an irrelevant instrument. Both estimators appear to converge to the true effect with $\sqrt{n}$-rate, as the root mean squared error (RMSE) is roughly cut by half when quadrupling the sample size. The average standard error across simulations (meanSE) based on the asymptotic variance approximation comes close to the respective estimator's standard deviation (sd). Finally, the coverage rate (coverage), i.e.\ the share of simulations in which the 95% confidence interval includes the true effect, is only slightly below the nominal level of 95%.

In a second simulation design, we set $\gamma=1$ and $\sigma^2_{U,V}=\left(

\right)$, such that selection is nonignorable, i.e.\ related to unobservables as discussed in Section \ref{ivident}, due to a strong correlation of $U$ and $V$. Table (ref) presents the results. DML MAR is no longer unbiased, while the bias of DML IV appears to approach zero as the sample size increases, at the price of somewhat higher standard deviation than DML MAR. However, DML IV dominates DML MAR under either sample size in terms of having a lower RMSE and has thus a more favorable bias-variance trade-off in the scenario considered. While coverage is quite satisfactory for DML IV, the 95% confidence interval mostly fails to include the true effect in the case of DML MAR, in particular under the larger sample size.

Application

As an empirical illustration, we apply our method to the Job Corps (JC) training program. The data come from the National Job Corps Study (NJCS), a randomized social experiment conducted in the mid-to-late 1990s in the United States to evaluate the effectiveness of JC on different labor market outcomes. The JC is the largest and most comprehensive job training program for disadvantaged youth in the US, in which participants are exposed to different types of academic and vocational instruction. The data set contains very detailed pre-treatment information about program participants, such as: expectations, motivations for applying to JC, age, gender, number of children at the moment of treatment assignment, occupation, household income, hourly wage, educational level, marital status, whether the individual was previously attending a school, JC program, or some other academic or vocational training, health status, past employment, types of crimes committed, family support for attending the training, and information on the mother and the father (e.g. education and employment). Furthermore, a large range of variables for instance related to labor market status, employment, income, and education are reassessed in several follow-up interviews after JC assignment.

SchBuGl01 and SchBuMc08 evaluate the impact of random program assignment on a wide range of labor market outcomes, showing positive effects on education, employment, and earnings in the longer run. Several studies focus on more specific program aspects of JC, e.g.\ on the effect of the length of exposure to training or of discrete sequences of training interventions on labor market outcomes like employment and earnings, as well as on JC's causal mechanisms, i.e.\ direct and indirect effects (operating via specific mediating variables) on labor market outcomes and health. FlGoNe12, for instance, acknowledge the existence of several types of instruction, which, along with the self-paced nature of the program, creates selective heterogeneity in the number of weeks the participants are exposed to vocational or academic training. Considering continuously distributed treatment doses while controlling for baseline covariates, they find a positive effect of weeks in training on earnings, however, with decreasing marginal returns as a function of weeks already accomplished. FlFl09 and Hu14 investigate the causal mechanisms underlying JC when considering work experience or employment as mediators, respectively, and find positive direct effects on earnings and general health, respectively, when invoking a selection-on-observables assumption. FlFl10 avoid the latter assumption by suggesting a partial identification strategy based on which they estimate bounds the on causal mechanisms of JC when considering the achievement of a GED, high school degree, or vocational degree as mediators. Under their strongest set of bounding assumptions, the results suggest a positive effect on labor market outcomes even net of the indirect mechanism via obtaining a degree. FrHu17 base their analysis of causal mechanisms on an instrumental variable approach and find a positive indirect effect of JC training on earnings through an increase in the number of hours worked.

Estimation approaches specifically dealing with truncated outcomes such as wages, which are only observed and defined for a selective subpopulation like those in employment, have also been considered. For instance, FrMePaRu12 and ZhRuMe09 consider the principal stratification approach of FrangakisRubin02 to evaluate the effects of JC on employment and wages among specific groups, e.g.\ those finding employment irrespective of training participation, rather among the full population. FrMePaRu12 simultaneously address the several identification issues related to noncompliance of training participation with JC assignment as well as missingness in wage outcomes due to survey non-response or non-employment based on a likelihood-based analysis using finite mixture models. They find positive effects on wages and evidence that the program should ideally have been designed differently for different subgroups of individuals depending on their personal characteristics.

Le05 considers a partial identification approach for bounding the effect of JC assignment on wage among those finding employment irrespective of the assignment. Rather than invoking MAR or IV assumptions for sample selection, this method merely relies on the monotonicity of employment in JC assignment (such that being randomized in never decreases the employment state), at the cost of giving up point identification. semenova2020 suggests a DML approach to tighten the bounds by controlling for covariates $X$ in a data-driven way. The results based on bounding generally point to a positive (intention-to-treat) effect of JC assignment on wages. Finally, BoHuLa2020 consider dynamic treatment evaluation based on DML, imposing a sequential conditional independence assumption to analyze discrete sequences of training in the first and second year after JC assignment. Controlling for both baseline characteristics and covariates measured after the first treatment period (i.e.\ one year after JC assignment) they find a positive effect of a sequence of vocational training on employment.

For our empirical analysis we similarly to FrHu17 consider female applicants to JC and aim at estimating the effects of academic or vocational training received in the first year of the program ($D$) on hourly wage ($Y$) in the short run measured in the last week of the first year or in the longer run, measured 4 years after random assignment to JC. Hourly wage is only observed conditional on employment ($S$) in the respective outcome period. Even though JC assignment is random, actual participation in training activities is likely selective and associated with individual factors, similarly to the selection into employment. As for example discussed in LECHNER2013111 and BiFiOsPau14, the previous labor market history and socio-economic characteristics are likely important confounders when assessing the impact of training interventions, which motivates our DML approach to account for a rich set of covariates in a data-driven way. To assess the short-run effects, we either assume MAR as discussed in Section (ref) and use our DML approach based on Theorem 1 to control for our all in all $355$ baseline covariates ($X$), or we impose the IV assumptions of Section (ref) to estimate the ATE based on Theorem 3, considering the number of young children in the household at JC assignment as instrument ($Z$) for employment. Even though numerous studies in labor economics consider children as instrument for employment, the presence of small children might arguably be associated with personal characteristics also affecting the wage and thus violate IV validity - a concern we aim to mitigate by including a rich set of individual pre-treatment characteristics in $X$ that is likely associated with both fertility and wages.

For assessing the longer run effects, we invoke the sequential conditional independence assumption of Section (ref) to apply DML based on Theorem 4. To this end, we additionally control for $619$ and $156$ post-treatment covariates ($M$) in the second and third year, respectively, after JC assignment, which include detailed information on the labor market participation after the first and prior to the second treatment. Appendix (ref) presents descriptive statistics for selected variables in $X$ and $M$. We also refer to BoHuLa2020 for a more detailed description of the pre-and post-treatment covariates used in our application and note that all numeric variables have been standardized to have a mean equal to $0$ and standard deviation equal to $0.5$ to facilitate the machine learning-based estimation of the nuisance parameters. For estimation, we apply the treatselDML and dyntreatDML commands of the causalweight package for R, using 3-fold cross-fitting and the random forest (with default options of the SuperLearner package) as machine learner. The latter is a nonparametric approach allowing for nonlinear associations between the outcome, treatment, and selection on the one hand and the covariates on the other hand.

Table (ref) reports the total number of females randomized into JC for whom participation in either vocational ($843$) or academic training ($830$) in the first year after program assignment is registered in the data. $200$ females did not participate in any JC training activity in the first year and serve as the control group in our analysis. Furthermore, $1698$ were randomized out of JC.

Table (ref) reports the ATE estimates for academic and vocational training based on our various DML approaches. The upper panel provides the short-run effects on hourly wages in the last week of the first year when assuming MAR. The point estimate of academic training is negative ($-0.683$ US \$), while that of vocational training is positive ($0.611$ US \$), but neither effect is statistically significant at any conventional level. The findings are very similar when considering the IV-based estimates shown in the intermediate panel, with negative and positive effects for the academic and vocational training, respectively, which are again not statistically significant. The lower panel provides the longer-run effects on hourly wages 4 years after assignment based on the sequential conditional independence assumption. While now both ATE estimates are positive, only the effect of vocational training, which amounts to an hourly increase of $0.567$ \$, is highly statistically significant. Our findings therefore suggest that JC-based education may facilitate human capital accumulation in a way that increases hourly wages after several years, in particular through vocational training, while there is no clear-cut evidence for short term effects.

Conclusion

In this paper, we discussed the evaluation of average treatment effects in the presence of sample selection or outcome attrition based on double machine learning. In terms of identifying assumptions, we imposed a selection-on-observables assumption on treatment assignment, which was combined with either selection-on-observables or instrumental variable assumptions concerning the outcome attrition/sample selection process. We also considered a sequential selection-on-observables assumption allowing for dynamic confounding such that covariates jointly affecting the outcome and sample selection may be affected by the treatment, which avoids exclusively relying on pre-treatment covariates. We proposed doubly robust score functions and formally showed the satisfaction of Neyman orthogonality, implying that estimators based on these score functions are robust to moderate (local) regularization biases in the machine learning-based estimation of the outcome, treatment, or sample selection models. Furthermore, we demonstrated the root-n consistency and asymptotic normality of our double machine learning approach to average treatment effect estimation under specific regularity conditions. We also provided an empirical illustration to the US Job Corps data, in which we assessed the effects of training on hourly wage one and four years after program assignment and found some statistical evidence for positive longer-run impacts. Our estimation procedure is available in the causalweight package for the statistical software R. \typeout