Doubly Robust Estimation of Direct and Indirect Quantile Treatment Effects with Machine Learning

Doubly Robust Estimation of Direct and Indirect Quantile Treatment Effects with Machine Learning

Introduction

Causal mediation analysis aims at understanding the mechanisms through which a treatment affects an outcome of interest. It disentangles the treatment effect into an indirect effect, which operates through a mediator, and a direct effect, which captures any causal effect not operating through the mediator. Such a decomposition of the total treatment effect permits learning the drivers of the effect, which may be helpful for improving the design of a policy or intervention. Causal mediation analysis typically focuses on the estimation of average indirect and direct effects, which may mask interesting effect heterogeneity across individuals. For this reason, several contributions focusing on total (rather than direct and indirect) effects consider quantile treatment effects (QTE) instead of average treatment effects (ATE). The QTE corresponds to the difference between the potential outcomes with and without treatment at a specific rank of the potential outcome distributions, but has so far received little attention in the causal mediation literature.

The main contribution of this paper is to propose doubly robust/debiased machine learning (DML) estimators of the direct and indirect QTE under a selection-on-observables (or sequential ignorability) assumption, implying that the treatment and the mediator are as good as random when controlling for observed covariates. The method computes the quantile of a potential outcome by inverting an DML estimate of its cumulative distributional function (c.d.f.). This approach makes use of the efficient score function of the c.d.f., into which models for the outcome, treatment, and mediator enter as plug-in or nuisance parameters. Relying on the efficient score function makes treatment effect estimation robust, i.e., first-order insensitive to (local) misspecifications of the nuisance parameters, a property known as Neyman1959-orthogonality. This permits estimating the nuisance parameters by machine learning (which generally introduces regularization bias) and still obtains root-n-consistent treatment effect estimators, given that certain regularity conditions hold. In addition, cross-fitting is applied to mitigate overfitting bias. Cross-fitting consists of estimating the nuisance parameter models and treatment effects in different subsets of the data and swapping the roles of the data to exploit the entire sample for treatment effect estimation, see CCDDHNR_2018. We then establish uniform consistency and asymptotic normality of the effect estimators.

For conducting statistical inference, we propose a multiplier bootstrap procedure and show the validity of the multiplier bootstrap. We also provide a simulation study to investigate the finite sample performance of our method. Finally, we apply our method to empirical data from the Job Corps Study to analyse the direct and indirect QTE of participation in a training program on earnings when considering general health as a mediator. The results point to positive direct effects of training across a large range of the earnings quantiles, while the indirect effects are generally close to zero and mostly statistically insignificant.

To more formally discuss the direct and indirect effects of interest, let $Y$ denote the outcome of interest, $D$ the binary treatment, $M$ the mediator, and $X$ a vector of pre-treatment covariates. Following Pearl00, we may represent causal relationships between $\left(Y,D,M,X\right)$ by means of a directed acyclic graph (DAG), as provided in Figure (ref). The causal arrows in the DAG imply that $\left(D,M,X\right)$ may affect $Y$, $\left(D,X\right)$ may affect $M$, and $X$ may affect $D$. We can therefore define the outcome as a function of the treatment and the mediator, $Y=Y\left(D, M\right)$, and the mediator as a function of the treatment, $M=M\left(D\right)$, while being agnostic about $X$. Furthermore, we make use of the potential outcome notation advocated by Neyman23 and Rubin74 to denote by $Y\left(d,m\right)$ the potential outcome if $D$ were set to a specific value $d\in\left\{ 0,1\right\}$ and $M$ were set to some value $m$ in the support of the mediator, while $M\left(d\right)$ denotes the potential mediator for $D=d$. Accordingly, $Y\left(d, M\left(d\right)\right)$ is the potential outcome if $D$ were set to $d$, implying that the mediator is not forced to take a specific value $m$, but corresponds to its potential value under $D=d$. Depending on the actual treatment and mediator values of an observation, $Y\left(d, M\left(d\right)\right)$, $Y\left(d,m\right)$, $M\left(d\right)$ is either observed or counterfactual. Furthermore, the potential outcome $Y\left(d, M\left(1-d\right)\right)$ is inherently counterfactual, as no observation can be observed in the opposite treatment states $d$ and $1-d$ at the same time.

Armed with this notation, we define the causal parameters of interest. The natural direct effect (NDE), which is for instance considered in RoGr92, Pearl01 and TS_2012, is based on a comparison of the potential outcomes when varying the treatment, but keeping the potential mediator fixed at treatment value $D=d$: $Y\left(1,M\left(d\right)\right)-Y\left(0,M\left(d\right)\right)$. The natural indirect effect (NIE) is based on a comparison of the potential outcomes when fixing the treatment to $D=d$, but varying the potential mediator according to the values it takes under treatment and non-treatment: $Y\left(d,M\left(1\right)\right)-Y\left(d,M\left(0\right)\right)$. It is worth noting that the NDE and the NIE defined upon opposite treatment states $d$ and $1-d$ sum up to the total effect (TE): $Y\left(1,M\left(1\right)\right)-Y\left(0,M\left(0\right)\right)$.

Previous methodological research on causal mediation predominantly focused on the estimation of averages of the aforementioned NDE, NIE and TE or of averages of related path-wise causal effects IKY_2010,TS_2012,HHL_2019,FHLLS_2022,Zhou_2022. We complement this literature by suggesting a method for estimating natural direct and indirect QTEs, which permits assessing the effects across the entire distribution of potential outcomes. The estimation of the total (rather than the direct or indirect) QTE has already been studied in multiple contributions AAI_2002,CH_2005,Firpo_2007,DH_2014,BCFH_2017,ALZ_2022,HLL_2022. Among the few studies considering QTEs in causal mediation is BVSC_2017, suggesting a two-stage quantile regression estimation to estimate the controlled direct QTE, i.e., $Y\left(1,m\right)-Y\left(0,m\right)$ at a specific rank, as well as a particular indirect QTE. The latter is based on first estimating the mediator at a specific rank and then including it in a quantile regression of the outcome, which generally differs from the natural indirect QTE considered in this paper. Furthermore, our approach is nonparametric and relies on results on semiparametric efficiency, very much in contrast to the parametric approach of BVSC_2017. HSS_2022 adapted the Changes-in-Changes (CiC) approach of AtheyImbens06 to estimate direct and indirect QTEs in subgroups defined in terms of how the mediator reacts to (or complies with) the treatment. The NDE and NIE investigated here differ from such subgroup-specific causal parameters and furthermore, our identification strategy relies on a selection-on-observables (or sequential ignorability) assumption rather than CiC.

The remainder of this study is organized as follows. Section (ref) introduces the natural direct and indirect QTE, the identifying assumptions, the effect estimators based on double/debiased machine learning, and the multiplier bootstrap procedure for inference. Section (ref) gives the theoretical results on the asymptotic behavior of our methods. Section (ref) presents a simulation study that investigates the finite sample properties of our method. Section (ref) provides an empirical application to data from the National Job Corps Study to assess the direct earnings effects of training across the earnings distribution, as well as the indirect effects operating via general health. Section (ref) concludes.

Methodology

Causal effects and Identifying Assumptions

