Multiply Robust Causal Mediation Analysis with Continuous Treatments

Multiply Robust Causal Mediation Analysis with Continuous Treatments

\def\spacingset#1{ {#1}} \spacingset{1}

}

{\bf Keywords:} Causal inference; Debiased estimator; Kernel smoothing; Semiparametric estimators.

\spacingset{1.7}

Introduction

Estimating the causal effect of a treatment, policy, or intervention on an outcome of interest is a fundamental task in various fields such as epidemiology, economics, medicine, and sociology. A common parameter of interest is the average causal effect (ACE), which has been extensively studied hernan2020causal, imbens2015causal. However, in addition to estimating the ACE, one can also be interested in the pathways and mechanisms through which the treatment affects the outcome of interest. Causal mediation analysis offers a precise and rigorous mathematical framework to answer such questions robins1992identifiability, tchetgen2012semiparametric, pearl2001direct, vanderweele2009marginal, goetgeluk2008estimation, imai2010identification, van2008direct, lange2011direct, lange2012simple.

Much of the literature on mediation analysis assumes that the treatment of interest is binary. However, interventions involving the dosage of a drug, and the duration or frequency of an activity are better described as continuous variables. In such cases, mediation effects are naturally represented by a multi-dimensional surface rather than a scalar parameter. This learning task is challenging if a priori shape constraints are not imposed on the surface. Additionally, the presence of continuous treatments complicates the estimation of nuisance parameters, making the estimation of causal parameters more challenging.

The challenges related to estimating ACE in the continuous treatment setting have been addressed in multiple work kennedy2017nonparametric, AI2021, imbens-continuous, kreif2015evaluation, 10.2307/2673642, su2019non,kallus2018policy, colangelo2020double, hill2011bayesian. A common method is based on outcome regression, which requires the correct specification of the relevant models, and hence machine learning methods such as Bayesian additive regression trees (BART) hill2011bayesian are often used. However, this inherits the rate of the outcome regression estimation, and complex machine learning methods tend to have a slower convergence rate than simple parametric methods wasserman2006all, tsybakov2009nonparametric. An alternative approach involves specifying a parametric form for the dose-response curve or projecting the true curve onto a parametric model, as presented in robins2000marginal, van1998locally, neugebauer2007nonparametric. However, these methods may suffer from bias when the dose-response curve is misspecified. In contrast to approaches involving parametric assumptions on the dose-response curve, kennedy2017nonparametric leverage semiparametric theory by utilizing a two-stage estimator that first constructs a doubly robust pseudo-outcome in the first stage, and then regresses the pseudo-outcome on the treatment in the second stage using non-parametric regression methods. colangelo2020double utilize double machine learning along with applying kernel smoothing to the augmented inverse propensity weighted (AIPW) score robins1995semiparametric. This allows a slower convergence rate of nuisance parameters, while still guaranteeing fast rates for the target parameter. However, these approaches are not investigated for mediation analysis in the presence of continuous treatments.

In this paper, we propose a kernel smoothing approach inspired by influence function-based estimators tsiatis2007semiparametric,newey1994asymptotic, bickel1993efficient,ichimura2015influence, tchetgen2012semiparametric to deal with continuous treatments for causal mediation analysis. We propose an estimator that, under mild regularity conditions, is consistent and asymptotically normal. Our work aims to extend the results for the continuous treatment ACE to the case of mediation analysis involving continuous treatments in the presence of high-dimensional nuisance functions. huber2020direct tackle this problem by weighting the observations with a generalized propensity score that involves two nuisance functions, which are the conditional density of treatment given covariates and the conditional density of treatment given mediators and covariates. In their proposed approach, the nuisance functions can be estimated either parametrically or non-parametrically. However, their estimator for the causal parameter is not robust with respect to the misspecification of the two nuisance functions and also inherits the rate of the nuisance function estimators, which could be slow. In contrast, we propose an approach motivated by influence functions and hence obtain many of the desirable properties of influence functions, namely allowing for slower estimation of nuisance functions, as well as robustness properties. Our work draws heavily on the existing causal mediation literature that discusses the identification and estimation of such effects pearl2001direct, imai2010identification, tchetgen2012semiparametric. Additionally, we utilize the cross-fitting strategy to relax the smoothness assumptions on the nuisance functions chernozhukov2018double.

The remainder of this paper is organized as follows. Section (ref) introduces the formal mediation framework, describes its identifying assumptions, and discusses an influence function-based estimator of mediation effects for binary treatments. Section (ref) extends the influence function-based approach to continuous treatment settings and describes the sample-splitting and smoothing procedures. In Section (ref), we provide our main results along with the required regularity conditions. Sections 5 and 6 present simulation and application results. We provide conclusions and discussion of future directions in Section 7.

Mediation Analysis

Let $A\in\mathcal{A}$ be the continuous treatment variable, $Y\in\mathcal{Y}$ be the outcome variable, and $M\in\mathcal{M}$ be a mediator variable that relays part of the causal effect of $A$ on $Y$. In addition, let $X\in\mathcal{X}$ denote the observed pre-treatment covariates in the system. See Figure (ref)$(c)$ for a graphical representation of the causal relationships between the variables. To describe the causal effect of the treatment on the outcome, we use the potential outcome framework rubin1974estimating. Let $Y^{(A=a)}$ be the random variable representing the potential outcome when the treatment is set to value $a$. Suppose we are interested in comparing the treatment values of $a$ and $a'$. A popular way to measure the causal effect of this change in treatment is to use the average causal effect (ACE), which captures the difference in the expected value of the potential outcome variables, that is \[ ACE=\mathbb{E}[Y^{(a)}-Y^{(a')}], \] where $\mathbb{E}[\cdot]$ denotes the population expectation operator.

The total ACE of the treatment $A$ on the outcome $Y$ can be partitioned into the part mediated by variable $M$ and the part directly affecting outcome $Y$ (see Figure (ref)). To formally define this partitioning, let $Y^{(a,m)}$ denote the potential outcome variable corresponding to the outcome when the treatment is set to value $a$ and the mediator is set to value $m$, and $M^{(a)}$ denote the mediator variable when the treatment is set to value $a$.

robins1992identifiability and pearl2001direct proposed the following partitioning of the ACE into the natural direct and indirect effects.

The natural direct effect (NDE) and natural indirect effect (NIE) can be described as follows. NDE captures the change in the expectation of the outcome if the value of the treatment variable is switched between the two arms of the experiment, while the mediator behaves as if the treatment has not changed. NIE captures the change in the expectation of the outcome if the value of the treatment variable is fixed, while the mediator behaves as if the treatment has been switched between the two arms of the experiment. In the following subsection, we discuss the estimation of NDE and NIE from observational data.

Estimating Natural Direct and Indirect Effects

