Relaxed doubly robust estimation in causal inference

Abstract

Causal inference plays a crucial role in biomedical studies and social sciences. Over the years, researchers have devised various methods to facilitate causal inference, particularly in observational studies. Among these methods, the doubly robust estimator distinguishes itself through a remarkable feature: it retains its consistency even when only one of the two components – either the propensity score model or the outcome mean model – is correctly specified, rather than demanding correctness in both simultaneously. In this paper, we focus on scenarios where semiparametric models are employed for both the propensity score and the outcome mean. Semiparametric models offer a valuable blend of interpretability akin to parametric models and the adaptability characteristic of nonparametric models. In this context, achieving correct model specification involves both accurately specifying the unknown function and consistently estimating the unknown parameter. We introduce a novel concept: the relaxed doubly robust estimator. It operates in a manner reminiscent of the traditional doubly robust estimator but with a reduced requirement for double robustness. In essence, it only mandates the consistent estimate of the unknown parameter, without requiring the correct specification of the unknown function. This means that it only necessitates a partially correct model specification. We conduct a thorough analysis to establish the double robustness and semiparametric efficiency of our proposed estimator. Furthermore, we bolster our findings with comprehensive simulation studies to illustrate the practical implications of our approach.

Keywords:

• Causal inference
• doubly robust
• relaxed doubly robust
• semiparametric model
• semiparametric efficiency
• model specification

1. Introduction

Causal inference plays a pivotal role across a variety of scientific disciplines, encompassing both randomized experiments and observational studies. While randomized experiments are the ideal gold standard for establishing causal relationships regarding treatment effects, their implementation can be infeasible and often challenging, due to ethical or logistical constraints in many real-world situations.

Alternatively, one can utilize observational data to draw causal inferences, provided that appropriate corrections are made to mitigate the bias arising from confounding factors resulting from nonrandomized treatment assignment. Given that all confounders affecting the relationship between treatment and outcome are observed, researchers have devised a wide array of techniques to effectively control for confounding, such as the weighting and matching methods based on propensity score (Rosenbaum & Rubin, 1983), outcome regression methods, and doubly robust (DR) methods based on the idea of augmenting the inverse probability weighting estimator (Bang & Robins, 2005; Lunceford & Davidian, 2004; Robins et al., 1994, 1995; Rotnitzky et al., 1998; Rotnitzky & Vansteelandt, 2014; Scharfstein et al., 1999). These methods have been well written into textbooks, such as Imbens and Rubin (2015).

Among all these methods, the doubly robust method has been popular because its consistency relies on the correct specification of either the propensity score model or the outcome mean model, not necessarily both. Further, if both propensity score model and outcome mean model are correctly specified, the doubly robust estimator becomes semiparametrically efficient (Bickel et al., 1993; Tsiatis, 2006). Various semiparametric models such as the single-index model have been popularly used in modelling the propensity score or the outcome mean, because semiparametric models typically possess the interpretability of parametric models and the flexibility of nonparametric models. Nonetheless, the advantageous property of double robustness may be compromised when semiparametric models are employed for either the propensity score or the outcome mean. In such cases, correct model specification necessitates the correct specification of the unknown function as well as the consistent estimation of the unknown parameter – two critical elements within the semiparametric model. In particular, the consistent estimation of the unknown parameter might not need the correct specification of the unknown function.

In this paper, we introduce a novel approach to doubly robust estimation specifically tailored for semiparametric models used in both the propensity score and the outcome mean. As mentioned earlier, the traditional doubly robust estimation hinges on the correct specification of the unknown function and the consistent estimation of the unknown parameter – two essential components that need to be accurately identified within either the propensity score or the outcome mean. In contrast, our newly proposed doubly robust estimation only demands the consistent estimation of the unknown parameter, reducing the requirement to a single piece of unknown information within the semiparametric model. Hence, we refer to our innovative method as the ‘relaxed doubly robust’ (RDR) estimation. Furthermore, when both parameters (one within the propensity score model and the other within the outcome mean model) are consistently estimated, our RDR estimator also attains the semiparametric efficiency bound.

The remaining of the paper is structured as follows. In Section 2, we first review the basics of the traditional DR estimation and point out its limitation when a semiparametric propensity score model or a semiparametric outcome mean model needs to be correctly specified. We then present the proposed RDR estimation in Section 3, that includes both algorithm in Section 3.1 and theory in Section 3.2. In Section 4, we conduct comprehensive simulation studies to illustrate our proposed estimator as well as its comparison with the traditional DR estimation. We conclude the paper with some discussions in Section 5.

2. Review of traditional DR estimation

We adopt the potential outcome framework (Neyman, 1923; Rubin, 1974) to briefly present the traditional DR estimation. We denote $\mathbf{X}$ a p-dimensional pre-treatment covariate. Suppose the treatment variable T is binary, $T\in \{0,1\}$, in that 1 stands for treatment and 0 control. For each level of the treatment T = t, under the standard stable unit treatment value assumption (SUTVA) (Rubin, 1980), we assume that there exists the potential outcome $Y_t$, representing the outcome had the subject, possibly contrary to the fact, been given the treatment t. We denote the observed outcome as Y that has the decomposition $Y=T{Y}_1+(1-T)Y_0$, based on the consistency assumption. We have independent and identically distributed data $\left\{\right(t_i,y_i,{\mathbf{x}}_i\left)\right\}$, $i=1,\dots ,n$, that are realizations of the triplet $(T,Y,\mathbf{X})$.

Throughout, we focus on estimating the average treatment effect (ATE), defined as $\mathrm{\Delta }=E\left(Y_1\right)-E\left(Y_0\right).$The ATE is one of the most popularly used causal estimands in a variety of scientific applications, and could generate important policy implications. The fundamental difficulty of estimating the ATE in causal inference is that one may observe at most one of $Y_1$ and $Y_0$ for every subject. We adopt the following two standard assumptions that are widely used in the causal inference literature. The first assumption is about the ignorability.

Assumption 2.1

Ignorability

We assume $(Y_1,Y_0)⫫T\mid \mathbf{X}$, which stands for conditional independence.

Assumption 2.1 requires $\mathbf{X}$ to include all potential confounders relevant to both treatment and outcome to be fully identified and completely observed. It indicates, the treatment assignment is independent of the potential outcomes $Y_1$ and $Y_0$ after conditioning on all pre-treatment confounders $\mathbf{X}$. Accordingly, we denote the propensity score model as $\pi \left(\mathbf{X}\right)\equiv \text{pr}(T=1\mid Y_1,Y_0,\mathbf{X})=\text{pr}(T=1\mid \mathbf{X}).$We also denote the outcome mean models as $\begin{array}{rl}{\mu }_1\left(\mathbf{X}\right)& \equiv E(Y_1\mid \mathbf{X})=E(Y_1\mid \mathbf{X},T=1),\text{ and }\\ {\mu }_0\left(\mathbf{X}\right)& \equiv E(Y_0\mid \mathbf{X})=E(Y_0\mid \mathbf{X},T=0).\end{array}$The second assumption pertains to the overlapping condition of the covariate distributions between the treated group and the controls.

Assumption 2.2

Overlap

There exist constants $c_1$ and $c_2$ such that $0<c_1\le \pi \left(\mathbf{X}\right)\le c_2<1$ almost surely.

Assumption 2.2 means that the covariates' distributions between the two groups are roughly somewhat similar to each other. If this assumption is not satisfied at some value of $\mathbf{X}$, the subjects with this value can only be in the treatment group or the control group, which would lead to the extrapolation at this value that makes the inference about the ATE inappropriate.