To define the direct and indirect QTEs of interest, let $Q_{Z}\left(\tau\right):=\inf\left\{ q\in\mathbb{R}:P\left(Z\leq q\right)\geq\tau\right\} $ denote the $\tau$th-quantile of a random variable $Z$, where $\tau\in(0,1)$. Furthermore, let $Q_{Z|V}\left(\tau\right):=\inf\left\{ q\in\mathbb{R}:P\left(Z\leq q|V\right)\geq\tau\right\}$ denote the $\tau$th-quantile of $Z$ conditional on another random variable (or a random vector) $V$, where $\tau\in(0,1)$. Let $F_{Z}\left(z\right)$ and $f_{Z}\left(z\right)$ denote cumulative distribution function (c.d.f.) and probability density or probability mass function (p.d.f. or p.m.f.) of $Z$ at $z$, and $F_{Z|V}\left(z|v\right)$ and $f_{Z|V}\left(z|v\right)$ denote the c.d.f. and p.d.f. (or p.m.f.) of $Z$ at $z$ conditional on $V=v$. We define the natural direct quantile treatment effect (NDQTE) at the $\tau$th-quantile as:

and the natural indirect quantile treatment effect (NIQTE) at the $\tau$th-quantile as:

The NDQTE in equation ((ref)) corresponds to the direct effect of the treatment when fixing the mediator at its potential value under non-treatment, $M(0)$. Alternatively, we may consider the NDQTE when conditioning on the potential mediator under treatment, $M(1)$:

Likewise, the NIQTE in equation ((ref)) is the indirect effect when varying the potential mediators but keeping the treatment fixed ad $D=1$, but we may also consider the indirect effect conditional on $D=0$:

If the effects in expressions ((ref)) and ((ref)) (or ((ref)) and ((ref))) are different, then this implies effect heterogeneity due to interaction effects between the treatment and the mediator. The sum of NDQTE (NDQTE') and NIQTE (NIQTE') yields the total quantile treatment effect (TQTE) at the $\tau$th-quantile, which includes all causal mechanisms through which the treatment affects the outcome:

We aim at estimating the quantile treatment effects ((ref)) to ((ref)).\footnote{We do not consider the controlled direct quantile treatment effect (CDQTE) at the $\tau$-quantile, $\text{CDQTE}\left(\tau\right):=Q_{Y\left(1,m\right)}\left(\tau\right)-Q_{Y\left(0,m\right)}\left(\tau\right)$, which can be identified under less stringent assumptions than required for the identification of natural effects. } To this end, we first need to estimate the $\tau$th-quantile of the relevant potential outcomes, by inverting estimates of the corresponding c.d.f.'s at the $\tau$th-quantile. Let $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$ denote the c.d.f. of the potential outcome $Y\left(d,M\left(d^{\prime}\right)\right)$ at value $a$. To identify $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$ in the data, we impose the following assumptions.

Assumption 1.1 implies the stable unit treatment value assumption (SUTVA), see Cox58 and Rubin80, stating the potential mediators and potential outcomes are only a function of an individual's own treatment and mediator states, respectively, which are well defined (ruling out multiple treatment or mediator versions). Assumptions 1.2 and 1.3 are sequential ignorability or selection-on-observables conditions IKY_2010 for causal mediation analysis. Assumption 1.2 states that conditional on $X$, the treatment variable $D$ is independent of the potential outcome $Y(d,m)$ and the potential mediator $M(d^{\prime})$. This assumption also implies that $Y\left(d,m\right)\perp D|M\left(d^{\prime}\right)=m^{\prime},X=x$. Assumption 1.3 requires that $Y(d,m)$ and $M(d^{\prime})$ are independent, too, conditional on $X$ and $D$. Even if treatment $D$ were random, this would not suffice to identify direct and indirect effects and for this reason, we need to impose an identifying assumption like Assumption 1.3 to tackle the endogeneity of the mediator. Assumption 1.4 is a common support condition, which says that the treatment is not deterministic in covariates $X$ and mediator $M$ such that for each covariate-mediator combination in the population, both treated and non-treated subjects exist.

Under Assumptions 1.1 to 1.4, we obtain the following identification result.

The proof of Proposition 1 is provided in the appendix. Under Proposition 1, we may estimate $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$ based on plug-in estimation of the nuisance parameters $F_{Y|D,M,X}\left(a|d,m,X\right)$ and $f_{M|D,X}\left(m|d^{\prime},X\right)$:

where

and $\hat{F}_{Y|D,M,X}\left(a|d,m,X_{i}\right)$ and $\hat{f}_{M|D,X}\left(m|d^{\prime},X_{i}\right)$ are estimates of $F_{Y|D,M,X}\left(a|d,m,X\right)$ and $f_{M|D,X}\left(m|d^{\prime},X_{i}\right)$. If $M$ is a continuous variable, we may avoid estimating the conditional density $f_{M|D,X}\left(m|d^{\prime},X\right)$, and use the following alternative estimator for estimating $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$:

where $\hat{E}\left[F_{Y|D,M,X}\left(a|d,M_{i},X_{i}\right)|d^{\prime},X_{i}\right]$ is an estimate of $E\left[F_{Y|D,M,X}\left(a|d,M_{i},X_{i}\right)|d^{\prime},X_{i}\right]$. For example, it might be based on a “regression-imputation” Zhou_2022, corresponding to the fitted value of a linear regression of $F_{Y|D,M,X}\left(a|d,X_{i},M_{i}\right)$ on $D_{i}$ and $X_{i}$ at $\left(d^{\prime},X_{i}\right)$.

The quality of estimators ((ref)) and ((ref)) crucially depends on the accuracy of nuisance parameter estimation. If the number of pretreatment covariates $X$ is small (low dimensional $X$) and the functional forms of the nuisance parameters are known, parametric methods can provide high-quality estimations on the nuisance parameters. In contrast, if $X$ is high dimensional and/or the nuisance parameters have complex forms, machine learning may be the preferred choice of estimation. However, applying ML directly to estimate expressions ((ref)) or ((ref)) may result in non-negligible bias induced by regularization and/or overfitting CCDDHNR_2018. Causal machine learning algorithms aim at avoiding such biases by applying ML estimation when making use of Neyman-orthogonal moment conditions, which imply that the estimation of causal parameters is first order insensitive to (regularization) bias in the nuisance parameters, and of cross-fitting, which avoids overfitting. One of these causal algorithms is double/debiased machine learning (DML) CCDDHNR_2018, which has been previously adapted to the estimation of average effects in causal mediation analysis FHLLS_2022, while this study extends it to the estimation of direct and indirect quantile treatment effects.

Let $Y_{a}=1\{Y\leq a\}$ be an indicator function which is one if outcome $Y$ is smaller than or equal to $a$ (and zero otherwise) and $W_{a}=\left(Y_{a}, D, M, X\right)$ be a vector of the observed variables. An estimator of the c.d.f. of the potential outcome that satisfies Neyman-orthogonality can be derived from the efficient influence function (EIF) of $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$:

where

for $\left(d,d^{\prime}\right)\in \{0,1\}^{2}$, $a\in \mathcal{A}$ and $v_{a}$ denoting the vector of nuisance parameters. Let $\theta_{d,d^{\prime},a}$ denote the value of $\theta$ that satisfies $E[\psi_{d,d^{\prime},a}^{\theta}\left(W_{a};v_{a}\right)]=0$:

