Causal mediation analysis with double machine learning

\doublespacing \pagestyle{plain}

Abstract: { This paper combines causal mediation analysis with double machine learning to control for observed confounders in a data-driven way under a selection-on-observables assumption in a high-dimensional setting. We consider the average indirect effect of a binary treatment operating through an intermediate variable (or mediator) on the causal path between the treatment and the outcome, as well as the unmediated direct effect. Estimation is based on efficient score functions, which possess a multiple robustness property w.r.t.\ misspecifications of the outcome, mediator, and treatment models. This property is key for selecting these models by double machine learning, which is combined with data splitting to prevent overfitting in the estimation of the effects of interest. We demonstrate that the direct and indirect effect estimators are asymptotically normal and root-$n$ consistent under specific regularity conditions and investigate the finite sample properties of the suggested methods in a simulation study when considering lasso as machine learner. We also provide an empirical application to the U.S.\ National Longitudinal Survey of Youth, assessing the indirect effect of health insurance coverage on general health operating via routine checkups as mediator, as well as the direct effect. We find a moderate short-term effect of health insurance coverage on general health which is, however, not mediated by routine checkups.}

{ Keywords: mediation, direct and indirect effects, causal mechanisms, double machine learning, efficient score.}

{ JEL classification: C21. \quad }