To estimate the natural direct and indirect effects, from the partitioning in Display (ref), it suffices to focus on estimating the parameter \[ \psi_0(a,a')=\mathbb{E}[Y^{(a,M^{(a')})}], \] for $a,a'\in\mathcal{A}$. Suppose i.i.d. data from a distribution $f$ on variables $O=\{A,X,M,Y\}$ are given. In general, the estimand $\psi_0$ is not identified from observational data, and identification assumptions are required to relate the distribution of the observational data to that of counterfactual variables. We require the following assumptions for the identification of $\psi_0$ from the observed distribution on variables, $f(O)$.

Under Assumption (ref), the estimand $\psi_0(a,a')$ can be identified from the observed distribution $f(O)$ using the following expression called the mediation formula, originally proposed in imai2010identification and pearl2001direct.

where $f_X$ is the marginal distribution of $X$.

Using Equation (ref), one can estimate the parameter of interest $\psi_0(a,a')$ by first estimating the nuisance functions $\mathbb{E}[Y|A,M,X]$ and $f_{M|A,X}$, and then using a plug-in estimator to estimate $\psi_0$ as follows \[ \frac{1}{n}\sum_{i=1}^n \int_\mathcal{M}\hat{\mathbb{E}}[Y_i|A=a,M=m,X_i]\hat{f}_{M|A,X}(m|A=a',X_i)dm. \] Unfortunately, this estimator is sensitive to bias in the estimation of nuisance functions. That is, misspecifying either of the nuisance functions induces bias in the estimation of the parameter of interest.