We can show that $\theta_{d,d^{\prime},a}=F_{Y(d,M(d^{\prime})}\left(a\right)$ if Assumptions 1.1 to 1.4 hold, see the appendix for a derivation of these results. Therefore, we may use the sample analogue of equation ((ref)) to estimate $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$. A similar strategy was previously used to derive the triply robust approach for estimating $E\left[Y\left(d,M\left(d^{\prime}\right)\right)\right]$ in TS_2012 and FHLLS_2022. If $d=d^{\prime}$, then the estimator of equation ((ref)) reduces to

where

for $d\in \{0,1\}$ and $a\in \mathcal{A}$. We may use the sample analogue of equation ((ref)) to estimate $F_{Y\left(d,M\left(d\right)\right)}\left(a\right)$. This is in analogy to the doubly robust approach for estimating $E\left[Y\left(d,M\left(d\right)\right)\right]$ in RRZ_1994 and Hahn_1998.

We note that by using the Bayes rule, we can rewrite equation ((ref)) alternatively as:

Therefore,

can also be used to construct an estimator of $F_{Y(d,M(d^{\prime})}\left(a\right)$. There are several differences between the estimators based on equations ((ref)) and ((ref)). Making use of equation ((ref)) requires estimating four nuisance parameters: $f_{D|X}\left(d|x\right)$, $f_{D|M,X}\left(d|m,x\right)$, $F_{Y|D,M,X}\left(a|d,m,x\right)$ and $g_{d,d^{\prime},a}\left(x\right)$. Since $D$ is binary, the first two nuisance parameters may for instance be estimated by a logit or probit model. The conditional c.d.f. $F_{Y|D,M,X}\left(a|d,m,x\right)$ can be estimated by distributional regression (DR) CFM_2013. $g_{d,d^{\prime},a}\left(x\right)$ might be estimated by regression imputation as outlined in equation ((ref)). However, the estimator based on equation ((ref)) requires only three nuisance parameter estimates of $f_{D|X}\left(d|x\right)$, $f_{M|D,X}\left(m|d,x\right)$ and $F_{Y|D,M,X}\left(a|d,m,x\right)$. We may estimate $g_{d,d^{\prime},a}\left(x\right)$ based on equation ((ref)) after having estimated $f_{M|D,X}\left(m|d,x\right)$ and $F_{Y|D,M,X}\left(a|d,m,x\right)$. The estimator based on equation ((ref)) appears particularly attractive if the mediator $M$ is discrete and takes a finite (and relatively small) number of values. However, if $M$ is continuous, estimation based on equation ((ref)) may appear more attractive, because it avoids estimating the conditional density $f_{M|D,X}\left(m|d,x\right)$ and the integral in equation ((ref)) to obtain an estimate of $g_{d,d^{\prime},a}\left(x\right)$.

The estimators based on equations ((ref)) and ((ref)) also differ in terms of their robustness to misspecification of the nuisance parameters. Let $\hat{\theta}_{d,d^{\prime},a}$ and $\hat{\theta}_{d,d^{\prime},a}^{\prime}$ denote estimators based on equations ((ref)) and ((ref)) and the respective estimators of the nuisance parameters. Applying the theorem of semiparametric efficiency in TS_2012 and Zhou_2022, we can show that under certain regularity conditions, the following results hold for estimating $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$ at outcome value $a$:

This implies that if two of the nuisance parameters entering equation ((ref)) are correctly specified, while also certain regularity conditions and Assumptions 1.1 to 1.4 hold, then $\hat{\theta}_{d,d^{\prime},a}^{\prime}$ is a consistent estimator of $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$ at outcome value $a$. In contrast, $\hat{\theta}_{d,d^{\prime},a}\overset{p.}{\longrightarrow}\theta_{d,d^{\prime},a}$ only holds if $f_{D|X}(d|x)$ is consistently estimated. If the latter holds and only one of the other three nuisance parameters in equation ((ref)) is misspecified, while certain regularity conditions and Assumptions 1.1 to 1.4 are satisfied, then $\hat{\theta}_{d,d^{\prime},a}$ remains a consistent estimator of $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$ at outcome value $a$. Finally, if all nuisance parameters are correctly specified and consistently estimated, while certain regularity conditions and Assumptions 1 to 4 also hold, then both $\hat{\theta}_{d,d^{\prime},a}$ and $\hat{\theta}_{d,d^{\prime},a}^{\prime}$ are semiparametrically efficient.

Improving Finite Sample Behavior

The estimate of the c.d.f. of $Y\left(d,M\left(d^{\prime}\right)\right)$ can be inverted at a specific rank $\tau$ to obtain an estimate of the $\tau$th quantile, which we denote by $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(\tau\right)$. Suppose $Y\left(d,M\left(d^{\prime}\right)\right)$ is continuous and let the grid of points used for the estimation be a non-decreasing sequence $\{a_{l}\}_{l=1}^{L}$, where $0<\underline{a}<a_{1}<a_{2}<\ldots<a_{L}<\bar{a}<\infty$. Let $\hat{p}_{l}$ denote an estimate of $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a_{l}\right)$ (e.g., the K-fold cross-fitting estimate, see Section 2.2). Note that $\hat{p}_{l}$ is not necessarily bounded away from 0 and 1 nor monotonically increasing in $a_{l}$, as required for a valid c.d.f. For this reason, we apply two additional constraints on the estimates $\hat{p}_{l}$, $l=1,\ldots,L$. The first one restricts their values to be within the range $\left[0,1\right]$. That is, we replace $\hat{p}_{l}$ with $\tilde{p}_{l}=\max\left\{ \min\left\{ \hat{p}_{l},1\right\} ,0\right\} $. Then we follow CFG_2010 and use the rearrangement operator to sort $\tilde{p}_{l}$ in non-decreasing order. Let $\left(\tilde{p}_{\left(1\right)},\tilde{p}_{\left(2\right)},\ldots,\tilde{p}_{\left(L\right)}\right)$ be the sorted sequence of $\tilde{p}_{l}$, $l=1,2,\ldots,L$. The sequence $\left(\tilde{p}_{\left(1\right)},\tilde{p}_{\left(2\right)},\ldots,\tilde{p}_{\left(L\right)}\right)$ is our final estimate of the c.d.f. of $Y\left(d,M\left(d^{\prime}\right)\right)$ at $\left(a_{1},a_{2},\ldots,a_{L}\right)$. We then fit a function for points $\left(\tilde{p}_{\left(l\right)},a_{l}\right)$, $l=1,2,\ldots,L$ with linear interpolation, and use the fitted function to calculate the value of $a$ at rank $\tau$ to estimate $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(\tau\right)$,\footnote{The monotonicity property is preserved under the linear interpolation.}, which permits estimating the quantile treatment effects ((ref)) to ((ref)). When $Y\left(d,M\left(d^{\prime}\right)\right)$ is discrete, we need not fit the function for points $\left(\tilde{p}_{\left(l\right)},a_{l}\right)$; we may obtain $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(\tau\right)$ by directly using the definition of the $\tau$th-quantile.

K-Fold Cross-Fitting

Neyman-orthogonality may mitigate regularization bias coming from machine learning-based estimation of the nuisance parameters in equations ((ref)) or ((ref)). To also safeguard against overfitting bias, we follow CCDDHNR_2018 and FHLLS_2022 and apply K-fold cross-fitting to estimate the nuisance parameters and the potential outcome distributions, $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$, in different parts of the data. To describe the approach, let $Y_{a,i}=1\{Y_{i}\leq a\}$ and $W_{a,i}=\left(Y_{a,i},D_{i},M_{i},X_{i}\right)$ denote the $i$th observation, $i=1,2,\ldots,n$. In the following, we use the estimator based on equation ((ref)) to illustrate K-fold cross-fitting.

We repeat steps 1 to 4 for a grid of points $a$, a non-decreasing sequence $\{a_{l}\}_{l=1}^{L}$, where $0<\underline{a}<a_{1}<a_{2}<\ldots<a_{L}<\bar{a}<\infty$, to construct the estimate of the c.d.f.\ profile of $Y\left(d,M\left(d^{\prime}\right)\right)$. In Section 3, we will establish the asymptotic properties uniformly over $a\in \cal{A}$ for the K-fold cross-fitting estimator of equation ((ref)).

Alternatively, we can also construct the K-fold cross-fitting estimator based on equation ((ref)):

where

Nuisance Parameter Estimation

In the case when $D$ and $M$ are binary variables, $f_{D|X}\left(d|x\right)$ and $f_{M|D,X}\left(m|d,x\right)$ can be estimated straightforwardly, for instance by logit or probit models. If $M$ is continuous, we might prefer $\hat{\theta}_{d,d^{\prime},a}$ to avoid estimating $f_{M|D,X}\left(m|d,x\right)$, while the required nuisance parameter $f_{D|M,X}\left(d|m,x\right)$ can be straightforwardly estimated if $D$ is binary. The estimation of any nuisance parameter may be based on machine learning, for example, the lasso or neural networks. For estimating $F_{Y|D,M,X}\left(a|d,m,x\right)$, we use distributional regression (DR). Conditional on $\left(D=d,M=m,X=x\right)$, the conditional c.d.f. of $Y$ may be written as

where $a$ is a constant and $a\in\mathcal{A}\subset\mathbb{R}$, with $\mathcal{A}$ being a countable subset of $\mathbb{R}$, and $Y_{a}=1\left\{ Y\leq a\right\} $ being an indicator function for the event $\left\{ Y\leq a\right\} $. Equation ((ref)) is the building block for constructing the DR estimator FP_1995, CFM_2013, which is based on using the binary dependent variable $Y_{a}$ to estimate $F_{Y|D,M,X}\left(a|d,m,x\right)$ by a regression approach that estimates $E\left[Y_{a}|d,m,x\right]$. For example, one may assume that the c.d.f. is linear in variables $\left(D,M,X\right)$, their higher-order terms and interaction terms, and estimate a linear probability model (LPM) by OLS. However, the LPM does not guarantee that the estimated $F_{Y|D,M,X}\left(a|d,m,x\right)$ will lie within the interval $\left[0,1\right]$. To overcome this difficulty, we may assume that \[ F_{Y|D,M,X}\left(a|d,m,x\right)=G_{a}\left(v_{a}\left(d,m,x\right)\right), \] where $v_{a}:\left(D,M,X\right)\mapsto\mathbb{R}$ and $G_{a}\left(.\right)$ is a link function which is non-decreasing and satisfies $G_{a}:\mathbb{R}\mapsto\left[0,1\right]$, $G_{a}\left(y\right)\rightarrow0$ if $y\rightarrow-\infty$ and $G_{r}\left(y\right)\rightarrow1$ if $y\rightarrow\infty$. The choices of $v_{a}\left(D,M,X\right)$ and the link function $G_{a}\left(.\right)$ are flexible. For example, $v_{a}\left(D,M,X\right)$ might be a neural network (see Section 3.1) or a transformation of $\left(D,M,X\right)$, which can vary with $a$. Depending on whether $Y$ is continuous or discrete, the link function $G_{a}\left(.\right)$ may be the logit, probit, linear, log-log, Gosset, the Cox proportional hazard function or the incomplete Gamma function. As pointed out by CFM_2013, for any given link function $G_{a}\left(.\right)$, we can approximate $F_{Y|D,M,X}\left(a|d,m,x\right)$ arbitrarily well if $v_{a}\left(D,M,X\right)$ is sufficiently flexible.

There are various ways to implement DR, and a popular choice is maximum likelihood estimation. Let $\hat{v}_{a}\left(D,M,X\right)$ denote the maximum likelihood estimator of $v_{a}\left(D,M,X\right)$, which is obtained by

where $\mathcal{V}_{a}$ is the parameter space of $\hat{v}_{a}\left(D,M,X\right)$. Then, $G\left(\hat{v}_{a}\left(D,M,X\right)\right)$ is an estimate of $F_{Y|D,M,X}\left(a|d,m,x\right)$. $v_{a}$ may be estimated by machine learning methods when the dimension of $X$ is large and/or the functional form of $v_{a}$ is complex.

Summarizing the Estimation Approach

Our estimation approach can be summarized as follows:

Multiplier Bootstrap

We propose a multiplier bootstrap for statistical inference. Let $\{\xi_i\}_{i=1}^n$ be a sequence of i.i.d.\ (pseudo) random variables, independent of the sample path $\{(Y_{a,i},M_i,D_i,X_i)\}_{i=1}^n$, with $E[\xi_i]=0$, $Var(\xi_i)=1$ and $E\left[\exp\left(\left|\xi_i\right|\right)\right]<\infty$. Let $\hat{v}_{k,a}$ denote a vector containing the K-fold cross-fitting estimates of the nuisance parameters, whose model parameters are estimated based on observations in the complement set, $W_{a,i}$ with ${i\in I^{c}_{k}}$. The proposed multiplier bootstrap estimator for $\hat{\theta}_{d,d^{\prime},a}$ in equation ((ref)) is given by:

The multiplier bootstrap estimator does not require re-estimating the nuisance parameters and re-calculating the causal parameters of interest in each bootstrap sample. This is particularly useful in our context, since our estimation approach is to be repeatedly applied across grid points $a$. After obtaining $\hat{\theta}_{d,d^{\prime},a}^{*}$ for all grid points, we use the procedures introduced in Section 2.1 to construct the bootstrap estimate of the c.d.f.\ of $Y(d,M(d^{\prime}))$ at these grid points and the bootstrap estimate of the $\tau$th quantile $Q_{Y(d,M(d^{\prime}))}^{*}(\tau)$. Section 3 establishes the uniform validity of the proposed multiplier bootstrap procedure.

Theoretical Results

Notation

In this section, we show that the proposed K-fold cross-fitting estimator in Section 2.2 is uniformly valid under certain conditions. We focus on establishing the theoretical properties of the estimator in equation ((ref)) and note that the assumptions and procedures required for demonstrating uniform validity of estimation based on equation ((ref)) are similar and omitted for this reason. To ease notation in our analysis, let $g_{1d}^{0}\left(X\right):=F_{D|X}^{0}\left(d|X\right)$, $g_{2d}^{0}\left(M,X\right):=F_{D|M,X}^{0}\left(d|M,X\right)$, $g_{3a}^{0}\left(D,M,X\right):=F_{Y|D,M,X}^{0}\left(a|D,M,X\right)$ and $g_{4ad}^{0}\left(D,X\right):=E\left[F_{Y|D,M,X}^{0}\left(a|d,M,X\right)|D,X\right]$ denote the nuisance parameters. Let $v_{a}^{0}$ denote the vector containing these true nuisance parameters and $\mathcal{G}_{a}$ be the set of all $v_{a}^{0}$. Let $F_{Y\left(d,M\left(d^{\prime}\right)\right)}^{0}\left(a\right)$ denote the true c.d.f.\ of the potential outcome $Y(d,M(d^{\prime}))$. The EIF of $F_{Y\left(d,M\left(d^{\prime}\right)\right)}^{0}\left(a\right)$ is $\psi_{d,d^{\prime},a}^{\theta}\left(W_{a};v_{a}^{0}\right)=\psi_{d,d^{\prime},a}\left(W_{a};v_{a}^{0}\right)-\theta$, where

Under Assumptions 1.1 to 1.4, it can be shown that

for all $a\in \mathcal{A}$ and $\left(d,d^{\prime}\right)\in\left\{ 0,1\right\} ^{2}$ (see proof of Theorem 1).

In the subsequent theoretical analysis, the expectation $E\left[.\right]$ is operated under the probability $P\in\mathcal{P}_{n}$. Let $N=n/K$ be the size in a fold or subsample, where $K$ is a fixed number. Let $E_{n}:=n^{-1}\sum_{i=1}^{n}\varsigma_{W_{i}}$ and $E_{N,k}:=N^{-1}\sum_{i\in I_{k}}\varsigma_{W_{i}}$ where $\varsigma_{w}$ is a probability distribution degenerating at $w$ and $I_{k}$ is a set of indices of observations in the $k$th subsample. Let $Z\rightsquigarrow Z^{\prime}$ denote a random variable $Z$ that weakly converges to a random variable $Z^{\prime}$. Let $\left\Vert x\right\Vert $ denote the $l^{1}$ norm and $\left\Vert x\right\Vert _{q}$ denote the $l^{q}$ norm, $q\geq2$ for a deterministic vector $x$. Let $\left\Vert X\right\Vert _{P,q}$ denote $\left(E[\left\Vert X\right\Vert ^{q}]\right)^{1/q}$ for a random vector $X$. The function $\psi_{d,d^{\prime},a}$ for identifying the parameter of interest and constructing the estimator is such that $\psi_{d,d^{\prime},a}\left(w,t\right):\mathcal{W}_{a}\times\mathcal{V}_{a}\longmapsto\mathbb{R}$, where $\left(d,d^{\prime}\right)\in\left\{ 0,1\right\} ^{2}$, $a\in\mathcal{A}\subset\mathbb{R}$, $\mathcal{W}_{a}\subset\mathbb{R}^{d_{w}}$ is a $d_{w}$ dimensional Borel set and $\mathcal{V}_{a}$ is a $G$ dimensional set of Borel measurable maps. Let $\boldsymbol{\psi}_{a}=\left(\psi_{1,1,a},\psi_{1,0,a}\psi_{0,1,a},\psi_{0,0,a}\right)$ and $\boldsymbol{\psi}_{a}:\mathcal{W}_{a}\times\mathcal{V}_{a}\longmapsto\mathbb{R}^{4}$. Let $v_{a,g}^{0}:\mathcal{U}_{a}\longmapsto\mathbb{R}$ denote the $g$th true nuisance parameter, where $\mathcal{U}_{a}\subseteq\mathcal{W}_{a}$ is a Borel set, and $v_{a}^{0}:=\left(v_{a,1}^{0},\ldots,v_{a,G}^{0}\right)\in\mathcal{V}_{a}$ denote the vector of these true nuisance parameters. Let $\hat{v}_{k,a,g}:\mathcal{U}_{a}\longmapsto\mathbb{R}$ denote an estimate of $v_{a,g}^{0}$, which is obtained by using the K-fold cross-fitting, such that the model parameters of the nuisance parameters are estimated based on observations in the complement set, $W_{a,j}$ with $j\in I_{k}^{c}$. Let $\hat{v}_{k,a}:=\left(\hat{v}_{k,a,1},\ldots,\hat{v}_{k,a,G}\right)$ denote the vector of these estimates. $v_{a}^{0}$ and $\hat{v}_{k,a}$ are both functions of $U_{a}\in\mathcal{U}_{a}$, a subvector of $W_{a}\in\mathcal{W}_{a}$. But to ease notation, we will write $v_{a}^{0}$ and $\hat{v}_{k,a}$ instead $v_{a}^{0}\left(U_{a}\right)$ and $\hat{v}_{k,a}\left(U_{a}\right)$.

$\boldsymbol{\psi}_{a}\left(W_{a},v\right)$ denotes $\boldsymbol{\psi}_{a}$ with elements $\psi_{d,d^{\prime},a}\left(W_{a};v\right)$, $\left(d,d^{\prime}\right)\in\left\{ 0,1\right\} ^{2}$. The parameter of interest is $F_{Y\left(d,M\left(d^{\prime}\right)\right)}^{0}\left(a\right)$, which can be identified by equation ((ref)), the expectation of $\psi_{d,d^{\prime},a}\left(W_{a};v\right)$ evaluated at the true nuisance parameters $v_{a}^{0}$. Let $\theta_{d,d^{\prime},a}^{0}:=E\left[\psi_{d,d^{\prime},a}\left(W_{a},v_{a}^{0}\right)\right]$. The proposed estimator of $\theta_{d,d^{\prime},a}^{0}$ is the K-fold cross-fitting estimator $\hat{\theta}_{d,d^{\prime},a}$ of ((ref)), in which $\hat{\theta}_{d,d^{\prime},a}^{\left(k\right)}=N^{-1}\sum_{i\in I_{k}}\psi_{d,d^{\prime},a}\left(W_{a,i};\hat{v}_{k,a}\right).$ In our case, the vector $\hat{v}_{k,a}$ contains estimates of the four true nuisance parameters $g_{1d}^{0}\left(X\right),g_{2d}^{0}\left(M,X\right),g_{3a}^{0}\left(D,M,X\right)$ and $g_{4ad}^{0}\left(D,X\right)$ when estimating the model parameters based on observations in the complement set, $W_{a,i}$ with $i\in I_{k}^{c}$. Let $\boldsymbol{\theta}_{a}^{0}$, $\hat{\boldsymbol{\theta}}_{a}$, $\hat{\boldsymbol{\theta}}_{a}^{\left(k\right)}$ and $\mathbf{F}^{0}(a)$ denote vectors containing $\theta_{d,d^{\prime},a}^{0}$, $\hat{\theta}_{d,d^{\prime},a}$, $\hat{\theta}_{d,d^{\prime},a}^{(k)}$ and $F_{Y\left(d,M\left(d^{\prime}\right)\right)}^{0}\left(a\right)$ for different $\left(d,d^{\prime}\right)\in\{0,1\}^{2}$. We note that $\hat{\boldsymbol{\theta}}_{a}=K^{-1}\sum_{k=1}^{K}\hat{\boldsymbol{\theta}}_{a}^{\left(k\right)}$ and $\boldsymbol{\theta}_{a}^{0}=E\left[\boldsymbol{\psi}_{a}\left(W_{a};v_{a}^{0}\right)\right]$, and if equation ((ref)) holds for all $a\in\mathcal{A}$ and $\left(d,d^{\prime}\right)\in \{0,1\}^{2}$, then $\mathbf{F}^{0}(a)=\boldsymbol{\theta}_{a}^{0}$.

Main Results

To establish the uniform validity of $\hat{\boldsymbol{\theta}}_{a}$ when estimating $\mathbf{F}^{0}(a)$, we impose the following conditions.

Note that in Assumption 2.4, $\hat{v}_{a}$ is by definition constructed based on observations in the complement set $\left(W_{a,i}\right)_{i\in I^{c}}$: $\hat{v}_{a}=\hat{v}_{a}\left(\left(W_{a,i}\right)_{i\in I^{c}}\right)$. When the random subset $I=I_{k}$, then $\hat{v}_{a}=\hat{v}_{k,a}$. Let $\mathbb{G}_{n}$ denote an empirical process $\mathbb{G}_{n}f\left(W\right)=\sqrt{n}\left(E_{n}f\left(W\right)-E\left[f\left(W\right)\right]\right)$, where $f$ is any $P\in\mathcal{P}_{n}$ integrable function on the set $\mathcal{W}$. Let $\mathbb{G}_{P}f\left(W\right)$ denote the limiting process of $\mathbb{G}_{n}f\left(W\right)$, which is a Gaussian process with zero mean and a finite covariance matrix $E\left[\left(f\left(W\right)-E\left[f\left(W\right)\right]\right)\left(f\left(W\right)-E\left[f\left(W\right)\right]\right)^{\top}\right]$ under probability $P$ (the $P$-Brownian bridge). Based on our notation and assumptions, we obtain the following result concerning the asymptotic behaviour of our estimator.

Next, we establish the uniform validity of the multiplier bootstrap under Assumptions 1 and 2. As in Assumption 2.4, let $\hat{v}_{a}$ denote the cross-fitting estimator of the nuisance parameters, whose model parameters are estimated based on observations in the complement set, $W_{a,i}$ with ${i\in I^{c}_{k}}$. Recall the multiplier bootstrap estimator in equation ((ref)) and the definition of the random variable $\xi$. By independence of $\xi$ and $W_{a}$, we have that \[ E\left[\xi\left(\psi_{d,d^{\prime},a}\left(W_{a};\hat{v}_{a}\right)-\hat{\theta}_{d,d^{\prime},a}\right)\right]=E\left[\xi\right]E\left[\psi_{d,d^{\prime},a}\left(W_{a};\hat{v}_{a}\right)-\hat{\theta}_{d,d^{\prime},a}\right]=0, \] and therefore,

Let $\hat{\boldsymbol{\theta}}_{a}^{*}$ denote a vector containing the multiplier bootstrap estimators $\hat{\theta}_{d,d^{\prime},a}^{*}$ for different $\left(d,d^{\prime}\right)\in \{0,1\}^{2}$. We write the previous result in vector form as \[ \sqrt{n}\left(\hat{\boldsymbol{\theta}}_{a}^{*}-\hat{\boldsymbol{\theta}}_{a}\right)=\mathbb{G}_{n}\xi\left(\boldsymbol{\psi}_{a}\left(W_{a};\hat{v}_{k,a}\right)-\hat{\boldsymbol{\theta}}_{a}\right), \] and let $Z_{n,P}^{*}:=\left(\mathbb{G}_{n}\xi\left(\boldsymbol{\psi}_{a}\left(W_{a};\hat{v}_{k,a}\right)-\hat{\boldsymbol{\theta}}_{a}\right)\right)_{a\in\mathcal{A}}$. Then we obtain the following result on the asymptotic behavior of the multiplier bootstrap.

Let $\phi_{\tau}\left(F_{X}\right):=\inf\left\{ a\in\mathbb{R}:F_{X}\left(a\right)\geq\tau\right\} $ be the $\tau$th quantile function of a random variable $X$ whose c.d.f. is $F_{X}$. The von Mises expansion of $\phi_{\tau}\left(F_{X}\right)$ (p.292 in Vaart_1998) is given by: \[ \phi_{\tau}\left(E_{n}\right)-\phi_{\tau}\left(E\right)=\frac{1}{\sqrt{n}}\phi_{\tau,E}^{\prime}\left(\mathbb{G}_{n}\right)+\ldots+\frac{1}{m!}\frac{1}{n^{m/2}}\phi_{\tau,E}^{\left(k\right)}\left(\mathbb{G}_{n}\right)+\ldots, \] where $\phi_{\tau,E}^{\prime}\left(.\right)$ is a linear derivative map and $\mathbb{G}_{n}$ denotes an empirical process: $\mathbb{G}_{n}f\left(W\right)=\sqrt{n}\left(E_{n}f\left(W\right)-E\left[f\left(W\right)\right]\right)$.

Let $\phi_{\boldsymbol{\theta}}^{\prime}:=\left(\phi_{\tau,\boldsymbol{\theta}}^{\prime}\right)_{\tau\in\mathcal{T}}$, where $\boldsymbol{\theta}=\left(\boldsymbol{\theta}_{a}\right)_{a\in\mathcal{A}}$. Let $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}^{0}\left(\tau\right):=\inf\left\{ a\in\mathbb{R}:F_{Y(d,M(d^{\prime}))}^{0}(a)\geq\tau\right\} $, $\hat{Q}_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(\tau\right):=\inf\left\{ a\in\mathbb{R}:\hat{\theta}_{d,d^{\prime},a}\geq\tau\right\} $ and $\hat{Q}_{Y\left(d,M\left(d^{\prime}\right)\right)}^{*}\left(\tau\right):=\inf\left\{a\in\mathbb{R}:\hat{\theta}_{d,d^{\prime},a}^{*}\geq\tau\right\}$. Let $\mathbf{Q}^{0}_{\tau}$, $\hat{\mathbf{Q}}_{\tau}$ and $\mathbf{\hat{Q}}^{*}_{\tau}$ denote the corresponding vectors containing $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}^{0}\left(\tau\right)$, $\hat{Q}_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(\tau\right)$ and $\hat{Q}_{Y\left(d,M\left(d^{\prime}\right)\right)}^{*}\left(\tau\right)$ over different $(d,d^{\prime})\in \{0,1\}^{2}$, respectively. We then obtain the following result of uniform validity for the estimation of quantiles, which can be proven by invoking the functional delta theorems (Theorems B.3 and B.4) of BCFH_2017.