{ {\tiny Addresses for correspondence: Helmut Farbmacher, Max Planck Society, Munich Center for the Economics of Aging, Amalienstr. 33, 80799 Munich, Germany; [email removed]. Martin Huber, University of Fribourg, Department of Economics, Bd.\ de P\'{e}rolles 90, 1700 Fribourg, Switzerland; [email removed]. Luk\'{a}\v{s} Laff\'{e}rs, Matej Bel University, Department of Mathematics, Tajovskeho 40, 97411 Bansk\'{a} Bystrica, Slovakia; [email removed]. Henrika Langen, University of Fribourg, Department of Economics, Bd.\ de P\'{e}rolles 90, 1700 Fribourg, Switzerland; [email removed]. Martin Spindler, University of Hamburg, Faculty of Business Administration, Moorweidenstr.\ 18, 20148 Hamburg, Germany; [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

Causal mediation analysis aims at decomposing the causal effect of a treatment on an outcome of interest into an indirect effect operating through a mediator (or intermediate outcome) and a direct effect comprising any causal mechanisms not operating through that mediator. Even if the treatment is random, direct and indirect effects are generally not identified by naively controlling for the mediator without accounting for its likely endogeneity, see RoGr92. While much of the earlier literature either neglected endogeneity issues or relied on restrictive linear models, see for instance Co57, JuKe81, and BaKe86, more recent contributions consider more general identification approaches using the potential outcome framework. Some of the numerous examples are RoGr92, Pearl01, Robins2003, PeSiva06, VanderWeele09, ImKeYa10, Hong10, AlbertNelson2011, ImYa2011, TchetgenTchetgenShpitser2011, VansteelandtBekaertLange2012, and Huber2012. Using the denomination of Pearl01, the literature distinguishes between natural direct and indirect effects, where mediators are set to their potential values `naturally' occurring under a specific treatment assignment, and the controlled direct effect, where the mediator is set to a `prescribed' value.

The vast majority of identification strategies relies on selection-on-observable-type assumptions implying that the treatment and the mediator are conditionally exogenous when controlling for observed covariates. Empirical examples in economics and policy evaluation include FlFl09, HeckmanPintoSavelyev2013, Keeleetal2015, ContiHeckmanPinto2016, Huber2015, HuberLechnerMellace2017, BellaniBia2018, BijwaardJones2018, and HuberLechnerStrittmatter2018. Such studies typically rely on the (implicit) assumption that the covariates to be controlled for can be unambiguously preselected by the researcher, for instance based on institutional knowledge or theoretical considerations. This assumes away uncertainty related to model selection w.r.t.\ covariates to be included and entails incorrect inference under the common practice of choosing and refining the choice of covariates based on their predictive power.

To improve upon this practice, this paper combines causal mediation analysis based on efficient score functions, see TchetgenTchetgenShpitser2011, with double machine learning as outlined in Chetal2018 for a data-driven control of observed confounders to obtain valid inference under specific regularity conditions. In particular, one important condition is that the number of important confounders (that make the selection-on-observables assumptions to hold approximately) is not too large relative to the sample size. However, the set of these important confounders need not be known a priori and the set of potential confounders can be even larger than the sample size.\footnote{Different from conventional semiparametric methods, the double machine learning framework does not require the set of potential confounders to be restricted by Donsker conditions, but permits the set to be unbounded and to grow with the sample size.} This is particularly useful in high dimensional data with a vast number of covariates that could potentially serve as control variables, which can render researcher-based covariate selection complicated if not infeasible. We demonstrate root-$n$ consistency and asymptotic normality of the proposed effect estimators under specific regularity conditions by verifying that the general framework of Chetal2018 for well-behaved double machine learning is satisfied in our context.

TchetgenTchetgenShpitser2011 suggest estimating natural direct and indirect effects based on the efficient score functions of the potential outcomes, which requires plug-in estimates for the conditional mean outcome, mediator density, and treatment probability. Analogous to doubly robust estimation of average treatment effects, see Robins+94 and RoRo95, the resulting estimators are semiparametrically efficient if all models of the plug-in estimates are correctly specified and remain consistent even if one model is misspecified. Our first contribution is to show that the efficient score function of TchetgenTchetgenShpitser2011 satisfies the so-called Neyman1959 orthogonality discussed in Chetal2018, which makes the estimation of direct and indirect effects rather insensitive to (local) estimation errors in the plug-in estimates. Second, we show that by an application of Bayes' Theorem, the score function of TchetgenTchetgenShpitser2011 can be transformed in a way that avoids estimation of the conditional mediator density and show it to be Neyman orthogonal. This appears particularly useful when the mediator is a vector of variables and/or continuous. Third, we establish the score function required for estimating the controlled direct effect along with Neyman orthgonality.

Neyman orthgonality is key for the fruitful application of double machine learning, ensuring robustness in the estimation of the nuisance parameters which is crucial when applying modern machine learning methods. Random sample splitting -- to estimate the parameters of the plug-in models in one part of the data, while predicting the score function and estimating the direct and indirect effects in the other part -- avoids overfitting the plug-in models (e.g.\ by controlling for too many covariates). It increases the variance by only using part of the data for effect estimation. This is avoided by cross-fitting which consists of swapping the roles of the data parts for estimating the plug-in models and the treatment effects to ultimately average over the effects estimates in either part. When combining efficient score-based effect estimation with sample splitting, $n^{-1/2}$-convergence of treatment effect estimation can be obtained under a substantially slower convergence of $n^{-1/4}$ for the plug-in estimates, see Chetal2018. Under specific regularity conditions, this convergence rate can attained by various machine learning algorithms including lasso regression, see Tibshirani96.

We investigate the estimators' finite sample behavior based on the score function of TchetgenTchetgenShpitser2011 and the alternative score suggested in this paper when using post-lasso regression as machine learner for the plug-in estimates. Furthermore, we apply our method to data from the National Longitudinal Survey of Youth 1997 (NLSY97), where a large set of potential control variables is available. We disentangle the short-term effect of health insurance coverage on general health into an indirect effect which operates via the incidence of a routine checkup in the last year and a direct effect covering any other causal mechanisms. While we find a moderate health-improving direct effect, the indirect effect is very close to zero. We therefore do not find evidence that health insurance coverage affects general health through routine checkups in the short run.

We note that basing estimation on efficient score functions is not the only framework satisfying the previously mentioned robustness w.r.t.\ estimation errors in plug-in parameters. This property is also satisfied by the targeted maximum likelihood estimation (TMLE) framework by vanderLaanRubin2006, see the discussion in Diaz2020. TMLE relies on iteratively updating (or robustifying) an initial estimate of the parameter of interest based on regression steps that involve models for the plug-in parameters. ZhengvanderLaan2012 have developed an estimation approach for natural direct and indirect effects using TMLE, where the plug-in parameters might by estimated by machine learners, e.g.\ the super learner, an ensemble method suggested by vanderLaanetal2007. This iterative estimation approach is therefore an alternative to the double machine learning-based approach suggested in this paper, for which we demonstrate $n^{-1/2}$-consistency under specific conditions.

This paper proceeds as follows. Section (ref) introduces the concepts of direct and indirect effect identification in the potential outcome framework. In Section (ref), we present the identifying assumptions and discuss identification based on efficient score functions. Section (ref) proposes an estimation procedure based on double machine learning and shows root-$n$ consistency and asymptotic normality under specific conditions. Section (ref) provides a simulation study. Section (ref) presents an empirical application to data from the NLSY97. Section (ref) concludes.

Definition of direct and indirect effects

We aim at decomposing the average treatment effect (ATE) of a binary treatment, denoted by $D$, on an outcome of interest, $Y$, into an indirect effect operating through a discrete mediator, $M$, and a direct effect that comprises any causal mechanisms other than through $M$. We use the potential outcome framework, see for instance Rubin74, to define the direct and indirect effects of interest, see also TenHaveetal2007 and Albert2008 for further examples in the context of mediation. $M(d)$ denotes the potential mediator under treatment value $d$ $\in$ $\{0,1\}$, while $Y(d,m)$ denotes the potential outcome as a function of both the treatment and some value $m$ of the mediator $M$.\footnote{Throughout this paper, capital letters denote random variables and small letters specific values of random variables.} The observed outcome and mediator correspond to the respective potential variables associated with the actual treatment assignment, i.e.\ $Y=D\cdot Y(1,M(1)) + (1-D)\cdot Y(0,M(0))$ and $M=D\cdot M(1) + (1-D)\cdot M(0)$, implying that any other potential outcomes or mediators are a priori (i.e.\ without further statistical assumptions) unknown.

We denote the ATE by $\Delta=E[Y(1,M(1))-Y(0,M(0))]$, which comprises both direct and indirect effects. To decompose the latter, note that the average direct effect, denoted by $\theta (d)$, equals the difference in mean potential outcomes when switching the treatment while keeping the potential mediator fixed, which blocks the causal mechanism via $M$:

The (average) indirect effect, $\delta (d)$, equals the difference in mean potential outcomes when switching the potential mediator values while keeping the treatment fixed to block the direct effect.

RoGr92 and Robins2003 referred to these parameters as pure/total direct and indirect effects, FlFl09 as net and mechanism average treatment effects, and Pearl01 as natural direct and indirect effects, which is the denomination used in the remainder of this paper.

The ATE is the sum of the natural direct and indirect effects defined upon opposite treatment states $d$, which can be easily seen from adding and subtracting the counterfactual outcomes $E[Y(0,M(1))]$ and $E[Y(1,M(0))]$:

The distinction between $\theta(1)$ and $\theta(0)$ as well as $\delta(1)$ and $\delta(0)$ hints to the possibility of heterogeneous effects across treatment states $d$ due to interaction effects between $D$ and $M$. For instance, the direct effect of health insurance coverage ($D$) on general health ($Y$) might depend on whether or not a person underwent routine check-ups ($M$). We note that a different approach to dealing with the interaction effects between $D$ and $M$ is a three-way decomposition of the ATE into the pure direct effect ($\theta(0)$), the pure indirect effect ($\delta(0)$) and the mediated interaction effect, see VanderWeele2013.

The so-called controlled direct effect, denoted by $\gamma (m)$, is a further parameter that received much attention in the mediation literature. It corresponds to the difference in mean potential outcomes when switching the treatment and fixing the mediator at some value $m$:

In contrast to $\theta (d)$, which is conditional on the potential mediator value `naturally' realized for treatment $d$ which may differ across subjects, $\gamma (m)$ is conditional on enforcing the same mediator state in the entire population. The two parameters are only equivalent in the absence of an interaction between $D$ and $M$. Whether the natural or controlled direct effect is more relevant depends on the feasibility and desirability to intervene on or prescribe the mediator, see Pearl01 for a discussion of the `descriptive' and `prescriptive' natures of natural and controlled effects. There is no indirect effect parameter matching the controlled direct effect, implying that the difference between the total effect and the controlled direct effect does in general not correspond to the indirect effect, unless there is no interaction between $D$ and $M$, see e.g.\ Kaufmanetal2004.

Assumptions and identification

Our identification strategy is based on the assumption that confounding of the treatment-outcome, treatment-mediator, and mediator-outcome relations can be controlled for by conditioning on observed covariates, denoted by $X$. The latter must not contain variables that are influenced by the treatment, such that $X$ is typically evaluated prior to treatment assignment. Figure (ref) provides a graphical illustration using a directed acyclic graph, with arrows representing causal effects. Each of $D$, $M$, 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 two or all three elements $(D,M,Y)$ conditional on $X$.

Formally, the first assumption invokes conditional independence of the treatment and potential mediators or outcomes given $X$. This restriction has been referred to as conditional independence, selection on observables, or exogeneity in the treatment evaluation literature, see e.g.\ Imbens03. This rules out confounders jointly affecting the treatment on the one hand and the mediator and/or the outcome on the other hand conditional on $X$. In non-experimental data, the plausibility of this assumption critically hinges on the richness of $X$. \newline Assumption 1 (conditional independence of the treatment):\newline $\{Y(d',m), M(d)\} \bot D | X$ for all $d',d \in \{0,1\}$ and $m$ in the support of $M$, \newline where `$\bot$' denotes statistical independence. The second assumption requires the mediator to be conditionally independent of the potential outcomes given the treatment and the covariates. \newline Assumption 2 (conditional independence of the mediator):\newline $ Y(d',m) \bot M | D=d, X=x $ for all $d',d \in \{0,1\}$ and $m,x$ in the support of $M,X$.\newline Assumption 2 rules out confounders jointly affecting the mediator and the outcome conditional on $D$ and $X$. If $X$ is pre-treatment (as is common to avoid controlling for variables potentially affected by the treatment), this implies the absence of post-treatment confounders of the mediator-outcome relation. Such a restriction needs to be rigorously scrutinized and appears for instance less plausible if the time window between the measurement of the treatment and the mediator is large in a world of time-varying variables.

The third assumption imposes common support on the conditional treatment probability across treatment states.\newline Assumption 3 (common support):\newline $\Pr(D=d| M=m, X=x)>0$ for all $d \in \{0,1\}$ and $m,x$ in the support of $M,X$. \newline The common support assumption, also known as positivity or covariate overlap assumption, restricts the conditional probability to be or not be treated given $M,X$, henceforth referred to as propensity score, to be larger than zero. It implies the weaker condition that $\Pr(D=d|X=x)>0$ such that the treatment must not be deterministic in $X$, otherwise no comparable units in terms of $X$ are available across treatment states. By Bayes' theorem, Assumption 3 also implies that $\Pr(M=m| D=d, X=x)>0$ if $M$ is discrete or that the conditional density of $M$ given $D,X$ is larger than zero if $M$ is continuous. Conditional on $X$, the mediator state must not be deterministic in the treatment, otherwise no comparable units in terms of the treatment are available across mediator states. Assumptions 1 to 3 are standard in the causal mediation literature, see for instance ImKeYa10, TchetgenTchetgenShpitser2011, VansteelandtBekaertLange2012, and Huber2012, or also Pearl01, PeSiva06, and Hong10, for closely related restrictions.

TchetgenTchetgenShpitser2011 discuss identification of the counterfactual $E[Y(d, M(1-d))]$ based on the efficient score function:

where $ f(M|D,X) $ denotes the conditional density of $ M $ given $ D $ and $ X $ (if $M$ is discrete, this is a conditional probability and integrals need to be replaced by sums), $ p_d(X) = \Pr(D=d|X) $ the probability of treatment $ D=d $ given $ X $, and $ \mu(D,M,X)=E(Y|D,M,X)$ the conditional expectation of outcome $ Y $ given $ D$, $ M $, and $ X $. (ref) satisfies a multiple robustness property in the sense that estimation remains consistent even if one out of the three models for the plug-in parameters $ f(M|D,X) $, $p_d(X)$, and $\mu(D,M,X)$ is misspecified.

To derive an alternative expression for identification, note that by Bayes' Law,

where $f(M|X)$ is the conditional distribution of $M$ given $X$ and $ p_d(X,M) = \Pr(D=d|X,M) $. Furthermore,

Therefore, an alternative and also multiply robust representation of (ref) is

Similarly as the approaches based on inverse probability weighting (rather than efficient scores) in Huber2012 and TchetgenTchetgen2013, (ref) avoids conditional mediator densities, which appears attractive if $M$ is continuous and/or multidimensional. On the other hand, it requires the estimation of an additional parameter, namely the nested conditional mean $E[ \mu(d,M,X) | D= 1-d, X ]$, as similarly found in Milesetal2020, who suggest a multiply robust score function for assessing path-specific effects. Alternatively to rearranging the score function by TchetgenTchetgenShpitser2011 as outlined above, ratios of conditional densities as for instance appearing in the first component of (ref) might be treated as additional nuisance parameter and estimated directly via density-ratio estimation, see e.g. Sugiyamaetal2010 for density-ratio estimation in high-dimensional settings. Such methods based on directly estimating the density ratio without going through estimating the densities in numerator and denominator separately are shown in several studies to compare favorably with estimating the densities separately, see e.g. Kanamorietal2012.

Efficient score-based identification of $E[Y(d, M(d))]$ under $ Y(d,m) \bot \{D,M\} | X=x$ (see Assumption 1) has been established in the literature on doubly robust ATE estimation, see for instance Robins+94 and Ha98:

where $ \mu(D,X)=E(Y|D,M(D),X)=E(Y|D,X)$ is the conditional expectation of outcome $ Y $ given $ D$ and $ X $.

For identifying the controlled direct effect, we now assume that $M$ is discrete (while this need not be the case in the context of natural direct and indirect effects) such that for all $m$ in the support of $M$, it must hold that $\Pr(M=m)>0$. As Assumptions 1 and 2 imply $ Y(d,m) \bot \{D,M\} | X=x$, doubly robust identification of the potential outcome $E[Y(d,m)]$, which is required for the controlled direct effect, follows from replacing $I\{D=d\}$ and $p_d(X)$ in (ref) by $I\{D=d, M=m\}= I\{M=m\}\cdot I\{D=d\}$ and $\Pr(D=d,M=m|X)=f(m|d,X)\cdot p_d(X)$:

Estimation of the counterfactual with K-fold Cross-Fitting

We subsequently propose an estimation strategy for the counterfactual $E[Y(d, M(1-d))]$ with $d \in \{0,1\}$ based on the efficient score function by TchetgenTchetgenShpitser2011 provided in (ref) and show its root-$n$ consistency under specific regularity conditions. To this end, let $\mathcal{W} = \{W_i|1\leq i \leq N\}$ with $W_i = (Y_i, M_i, D_i, X_i)$ for $i=1,\ldots, n$ 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,M,X),\hat{f}(M|D,X), \hat{p_d}(X) \}$ and the true nuisance parameters by $\eta_0 = \{\mu_0(D,M,X), f_0(M|D,X), p_{d0}(X) \}$. Finally, $\Psi_{d0}=E[Y(d, M(1-d))]$ denotes the true counterfactual.

We suggest estimating $\Psi_{d0}$ using the following algorithm that combines orthogonal score estimation with sample splitting and is root-$n$ consistent under conditions outlined further below. \newline Algorithm 1: Estimation of $E[Y(d, M(1-d))]$ based on equation (ref)

Algorithm 1 can be adapted to estimate the counterfactuals required for the controlled direct effect, see (ref). To this end, denote by $\Psi_{dm0}=E[Y(d, m)]$ the true counterfactual of interest, which is estimated by replacing $\psi_{d}$ and $\Psi_d$ by $\psi_{dm}$ and $\Psi_{dm0}$, respectively, everywhere in Algorithm 1.

In order to achieve root-$n$ consistency for counterfactual estimation, we make specific assumptions about the prediction qualities of the machine learners for our plug-in estimates of the nuisance parameters. Closely following Chetal2018, we to this end 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.$ Also, let $c, \epsilon, C, \underline{f},\overline{f}$ and $q$ be positive constants such that $q>2,$ and let $K \geq 2$ be a fixed integer. Furthermore, for any random vector $Z = (Z_1,...,Z_l)$, let $\left\| Z \right\|_{q} = \max_{1\leq j \leq l}\left\| Z_l \right\|_{q},$ where $\left \| Z_l \right\|_{q} = \left( E\left[ \left| Z_l \right|^q \right] \right)^{\frac{1}{q}}$. For the sake of easing notation, we assume that $n/K$ is an integer. For brevity, we omit the dependence of probability $\Pr_P(\cdot),$ expectation $E_P(\cdot),$ and norm $\left\| \cdot \right\|_{P,q}$ on the probability measure $P$. \newline Assumption 4 (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,M,X)$ for $d \in \{0,1\}$:

For demonstrating root-$n$ consistency of the proposed estimation strategy for the counterfactual, we heavily draw from Chetal2018. We show that our estimation strategy satisfies the requirements for their double machine learning framework by first verifying the satisfaction of a specific moment condition as well as linearity and Neyman orthogonality of the score (see Appendix (ref)). Then, as e.g.\ $\psi_d(W, \eta, \Psi_{d0})$ is smooth in $ (\eta,\Psi_{d0}) $, the plug-in estimators must converge with rate $ n^{-1/4} $ in order to achieve $ n^{-1/2} $-convergence for the estimation of $ \hat \Psi_d $. This convergence rate of $ n^{-1/4} $ is achievable for many commonly used machine learners such as lasso, random forest, boosting and neural nets. The rates for $L_2$-boosting were, for instance, derived in LuoSpindler2016. \newline

Theorem 1\\ Under Assumptions 1-4, it holds for estimating $E[Y(d, M(1-d))]$, $E[Y(d, m)]$ 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]$. \\ $\sqrt{n} \Big(\hat \Psi_{dm} - \Psi_{dm0} \Big) \rightarrow N(0,\sigma^2_{\psi_{dm}})$, where $\sigma^2_{\psi_d}= E[(\psi_d-\Psi_{dm0})^2]$.\\ The proof is provided in Appendix (ref). \newline

Analogous results follow for the estimation of $\Lambda=E[Y(d, M(d))]$ when replacing $\hat \psi_d$ in the algorithm above by an estimate of score function $\alpha_d$ from (ref),

where $\hat{\mu}^k(d,x)$ is an estimate of $\mu(d,x)$. This approach has been discussed in literature on ATE estimation based on double machine learning, see for instance Bellonietal2017 and Chetal2018. Denoting by $\hat \Lambda$ the estimate of $\Lambda$, it follows under Assumptions 1-4 that $\sqrt{n} \Big(\hat \Lambda_d - \Lambda_d \Big) \rightarrow N(0,\sigma^2_{\alpha_d})$, where $\sigma^2_{\alpha_d}= E[(\alpha_d-\Lambda_d)^2]$. Therefore, root-$n$-consistent estimates of the total as well as the direct and indirect effects are obtained as difference of the estimated potential outcomes, which we denote by $\hat \Delta$, $\hat \theta(d)$, and $\hat \delta(d)$. That is, $\hat \Delta=\hat \Lambda_1-\hat \Lambda_0$, $\hat \theta(1)=\hat \Lambda_1-\hat \Psi_0$, $\hat \theta(0)=\hat \Psi_1 - \hat \Lambda_0$, $\hat \delta(1)=\hat \Lambda_1-\hat \Psi_1$, and $\hat \delta(0)=\hat \Psi_0 - \hat \Lambda_0$.

Naturally, the asymptotic variance of any effect is obtained based on the variance of the difference in the score functions of the potential outcomes required for the respective effect. For instance, the asymptotic variance of $\hat \theta(1)$ is given by $Var(\hat \theta(1))=Var(\alpha_1-\psi_0)/n=(\sigma^2_{\alpha_1}+\sigma^2_{\psi_0}-2Cov(\alpha_1,\psi_0))/n$.

Chetal2018 show that under Assumptions 1-4, $\hat \sigma^2_{\psi_d}$ can be estimated as:

The asymptotic variance of $\alpha_d$ can be estimated accordingly, with $\psi_d$ and $\hat \Psi_{d0}$ substituted by $\alpha_d$ and $\hat \Lambda_{d0}$.

We subsequently discuss estimation based on the score function $\psi_d^*$ in expression (ref). We note that in this case, one needs to estimate a nested nuisance parameter $E\Big[ \mu(d,M,X) \Big| D= 1-d, X \Big]$. To avoid overfitting, the models for $\mu(d,M,X)$ and $E\Big[ \mu(d,M,X) \Big| D= 1-d, X \Big]$ are estimated in different subsamples. The plug-in estimates for the conditional mean outcome, mediator density and treatment probability are referred to by $\hat{\eta}^* = \{ \hat \mu(D,M,X), \hat \omega(D,M,X), \hat p(M,X), \hat p_{d}(X) \}$ and the true nuisance parameters by $\eta_0^* = \{\mu_0(D,M,X), \omega(D,M,X), p_{d0}(M,X), p_{d0}(X) \}$. \newline Algorithm 2: Estimation of $E[Y(d, M(1-d))]$ based on equation (ref)

Also this approach can be shown to be root-$n$-consistent under specific regularity conditions outlined below.\newline Assumption 5 (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,M,X)$ for all $d \in \{0,1\}$:

Theorem 2\\ Under Assumptions 1-3 and 5, it holds for estimating $E[Y(d, M(1-d))]$ based on Algorithm 2: \\ $\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).

Simulation study

This section provides a simulation study to investigate the finite sample behaviour of the proposed methods based on the following data generating process:

Outcome $Y$ is a function of the observed variables $D,M,X$, including an interaction between the mediator and the treatment, and an unobserved term $U$. The binary mediator $M$ is a function of $D,X$ and the unobservable $V$, while the binary treatment $D$ is determined by $X$ and the unobservable $W$. $X$ is a vector of covariates of dimension $p$, which is drawn from a multivariate normal distribution with zero mean and covariance matrix $\Sigma$. The latter is defined based on setting the covariance of the $i$th and $j$th covariate in $X$ to $\Sigma_{ij}=0.5^{|i-j|}$.\footnote{The results presented below are hardly affected when setting $\Sigma$ to the identity matrix (zero correlation across $X$).} Coefficients $\beta$ gauge the impact of $X$ on $Y$, $M$, and $D$, respectively, and thus, the strength of confounding. $ U, V, W$ are random and standard normally distributed scalar unobservables. We consider two sample sizes of $n=1000,4000$ and run $1000$ simulations per data generating process.

We investigate the performance of effect estimation based on (i) Theorem 1 using the identification result in expression (ref) derived by TchetgenTchetgenShpitser2011 as well as (ii) Theorem 2 using the modified score function in expression (ref) which avoids conditional mediator densities. The nuisance parameters are estimated by post-lasso regression based on the `hdm' package by Rhdm for the statistical software `R' with its default options, using logit specifications for $p_d(X)$, $p_d(M,X)$, and $f(M|D,X)$ and linear specifications for $\mu(D,M,X)$ and $E[ \mu(d,M,X)| D=1-d, X ]$. The estimation of direct and indirect effects is based on 3-fold cross-fitting. For all methods investigated, we drop observations whose (products of) estimated conditional probabilities in the denominator of any potential outcome expression are close to zero, namely smaller than a trimming threshold of $0.05$ (or 5%). Our estimation procedure is available in the `causalweight' package for `R' by BodoryHuber2018.

In our first simulation design, we set $p=200$ and the $i$th element in the coefficient vector $\beta$ to $0.3/i^2$ for $i=1,...,p$, meaning a quadratic decay of covariate importance in terms of confounding. This specification implies that the $R^2$ of $X$ when predicting $Y$ amounts to 0.22 in large samples, while the Nagelkerke1991 pseudo-$R^2$ of $X$ when predicting $D$ and $M$ by probit models amounts to 0.10 and 0.13, respectively. The left panel of Table (ref) reports the results for either sample size. For $n=1000$, double machine learning based on Theorem 2 on average exhibits a slightly lower absolute bias (`abias') and standard deviation (`sd') than estimation based on Theorem 1. The behavior of both approaches improves when increasing sample size to $n=4000$, as the absolute bias is very close to zero for any effect estimate and standard deviation is roughly cut by half. Under the larger sample size, differences in terms of root mean squared error (`rmse') between estimation based on Theorems 1 and 2 are very close to zero. By and large, the results suggest that the estimators converge to the true effects at root-$n$ rate.

In our second simulation, confounding is increased by setting $\beta$ to $0.5/i^2$ for $i=1,...,p$. This specification implies that the $R^2$ of $X$ when predicting $Y$ amounts to 0.42, while the Nagelkerke1991 pseudo-$R^2$ of $X$ when predicting $D$ and $M$ amounts to 0.23 and 0.28, respectively. The results are displayed in the right panel of Table (ref). Again, estimation based on Theorem 2 slightly dominates in terms of having a smaller absolute bias and standard deviation, in particular for $n=1000$. However, in other settings, the two methods might compare differently in terms of finite sample performance. Both methods based on Theorems 1 and 2, respectively, appear to converge to the true effects at root-$n$ rate, and differences in terms of root mean squared errors are minor for $n=4000$.

Appendix (ref) reports the simulation results (namely the absolute bias, standard deviation, and root mean squared error) for the standard errors obtained by an asymptotic approximation based on the estimated variance of the score functions. The results suggest that the asymptotic standard errors decently estimate the actual standard deviation of the point estimators.

Application

In this section, we apply our method to data from the National Longitudinal Survey of Youth 1997 (NLSY97), a survey following a U.S. nationally representative sample of 8,984 individuals born in the years 1980-84. Since 1997, the participants have been interviewed on a wide range of demographic, socioeconomic, and health-related topics in a one- to two-year circle. We investigate the causal effect of health insurance coverage ($D$) on general health ($Y$) and decompose it into an indirect pathway via the incidence of a regular medical checkup ($M$) and a direct effect entailing any other causal mechanisms. Whether or not an individual undergoes routine checkups appears to be an interesting mediator, as it is likely to be affected by health insurance coverage and may itself have an impact on the individual's health, because checkups can help identifying medical conditions before they get serious to prevent them from affecting a person's general health state.

The effect of health insurance coverage on self-reported health has been investigated in different countries with no compulsory medical insurance and no publicly provided universal health coverage, see for example simon2017impact, sommers2017three, baicker2013oregon, yoruk2016health and cardella2014effect for the U.S.\ and king2009public for Mexico). Most of these studies find a significant positive effect of insurance coverage on self-reported health. The impact of insurance coverage on the utilization of preventive care measures, particularly routine checkups like cancer, diabetes and cardiovascular screenings, is also extensively covered in public health literature. Most studies find that health insurance coverage increases the odds of attending routine checkups. While some contributions include selected demographic, socioeconomic and health-related control variables to account for the endogeneity of health insurance status (see e.g. faulkner1997effect, press2014insurance, burstin1998effect, fowler2007risk), others exploit natural experiments: simon2017impact estimate a difference-in-differences model comparing states which did and did not expand Medicaid to low-income adults in 2005, while baicker2013oregon exploit that the state of Oregon expanded Medicaid based on lottery drawings from a waiting list. The results of both studies suggest that the Medicaid expansions increased use of certain forms of preventive care. In a study on Mexican adults, pagan2007health use self-employment and commission pay as instruments for insurance coverage and also find a more frequent use of some types of preventive care by individuals with health insurance coverage.

While the bulk of studies investigating checkups focus on one particular type of screening (rather than general health checkups), see maciosek2010greater for a literature review, several experimental contributions also assess general health checkups. For instance, rasmussen2007preventive conduct an experiment with individuals aged 30 to 49 in Denmark by randomly offering a set of health screenings, including advice on healthy living and find a significant positive effect on life expectation. In a study on Japan’s elderly population, nakanishi1996preventive find a significantly negative correlation between the rate of attendance at health check-ups and hospital admission rates. Despite the effects of health insurance coverage and routine checkups being extensively covered in the public health literature, the indirect effect of insurance on general health operating via routine checkups as mediator has to the best of our knowledge not yet been investigated. A further distinction to most previous studies is that we consider comparably young individuals with an average age below 30. For this population, the relative importance of different health screenings might differ from that for other age groups. We also point out that our application focuses on short-term health effects.

We consider a binary indicator for health insurance coverage, equal to one if an individual reports to have any kind of health insurance when interviewed in 2006 and zero otherwise. The outcome, self-reported general health, is obtained from 2008 interview and measured with an ordinal variable, taking on the values `excellent', `very good', `good', `fair' and `poor'. In the 2007 interview, participants were asked whether they have gone for routine checkups since the 2006 interview. This information serves as binary mediator, measured post-treatment but pre-outcome.

To ensure that the control variables ($X$) are not influenced by the treatment, they come from the pre-treatment 2005 and earlier interview rounds. They cover demographic characteristics, family background and quality of the home environment during youth, education and training, labor market status, income and work experience, marital status and fertility, household characteristics, received monetary transfers, attitudes and expectations, state of physical and mental health as well as health-related behavior regarding e.g. nutrition and physical activity. For some variables, we only consider measurements from 2005 or from the initial interview round covering demographics and family related topics. For other variables we include measurements from both the indiviuals' youth and 2005 in order to capture their social, emotional and physical development. Treatment and mediator state in the pre-treatment period (2005) are also considered as potential control variables. Item non-response in control variables is dealt with by including missing dummies for each control variable and setting the respective missing values to zero. In total, we end up with a set of 770 control variables, 601 of which are dummy variables (incl. 252 dummies for missing values).

After excluding 1923 observations with either mediator or treatment status missing, we remain with 7,061 observations. Table (ref) presents some descriptive statistics for a selection of control variables. It shows that the group of individuals with and without health insurance coverage differ substantially. There are significant differences with respect to most of the control variables listed in the table. Females are significantly more likely to have health insurance coverage. Education and household income also show a significant positive correlation with health insurance coverage while the number of household members for example is negatively correlated with insurance coverage. Regarding the mediator, we find a similar pattern as for the treatment. With respect to many of the considered variables, the group of individuals who went for medical checkup differs substantially from those who did not. Further, we see that the correlation between many control variables and the treatment appear to have the same sign as that with the mediator.

In order to assess the direct and indirect effect of health insurance coverage on general health, we consider estimation based on Theorem 1 and expression (ref) derived by TchetgenTchetgenShpitser2011 well as (ii) Theorem 2 and expression (ref). We estimate the nuisance parameters and treatment effects in the same way as outlined in Section (ref) (i.e.\ post-lasso regression for modeling the nuisance parameters and 3-fold cross fitting for effect estimation) after augmenting the set of covariates with 3101 interaction and higher order terms. The trimming threshold for discarding observations with too extreme propensity scores is set to $0.02$ (2%), such that 893 and 136 observations are dropped when basing estimation on Theorems 1 and 2, respectively.

Table (ref) provides the estimated effects along with the standard error (`se') and p-value (`p-val') and also provides the estimated mean potential outcome under non-treatment for comparison (`$\hat{E}[Y(0,M(0))]$'). The ATEs of health insurance coverage on general health in the year 2008 (columns 2 and 8), estimated based on Theorems 1 or 2, are statistically significant at the 5% and 10% levels, respectively. As the outcome is measured on an ordinal scale ranging from `excellent' to `poor', the negative ATEs suggest a short-term health improving effect of health coverage. The direct effects under treatment (columns 3 and 9) under non-treatment (columns 4 and 10) are very similar to the ATEs and statistically significant (at least) at the 10% level in 3 out ouf 4 cases. In contrast, the indirect effects under treatment (columns 5 and 11) and non-treatment (columns 6 and 12) are generally close to zero and not statistically significant at the 10% level in 3 out of 4 cases. Thus, health insurance coverage does not seem to importantly affect general health of young adults in the U.S.\ through routine checkups in the short run, but rather through other mechanisms.

Conclusion

In this paper, we combined causal mediation analysis with double machine learning under selection-on-observables assumptions which avoids adhoc pre-selection of control variables. Thus, this approach appears particularly fruitful in high-dimensional data with many potential control variables. We proposed estimators for natural direct and indirect effects as well as the controlled direct effect exploiting efficient score functions, sample splitting, and machine learning-based plug-in estimates for conditional outcome means, mediator densities, and/or treatment propensity scores. We demonstrated the root-$n$ consistency and asymptotic normality of the effect estimators under specific regularity conditions. Furthermore, we investigated the finite sample behavior of the proposed estimators in a simulation study and found the performance to be decent in samples with several thousand observations. Finally, we applied our method to data from the U.S.\ National Longitudinal Survey of Youth 1997 and found a moderate short-term effect of health insurance coverage on general health, which was, however, not mediated by routine checkups. The estimators considered in the simulation study and the application are available in the `causalweight' package for the statistical software `R'.