As an alternative approach, in the case of binary treatment, that is, $\mathcal{A}=\{0,1\}$, tchetgen2012semiparametric developed a general semiparametric framework for inference. They derived the efficient influence function for $\psi_0(a,a')$ as

where $\lambda(a,X) := 1/{f(a|X)}$, $\alpha(a,M,X) := f(M| a, X)$, and $\gamma(X,M,a) := \mathbb{E}[Y|A=a,M,X]$ are the nuisance functions, $a,a'\in\{0,1\}$, $I(\cdot)$ denotes the indicator function, and

\[ \eta(a, a', X) = \int_\mathcal{M} \gamma(X, m, a) \alpha(a', m,X) dm. \] Note that $IF_{\psi_0}$ is a function of three nuisance functions: $\lambda$, $ \alpha$, and $\gamma$. tchetgen2012semiparametric showed that the estimator based on this influence function has the triple robustness property, that is, it is consistent even if the model for one (but not more than one) nuisance function is misspecified. Formally, let

The estimator for $\psi_0$ based on the influence function $IF_{\psi_0}$ defined as

is consistent when the truth lies in the submodel union $\mathfrak{M}_{ym}\cup\mathfrak{M}_{ya}\cup\mathfrak{M}_{ma}$ and the estimators of nuisance functions in the correctly specified submodel(s) are consistent, where $$\hat{\eta}(a, a', X) = \int_\mathcal{M} \hat{\gamma}(X,m,a) \hat{\alpha}(a',m, X)dm.$$ Inspired by this result, in the following section, we propose a kernel-based estimator for meditation effects in settings with continuous treatment variables, while preserving multiple robustness and allowing for the nuisance parameters to be estimated at a slower rate than the parameter of interest.

Continuous Mediation Analysis

In the case of continuous treatments, the parameter of interest, $\psi_0(a,a')$, is no longer regular bickel1993efficient,tsiatis2007semiparametric. Therefore, the method of tchetgen2012semiparametric cannot be applied directly. However, their estimator can be modified to be suitable for inference in the case of continuous treatments, while still obtaining desirable properties such as asymptotic normality, robustness to misspecification of nuisance functions, and valid inference while allowing for the nuisance parameters to be estimated at a slower rate than the parameter of interest. Specifically, we modify $\hat{\psi}^{TTS}(a,a')$ by employing a kernel smoothing technique, wherein the indicators in the calculation of $\hat{\psi}^{TTS}(a, a')$ are replaced by kernel-based weights. The weights are functions of treatment values falling within a neighborhood (defined by the bandwidth parameter $h$) of $a$ and $a'$. This modification introduces several challenges in the inference, which we will present and address in Section (ref).

Let $d_A$ denote the dimension of the treatment variable, and let \[ K_h(a):=\frac{1}{h^{d_A}}\prod_{j=1}^{d_A}k\Big(\frac{a_j}{h}\Big), \] where $k(\cdot)$ is a kernel function, and $h$ denotes the bandwidth parameter. We propose to use the following modification of the efficient influence function in Equation (ref) for any $a$ and $a' \in \mathcal{A}$:

To derive desired results on consistency, asymptotic normality, and multiple robustness, we require the kernel $k(\cdot)$ to satisfy the following conditions.

These assumptions are satisfied by common kernels such as the Gaussian kernel and Epanechnikov kernel. Note that in the moment function in Equation (ref), the nuisances are not functions of the parameter of interest $\psi_0$. Therefore, having estimators for nuisance functions suffices for obtaining an estimator for $\psi_0$. Next, we describe the estimation procedure for utilizing Equation (ref) to estimate $\psi_0(a, a')$.

Estimation Procedure. We use the cross-fitting estimation approach of chernozhukov2018double for separating the estimation of the nuisance functions from the parameter of interest. This approach is beneficial since weaker smoothness requirements are needed for the estimation of nuisance functions. In the cross-fitting approach, we partition the sample indices into $L$ folds $\{I_1,...,I_L\}$ of roughly equal size. Data from the $\ell$-th fold is denoted by $O_{I_\ell}$, and the data in the rest of the folds is denoted by $O^c_{I_\ell}$. For $\ell\in\{1,...,L\}$, we estimate the nuisance functions $\alpha$, $ \lambda$, $\gamma$ by $\hat{\alpha}_\ell,\hat{\lambda}_\ell,\hat{\gamma}_\ell$ based on data $O^c_{I_{\ell}}$. For all $\ell$, let $\hat{\psi}_\ell$ be the estimator for $\psi_0$ obtained by solving \[ \frac{1}{|I_\ell|}\sum_{i\in I_\ell}m(O_i;\hat{\alpha}_\ell,\hat{\lambda}_\ell,\hat{\gamma}_\ell,\hat{\psi}_\ell(a,a'))=0. \] Our proposed estimator for $\psi_0$ is

where MR stands for multiply robust. In the next section, we present the asymptotic properties of our proposed estimator, along with the required regularity conditions.

Asymptotic Analysis

In this section, we provide asymptotic properties of our proposed estimator $\hat{\psi}^{MR}(a,a')$ in Equation (ref). We start by stating the required regularity conditions.

In addition to the regularity conditions, we require the following conditions regarding the convergence of the estimators of the nuisance functions.

Similar to influence function-based estimators, in Assumption (ref), we do not require individual nuisance estimators to satisfy convergence rate conditions. However, in our proposed method, we have requirements on the convergence rate of the product for the nuisance estimators as follows.

As seen in Assumption (ref), our requirements on the convergence rate of nuisance function estimators are on the product of the error rates, rather than on the individual nuisance function estimators. Therefore, if one of the estimators converges at a slow rate, the other estimator can compensate. This is a desirable property when working with non-parametric estimators since they typically have slow rates of convergence.

In the case of binary treatment variables, the combination of Assumptions (ref) and (ref) can lead to asymptotic normality, which is used to construct Wald-style confidence intervals. However, when the treatments are continuous, the Central Limit Theorem (CLT) cannot be directly applied to our proposed method because the bandwidth $h$ varies as a function of the sample size $n$, implying that the distribution of equation (ref) changes with $n$. Instead, we impose additional assumptions stated below to satisfy the Lyapunov's condition for CLT and achieve asymptotic normality.

Having stated the assumptions in our setting, we now provide the following result regarding the asymptotic behavior of the proposed estimator in Equation (ref).

All proofs are provided in the Supplementary Material. Theorem (ref) provides results on the point-wise convergence of $\hat{\psi}^{MR}(a,a')$ and establishes the asymptotic normality of our estimator. Additionally, $\hat{\psi}^{MR}(a,a')$ has a multiple robustness property analogous to the $\hat{\psi}^{TTS}(a,a')$, formally stated in Lemma (ref).

While Theorem (ref) and Lemma (ref) establish properties of $\hat{\psi}^{MR}(a,a')$ that are desirable for point estimation, uncertainty quantification through the calculation of valid confidence intervals requires the estimation of $V(a, a')$ and $B(a, a')$. However, these are hard to estimate well because of their complicated analytical form. Nevertheless, by choosing an undersmoothing bandwidth $h$ that satisfies $\sqrt{nh^{d_A + 4} } \rightarrow 0$, valid confidence intervals can still be constructed without estimating $B(a, a')$. In this case, we only need to focus on estimating $V(a, a')$. Naturally, an estimator for $V(a, a')$ can be constructed as follows. \[\widehat{V}(a, a') = h^{d_A}\frac{1}{L}\sum_{\ell = 1}^L \frac{1}{|I_\ell|}\sum_{i \in I_\ell} m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}_\ell(a,a')). \]

We present the additional assumptions necessary for the consistency of $\widehat{V}(a, a')$ below.

Next, we prove that $\widehat{V}(a, a')$ is a consistent estimator of $V(a, a')$ in Lemma (ref).

Using the result of Lemma (ref), $\widehat{V}(a, a')$ can be used to construct asymptotically valid confidence intervals as follows. Choose an undersmoothing bandwidth $h$ that satisfies $\sqrt{nh^{d_A + 4}} \rightarrow 0$, so that $\sqrt{nh^{d_A}}h^2 B(a, a')$ is asymptotically negligible. Then, the $(1 - \alpha)$ confidence interval is given as

where $\Phi$ is the CDF of $\mathcal{N}(0,1)$. However, there is little guidance about how to undersmooth in practice and it is mostly a technical device kennedy2017nonparametric. In simulation and application, we express uncertainty about the estimator $\hat{\psi}(a,a')$ using a bandwidth chosen under the Silverman rule of thumb van2000asymptotic, silverman2018density using pointwise confidence intervals (ref). As stated by wasserman2006all, such an approach does not necessitate the choice of an undersmoothing bandwidth, but the constructed confidence intervals center around the parameter of interest $\psi_0(a,a')$ with an additional asymptotic bias.

Simulation Study

We now present a simulation study to demonstrate that the proposed estimator is consistent and multiply robust, and a sensitivity analysis to assess the uncertainty of the proposed estimator under different bandwidths and sample sizes. The data-generating process is as follows:

The parameter of interest is $\psi_0(a, a')$ at $a = 4.5$ and $a' = 6$. Under the described simulation setting, the true parameter value is $9.1$, calculated based on Monte Carlo approximation of $\psi_0(a, a') = \int_\mathcal{X}\eta(a,a',X=x)dx$. To demonstrate the multiple robustness property, we consider various types of model misspecification in Table (ref), where we also compare the proposed estimator $\hat{\psi}^{MR}(a, a')$ to the estimator in huber2020direct, $\hat{\psi}^H(a, a')$, and the estimator without bias correction, $\hat{\psi}^{\eta}(a,a')$. To ensure comparability, we calculate the three estimators under the cross-fitting approach and define $\hat{\psi}^H(a, a')$ and $\hat{\psi}^{\eta}(a, a'))$ as $$\hat{\psi}^H(a, a') = \frac{1}{L}\sum^L_{\ell = 1}\hat{\psi}^H_\ell(a, a')\quad\text{and}\quad \hat{\psi}^\eta(a, a') = \frac{1}{L}\sum^L_{\ell = 1}\hat{\psi}^\eta_\ell(a, a'),$$ where

$\hat{\eta}(a,a',X) = \int_\mathcal{M} \hat{\gamma}(X, m, a) \hat{\alpha}(a',m,X)dm$, and we set $L=3$. We use 1000 simulation replicates for each of the sample sizes 2000, 5000, and 8000, and chose the kernel bandwidth using the Silverman rule of thumb silverman2018density under Gaussian kernels. The types of model misspecification considered include the scenario where all three models, $\mathbb{E}[Y|A,M,X]$, $f(M|A,X)$, and $f(A | X)$, are correctly specified, scenarios where only two out of the three models are correctly specified, and the scenario where all three models are misspecified. As shown in Table (ref), our proposed estimator has minimal or close to minimal bias for all scenarios except when all models are misspecified, demonstrating its theoretically proven multiple robustness property. Additionally, bias and the root mean square error (RMSE) across simulation replicates reduce as the sample size gets larger, showing the consistency of our estimator. The bias becomes significant when all models are misspecified for all sample sizes and considered estimators.

Table (ref) summarizes the sensitivity of the estimated $\hat{\psi}^{MR}(a,a')$ under correct model specifications to different sample sizes and kernel smoothing bandwidths, by reporting absolute average bias, average of $\hat{V}(a,a')^{1/2}$, and Monte Carlo coverage. We define coverage as the proportion of simulation replicates that include the true value $\psi_0(a, a')$ in the estimated 95$\%$ confidence interval. Recall that $\hat{V}(a, a')$ is the empirical variance of the estimating functions; we construct pointwise 95% confidence intervals at significance level $\alpha = 0.05$ via a Wald approach as in equation (ref). The Silverman kernel bandwidths for n = 2000, 5000, 8000, and under $L=3$ are approximately 0.28, 0.23, and 0.20. For the estimated mediation function $\hat{\psi}^{MR}(a,a')$, the presence of bias becomes apparent with reduced length of confidence intervals and decreasing coverage when the kernel bandwidths are larger than the Silverman-suggested optimal bandwidths, i.e., when the bandwidths are greater or equal to $0.3$. This pattern persists across different sample sizes. Coverage is theoretically guaranteed when the sample size $n$ goes to infinity and $\sqrt{nh^{d_A + 4}} \rightarrow 0$. For example, when $d_A = 1$, an undersmoothing bandwidth would satisfy the requirement. We can see from Table (ref) that when bandwidth equals 0.1 (undersmoothed), coverage of the proposed pointwise 95% confidence interval is indeed at least 0.95 for all sample sizes. In practice, choosing an undersmoothed kernel bandwidth can guarantee a relatively smaller bias in finite sample settings (for the price of conservative coverage). However, as previously discussed, there is no guidance on the choice of an undersmoothing bandwidth in practice.

Application

The proposed approach is applied to the Job Corps study huber2020direct,schochet2008does,schochet2001national. Study participants were enrolled between 16 and 24 years old and from low-income households. The program provides eight months or approximately 1,200 hours of training on average. We aim to study the effect of the duration of Job Corp training $(A)$ on the binary outcome of the occurrence of any criminal arrests in the fourth year following program participation $(Y)$, with the proportion of weeks employed in the second year being the mediator $(M)$. Our study design follows huber2020direct, who considered a similar causal mechanism and focused on the actual number of arrests in the fourth year as the outcome.

We consider a rich set of time-invariant socioeconomic variables as pretreatment confounders $X$, similar to the study in huber2020direct. Table (ref) presents summary statistics of the following variables: outcome $Y$, mediator $M$, treatment $A$, and confounders $X$. Missing values in confounders are addressed by including the indicators of missingness as covariates. Moreover, following previous work on this dataset huber2020direct,flores2012estimating, we apply our evaluation to the 4,000 individuals in the dataset who received training, i.e., with a training duration in the program strictly greater than zero. Table (ref) shows that on average, 5.1$\%$ of the participants had a history of imprisonment, and 23.75$\%$ had been arrested at least once before joining the study. Additionally, 8.7$\%$ of the individuals included in the study were arrested for criminal activities during the fourth year after study participation.

We investigate whether longer Job Corps training reduces criminal behavior by evaluating the natural direct effect for treatment durations at $a\in\{100,200,\ldots,2000\}$ hours versus just $a'=40$ hours. Following huber2020direct, we assume treatment to follow a log-normal distribution and parametric linear models for the outcome, mediator, and log-treatment. Let $\hat{f}(a|X)$ be a model-based propensity score estimator at treatment $a$. To improve the stability of our estimator, we use a H\'ajek-type stabilized weighted propensity score hernan2020causal, which is defined as $\hat{f}(a|X_i) \times \frac{1}{|I_{-\ell}|}\sum^{|I_{-\ell}|}_{j\in I_{-\ell}} \frac{K_h (A_j-a)}{\hat{f}(a|X_j)}$. This prevents a large influence from observations that have extremely small propensities $\hat{f}(a|X)$.

In Supplementary Material Section 5, we show that the proposed Hajek-type stabilized weighted propensity score is a consistent estimator of the conditional treatment probability when the model-based propensity estimator is consistent. Figure (ref) displays the mean and 95$\%$ confidence interval of natural direct effect over the range of values for $a$, under H\'ajek-type stabilized weighted propensities clipped at 0.01 ionides2008truncated and Gaussian kernels with bandwidth chosen using the Silverman-type rule of thumb silverman2018density. The confidence interval is obtained via Equation (ref). Figure (ref) reports the sensitivity analysis and demonstrates how the estimated mean and empirical standard deviation vary by different clipping thresholds and smoothing bandwidths.

Both schochet2008does and huber2020direct identified significant effects of the training program in reducing criminal arrests, especially when the training duration is over 1000 hours. Similar to huber2020direct, we observed a nonlinear effect of the training duration on the occurrence of arrests. However, our sensitivity analysis showed that a significant reduction in the likelihood of arrests (at 5$\%$ level) does not always occur after a certain number of hours in Job Corps. As shown in the sensitivity analysis, the significance and size of the natural direct effect also depend on two factors: the choice of bandwidth for smoothing the continuous treatment and the threshold for clipping the estimated propensity score. Regardless of the actual treatment effect level, we generally noticed a decline in the probability of experiencing at least one arrest in year 4 for those with training durations between approximately 500 and 1,300 hours. From summary descriptions in Table (ref), we know that the interquartile range (IQR) of the training durations in the observed data falls within the range of approximately 405 to 1,767 hours. Hence, it can be inferred that extended training has a tendency to diminish criminal behaviors among the majority of the study participants.

Conclusion

In this paper, we present an estimator for the direct and indirect components of the causal effect in the presence of continuous treatments, while allowing for complex data generating processes. Our estimator is inspired by influence function based estimation strategies from the literature on semiparametric statistics bickel1993efficient, tsiatis2007semiparametric. We provide results on the convergence rate and asymptotic normality of our proposed estimator, which allows for the construction of confidence intervals and hypothesis tests. Extending the point-wise results to uniform results and providing optimal rates is left for future work.

in the presence of high-dimensional nuisance parameters, the estimator may inherit the slow rates of the nuisance estimators, adversely affecting the inference of the target parameter.

Before we start with the proofs, we establish some lemmas that will help us with the proofs in the rest of the Supplementary Material.

Proof For Asymptotic Normality

The proof for asymptotic normality follows from an application of the Lyapunov Central Limit theorem to the terms $\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))$. We first prove the Lyapunov condition holds for $\delta = 1$, i.e. \[ \lim_{n \rightarrow \infty} \frac{1}{s^{3}_n} \sum_{i = 1}^n \mathbb{E}\left[\Big|\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) - \mu_{i} \Big|^3\right] = 0 \] Where $\mu_i$ equals $\mathbb{E}\left[\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\right]$ and $s^2_n = \sum_{i = 1}^n \sigma^2_i$ where $\sigma^2_i$ is the variance of of $\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))$. To prove the Lyaponuv condition holds, we first derive $\mu_i$ and $\sigma^2_i$.

Calculation for $B(a, a')$ and $\mu_i$

Given

Since $\mathbb{E}[\eta(a, a', X_i) - \psi_0(a,a')] = 0$, we focus on $\frac{K_h(A - a)f(M \mid A = a', X)}{f(M \mid A = a, X)f(a \mid X)}\{Y - \mathbb{E}[Y \mid X, M , A = a]\}+ \frac{K_h(A - a')}{f(a' \mid X)}\{\mathbb{E}[Y \mid X, M, A = a] - \eta(a, a', X)\}$.

We start by computing the expectation of each the individual terms one at a time.

Expectation Part 1 $$\frac{K_h(A - a)f(M \mid A = a', X)}{f(M \mid A = a, X)f(a \mid X)}\{Y - \mathbb{E}[Y \mid X, M , A = a]\}$$

From $\mathbb{E}\{\mathbb{E}[\gamma(X,M,a)|X,M]\} = \mathbb{E}\{\mathbb{E}[Y|X,M]\}$ and $\gamma(X,M,a) = \mathbb{E}(Y|X,M,A)$, expectation of the first term

The inner product further expands as follows,

for all $X,M$ in respective range. Inserting this back into the original expectation we get, {

}

Expectation Part 2 $$\frac{K_h(A - a')}{f(a' \mid X)}\{\gamma(X,M,a) - \eta(a, a', X)\}$$

The inner expectation can be written as

Plugging this back into the above expectation {

} from having the first term in this expectation equal to zero, which we prove below

Hence, letting

we have $\mathbb{E}\left[m(O_i;\alpha,\lambda,\gamma, \psi_0(a,a'))\right] = h^2 B(a,a')$. Additionally from this derivation $\mathbb{E}\left[\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\right] = O(\sqrt{\frac{h^{d_A + 4}}{n}})$. Next, we prove the properties of variance.

Calculation for $V(a, a')$ and $s^2_n$

From the definition of $s^2_n$, we have

Consequently, we calculate

Using the property that $var(X) = \mathbb{E}[X^2] - \mathbb{E}[X]^2$ and constant values do not contribute to the variance, the variance term above can be re-written as

Examining each of the terms above one by one, the first term can be expanded as

We analyze each of these terms part by part.

Variance Part 1

Because $0 < \int u^6 k(u) du < \infty$ from Assumption 2 (3), we also have boundedness of $\int u^6k^2(u)du$. The inner expectation can be written as

where $\Bar{a}_v,\Bar{a}_\gamma$, and $\Bar{a}_f$ are between $a$ and $a+uh$. Hence, part 1 of the variance

Variance Part 2

The inner expectation can be written as

the last equation is from

Hence, the part 2 of variance

Variance Part 3

This holds because we assume $\eta$ is bounded.

Variance Part 4

The inner expectation

Hence, the part 4 of variance

Variance Part 5

Applying the same expansion in Expectation Part 1, we can write the inner expectation as

Inserting this back into the full expectation, combined with the boundedness of $\eta$, $f(M \mid a, X)$ and $f(a \mid X)$, we get

Variance Part 6

Using the expansion from Part 2 of the expectation on $\mathbb{E}[K_h(A - a')| X,M]$, we get

Plugging this back into the full expectation, we get

And using the calculation for the bias,

Finally, putting the pieces of the variance together, we have

where the term converges to $V(a, a^\prime)$ as $h\rightarrow 0$ and

Having derived the bias and variance terms, we now prove the Lyapunov condition for $\delta = 1$.

Proof for Lyapunov Condition

We now prove the Lyapunov condition \[ \lim_{n \rightarrow \infty} \frac{1}{s^{3}_n} \sum_{i = 1}^n \mathbb{E}\left[\Big|\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) - \mu_{i} \Big|^3\right] = 0 \]

Note that $$\Big|\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) - \mu_{i} \Big| \leq \Big|\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\Big| + \Big| \mu_{i} \Big|$$ Since both sides are positive,

From the monotonicity of the expected value, we have

Since $\sum_{i = 1}^n\Big| \mu_{i} \Big|^3 = O(h^{(d_A + 4)3/2}n^{-1/2}) = o(1)$, $\sum_{i = 1}^n 3h^{d_A}n^{-1}\mathbb{E}\left[\Big| m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\Big|^2\right]\Big|\mu_i\Big| = O(\sqrt{\frac{h^{d_A + 4}}{n}}) = o(1)$, and $\sum_{i = 1}^n 3(h^{d_A}n^{-1})^{1/2}\mathbb{E}\left[\Big| m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\Big|\right]\Big|\mu_i\Big|^2 = o(1)$, it suffices to prove the following condition \[ \lim_{n \rightarrow \infty} \frac{1}{s^{3/2}_n} \sum_{i = 1}^n \mathbb{E}\left[\Big|\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) \Big|^3\right] = 0 \]

Following a similar derivation as in the proof for consistency of $\hat{V}(a, a')$, from the assumption $\mathbb{E}\{[Y-\gamma(X,M,a)]^3 | A = a', M = m, X = x\}$ over any $(a,a',m,x)\in \mathcal{A}\times\mathcal{A}\times \mathcal{M} \times \mathcal{X}$, along with $\int^{\infty}_{-\infty} k(u)^{c_1}k(u+\tilde{c})^{c_2} du < \infty$ and $\int^{\infty}_{-\infty} u^2k(u)^{c_1}k(u+\tilde{c})^{c_2} du < \infty$ for $\tilde{c} \in \mathcal{R}$ and $c_1+c_2 \in \{2,3\}$ for $c_1,c_2\in\{0,1,2,3\}$, we can bound $\mathbb{E}\left[\Big|m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) \Big|^3\right] = O(\frac{1}{h^{2d_A}})$. Hence \[ \sum_{i = 1}^n \mathbb{E}\left[\Big|\sqrt{nh^{d_A}}n^{-1} m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) \Big|^3\right] = O((nh^{d_A})^{-1/2}) = o(1) \]

Combining this with $s^2_n = V(a, a') + o(1)$ proves the Lyapunov condition. Hence, $$ \frac{1}{s_n} \sum_{i = 1}^n \left(\sqrt{\frac{h^{d_A}}{n}}m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a')) - \mu_i\right) \xrightarrow[]{d} \mathcal{N}(0, 1) $$

An application of Slutsky's theorem provides the desired result that $$\sqrt{nh^{d_A}} (\hat{\psi}^{MR}(a,a')- \psi_0(a,a') -B(a,a')) \xrightarrow[]{d} N(0, V(a,a'))$$

Proof for Multiple Robustness

Following a similar breakdown as that in Theorem 1, $\hat{\psi}^{MR}(a, a') - \psi(a, a')$ can be expanded as

Following the result in Theorem 1 on asymptotic normality by application of the Lyapunov CLT, the term $\sqrt{\frac{h^{d_A}}{n}}\sum^L_{\ell=1}\sum_{i\in I_\ell} \Big\{ m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\Big\} = O_p(1)$. Since $\frac{1}{\sqrt{nh^{d_A}}} = o_p(1)$, the following holds

\[\frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{\frac{h^{d_A}}{n}}\sum^L_{\ell=1}\sum_{i\in I_\ell} \Big\{ m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\Big\} = o_p(1). \]

The remainder of the proof demonstrates the remaining the remaining terms are $o_p(1)$, i.e. \[ \frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{\frac{h^{d_A}}{n}}\sum^L_{\ell=1}\sum_{i\in I_\ell}\Big\{m(O_i;\hat{\alpha}_\ell,\hat{\lambda}_\ell,\hat{\gamma}_\ell,\psi_0(a,a'))- m(O_i;\alpha,\lambda,\gamma,\psi_0(a,a'))\Big\}. \]

To see this, we first expand these terms identically as the proof of Theorem 1, and provide proofs for the convergence of the terms (CS1) - (CS6), (E1) - (E8) and (TR1) - (TR5) under the assumption that any two out of three nuisance models are correctly specified in the following sub-sections.

Proof for Terms (CS1)-(CS6)

All of these terms contain the product of two or more errors and can be treated similarly. We provide a detailed proof for (CS2), and a similar method can be followed for the rest of the terms.

For (CS2), write $\Delta_{i\ell}=K_h(A_i - a)\big[\hat{R}(M_i,X_i)-R(M_i,X_i)\big]\big[\hat{\gamma}(X_i,M_i,a)-\gamma(X_i,M_i,a)\big]\big[Y_i - \gamma(X_i,M_i,a)\big]$. Following Lemma (ref), it suffices to bound $\mathbb{E}\left[|\frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{\frac{h^{d_A}}{n}} \sum_{i \in I_\ell} \Delta_{i\ell} | \Big| O^c_{I_\ell}\right]$ as $o_p(1)$ in order to show that $$\frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{\frac{h^{d_A}}{n}} \sum_{i \in I_\ell}\Delta_{i\ell} = o_p(1).$$

First, from the triangle inequality, $\mathbb{E}\left[\left|\frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{\frac{h^{d_A}}{n}} \sum_{i \in I_\ell} \Delta_{i\ell} \right| \mid O^c_{I_\ell}\right] \leq \frac{1}{L}\mathbb{E}\left[\left|\Delta_{i\ell} \right| \mid O^c_{I_\ell}\right]$, and so it suffices to bound $\mathbb{E} \bigg[\big|\Delta_{i\ell} \big|\bigg| O^c_{I_\ell}\bigg]$.

{

} Next, Assumption 3.1 on the boundedness of $\gamma(X,M,a)$ and Assumption 3.3 on the boundedness of $var(Y_i|a,m,x)$, along with an application of Lemma (ref) on $f(a \mid M, X)$, we get {

} As long as either Assumption 4.2 or Assumption 4.3 hold, then combined with Assumption 3.2, (CS2) will be $o_p(1)$. A similar approach can be used to bound the remaining CS terms.

Proof for Terms (E1)-(E8)

Terms (E1)-(E8) are normalized terms of the form of a bias times a bounded quantity; they can all be treated similarly. We only provide the proof of the convergence in probability to zero for the term (E2). (E2) is given as

To prove this, we set $\hat{\Delta}_{i\ell}$ as (E2). By construction, $O^c_{I_\ell}$ and $O_i$ are independent, $i\in I_\ell$, and consequently $\mathbb{E}\left[\hat{\Delta}_{i\ell}|O^c_{I_\ell}\right]=0$ and $\mathbb{E}\left[\hat{\Delta}_{i\ell}\hat{\Delta}_{j\ell}|O^c_{I_\ell}\right]=0$ for $i,j \in I_\ell$ and all $a', a\in \mathcal{A}_0$. Next we note that

Where (a) follows from Assumption 3.1 on the boundedness of $f(a \mid M, X)$, along with Assumption 3.1 and Assumption 3.3 combined with the derivation provided below

Next, (b) follows from Assumption 2.4, and finally, (c) follows Assumption 3.2 along with Assumption 4.1. Then $\mathbb{E} \bigg[\left(\frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{h^{d_A}/n} \sum_{i\in I_\ell} \hat{\Delta}_{i\ell}\right)^2\bigg| O^c_{I_\ell}\bigg] =\frac{1}{n^2} \sum_{i\in I_\ell} \mathbb{E}\left[\hat{\Delta}^2_{i\ell}|O^c_{I_\ell}\right] = O(\frac{1}{n})\mathbb{E}\left[\hat{\Delta}^2_{i\ell}|O^c_{I_\ell}\right] = O_p(\frac{1}{nh^{d_A}}) = o_p(1)$.

Applying Lemma 1 to the above gives $\frac{1}{\sqrt{nh^{d_A}}} \times \sqrt{h^{d_A}/n} \sum^L_{l=1}\sum_{i\in I_\ell} \hat{\Delta}_{i\ell}\xrightarrow[]{P} 0$.

Proof for Terms (TR1)-(TR5)

The proofs of the convergence in probability to zero for the terms (TR1)-(TR5) follows a similar outline as Theorem 1, and we prove convergence on a case by case below.

Proof for Terms TR1 and TR2

Terms (TR1) and (TR2) are similar; we only provide the proof of the convergence in probability to zero for the term (TR2).

To bound TR2, first set $\hat{\Delta}_{i\ell}=K_h(A_i - a) \left[\hat{\lambda}(a, X_i) - \lambda(a,X_i) \right]R(M_i, X_i)\left[Y_i - \gamma(X_i,M_i,a)\right]$. Bounding (TR2) amounts to showing $\mathbb{E}[\hat{\Delta}_{i\ell}|O^c_{I_\ell}] = o_p(1)$.

where the equalities follow identically as in the proof of Theorem 1, and the final equality follows from $h \rightarrow 0$ along with the boundedness assumptions in Assumption 3.

Proof for TR3

For Term (TR3), we have

where (b) follows from the definition of $\lambda_{a^\prime}(X_i)$, $h\rightarrow 0$, Assumption 3 (boundedness of $\lambda$, $\hat{\eta}$ and $\eta$).

Proof For TR4

Demonstrating the bound for (TR4), we have

TR-4-1 can be written as

An application of Lemma (ref) to TR-4-1 gives

A similar approach applied to TR-4-2 gives

Now, the first terms of TR-4-1 and TR-4-2 cancel out with each other, with an identical proof to that used in the proof of Theorem 1.

Consequently only the remaining terms must be bounded. \[ \int \big(\hat{\gamma}(X_i,M_i,a)-\gamma(X_i,M_i,a)\big)\lambda(a', X_i) O(h^2)f(M, X) dM_i dX_i = o_p(1) \] The second term in TR-4-1 and TR-4-2 can be bounded from $h \rightarrow 0$, combined with the boundedness assumptions in Assumption 3.2.

Proof For TR5

Finally, for term (TR5), we note that

Where $(a)$ follows from an application of Lemma (ref), (b) follows from the definition of $\eta$, and (c) follows from an application of Cauchy-Schwartz combined with the consistency of $\hat{\lambda}$.

Proof for Consistency of $\widehat{V}(a, a')$

Recall that

To prove consistency of $\widehat{V}(a, a')$, we first prove propositions (I), (II) and (III), which together prove the desired result.

(I) $h^{d_A}n^{-1} \sum_{i\in I_\ell} m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a')) - V(a,a') = o_p(1)$

To simplify notation, denote $m(O_i;\alpha,\lambda,\gamma,\psi(a,a'))$ as $m_i(a,a')$. From the proof of Theorem 1, we have $h^{d_A} \mathbb{E}[m^2_i(a,a')] = V(a,a') + o_p(1)$.

We write

Then,

We only need to investigate the terms $\mathbb{E}(U_1^{c_1} U_2^{c_2}U_3^{c_3})$ for any $c_1 \ge 0$, $c_2 \ge 0$ and $c_3 \ge 0$ with $c_1 + c_2 + c_3 = 4$. To be specific, dropping the terms with power index being zero, we will be studying $\mathbb{E}(U_1^{c_1})$, $\mathbb{E}(U_2^{c_2})$, $\mathbb{E}(U_1^{c_1} U_2^{c_2})$, $\mathbb{E}(U_1^{c_1} U_3^{c_3})$, $\mathbb{E}(U_2^{c_2} U_3^{c_3})$, and $\mathbb{E}(U_1^{c_1} U_2^{c_2}U_3^{c_3})$ for positive $c_1, c_2,$ and $c_3$.

Combining all the terms, we obtain $\mathbb{E} (m^4_i) = O(h^{-3d_A})$. Then by Markov inequality, for any $\epsilon > 0$,

where the equality in the last row comes from $var( m^2_i) = O(\mathbb{E} (m^4_i)) = O(h^{-3d_A})$.

(II): $h^{d_A}|I_\ell|^{-1}\sum_{i\in I_\ell}\mathbb{E}[ m^2(O_i;\hat{\alpha}_\ell,\hat{\lambda}_\ell,\hat{\gamma}_\ell,\hat{\psi}_\ell(a,a')) - m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a')) \mid O^c_{I_\ell}] = o_p(1)$

For simplicity in notation, we ignore the subscripts $\ell$ below for nuisance parameters estimated from $O^c_{I_\ell}$. First, we analyze $h^{d_A}\mathbb{E}[m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a'))\mid O^c_{I_\ell}]$ as follows. We write

Denote $m(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a'))$ as $\hat{m}_i$. Then,

Combining all the six parts, we have $h^{d_A}\mathbb{E}[m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a'))]$ equal to

Next, $h^{d_A}\mathbb{E}[m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a'))]$ can be written as

Define $\int \tilde{k}(u)^2 du = R^2_{d_A}$,

Then $h^{d_A}\mathbb{E}[m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a'))]- h^{d_A}\mathbb{E}[m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a'))] = \omega_1 + \omega_2 + o_p(1).$ First, we focus on simplifying $\omega_2$, which equals

From expressing $\frac{1}{\hat{f}(a'|X)}\Big(\hat{\gamma}(X,M,a) - \hat{\eta}(a, a', X)\Big)$ as

there is

We show each of these terms are $o_p(1)$ as follows. Because $f^2(a'| X)$ is bounded away from $0$ based on Assumption 3 (ii) and the consistency of $\hat{\gamma}$ from Assumption 4(iii), there is $$\mathbb{E}\bigg\{ \frac{1}{f^2(a'|X)}\Big[\hat{\gamma}(X,M,a) - \gamma(X,M,a)\Big]^2\Big| O^c_{I_\ell} \Big\} = o_p(1).$$ Under a similar argument and from Assumption 4(iv), $$\mathbb{E}\bigg\{ \frac{1}{f^2(a'|X)}\Big[\eta(a, a', X) - \hat{\eta}(a, a', X)\Big]^2\Big| O^c_{I_\ell} \Big\} = o_p(1).$$ Based on the boundedness of nuisance estimators from Assumption 3(ii) and the consistency of $\hat{f}(a'|X)$ from Assumption 4(i), there is $$\mathbb{E}\bigg\{ \Big[\hat{\gamma}(X,M,a) - \hat{\eta}(a, a', X)\Big]^2\Big[\frac{1}{\hat{f}(a'|X)} - \frac{1}{f(a'|X)}\Big]^2\Big| O^c_{I_\ell} \Big\} = o_p(1).$$ Each of the remaining cross terms is a product of a term that is $o_p(1)$ from the estimator's consistency and a term that is bounded. Hence, we have $\omega_2 = o_p(1)$.

Next, we employ a similar derivation to simplify $\omega_1$.

Note that

Hence,

After further expansions, we can show that the squared terms contain a component that is bounded based on Assumption 3 and another component that is $o_p(1)$ from Assumption 4. The $(Y - \gamma(X,M,a))^2$ in some squared terms is integrated out as a bounded component due to $var(Y|X,M,a)$ being bounded as assumed in Assumption 3(3). For interaction terms, those containing $(Y - \gamma(X,M,a))$ equals zero because $\int (Y - \gamma(X,M,a))f(Y|X,M,a)dY = 0$. All of the interaction terms contain a bounded component and a $o_p(1)$ component. Consequently, $\omega_1 = o_p(1)$, leading to $h^{d_A}\mathbb{E}[m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a'))]- h^{d_A}\mathbb{E}[m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a'))] = o_p(1).$

(III) $h^{d_A}|I_\ell|^{-1} \sum_{i\in I_\ell} \Delta_{i\ell} = o_p(1)$, where $$ \Delta_{i\ell} = m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a')) - m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a')) - \mathbb{E}\Big\{m^2(O_i;\hat{\alpha},\hat{\lambda},\hat{\gamma},\hat{\psi}(a,a')) - m^2(O_i;\alpha,\lambda,\gamma,\psi(a,a'))\Big| O^c_{I_\ell}\Big\}. $$

By Lemma (ref), it suffices to bound $\mathbb{E}\left[\Big(h^{d_A}|I_\ell|^{-1} \sum_{i\in I_\ell} \Delta_{i\ell} \Big)^2 \Big| O^c_{I_\ell}\right] =h^{2d_A}|I_\ell|^{-1} \mathbb{E}\left[\Delta^2_{i\ell} \Big| O^c_{I_\ell}\right] $ as $o_p(1)$. Note that $\mathbb{E}[\Delta_{i\ell}] = 0$ and interaction terms are zero due to conditional independence. We start with analyzing $\mathbb{E}[\Delta^2_{i\ell} \mid O^c_{I_\ell}]$ as follows. For simplicity of notation, we adopt the notation definitions in parts (I) and (II), ignoring the subscripts $\ell$ for nuisance estimators. We have

From (II), we know that $\mathbb{E} \left(\hat{m}^2_{i} - m^2_i \mid O^c_{I_\ell}\right)^2 = o_p( h^{-2d_A})$. To bound $\mathbb{E}\bigg[\big(\hat{m}^2_{i} - m^2_i \big)^2\Big| O^c_{I_\ell}\bigg]$, by $\hat{m}_{i} = \hat{U}_1 + \hat{U}_2 + \hat{U}_3$ and $m_i = U_1 + U_2 + U_3$, we can rewrite the term as

where $\bar{c} = (c_1, \ldots, c_6)$ and $\mathcal{W}$ represents the possible combinations of $\bar{c}$ from the decomposition. We will prove that $\mathbb{E}\bigg\{\hat{U}^{c_1}_1\hat{U}^{c_2}_2\hat{U}^{c_3}_3 U^{c_4}_1U^{c_5}_2U^{c_6}_3 \Big|O^c_{I_\ell}\bigg\} = O(h^{-(c_1 + c_2 + c_4 + c_5 - 1)d_A})$. Note that $$\mathbb{E}\bigg\{\hat{U}^{c_1}_1\hat{U}^{c_2}_2\hat{U}^{c_3}_3 U^{c_4}_1U^{c_5}_2U^{c_6}_3 \Big| O^c_{I_\ell}\bigg\} = \iint \hat{U}^{c_1}_1\hat{U}^{c_2}_2\hat{U}^{c_3}_3 U^{c_4}_1U^{c_5}_2U^{c_6}_3 f(Y, A, M, X \mid O^c_{I_\ell})dYdMdAdX.$$ By the boundedness of nuisance parameters and their estimates (Assumption 3(ii)), the above term equals { $$ O\Big(\iint K_h(A - a)^{c_1 + c_4}K_h(A - a')^{c_2 + c_5}\Big|[Y - \hat{\gamma}(X,M,a)]^{c_{1}}[Y - \gamma(X,M,a)]^{c_{4}}\Big| f(Y, A, M, X \mid O^c_{I_\ell})dYdMdAdX \Big).$$} The possible combinations of $c_1, c_4$ in $\bar{c}$ are $\{(c_1, c_4): (1, 1), (2, 0), (0, 2), (2, 1), (1, 2), (3, 0), (0, 3)\}$. Similar to the derivation in part (I), we will prove that the rate is $O(h^{-(c_1 + c_2 + c_4 + c_5 - 1)d_A})$ case-by-case. For the terms with $c_1 = 0$, the boundedness of $\mathbb{E}[|Y - \gamma|^4 \mid X, M, A]$ from Assumption 7 provides the boundedness of lower moments by separately considering the regions on which $|Y-\gamma|^{c_4}$ is $\ge$ or $< 1$. Next, we prove for the remaining terms.

The remaining terms to bound are $\mathbb{E}[\hat{U}^2_1-U^2_1)^2\mid O^c_{I_\ell}]$, $ \mathbb{E}[(\hat{U}^2_2-U^2_2)^2 \mid O^c_{I_\ell}]$, and $ \mathbb{E}[(\hat{U}^2_3-U^2_3)^2 \mid O^c_{I_\ell}]$. First, $\mathbb{E}[(\hat{U}^2_3-U^2_3)^2 \mid O^c_{I_\ell}]$ can be bounded from Assumption 3(ii). Next, we demonstrate the boundedness of $\mathbb{E}[(\hat{U}^2_2-U^2_2)^2 \mid O^c_{I_\ell}]$; a similar derivation applies to $\mathbb{E}[\hat{U}^2_1-U^2_1)^2\mid O^c_{I_\ell}]$. To start with, we re-express the term

From Assumption 7,

Hence,

From the expansion of $\omega_2$ in proving term 6 of the part (II), we can express $\frac{1}{\hat{f}(a'|X)^2}\big[\hat{\gamma}(X,M,a) - \hat{\eta}(a, a', X)\big]^2 - \frac{1}{f(a'|X)^2}\big[\gamma(X,M,a) - \eta(a, a', X)\big]^2$ as a summation of 9 components, i.e.

For the multiplication of any two of the nine components chosen with replacement, the corresponding conditional expectation $\mathbb{E}(\cdot|O^c_{I_\ell})$ is a construct of a subcomponent that is $o_p(1)$ from the consistency of nuisance parameters multiplied by other subcomponents that are bounded from Assumption 3. As a consequence, we obtained that $\mathbb{E}[(\hat{U}^2_2-U^2_2)^2 \mid O^c_{I_\ell}] = o_p(h^{-3d_A})$. A similar argument can be used to prove $\mathbb{E}[\hat{U}^2_1-U^2_1)^2\mid O^c_{I_\ell}] = o_p(h^{-3d_A})$ by utilizing the boundedness of $\mathbb{E}[(Y - \gamma)^4 \mid X, M, A]$ from Assumption 7 (i).

Because $ c_1+c_2+c_4+c_5 \le 4$, $O(1) \le O(h^{-(c_1 + c_2 + c_4 + c_5 - 1)d_A}) \le O(h^{-3d_A})$. We conclude that $\mathbb{E}\left[\Delta^2_{i\ell} \Big| O^c_{I_\ell}\right] = O(h^{-3d_A})$ and $$\mathbb{E}\left[\Big(h^{d_A}|I_\ell|^{-1} \sum_{i\in I_\ell} \Delta_{i\ell} \Big)^2 \Big| O^c_{I_\ell}\right] =h^{2d_A}|I_\ell|^{-1} \mathbb{E}\left[\Delta^2_{i\ell} \Big| O^c_{I_\ell}\right] =O\left([nh^{d_A}]^{-1}\right) = o_p(1).$$

Consistency of Hajek-type Propensity Estimator in Cross Validation

Given a consistent estimator of propensity score at treatment value $a$ for person $i$ in cross validation fold $I_\ell$, $\hat{f}(a|X_i)$, we define the corresponding Hajek-type stabilized weighted propensity score as follows, $$\hat{f}(a|X_i) \times \frac{1}{|I_{-\ell}|}\sum^{|I_{-\ell}|}_{j\in I_{-\ell}} \frac{K_h (A_j-a)}{\hat{f}(a|X_j)}. $$

The goal here is to prove that $$\lim_{|I_{-\ell}| \rightarrow \infty}\hat{f}(a|X_i) \times \frac{1}{|I_{-\ell}|}\sum^{|I_{-\ell}|}_{j\in I_{-\ell}} \frac{K_h (A_j-a)}{\hat{f}(a|X_j)}\stackrel{p}{=} f(a|X_i). $$