Simulation

Simulation Design

This section presents a simulation study to examine the finite sample performance of the proposed DML estimators in equations ((ref)) and ((ref)). We consider the following data-generating process for the observed covariates $X=\left(X_{1},X_{2},X_{3}\right)$, where

with $V_{1}$, $V_{2}$ and $V_{3}$ being i.i.d.\ random variables following a chi-squared distribution with 1 degree of freedom. The binary treatment variable $D$ is generated based on the following model:

For the binary mediator $M$, the data-generating process is \[ M=1\left\{ -0.070+0.710D+X^{\top}\left(-0.054,-0.482,0.299\right)+\varepsilon_{M}>0\right\} . \] The model for outcome $Y$ is defined as follows:

The error terms $\left(\varepsilon_{Y},\varepsilon_{D},\varepsilon_{M}\right)$ are mutually independent standard normal random variables and also independent of $X$. Analytically computing the unconditional c.d.f. and quantiles of the potential outcome $Y\left(d,M\left(d^{\prime}\right)\right)$ is difficult for the data generating process considered. For this reason, we use a Monte Carlo simulation to approximate the true values of the c.d.f. and the quantiles. We draw 40 million observations of $\left(\varepsilon_{Y},\varepsilon_{M},V_{1},V_{2},V_{3}\right)$ from their respective true distributions, and for each observation, we calculate the corresponding potential outcomes $Y\left(d,M\left(d^{\prime}\right)\right)$ for $(d,d^{\prime})\in \{0,1\}^{2}$. In the next step, we evaluate profiles of the empirical c.d.f.'s and quantiles of the 40 million sampled potential outcomes and use the evaluated profiles as approximations to their true profiles. In Figure (ref) in the appendix, the upper panel shows the approximate true profiles of the c.d.f.'s $F_{Y\left(d,M\left(d^{\prime}\right)\right)}\left(a\right)$, while the lower panel provides the approximate true profiles of quantiles $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}(\tau)$ for $(d,d^{\prime})\in \{0,1\}^{2}$. Figure (ref) in the appendix depicts the approximate true profiles of NDQTE (NDQTE'), NIQTE (NIQTE') and TQTE across quantiles.