Under Assumptions 2.1 and 2.2, the ATE can be identified and then estimated via a variety of well known methods (Imbens, 2004; Rosenbaum, 2002). These methods would need either an estimate of the propensity score model, say, $\hat{\pi }\left(\mathbf{x}\right)$, that correctly specifies $\pi \left(\mathbf{x}\right)$, or an estimate of the outcome mean models, say, ${\hat{\mu }}_1\left(\mathbf{x}\right)$ and ${\hat{\mu }}_0\left(\mathbf{x}\right)$, that respectively specifies ${\mu }_1\left(\mathbf{x}\right)$ and ${\mu }_0\left(\mathbf{x}\right)$ correctly, or both of them. These methods would also induce different estimation efficiencies. In particular, the efficient influence function (EIF) (Tsiatis, 2006) for estimating the ATE Δ is $\varphi =\frac{t}{\pi \left(\mathbf{x}\right)}\{y_1-{\mu }_1(\mathbf{x}\left)\right\}-\frac{1-t}{1-\pi \left(\mathbf{x}\right)}\{y_0-{\mu }_0(\mathbf{x}\left)\right\}+{\mu }_1\left(\mathbf{x}\right)-{\mu }_0\left(\mathbf{x}\right)-\mathrm{\Delta }.$Accordingly, the traditional DR estimator is ${\hat{\mathrm{\Delta }}}_{\mathrm{dr}}=\frac{1}{n}\sum _{i=1}^n\left[\frac{t_i}{\hat{\pi }\left({\mathbf{x}}_i\right)}\right\{y_{1i}-{\hat{\mu }}_1\left({\mathbf{x}}_i\right)\}-\frac{1-t_i}{1-\hat{\pi }\left({\mathbf{x}}_i\right)}\{y_{0i}-{\hat{\mu }}_0\left({\mathbf{x}}_i\right)\}+{\hat{\mu }}_1({\mathbf{x}}_i)-{\hat{\mu }}_0({\mathbf{x}}_i\left)\right].$It is called DR since it only needs that $\hat{\pi }\left(\mathbf{x}\right)$ correctly specifies $\pi \left(\mathbf{x}\right)$ or that ${\hat{\mu }}_t\left(\mathbf{x}\right)$ correctly specifies ${\mu }_t\left(\mathbf{x}\right)$ for t = 1, 0, but not necessarily both. If all of these nuisance functions are correctly specified, the asymptotic variance of this estimator achieves the semiparametric efficiency bound $E\left({\varphi }^2\right)$ and this traditional DR estimator is called semiparametrically efficient. Among all regular and asymptotically linear estimators, $E\left({\varphi }^2\right)$ is the minimum possible asymptotic variance that one can attain.

2.1. Limitation of traditional DR estimation

In reality, practitioners may choose to model the nuisance functions $\pi \left(\mathbf{x}\right)$, ${\mu }_1\left(\mathbf{x}\right)$ and ${\mu }_0\left(\mathbf{x}\right)$ differently based on their own preferences. For example, one may impose parametric restrictions, say, $\pi (\mathbf{x},\mathit{\alpha })$, ${\mu }_1(\mathbf{x},\text{β})$ and ${\mu }_0(\mathbf{x},\mathit{\gamma })$, where only the parameters $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ are unknown. On the other hand, one may choose to impose no parametric restrictions and regard all of them $\pi \left(\mathbf{x}\right)$, ${\mu }_1\left(\mathbf{x}\right)$ and ${\mu }_0\left(\mathbf{x}\right)$ as nonparametric functions.

As another alternative, semiparametric models such as single index models generally possess the interpretability of parametric models as well as the flexibility of nonparametric models, and thus have been widely used in practice to model either the propensity score or the outcome mean or the both. Throughout, we consider single index models $\pi \left({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}\right)$, ${\mu }_1\left({\text{β}}^{\mathrm{\top }}\mathbf{x}\right)$ and ${\mu }_0\left({\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}\right)$, where $\pi (\cdot )$, ${\mu }_1(\cdot )$ and ${\mu }_0(\cdot )$ are unknown functions, $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ are unknown p-dimensional vectors whose first elements are all one's. Under such a scenario, the traditional DR estimator becomes ${\hat{\mathrm{\Delta }}}_{\mathrm{dr}}=\frac{1}{n}\sum _{i=1}^n\left[\frac{t_i}{\hat{\pi }\left({\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i\right)}\right\{y_{1i}-{\hat{\mu }}_1\left({\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i\right)\}-\frac{1-t_i}{1-\hat{\pi }\left({\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i\right)}\{y_{0i}-{\hat{\mu }}_0\left({\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i\right)\}+{\hat{\mu }}_1({\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i)-{\hat{\mu }}_0({\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i\left)\right].$Clearly, the double robustness then refers to either the correct specification of $\pi (\cdot )$ as well as the consistent estimate of $\mathit{\alpha }$, or the correct specifications of ${\mu }_1(\cdot )$ and ${\mu }_0(\cdot )$ as well as the consistent estimates of $\text{β}$ and $\mathit{\gamma }$. This might not be appealing for practical use since the correct model specification includes both correct specification of functional forms and consistent estimate of unknown parameters.

3. Proposed RDR estimation

We propose a relaxed doubly robust (RDR) estimator for estimating ATE under the scenario that both propensity score and outcome mean are single index models. We term it RDR since its double robustness refers to the consistent estimate of the parameter $\mathit{\alpha }$ or of the parameters $\text{β}$ and $\mathit{\gamma }$ only. Further, when all of the parameters $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ are consistently estimated, the RDR estimator would also be semiparametrically efficient. The correct specification of functional forms $\pi (\cdot )$, ${\mu }_1(\cdot )$ and ${\mu }_0(\cdot )$ would never be required.