In our simulation design, the nuisance parameters have the following functional forms:

where $\Phi\left(.\right)$ is the c.d.f.\ of a standard normal random variable. The vector of parameters satisfies $\left(\beta_{0,a},\alpha_{0,a},\alpha_{1,a},\beta_{1,a},\boldsymbol{\beta}_{2,a}^{\top}\right)=a\times\left(\beta_{0},\alpha_{0},\alpha_{1},\beta_{1},\boldsymbol{\beta}_{2}^{\top}\right)$. When running the simulations, we also include a set of auxiliary (exogenous) variables: $X^{aug}:=\left(X_{1}^{aug},X_{2}^{aug},\ldots,X_{J}^{aug}\right)$, $X_{j}^{aug}=U_{1j}\left(V_{1}+V_{2}+V_{3}\right)+U_{2j}\left(Z_{j}\right)^{2}$, $j=1,\ldots,J$, where $U_{1j}\sim i.i.d.U\left(0,0.2\right)$, $U_{2j}\sim i.i.d.U\left(0.8,1\right)$. $\mathbf{Z}:=\left(Z_{1},Z_{2},\ldots,Z_{J}\right)$ follow a multivariate normal distribution with a mean vector $\mathbf{0}$ and a covariance matrix with elements $0.5^{\left|j-l\right|},$ $j,l=1,\ldots,J$. $V_{1},V_{2},V_{3},U_{1j},U_{2j}$, $\mathbf{Z}$ and the error terms $\left(\varepsilon_{Y},\varepsilon_{D},\varepsilon_{M}\right)$ are mutually independent. Depending on the realized values of $U_{1j}$ and $U_{2j}$, the correlations $cor\left(X_{1},X_{j}^{a}\right)$, $cor\left(X_{2},X_{j}^{a}\right)$, $cor\left(X_{3},X_{j}^{a}\right)$ vary, and on average they amount to 0.139, 0.147 and 0.135, respectively.

The Post Lasso Estimator

Let $W_{a,i}^{aug}=\left(Y_{i},D_{i},M_{i},X_{i},X_{i}^{aug}\right)$ denote the $i$th observation of the simulated data. When applying K-fold cross-fitting to $W_{a,i}^{aug}$, we estimate the models of the nuisance parameters based on post-lasso regression: we first estimate the models by lasso regression and then re-estimate the models without (lasso) penalization when including only those regressors with non-zero coefficients in the respective previous lasso steps. We denote the lasso estimator of the coefficients of $F_{Y|D,M,X}\left(a|D,M,X\right)$ in K-fold cross-fitting by

where $p$ is the number of covariates, $\left|I_{k}^{c}\right|$ is the number of observations in the complement set $I_{k}^{c}$, $\gamma$ is the penalty parameter, $\left\Vert .\right\Vert _{1}$ denotes the $l_{1}$ norm and $\hat{\varPsi}$ is a diagonal matrix of penalty loadings. Here, the loss function $L\left(.\right)$ corresponds to that in equation ((ref)) and the link function $G_{a}\left(.\right)$ is $\Phi\left(.\right)$. When solving the lasso estimation problem of expression ((ref)), we only impose a penalty on $\left(X,X^{aug}\right)$ and therefore, the first four diagonal elements (for the intercept term, $D$, $M$ and $DM$) of $\hat{\varPsi}$ are ones, while the remaining diagonal elements are zeros. The value of the penalty parameter is determined based on the procedure outlined in BCFH_2017. The other nuisance parameters $f_{D|X}\left(d|X\right)$ and $f_{M|D,X}\left(m|D,X\right)$ are estimated in an analogous way and we denote by $\hat{\boldsymbol{\gamma}}_{D}$ and $\hat{\boldsymbol{\gamma}}_{M}$ their corresponding lasso estimators. Let $\tilde{\Xi}$ denote the union of variables in $\left(X,X^{aug}\right)$ with non-zero lasso coefficient estimates in one or several lasso regressions of the three nuisance parameters, with $\tilde{\Xi}\subseteq\text{supp}\left(\hat{\boldsymbol{\gamma}}_{Y,a}\right)\cup\text{supp}\left(\hat{\boldsymbol{\gamma}}_{D}\right)\cup\text{supp}\left(\hat{\boldsymbol{\gamma}}_{M}\right)$. The post-lasso estimator of $F_{Y|D,M,X}\left(a|D,M,X\right)$ is defined as