Our key idea is to rewrite the definition of Δ as $\mathrm{\Delta }=E\left\{E\right(Y_1\mid {\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x},{\text{β}}^{\mathrm{\top }}\mathbf{x}\left)\right\}-E\left\{E\right(Y_0\mid {\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x},{\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}\left)\right\},$and then propose the RDR estimator as ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}\equiv {\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}(\hat{\mathit{\alpha }},\widehat{\text{β}},\hat{\mathit{\gamma }})=\frac{1}{n}\sum _{i=1}^n\left\{\hat{E}\right(Y_1\mid {\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i)-\hat{E}(Y_0\mid {\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i\left)\right\},$where the estimates of parameters $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ can be obtained through some off-the-shelf dimension reduction methods, and the estimate of the conditional expectation $\hat{E}(\cdot \mid \cdot )$ can be obtained from the standard kernel regression (Cheng, 1994; Hu et al., 2012), which will be detailed below.

3.1. Algorithm

The algorithm entails the following six steps.

1. Estimate the coefficients $\mathit{\alpha }$, $\text{β}$, and $\mathit{\gamma }$ from the regression models of regressing T on $\mathbf{X}$, regressing $Y_1$ on $\mathbf{X}$ given T = 1, regressing $Y_0$ on $\mathbf{X}$ given T = 0, using some dimension reduction technique such as the sliced inverse regression (SIR). These estimates are denoted by $\hat{\mathit{\alpha }}$, $\widehat{\text{β}}$, and $\hat{\mathit{\gamma }}$, respectively.

2. Compute ${\hat{\mathbf{S}}}_i=({\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i)$ and ${\hat{\mathbf{M}}}_i=({\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i)$ for $i=1,2,\dots ,n$. Now the original data $(\mathbf{X},T,Y)$ has been transformed to $(\hat{\mathbf{S}},T,Y)$ and $(\hat{\mathbf{M}},T,Y)$.

3. For $i=1,2,\dots ,n$, estimate $E\left(Y_1|{\hat{\mathbf{S}}}_i\right)$ by a Nadaraya-Watson kernel regression of $Y_1$ on $\hat{\mathbf{S}}$, $\hat{E}\left(Y_1|{\hat{\mathbf{S}}}_i\right)=\hat{E}(Y_1|{\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i)=\sum _{j=1}^n\frac{T_j{K}_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j)}{\sum _{k=1}^n{T}_k{K}_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_k)}{Y}_{1j},$where $K(\cdot )$ is a bivariate kernel function, $K_H\left(\mathit{u}\right)=\frac{K\left(H^{-1}\mathit{u}\right)}{det\left(H\right)h_n}$ with $\mathit{u}=(u_1,u_2{)}^{\mathrm{\top }}$ and $H={\hat{\mathrm{\Sigma }}}_1^{1/2}$ as the square root of the 2-by-2 sample covariance matrix of $\mathbf{S}$. Here $h_n$ is the bandwidth, and we follow Hu et al. (2012) to choose the optimal bandwidth $h_n\sim n^{-\frac{1}{3}}$.

4. Similarly, for $i=1,2,\dots ,n$, estimate $E\left(Y_0|{\hat{\mathbf{M}}}_i\right)$ by a Nadaraya-Watson kernel regression of $Y_0$ on $\hat{\mathbf{M}}$, $\hat{E}\left(Y_0|{\hat{\mathbf{M}}}_i\right)=\hat{E}(Y_0|{\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i)=\sum _{j=1}^n\frac{(1-T_j)K_H({\hat{\mathbf{M}}}_i-{\hat{\mathbf{M}}}_j)}{\sum _{k=1}^n(1-T_k)K_H({\hat{\mathbf{M}}}_i-{\hat{\mathbf{M}}}_k)}{Y}_{0j},$where $K(\cdot )$ is a bivariate kernel function, $K_H\left(\mathit{u}\right)=\frac{K\left(H^{-1}\mathit{u}\right)}{det\left(H\right)h_n}$ with $\mathit{u}=(u_1,u_2{)}^{\mathrm{\top }}$ and $H={\hat{\mathrm{\Sigma }}}_0^{1/2}$ as the square root of the 2-by-2 sample covariance matrix of $\mathbf{M}$.

5. Estimate $E\left(Y_1\right)$ and $E\left(Y_0\right)$ by the sample average of $\hat{E}(Y_1|{\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i)$ and $\hat{E}(Y_0|{\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i)$. That is, ${\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)=\frac{1}{n}\sum _{i=1}^n\hat{E}(Y_1|{\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i),\text{ and }{\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)=\frac{1}{n}\sum _{i=1}^n\hat{E}(Y_0|{\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i,{\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i).$

6. The proposed RDR estimator is ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}={\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)-{\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right).$

3.2. Theory

To facilitate the presentation of the theory, we define the following conditions.

Condition (i): $\mathit{\alpha }$ is consistently estimated and its estimate $\hat{\mathit{\alpha }}$ satisfies $\hat{\mathit{\alpha }}-\mathit{\alpha }=O_p\left(n^{-1/2}\right)$.

Condition (ii): $\text{β}$ and $\mathit{\gamma }$ are consistently estimated and their estimates $\widehat{\text{β}}$ and $\hat{\mathit{\gamma }}$ satisfy $\widehat{\text{β}}-\text{β}=O_p\left(n^{-1/2}\right)$ and $\hat{\mathit{\gamma }}-\mathit{\gamma }=O_p\left(n^{-1/2}\right)$.

We will defer some discussions on how to consistently estimate these parameters $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ to Section 5. In general, interested readers should refer to the literature on semiparametric modelling techniques, including the single-index models and the dimension reduction techniques.

Other than Assumptions 2.1 and 2.2, our theory also requires the following regularity conditions.

1. The kernel function in $K_H(\cdot )$ satisfies $\int \mathit{u}K\left(\mathit{u}\right)\mathrm{d}\mathit{u}=\mathbf{0}$, $\int \mathit{u}{\mathit{u}}^{\mathrm{\top }}K\left(\mathit{u}\right)\mathrm{d}\mathit{u}={\gamma }_K\mathbf{I}$ for ${\gamma }_K\le \mathrm{\infty }$, and $\int K^2\left(\mathit{u}\right)\mathrm{d}\mathit{u}={\tau }_K$ for ${\tau }_K\le \mathrm{\infty }$.

2. The density of $\mathbf{X}$ is bounded away from 0.

3. As $n\to \mathrm{\infty }$, the bandwidth matrix $H\to 0$ such that $\sqrt{\mathrm{tr}\left(H{H}^{\mathrm{\top }}\right)}\to 0$ and $n\det \left(H\right)\to \mathrm{\infty }$, where $\mathrm{tr}$ is the trace of a matrix.

4. $E\{|\mathbf{X}-E(\mathbf{X}\left){|}^k\right\}<\mathrm{\infty }$ for some k>6.

Theorem 3.1

Relaxed Double Robustness

Given Assumptions 2.1 and 2.2 and the above regularity conditions, if either Condition (i) or Condition (ii) is satisfied, then the proposed estimator ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$ is asymptotically consistent; that is, ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}⟹p\mathrm{\Delta }$ as $n\to \mathrm{\infty }$.

Theorem 3.2

Semiparametric Efficiency

Given Assumptions 2.1 and 2.2 and the above regularity conditions, if both Condition (i) and Condition (ii) are satisfied, then the proposed estimator ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$ achieves the semiparametric efficiency; that is $\sqrt{n}({\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}-\mathrm{\Delta })⟹dN\{0,E({\varphi }^2\left)\right\},$where $\varphi =\frac{t}{\pi \left({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}\right)}\{y_1-{\mu }_1({\text{β}}^{\mathrm{\top }}\mathbf{x}\left)\right\}-\frac{1-t}{1-\pi \left({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}\right)}\{y_0-{\mu }_0({\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}\left)\right\}+{\mu }_1\left({\text{β}}^{\mathrm{\top }}\mathbf{x}\right)-{\mu }_0\left({\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}\right)-\mathrm{\Delta }$ is the EIF for estimating Δ under the scenario that both propensity score and outcome mean are single index models.

The proofs of Theorems 3.1 and 3.2 are contained in the Appendix.

4. Simulation studies

In this section, we conduct comprehensive simulation studies to illustrate the proposed RDR estimator, as well as its comparison with the oracle estimator (Oracle, using all simulated $Y_1$ and $Y_0$ data), the regression based estimator (REG), the IPW estimator (IPW) and the traditional DR estimator (DR).

We consider the following three settings pertaining to the misspecification of the nuisance functional forms.

• Setting 1: neither $\pi (\cdot )$ nor ${\mu }_t(\cdot )$, t = 1, 0 is correctly specified;

• Setting 2: only $\pi (\cdot )$ is correctly specified;

• Setting 3: only ${\mu }_t(\cdot )$, t = 1, 0 is correctly specified.

We expect that, in every setting, the RDR would perform similarly as the Oracle. Because of the constraints on the correct model specification, we expect that all three methods REG, IPW and DR would be biased in Setting 1, REG would be biased in Setting 2, and IPW would be biased in Setting 3.

In each setting, we consider four cases with sample size n = 1000 and n = 2000 blended with independent and correlated covariates: $\begin{array}{rl}& \mathbf{X}\sim N(\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\left[\begin{array}{cccc}0.5& & & \\ & 0.5& & \\ & & 0.5& \\ & & & 0.5\end{array}\right]),\text{ and}\\ & \mathbf{X}\sim N(\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\left[\begin{array}{cccc}0.5& {0.5}^2& {0.5}^3& {0.5}^4\\ {0.5}^2& 0.5& {0.5}^2& {0.5}^3\\ {0.5}^3& {0.5}^2& 0.5& {0.5}^2\\ {0.5}^4& {0.5}^3& {0.5}^2& 0.5\end{array}\right]).\end{array}$For Setting 1, we consider the nuisance models $\begin{array}{rl}& \pi \left(\mathbf{x}\right)=\text{probit}\left(2\right({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}{)}^2-1.3|{\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}|),\\ & {\mu }_1\left(\mathbf{x}\right)=\frac{{\text{β}}^{\mathrm{\top }}\mathbf{x}}{0.5+({\text{β}}^{\mathrm{\top }}\mathbf{x}+1.5{)}^2},\\ & {\mu }_0\left(\mathbf{x}\right)=\frac{{\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}}{0.5+({\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}+1{)}^2},\end{array}$where $\mathit{\alpha }=(1,-0.5,1,-0.5{)}^{\mathrm{\top }}$, $\text{β}=(-1,-1,-0.5,1{)}^{\mathrm{\top }}$ and $\mathit{\gamma }=(0.5,1,1,-1{)}^{\mathrm{\top }}$. For Setting 2, we consider $\pi \left(\mathbf{x}\right)=\text{logit}\left({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}\right),$and the same outcome means models as Setting 1. For Setting 3, we consider $\begin{array}{rl}& \pi \left(\mathbf{x}\right)=\text{probit}\left(\right({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}{)}^2-|{\mathit{\alpha }}^{\mathrm{\top }}\mathbf{x}|),\\ & {\mu }_1\left(\mathbf{x}\right)={\text{β}}^{\mathrm{\top }}\mathbf{x},\\ & {\mu }_0\left(\mathbf{x}\right)={\mathit{\gamma }}^{\mathrm{\top }}\mathbf{x}.\end{array}$For every specific case, we use the sliced inverse regression (SIR) method to obtain $\hat{\mathit{\alpha }}$, $\widehat{\text{β}}$ and $\hat{\mathit{\gamma }}$ for the proposed RDR estimator. For the traditional REG, IPW and DR methods, we simply fit logistic regression for $\pi \left(\mathbf{x}\right)$ and linear regression for ${\mu }_t\left(\mathbf{x}\right)$, t = 1, 0. We conduct 500 Monte Carlo simulations and report the empirical bias and empirical mean squared error in Tables 1, 2 and 3, respectively.

Table 1. Simulation results for setting 1.

			Bias($\times {10}^2$)					Mean Squared Error($\times {10}^2$)
Scenario	Sample Size	Estimand	Oracle	REG	IPW	DR	RDR	Oracle	REG	IPW	DR	RDR
		$E\left(Y_1\right)$	0.07	−2.83	−2.89	−2.88	0.36	0.05	0.15	0.15	0.15	0.05
		$E\left(Y_0\right)$	−0.01	2.82	2.81	2.79	0.22	0.05	0.15	0.15	0.15	0.06
	1000	ATE	0.08	−5.65	−5.70	−5.67	0.14	0.14	0.53	0.52	0.53	0.15
		$E\left(Y_1\right)$	0.04	−2.82	−2.86	−2.84	0.22	0.02	0.11	0.12	0.12	0.03
		$E\left(Y_0\right)$	0.08	2.86	2.86	2.84	0.20	0.03	0.12	0.12	0.12	0.03
independent	2000	ATE	−0.04	−5.68	−5.72	−5.68	0.02	0.07	0.42	0.42	0.42	0.08
		$E\left(Y_1\right)$	−0.11	−4.68	−4.78	−4.74	0.10	0.05	0.31	0.31	0.31	0.05
		$E\left(Y_0\right)$	−0.19	2.76	2.78	2.73	−0.14	0.05	0.15	0.15	0.15	0.05
	1000	ATE	0.09	−7.44	−7.56	−7.47	0.24	0.13	0.78	0.78	0.79	0.14
		$E\left(Y_1\right)$	−0.14	−4.53	−4.58	−4.57	0.05	0.03	0.25	0.25	0.25	0.03
		$E\left(Y_0\right)$	0.00	3.08	3.09	3.06	0.04	0.03	0.13	0.13	0.13	0.03
correlated	2000	ATE	−0.14	−7.61	−7.67	−7.63	0.01	0.08	0.69	0.69	0.69	0.08

Table 2. Simulation results for setting 2.

			Bias($\times {10}^2$)					Mean Squared Error($\times {10}^2$)
Scenario	Sample Size	Estimand	Oracle	REG	IPW	AIPW	RDR	Oracle	REG	IPW	AIPW	RDR
		$E\left(Y_1\right)$	0.07	1.64	0.07	0.14	0.23	0.05	0.09	0.06	0.07	0.05
		$E\left(Y_0\right)$	−0.01	1.10	0.05	0.06	0.03	0.05	0.09	0.07	0.08	0.06
	1000	ATE	0.08	0.55	0.02	0.08	0.20	0.14	0.21	0.17	0.22	0.15
		$E\left(Y_1\right)$	0.04	1.50	0.06	0.05	0.15	0.03	0.06	0.03	0.03	0.03
		$E\left(Y_0\right)$	0.08	1.18	0.13	0.18	0.10	0.03	0.05	0.03	0.04	0.03
independent	2000	ATE	−0.04	0.32	−0.07	−0.13	0.05	0.07	0.10	0.09	0.10	0.08
		$E\left(Y_1\right)$	−0.11	1.84	−0.08	0.02	0.03	0.06	0.11	0.07	0.08	0.06
		$E\left(Y_0\right)$	−0.07	1.12	−0.06	−0.08	−0.03	0.06	0.09	0.06	0.08	0.06
	1000	ATE	−0.04	0.72	−0.02	0.10	0.06	0.16	0.22	0.19	0.23	0.16
		$E\left(Y_1\right)$	−0.14	1.68	−0.11	−0.10	−0.03	0.03	0.07	0.04	0.04	0.03
		$E\left(Y_0\right)$	0.00	1.22	0.04	0.07	0.02	0.03	0.05	0.03	0.04	0.03
correlated	2000	ATE	−0.14	0.46	−0.15	−0.18	−0.06	0.08	0.12	0.10	0.11	0.08

Table 3. Simulation results for setting 3.

			Bias($\times {10}^2$)					Mean Squared Error($\times {10}^2$)
Scenario	Sample Size	Estimand	Oracle	REG	IPW	AIPW	RDR	Oracle	REG	IPW	AIPW	RDR
		$E\left(Y_1\right)$	0.09	0.10	2.52	0.09	−0.07	0.05	0.05	0.12	0.05	0.05
		$E\left(Y_0\right)$	−0.01	0.01	−4.85	0.01	−1.09	0.05	0.05	0.29	0.05	0.06
	1000	ATE	0.10	0.09	7.37	0.09	1.02	0.19	0.19	0.74	0.19	0.20
		$E\left(Y_1\right)$	−0.04	−0.04	2.39	−0.04	−0.13	0.03	0.03	0.09	0.03	0.03
		$E\left(Y_0\right)$	0.04	0.04	−4.81	0.04	−0.84	0.03	0.03	0.26	0.03	0.03
independent	2000	ATE	−0.08	−0.08	7.20	−0.08	0.71	0.10	0.10	0.63	0.10	0.11
		$E\left(Y_1\right)$	0.01	0.02	2.13	0.02	−0.23	0.06	0.06	0.11	0.06	0.06
		$E\left(Y_0\right)$	0.02	0.04	−2.95	0.04	−0.74	0.05	0.05	0.14	0.05	0.06
	1000	ATE	−0.01	−0.02	5.08	−0.02	0.51	0.22	0.21	0.47	0.21	0.22
		$E\left(Y_1\right)$	−0.10	−0.09	2.03	−0.09	−0.21	0.03	0.03	0.08	0.03	0.03
		$E\left(Y_0\right)$	0.07	0.07	−2.88	0.07	−0.52	0.03	0.03	0.11	0.03	0.03
correlated	2000	ATE	−0.16	−0.15	4.91	−0.15	0.31	0.11	0.11	0.36	0.11	0.11

The results of the simulation Setting 1 are shown in Table 1. We could see that when both models are misspecified, the relaxed doubly robust estimator has the smallest biases for $E\left(Y_1\right)$, $E\left(Y_0\right)$, and Δ among all these estimators. However, the regression, IPW, and DR estimators are significantly biased for these estimands. Also, the mean squared errors of the RDR estimator are quite close to that of the oracle estimator under different scenarios, but other estimators are not.

From Table 2, we could check the numerical results for our simulation Setting 2. When the regression models are misspecified, the regression estimator has significantly larger biases than other estimators while the IPW and DR estimators have small biases, which echoes with our expectation. Our RDR estimator achieves relatively smaller biases and MSEs than the regression estimators, and its MSEs are quite close to that of the oracle estimator, which means that the performance of the RDR estimator is quite good under this condition.

Finally, we can find the similar results in Table 3 for simulation Setting 3. When regression models are correctly specified, the regression and DR estimators both achieve relatively small biases and mean square errors. The IPW estimator has large biases and MSEs due to the model misspecification. Our RDR estimator achieves small biases and mean square errors that are close to the oracle estimator, which means that it has good performance when the propensity score model is misspecified.

5. Discussion

In this paper, we propose a relaxed doubly robust (RDR) estimator for estimating the average treatment effect under the scenario that both propensity score and outcome mean are single index models. In contrast to the traditional doubly robust (DR) estimation, our new estimator only necessitates partial correctness in model specification. It only requires consistent estimation of the unknown parameter in either the propensity score model or the outcome mean model, or both, to achieve double robustness and semiparametric efficiency, respectively. The correct specification of unknown functions in these single-index models is no longer a need. This is why we term it ‘relaxed double robustness’.

To implement the proposed RDR estimator, one needs to consistently estimate $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ in a semiparametric modelling framework. This is a stand-alone problem that was fortunately well studied in the literature. One can either go with standard semiparametric methods in single-index models (Ichimura, 1993; Klein & Spady, 1993) or choose to estimate these unknown parameters in the context of dimension reduction (Cook, 2009; Ma & Zhu, 2012, 2013). In our numerical studies, we choose the standard sliced inverse regression (Li, 1991) in dimension reduction to estimate these parameters $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$.

Finally, as the Associate Editor suggested, we point out that different relaxed doubly robust estimations also exist. For instance, after estimating the unknown parameters $\mathit{\alpha }$, $\text{β}$ and $\mathit{\gamma }$ same as what we propose, one may opt to using some one-dimensional nonparametric regression method such as kernels to estimate $\hat{\pi }(\cdot )$, ${\hat{\mu }}_1(\cdot )$ and ${\hat{\mu }}_0(\cdot )$ respectively. It will ultimately lead to the estimator $\frac{1}{n}\sum _{i=1}^n\left[\frac{t_i}{\hat{\pi }\left({\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i\right)}\right\{y_{1i}-{\hat{\mu }}_1\left({\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i\right)\}-\frac{1-t_i}{1-\hat{\pi }\left({\hat{\mathit{\alpha }}}^{\mathrm{\top }}{\mathbf{x}}_i\right)}\{y_{0i}-{\hat{\mu }}_0\left({\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i\right)\}+{\hat{\mu }}_1({\widehat{\text{β}}}^{\mathrm{\top }}{\mathbf{x}}_i)-{\hat{\mu }}_0({\hat{\mathit{\gamma }}}^{\mathrm{\top }}{\mathbf{x}}_i\left)\right],$which is also relaxed doubly robust. This estimator and the proposed estimator have different ideas. The proposed RDR estimator achieves the double robustness by fitting two augmented outcome mean models $E(Y_1\mid {\hat{\mathit{\alpha }}}^{\mathrm{\top }}\mathbf{x},{\widehat{\text{β}}}^{\mathrm{\top }}\mathbf{x})$ and $E(Y_0\mid {\hat{\mathit{\alpha }}}^{\mathrm{\top }}\mathbf{x},{\hat{\mathit{\gamma }}}^{\mathrm{\top }}\mathbf{x})$; see the details in steps (iii)-(v) of the algorithm in Section 3.1. To the contrary, the aforementioned estimator achieves the double robustness because it was directly implemented according to the efficient influence function.

Acknowledgements

The authors would like to thank the Editor, an Associate Editor, and one reviewer for their insightful comments which have helped improve the manuscript substantially.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The research is supported in part by US National Science Foundation [Grant Numbers DMS 1953526, 2122074, 2310942], US National Institutes of Health [Grant Number R01DC021431] and the American Family Funding Initiative of UW-Madison.

Appendix

Proof

Proof of Theorem 3.1

First, we can prove that $\hat{E}\left(Y_1|\mathbf{S}\right)$, i.e., $\hat{E}(Y_1|{\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X},{\text{β}}^{\mathrm{\top }}\mathbf{X})$, is consistent for $E\left(Y_1|\mathbf{S}\right)$ if either of the following conditions is satisfied.

• Condition (a): $\mathit{\alpha }$ is correctly specified, where $\mathit{\alpha }$ is the true coefficient for the propensity score model;

• Condition (b): $\text{β}$ is correctly specified, where $\text{β}$ is the true coefficient for the linear model of $Y_1$ on $\mathbf{X}$.

If $\mathit{\alpha }$ is correctly specified, the propensity score $\pi \left(\mathbf{X}\right)$ should be some function of ${\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X}$ with probability 1, which means that $\pi \left(\mathbf{X}\right)$ could also be written as a function of $\mathbf{S}=({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X},{\text{β}}^{\mathrm{\top }}\mathbf{X})$. Denote this function by f, and we have $\pi \left(\mathbf{X}\right)=f\left(\mathbf{S}\right)\text{ with probability 1}.$If $\text{β}$ is correctly specified, $E\left(Y_1|\mathbf{X}\right)$ should be some function of ${\text{β}}^{\mathrm{\top }}\mathbf{X}$ with probability 1, which means that $E\left(Y_1|\mathbf{X}\right)$ could also be written as a function of $\mathbf{S}=({\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X},{\text{β}}^{\mathrm{\top }}\mathbf{X})$. Denote this function by g, and we have $E\left(Y_1|\mathbf{X}\right)=g\left(\mathbf{S}\right)\text{ with probability 1}.$Note that the following equation holds if $\text{β}$ is correctly specified (A1) $E\left(Y_1|\mathbf{S}\right)=E\left\{E\right(Y_1|\mathbf{X}\left)|\mathbf{S}\right\}=g\left(\mathbf{S}\right)=E\left(Y_1|\mathbf{X}\right).$(A1) For the estimation in the algorithm, we estimate $E(Y_1|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\text{β}}^{\mathrm{\top }}{\mathbf{x}}_i)$ by a Nadaraya-Watson kernel regression of $Y_1$ on $\hat{\mathbf{S}}$. $\begin{array}{rl}\hat{E}(Y_1|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\text{β}}^{\mathrm{\top }}{\mathbf{x}}_i)& =\sum _{j=1}^n\frac{T_j{K}_H({\mathbf{S}}_i-{\mathbf{S}}_j)}{\sum _{k=1}^n{T}_k{K}_H({\mathbf{S}}_i-{\mathbf{S}}_k)}{Y}_{1j}\\ & =\frac{n^{-1}\sum _{j=1}^n{T}_j{Y}_{1j}{K}_H({\mathbf{S}}_i-{\mathbf{S}}_j)}{n^{-1}\sum _{k=1}^n{T}_k{K}_H({\mathbf{S}}_i-{\mathbf{S}}_k)}.\end{array}$Under the regularity conditions, we can show that $\hat{E}(Y_1|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\text{β}}^{\mathrm{\top }}{\mathbf{x}}_i)$ converges in probability to $\frac{E\left(T{Y}_1|{\mathbf{S}}_i\right)}{E\left(T|{\mathbf{S}}_i\right)}=\frac{E\left\{E\right(T{Y}_1|X_i\left)|{\mathbf{S}}_i\right\}}{E\left\{E\right(T|X_i\left)|{\mathbf{S}}_i\right\}}=\frac{E\left\{\pi \right(X_i\left)E\right(Y_1|X_i\left)|{\mathbf{S}}_i\right\}}{E\left\{\pi \right(X_i\left)|{\mathbf{S}}_i\right\}},$which can be shown equal to $E\left(Y_1|{\mathbf{S}}_i\right)$ under either Condition (a) or Condition (b).

When Condition (a) is satisfied, i.e., ${\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X}$ is correctly specified, we have shown that $\pi \left({\mathbf{X}}_i\right)=f\left({\mathbf{S}}_i\right)$ with probability 1, whereby we can further simplify the limit above. $\begin{array}{rl}& \frac{E\left\{\pi \right(X_i\left)E\right(Y_1|X_i\left)|{\mathbf{S}}_i\right\}}{E\left\{\pi \right(X_i\left)|{\mathbf{S}}_i\right\}}\\ & =\frac{f\left({\mathbf{S}}_i\right)E\left\{E\right(Y_1|X_i\left)|{\mathbf{S}}_i\right\}}{f\left({\mathbf{S}}_i\right)}\\ & =E\left\{E\right(Y_1|X_i\left)|{\mathbf{S}}_i\right\}\\ & =E\left(Y_1|{\mathbf{S}}_i\right)\text{from the law of total expectation}.\end{array}$When Condition (b) is satisfied, i.e., ${\text{β}}^{\mathrm{\top }}\mathbf{X}$ is correctly specified, we have shown that $E\left(Y_1|{\mathbf{X}}_i\right)=g\left({\mathbf{S}}_i\right)$ with probability 1. The limit could be further simplified as follows. $\begin{array}{rl}& \frac{E\left\{\pi \right(X_i\left)E\right(Y_1|X_i\left)|{\mathbf{S}}_i\right\}}{E\left\{\pi \right(X_i\left)|{\mathbf{S}}_i\right\}}\\ & =\frac{g\left({\mathbf{S}}_i\right)E\left\{\pi \right(X_i\left)|{\mathbf{S}}_i\right\}}{E\left\{\pi \right(X_i\left)|{\mathbf{S}}_i\right\}}\\ & =g\left({\mathbf{S}}_i\right)\\ & =E\left(Y_1|{\mathbf{S}}_i\right)\text{from Equation (A1)}.\end{array}$Thus, $\hat{E}(Y_1|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\text{β}}^{\mathrm{\top }}{\mathbf{x}}_i)$ is consistent for $E(Y_1|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\text{β}}^{\mathrm{\top }}{\mathbf{x}}_i)$, i.e., $\hat{E}\left(Y_1|{\mathbf{S}}_i\right)$ is consistent for $E\left(Y_1|{\mathbf{S}}_i\right)$ when either ${\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X}$ or ${\text{β}}^{\mathrm{\top }}\mathbf{X}$ is correctly specified.

Similarly, we could also show that $\hat{E}(Y_0|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\mathit{\gamma }}^{\mathrm{\top }}{\mathbf{x}}_i)$ is consistent for $E(Y_0|{\mathit{\alpha }}^{\mathrm{\top }}{\mathbf{x}}_i,{\mathit{\gamma }}^{\mathrm{\top }}{\mathbf{x}}_i)$, i.e., $\hat{E}\left(Y_0|{\mathbf{M}}_i\right)$ is consistent for $E\left(Y_0|{\mathbf{M}}_i\right)$ when either ${\mathit{\alpha }}^{\mathrm{\top }}\mathbf{X}$ or ${\mathit{\gamma }}^{\mathrm{\top }}\mathbf{X}$ is correctly specified.

Combining these two consistencies, we know that ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)-{\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ is a consistent estimator for Δ when either Condition (1) or Condition (2) is satisfied.

Note that the proposed estimator ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$ is ${\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)-{\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$ instead of ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)-{\hat{E}}_{\mathbf{M}}\left(Y_0\right)$. Next, we want to show that ${\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)$ and ${\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$ are asymptotically equivalent to ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ and ${\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ respectively under certain conditions, whereby the proposed estimator ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$ is asymptotically equivalent to ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)-{\hat{E}}_{\mathbf{M}}\left(Y_0\right)$, and has the same double robustness.

Define $\mathit{\zeta }\mathbf{=}\mathbf{(}{\mathit{\alpha }}^{\mathrm{\top }}\mathbf{,}{\text{β}}^{\mathrm{\top }}{\mathbf{)}}^{\mathrm{\top }}$ and $\hat{\mathit{\zeta }}=({\hat{\mathit{\alpha }}}^{\mathrm{\top }},{\widehat{\text{β}}}^{\mathrm{\top }}{)}^{\mathrm{\top }}$. Thus, $\mathbf{S}=b({\mathbf{X}}_i;\mathit{\zeta })$ and $\hat{\mathbf{S}}=b({\mathbf{X}}_i;\hat{\mathit{\zeta }})$, where b is a continuous function of $\mathbf{X}$.

Define $\phi (\mathbf{X};\mathit{\zeta })=\mathrm{\partial b}(\mathbf{X};\mathit{\zeta })/\mathrm{\partial }\mathit{\zeta }-E\left\{\mathrm{\partial b}\right(\mathbf{X};\mathit{\zeta }\left)/\mathrm{\partial }\mathit{\zeta }\right\}=\left[\begin{array}{cc}\mathbf{X}-E\left(\mathbf{X}\right)& \\ & \mathbf{X}-E\left(\mathbf{X}\right)\end{array}\right].$We want to show that ${\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)$ is asymptotically equivalent to ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ if $E\{|\mathbf{X}-E(\mathbf{X}\left){|}^k\right\}<\mathrm{\infty }$ for some k>6 and $\hat{\mathit{\zeta }}-\mathit{\zeta }=O_p\left(n^{-1/2}\right)$ are satisfied.

To prove the claim above, it is enough to prove that $K_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j)$ is asymptotically equivalent to $K_H({\mathbf{S}}_i-{\mathbf{S}}_j)$ under the given conditions.

We can first simplify the bivariate kernel function $K_H\left(\mathit{u}\right)$ by normalization. Note that if we use ${\mathbf{S}}^{\ast }={\hat{\mathrm{\Sigma }}}_1^{-1/2}\mathbf{S}$, ${\mathbf{S}}^{\ast }$ can also preserve the same property as $\mathbf{S}$ since there exists a continuous function g such that $f\left(\mathbf{S}\right)=f\left({\hat{\mathrm{\Sigma }}}_1^{1/2}{\mathbf{S}}^{\ast }\right)=g\left({\mathbf{S}}^{\ast }\right)$. After the linear transformation, ${\mathbf{S}}^{\ast }$ has the identity covariance, which means that we can simplify the bivariate kernel function to $K_H\left(\mathit{u}\right)=h_n^{-2}K\left(\frac{\mathit{u}}{h_n}\right).$Thus, we can rewrite the bivariate kernel functions as $K_H({\mathbf{S}}_i-{\mathbf{S}}_j)=h_n^{-2}K\left(\frac{{\mathbf{S}}_i-{\mathbf{S}}_j}{h_n}\right)$and $K_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j)=h_n^{-2}K\left(\frac{{\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j}{h_n}\right),$respectively.

Notice that $K_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j)=K_H\left(b\right({\mathbf{X}}_i;\hat{\mathit{\zeta }})-b({\mathbf{X}}_j;\hat{\mathit{\zeta }}\left)\right)$ is a function of $\hat{\mathit{\zeta }}$. By Taylor's Expansion, it could be further written as $\begin{array}{rl}{K}_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j)& =h_n^{-2}[K\left(\frac{{\mathbf{S}}_i-{\mathbf{S}}_j}{h_n}\right)+K\left(\frac{\left\{\phi \right({\mathbf{X}}_i;{\mathit{\zeta }}^{\mathrm{\prime }})-\phi ({\mathbf{X}}_j;{\mathit{\zeta }}^{\mathrm{\prime }}\left){\}}^{\mathrm{\top }}\right(\hat{\mathit{\zeta }}-\mathit{\zeta })}{h_n}\right)]\\ & =h_n^{-2}K(\frac{{\mathbf{S}}_i-{\mathbf{S}}_j}{h_n}+\frac{\left\{\phi \right({\mathbf{X}}_i;{\mathit{\zeta }}^{\mathrm{\prime }})-\phi ({\mathbf{X}}_j;{\mathit{\zeta }}^{\mathrm{\prime }}\left){\}}^{\mathrm{\top }}\right(\hat{\mathit{\zeta }}-\mathit{\zeta })}{h_n}),\end{array}$where ${\mathit{\zeta }}^{\mathrm{\prime }}$ is on the line segment between $\hat{\mathit{\zeta }}$ and $\mathit{\zeta }$. Since $E\left\{|\phi \right(\mathbf{X};{\mathit{\zeta }}^{\mathrm{\prime }}\left){|}^k\right\}<\mathrm{\infty }$, we have $|\phi (\mathbf{X};\mathit{\zeta })|\le n^{1/k}$ almost surely by the Borel-Cantelli lemma. Thus, $|\phi ({\mathbf{X}}_j;{\mathit{\zeta }}^{\mathrm{\prime }})-\phi ({\mathbf{X}}_i;{\mathit{\zeta }}^{\mathrm{\prime }})|\le 2n^{1/k}$ almost surely by the triangle inequality. Combined with the optimal bandwidth $h_n\sim n^{-1/3}$ and the condition that $\hat{\mathit{\zeta }}-\mathit{\zeta }=O_p\left(n^{-1/2}\right)$, $K\left(\frac{\left\{\phi \right({\mathbf{X}}_i;{\mathit{\zeta }}^{\mathrm{\prime }})-\phi ({\mathbf{X}}_j;{\mathit{\zeta }}^{\mathrm{\prime }}\left){\}}^{\mathrm{\top }}\right(\hat{\mathit{\zeta }}-\mathit{\zeta })}{h_n}\right)$ is $O\left(n^{1/k+1/3-1/2}\right)$. Since k>6, this term goes to zero, which means that $K_H({\hat{\mathbf{S}}}_i-{\hat{\mathbf{S}}}_j)$ is asymptotically equivalent to $K_H({\mathbf{S}}_i-{\mathbf{S}}_j)$.

Thus, that the kernel regression estimator ${\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)$ is consistent for ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ if $E\{|\mathbf{X}-E(\mathbf{X}\left){|}^k\right\}<\mathrm{\infty }$ for some k>6 and $\hat{\mathit{\zeta }}-\mathit{\zeta }=O_p\left(n^{-1/2}\right)$ are satisfied.

Similarly, to prove the consistency of ${\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$, we can define ${\mathit{\zeta }}^{\mathbf{\ast }}\mathbf{=}\mathbf{(}{\mathit{\alpha }}^{\mathrm{\top }}\mathbf{,}{\mathit{\gamma }}^{\mathrm{\top }}{\mathbf{)}}^{\mathrm{\top }}$ and ${\hat{\mathit{\zeta }}}^{\ast }=({\hat{\mathit{\alpha }}}^{\mathrm{\top }},{\hat{\mathit{\gamma }}}^{\mathrm{\top }}{)}^{\mathrm{\top }}$. We can use a similar argument to show that ${\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$ is asymptotically equivalent to ${\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ if $E\{|\mathbf{X}-E(\mathbf{X}\left){|}^k\right\}<\mathrm{\infty }$ for some k>6 and ${\hat{\mathit{\zeta }}}^{\ast }-{\mathit{\zeta }}^{\ast }=O_p\left(n^{-1/2}\right)$ are satisfied.

Hence, ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}={\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)-{\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$ is asymptotically equivalent to ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)-{\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ under the given conditions, which indicates that ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$ is asymptotically consistent if either Condition (i) or Condition (ii) is satisfied.

Proof

Proof of Theorem 3.2

To prove the optimal efficiency for ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$, we need to prove the optimal efficiency for ${\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)$ and ${\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$.

We can first show that if $\mathit{\alpha }$ is correctly specified, ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ has an asymptotic normal distribution $n^{1/2}\left({\hat{E}}_{\mathbf{S}}\right(Y_1)-E(Y_1\left)\right)\to N(0,\text{var}\{E\left(Y_1|\mathbf{S}\right)\}+E\left\{\frac{\text{var}\left(Y_1|\mathbf{S}\right)}{\pi \left(\mathbf{X}\right)}\right\}).$And then, if $\text{β}$ is additionally satisfied, we can further prove that ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ achieves optimal efficiency.

To prove that $n^{1/2}\left({\hat{E}}_{\mathbf{S}}\right(Y_1)-E(Y_1\left)\right)$ has such asymptotic distribution, we could split it into three parts, $\begin{array}{rl}{n}^{1/2}\left({\hat{E}}_{\mathbf{S}}\right(Y_1)-E(Y_1\left)\right)& =n^{1/2}\left[\frac{1}{n}\sum _{i=1}^n\left\{\hat{E}\right(Y_1|{\mathbf{S}}_i)-E(Y_1\left)\right\}\right]\\ & =n^{1/2}{A}_n+n^{1/2}{B}_n+n^{1/2}{C}_n,\end{array}$where $\begin{array}{rl}{A}_n& =\frac{1}{n}\sum _{i=1}^n\left\{E\right(Y_1|{\mathbf{S}}_i)-E(Y_1\left)\right\},\\ B_n& =\frac{1}{n}\sum _{i=1}^{n}E\left\{\hat{E}\right(Y_1|{\mathbf{S}}_i)-E(Y_1|{\mathbf{S}}_i\left)|{\mathrm{\Theta }}_i\right\},\\ C_n& =\frac{1}{n}\sum _{i=1}^n\left[\hat{E}\right(Y_1|{\mathbf{S}}_i)-E(Y_1|{\mathbf{S}}_i)-E\left\{\hat{E}\right(Y_1|{\mathbf{S}}_i)-E(Y_1|{\mathbf{S}}_i\left)|{\mathrm{\Theta }}_i\right\}],\\ {\mathrm{\Theta }}_i& =\left\{\right({\mathbf{X}}_j,Y_{1j},T_j)|j\ne i\}.\end{array}$For $A_n$, we can show that $n^{1/2}{A}_n$ asymptotically converges to $N(0,E(Y_1|{\mathbf{S}}_i\left)\right)$ due to the central limit theorem. For $B_n$, we can first rewrite $\hat{E}\left(Y_1|{\mathbf{S}}_i\right)$ as $\begin{array}{rl}\hat{E}\left(Y_1|{\mathbf{S}}_i\right)& =\sum _{j=1}^n\frac{Y_{1j}{T}_j{K}_H({\mathbf{S}}_i-{\mathbf{S}}_j)}{\sum _{k=1}^n{T}_k{K}_H({\mathbf{S}}_i-{\mathbf{S}}_k)}\\ & =\frac{n^{-1}\sum _{j=1}^n{T}_j{Y}_{1j}{K}_H({\mathbf{S}}_i-{\mathbf{S}}_j)}{n^{-1}\sum _{k=1}^n{T}_k{K}_H({\mathbf{S}}_i-{\mathbf{S}}_k)}.\end{array}$When $\mathit{\alpha }$ is correctly specified, $\pi \left(\mathbf{X}\right)=f\left(\mathbf{S}\right)$ with probability 1. The denominator could be further written as $n^{-1}\sum _{k=1}^n{T}_k{K}_H({\mathbf{S}}_i-{\mathbf{S}}_k)=f\left({\mathbf{S}}_i\right)+o_p\left(1\right)$. Hence, $B_n$ could be rewritten as $B_n=\frac{1}{n}\sum _{j=1}^n{T}_{j}E\left[\frac{K_H({\mathbf{S}}_i-{\mathbf{S}}_j)\{Y_{1}j-E(Y_1|{\mathbf{S}}_i\left)\right\}}{f\left({\mathbf{S}}_i\right)}\right|{\mathrm{\Theta }}_i]\{1+o_p(1\left)\right\}.$We can show that $\sqrt{n}(B_n-B_n^{\ast })=o_p\left(1\right)$ with $B_n^{\ast }=\frac{1}{n}\sum _{j=1}^n{T}_j\frac{Y_{1j}-E\left(Y_1|{\mathbf{S}}_i\right)}{f\left({\mathbf{S}}_j\right)}$by Theorem 2.1 from Cheng (1994). Thus, by the central limit, we can show that $n^{1/2}{B}_n$ asymptotically converges to $N(0,E\left\{\frac{\text{var}\left(Y|\mathbf{S}\right)}{f\left(\mathbf{S}\right)}\right\})$. For $C_n$, we know that $E\left(C_n\right)=0$ and $\mathrm{nE}\left(C_n^2\right)\le E\left[\right\{\hat{E}\left(Y_1|{\mathbf{S}}_i\right)-E\left(Y_1|{\mathbf{S}}_i\right){\}}^2]=O\left(\right\{\text{tr}\left(H{H}^{\mathrm{\top }}\right){\}}^2+\left\{n\det \right(H\left){\}}^{-1}\right).$Thus, $C_n=o_p\left(n^{-1/2}\right)$.

Combining all three parts together, we know that we have the following asymptotic distribution: $n^{1/2}\left({\hat{E}}_{\mathbf{S}}\right(Y_1)-E(Y_1\left)\right)\to N(0,\text{var}\{E\left(Y_1|S\right)\}+E\left\{\frac{\text{var}\left(Y_1|S\right)}{\pi \left(\mathbf{X}\right)}\right\}).$Note that the asymptotic variance $n\left[\text{var}\right\{{\hat{E}}_S\left(Y_1\right)\left\}\right]$ could be rewritten as $\begin{array}{rl}n\left[\text{var}\right\{{\hat{E}}_S\left(Y_1\right)\left\}\right]& =\text{var}\left\{E\right(Y_1|\mathbf{S}\left)\right\}+E\left\{\frac{\text{var}\left(Y|\mathbf{S}\right)}{\pi \left(\mathbf{X}\right)}\right\}\\ & =\text{var}\left(Y_1\right)-E\left[\text{var}\right(Y_1|\mathbf{S}\left)\right]+E\left\{\frac{\text{var}\left(Y_1|\mathbf{S}\right)}{\pi \left(\mathbf{X}\right)}\right\}\\ & =\text{var}\left(Y_1\right)+E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{S}\right)\right].\end{array}$Additionally, if $\text{β}$ is correctly specified, we have the equation $E\left(Y|\mathbf{X}\right)=E\left(Y|\mathbf{S}\right)$ according to the argument in the proof of Theorem 3.1. From this equation, $\text{var}\left(Y_1|\mathbf{S}\right)$ could be further rewritten as $E\left\{\text{var}\right(Y_1|\mathbf{X}\left)|\mathbf{S}\right\}$: $\begin{array}{rl}\text{var}\left(Y_1|\mathbf{S}\right)& =E\left[\right\{Y_1-E\left(Y_1|\mathbf{S}\right){\}}^2|\mathbf{S}]\\ & =E\left(E\left[\right\{Y_1-E\left(Y_1|\mathbf{S}\right){\}}^2|\mathbf{X}]|\mathbf{S}\right)\\ & =E\left\{\text{var}\right(Y_1|\mathbf{X}\left)|\mathbf{S}\right\}.\end{array}$Thus, the second term $E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{S}\right)\right]$ could be rewritten according to the equation $E\left(Y|\mathbf{X}\right)=E\left(Y|\mathbf{S}\right)$ and the condition that $\pi \left(\mathbf{X}\right)=f\left(\mathbf{S}\right)$ with probability 1. $\begin{array}{rl}E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{S}\right)\right]& =E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}E\left\{\text{var}\right(Y_1|\mathbf{X}\left)|\mathbf{S}\right\}\right]\\ & =E\left[\right\{{f\left(\mathbf{S}\right)}^{-1}-1\left\}E\left\{\text{var}\right(Y_1|\mathbf{X}\left)|\mathbf{S}\right\}\right]\\ & =E\left[E\left\{\right({f\left(\mathbf{S}\right)}^{-1}-1\left)\text{var}\right(Y_1|\mathbf{X}\left)|\mathbf{S}\right\}\right]\\ & =E\left[\right\{{f\left(\mathbf{S}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{X}\right)\right]\\ & =E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{X}\right)\right].\end{array}$Combining all the equations above, we know that the asymptotic variance could be transformed into the following expression, $n\left[\text{var}\right\{{\hat{E}}_S\left(Y_1\right)\left\}\right]=\text{var}\left(Y_1\right)+E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{S}\right)\right]=\text{var}\left(Y_1\right)+E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{X}\right)\right].$According to Hahn (1998), $\text{var}\left(Y_1\right)+E\left[\right\{{\pi \left(\mathbf{X}\right)}^{-1}-1\left\}\text{var}\left(Y_1|\mathbf{X}\right)\right]$ is the optimal efficient variance that a double robust estimator can achieve. It means that when $\mathit{\alpha }$ and $\text{β}$ are both correctly specified, ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ achieves optimal efficiency.

We can use a similar argument to prove that when $\mathit{\alpha }$ and $\mathit{\gamma }$ are both correctly specified, ${\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ achieves optimal efficiency.

Note that the optimal efficiency for ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)$ and ${\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ indicates that ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)-{\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ achieves optimal efficiency for Δ when both Condition (1) and Condition (2) are satisfied.

According to the proof of Theorem 3.1, we know that the proposed estimator ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}={\hat{E}}_{\hat{\mathbf{S}}}\left(Y_1\right)-{\hat{E}}_{\hat{\mathbf{M}}}\left(Y_0\right)$ is asymptotically equivalent to ${\hat{E}}_{\mathbf{S}}\left(Y_1\right)-{\hat{E}}_{\mathbf{M}}\left(Y_0\right)$ under the given conditions. It implies that when both Condition (1) and Condition (2) are satisfied, ${\hat{\mathrm{\Delta }}}_{\mathrm{rdr}}$ also achieves optimal efficiency for Δ.