The post-lasso estimators of $f_{D|X}\left(d|X\right)$ and $f_{M|D,X}\left(m|D,X\right)$ are obtained analogously. Based on the post-lasso estimates of the coefficients, we estimate the nuisance parameters among observations $W_{a,i}^{aug}$, $i\in I_{k}$, and use them to compute estimators ((ref)) or ((ref)).

When using the estimator based on equation ((ref)), we approximate the nuisance parameters $f_{D|M,X}\left(d|M,X\right)$ by a probit model $\Phi\left(\lambda_{2}+\lambda_{3}M+X^{\top}\boldsymbol{\lambda}_{4}\right)$. The post-lasso approach for estimating $f_{D|M,X}\left(d|M,X\right)$ is the same as before. To estimate $E\left[F_{Y|D,M,X}\left(a|D,M,X\right)|d^{\prime},X\right]$, we approximate $E\left[F_{Y|D,M,X}\left(a|D,M,X\right)|D,X\right]$ by a linear model $\beta_{3}+\beta_{4}D+X^{\top}\boldsymbol{\beta}_{5}$. We calculate post-lasso estimates of $F_{Y|D,M,X}\left(a|D,M,X\right)$ among observations in the complement set, $W_{a,i}^{aug}$ with $i\in I_{k}^{c}$, and estimate $\left(\beta_{3},\beta_{4},\boldsymbol{\beta}_{5}\right)$ by linearly regressing these estimates on $D$ and those covariates previously selected for computing the post-lasso estimate of $F_{Y|D,M,X}\left(a|D,M,X\right)$. We then use the linear regression coefficients coming from the complement set to make cross-fitted predictions among observations with indices $i\in I_{k}$ and $D_{i}=d^{\prime}$, which serve as estimates of $E\left[F_{Y|D,M,X}\left(a|D,M,X\right)|d^{\prime},X\right]$.

Simulation Results

To evaluate the performance of the proposed DML estimators of the c.d.f.'s of the potential outcomes across grid values $a$, we calculate integrated mean squared error (IMSE) and integrated Anderson--Darling weighted MSE (IWMSE) for each simulation:

To assess the performance of the estimators of the quantiles of the potential outcomes, $Q_{Y\left(d,M\left(d^{\prime}\right)\right)}$, we compute the integrated absolute error (IAE) across ranks $\tau$:

where $\mathcal{T}$ is the grid of ranks, which we set to $(0.05,0.06,\ldots,0.95)$. Furthermore, we calculate the IAE for the estimators of the quantile treatment effects defined in equations ((ref)) to ((ref)). In the simulations, we set $K=3$ for 3-fold cross-fitting. The number of auxiliary variables is $J=250$ and we consider sample sizes 2,500, 5,000 and 10,000 observations in the simulations. The reported performance measures are averages for $\left(d,d^{\prime}\right)\in\{0,1\}^{2}$ over 1,000 simulations.

Table (ref) reports the results for the IMSE and IWMSE (scaled by 1,000), Table (ref) those for the IAE. All the performance measures behave rather favorably. As the sample size increases, the performance measures (and thus, estimation errors) decline sharply. However, for different combinations of $(d,d^{\prime})$, the levels of the performance measures are different, especially when the sample size is small. Estimation errors are significantly larger if $(d,d^{\prime}) = (0,1)$ and $(0,0)$, rather than $(d,d^{\prime})=(1,1)$ and $(1,0)$. This is also reflected by the performance measures of NDQTE (NDQTE') and TQTE, which point to higher errors than those of NIQTE (NIQTE'). When comparing the performance measures of the two estimators based on equations ((ref)) and ((ref)), we find some differences in their levels when the sample size is small. However, the differences vanish as the sample size increases, which suggests that the two estimators perform equally well asymptotically in the simulation design considered.

Empirical Application

The Job Corps Data

We apply the proposed estimators of natural direct and indirect quantile treatment effects to data from the National Job Corps Study, in order to evaluate the impact of the Job Corps (JC) training program on earnings of young individuals with disadvantaged backgrounds. JC is the largest and most comprehensive job training program for disadvantaged youth in the US. It provides participants with vocational training and/or classroom education, housing, and board over an average duration of 8 months. Participants also receive health education as well as health and dental care. SchBuGl01 and SchBuMc08 assess the average effects of random assignment to JC on several labor market outcomes and find it to increase education, employment, and earnings in the longer run. Other contributions evaluate more specific aspects or components of JC, like the average effect of the time spent in training or of particular training sequences on employment and earnings, see e.g.\ FlGoNe12 and BoHuLa2020.

Furthermore, several studies conduct mediation analyses to assess the average direct and indirect effects of program participation. FlFl09 and Hu14 consider work experience or employment as mediators, respectively, and find positive direct effects of JC on earnings and general health, respectively, when invoking a selection-on-observables assumption. FlFl10 avoid the latter assumption based on a partial identification approach based on which they compute upper and lower bounds for the 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, they find a positive direct effect on labor market outcomes, net of the indirect mechanism via obtaining a degree. FrHu17 base their mediation analysis on separate instrumental variables for the treatment and the mediator and find a positive indirect effect of training on earnings through an increase in the number of hours worked. We contribute to the causal mediation literature on the effectiveness of the JC program by considering quantile treatment effects across different ranks of the potential outcome distributions, which provides more insights on effect heterogeneity than the evaluation of average effects.

For our empirical analysis, we consider the JC data provided in the causalweight package by BH_2022 for the statistical software R, which is a processed data set with 9,240 observations that contains a subset of the variables available in the original National Job Corps Study. Our outcome of interest is weekly earnings in the third year after the assignment (the variable earny3 in the JC data frame), while the treatment is a binary indicator for participation in any (classroom-based or vocational) training in the first year after program assignment (trainy1). We aim at assessing whether training directly affects the earnings outcome, and whether it also has an indirect effect by affecting health. For this reason, we consider general health one year after program assignment (health12) as mediator, a categorical variable ranging from 1 (excellent health) to 4 (poor health). The motivation is that participation in training aimed at increasing human capital and labor market perspectives may have an impact on mental health, which in turn may affect labor market success. Furthermore, JC might also affect physical health through health education and health/dental care, which can influence labor market success, too. For this reason, we aim at disentangling the direct earnings effect of training and its indirect effect operating via health.

The data set also contains 28 pre-treatment covariates, which include socio-economic information such as a study participant's gender, age, ethnicity, (own) education and parents' education, mother tongue, marital status, household size, previous employment, earnings and welfare receipt, health status, smoking behavior, alcohol consumption, and whether a study participant has at least one child. We assume that sequential ignorability of the treatment and the mediator holds conditional on these observed characteristics, implying that the permit controlling for any factors jointly affecting training participation and the earnings outcome, training participation and health 12 months after assignment, or health and earnings. To make lasso-based estimation of the nuisance parameters in our DML approach more flexible, we create interaction terms between all of the 28 covariates and squared terms for any non-binary covariates. This entails a total of 412 control variables that include both the original covariates and the higher order/interaction terms which we include in our DML approach. Table (ref) provides summary statistics for the outcome, the treatment, the mediator and the covariates.

Effect Estimates

Before considering quantile treatment effects, we first estimate the average direct and indirect effects by a K-fold cross-fitting estimator based on Theorem 2 in FHLLS_2022, as implemented in the causalweight package for R. Table (ref) reports the estimated average total effect (TE) of training, the average natural direct effects (NDE and NDE') and the average natural indirect effects (NIE and NIE') operating via general health. The TE estimate (Effect) suggests that participation in JC increases average weekly earnings in the third year by roughly 16 to 17 USD. As the estimated mean potential outcome under non-treatment amounts to approximately 161 USD, the program increases weekly earnings by roughly 10% according to our estimate. The TE is highly statistically significant as the standard error (Sdt.err) of 3.740 is rather low relative to the effect estimate, such that p-value that is close to zero.

The total effect seems to be predominantly driven by the direct impact of training on earnings, as both NDE and NDE' are of similar magnitude as TE and highly statistically significant. In contrast, the indirect effect under non-treatment ($d=0$), NIE', is close to zero and insignificant, while that under treatment ($d=1$), NIE, amounts to -0.403 USD and is statistically significant at the 5% level. Bearing in mind that the health mediator is inversely coded (a smaller value implies better health), this negative estimate suggests a positive average indirect effect of training participation on earnings under treatment, which is, however, rather modest. Furthermore, the effect heterogeneity across NIE and NIE' points to moderate interaction effects of the treatment and the mediator: the impact of health on earnings appears to be somewhat more important under training than without training.

The average effects might mask interesting effect heterogeneity across ranks of the earnings distribution. For this reason, we estimate the total quantile treatment effect (TQTE), natural direct quantile treatment effects (NDQTE and NDQTE') and natural indirect quantile treatment effects (NIQTE and NIQTE') across ranks ($\tau$) 0.2 to 0.9. To this end, we invert our K-fold cross-fitting estimator $\hat{\theta}_{d,d^{\prime},a}$ of equation ((ref)) and estimate the nuisance parameters by post-lasso regression as outlined in Section 4.2. Figures (ref) and (ref) depict the estimates of the causal effects (on the y-axis) across $\tau$ (on the x-axis), which correspond to the solid lines in the respective graphs. The dashed lines provide the 95% confidence intervals based on the multiplier bootstrap introduced in Section 2.5.

The quantile treatment effects are by and large in line with the average treatment effects. TQTE, NDQTE and NDQTE' are statistically significantly positive at the 5% across almost all ranks $\tau$ considered and generally quite similar to each other. In contrast, all of the NIQTE estimates (the indirect effects under $d=1$) are relatively close to zero and statistically insignificant. The majority of the NIQTE' estimates (the indirect effects under $d=0$) are not statistically significantly different from zero either. However, several of the negative effects measured at lower ranks (roughly between the 0.2th and 0.4th quantiles) are marginally statistically significant and point to an earnings-increasing indirect effect under non-treatment (due to inverse coding of the health mediator). This potentially interesting pattern is averaged out when considering NIE' (the average indirect effect under $d=0$), which we found to be virtually zero and insignificant, see Table (ref). Finally, the non-monotonic shape of the point estimates of TQTE, NDQTE and NDQTE' across ranks $\tau$ suggests heterogeneous effects at different quantiles of the potential earnings distributions. At the same time, the width of the confidence intervals suggests that the null hypothesis of homogeneous effects cannot be rejected for most of the quantiles considered.

Conclusion

We proposed a DML approach for estimating natural direct and indirect quantile treatment effects under a sequential ignorability assumption. The method relies on the efficient score functions of the potential outcomes' cumulative distributional functions, which are inverted to compute the quantiles as well as the treatment effects (as the differences in potential outcomes at those quantiles). The robustness property of the efficient score functions permits estimating the nuisance parameters (outcome, treatment, and mediator models) by machine learning and cross-fitting avoids overfitting bias. We demonstrated that our quantile treatment effect estimators are root-n-consistent and asymptotically normal. Furthermore, we suggested a multiplier bootstrap and demonstrated its consistency for uniform statistical inference. We also investigated the finite sample performance of our estimators by means of a simulation study. Finally, we applied our method to data from the National Job Corp Study to evaluate the direct earnings effects of training across the earnings distribution, as well as the indirect effects operating via general health. We found positive and statistically significant direct effects across a large range of the earnings quantiles, while the indirect effects were generally close to zero and mostly statistically insignificant.

{ \setcounter{equation}{0}