Multiply robust estimation for average treatment effect among treated

Abstract

We propose a multiply robust estimator for the Average Treatment Effect Among the Treated (ATT). The proposed estimation procedure can simultaneously accommodate multiple working models for both the propensity score and the conditional mean of the counterfactual outcome given covariates. In addition, it can explicitly balance a set of user-specified moments of the covariate distributions between the treatment groups. The resulting estimator is consistent if any working model is correctly specified. With the data generating process typically unknown for observational studies, the proposed method provides substantial robustness against possible model misspecifications compared to existing estimators of the ATT. Simulation results show the excellent finite sample performance of the proposed estimator.

Keywords:

• Causal inference
• covariate balancing
• empirical likelihood
• multiple robustness
• propensity score

1. Introduction

In causal inference, since we only observe what happens to an individual under the treatment condition they actually receive, it is generally impossible to estimate the causal effects for individuals. The causal effects one typically considers involve summary statistics of the individual effects across populations or sub-populations of interest. Two widely considered causal summaries are the average treatment effect (ATE) and the average treatment effect among the treated (ATT). The ATE of a treatment relative to a control is the comparison of the mean outcome which had the entire population been treated versus had the entire population been the control. The ATT is the comparison of the mean outcome under treatment among those who are treated with the mean outcome which had the treated subjects received control instead. A study can estimate both ATE and ATT, but one or the other may be better suited for a particular situation. The ATE may be of more interest if each treatment can potentially be offered to every member of the population. Conversely, if the research question focuses on the effectiveness of an alternative treatment were it to replace the standard treatment, and then the ATT may be of more interest because it measures the relative effectiveness of the two treatment options on the population that is receiving the standard treatment.

There has been a large literature on estimating the ATE and the ATT for observational studies, where a typical consideration is confounding in the sense that individual characteristics are related to both the treatment assignment and the outcome of interest. Propensity score (PS) based methods, where the PS is the probability of receiving a treatment given covariates (Rosenbaum & Rubin, 1983), are commonly used to deal with confounding and to achieve balance of covariate distributions between different treatment groups. Such methods include PS matching (e.g., Abadie & Imbens, 2006; Rosenbaum & Rubin, 1985) and weighting (e.g., Hirano et al., 2003). Refer to Imbens and Rubin (2015) and Hernán and Robins (2018) for more details. Parametric modelling of the PS is common, especially when the dimension of covariates is moderate to large. Misspecification of the PS model is usually a major concern as it may lead to substantial estimation bias.

To mitigate the impact of misspecification of the PS model, substantial interests have been given to doubly robust estimators, which involve models for both the PS and the conditional mean of the counterfactual outcome, or outcome regression (OR), and remain consistent if either model is correctly specified. The original doubly robust estimator was constructed through augmented inverse probability weighting (AIPW) in the missing data context (Robins et al., 1994). Since then, a large number of doubly robust estimators have been proposed in both missing data and causal inference settings (e.g., Bang & Robins, 2005; Cao et al., 2009; Han, 2012; Kang & Schafer, 2007; Qin et al., 2008; Qin & Zhang, 2007; Rotnitzky et al., 2012; Tan, 2010; van der Laan & Gruber, 2010). When both the PS model and the OR model are wrong, these estimators are in general no longer consistent.

As an improvement over double robustness, in the missing data context, Han and Wang (2013) proposed a multiply robust estimation procedure that can simultaneously accommodate multiple working models for both PS and OR, and the resulting estimators are consistent if any of these models is correctly specified. Since the data generating process for observational studies is typically unknown, it is common that several candidate models all seem reasonable yet none rules out the possibility of others, especially when the dimension of covariates is moderate to large. In such a case, the multiply robust method in Han and Wang (2013) provides a useful tool for data analysis with more protection on estimation consistency. Such a method has attracted considerable interest in both missing data and causal inference research (e.g., Chan & Yam, 2014; Chan et al., 2016; S. Chen & Haziza, 2017; Duan & Yin, 2017; Han, 2014a, 2014b, 2016a, 2016b; Han et al., 2019; Li et al., 2020). Especially, Wang (2019) extended the multiply robust method to the estimation of the ATE. It is worth pointing out that the term “ multiply robust” has also been used by other authors in different settings with different meanings. For example, in Molina et al. (2017), it refers to estimation consistency when various combinations of the components of a factorized likelihood are correctly modelled, while in Wang and Tchetgen Tchetgen (2018) and Shi et al. (2020) it refers to estimation consistency being achieved across the union of three different observed data models. We use “ multiply robust” to refer to the property that estimation consistency is achieved if one of the multiple working models for the same quantity is correctly specified.

Despite the success in dealing with missing data problems and in estimating the ATE, multiply robust estimators have not been developed for estimating the ATT. In this paper, we construct such a desirable estimator for ATT. In addition to being multiply robust, the proposed estimator can easily achieve certain level of balancing of covariate distributions between treatment groups. Covariate balancing is a highly desired property for causal effect estimation. For the estimation of ATT, Hainmueller (2012) proposed the entropy balancing (EB) method by imposing a set of balance constraints so that certain moments of covariate distributions for different treatment groups match exactly. Zhao and Percival (2017) found that EB implicitly fits a logistic linear regression model for the PS and a linear regression model for the OR, and when either model is correctly specified the EB estimator is consistent. The estimator we propose preserves the same balance of covariate distributions as EB, and in addition it accounts for multiple models for PS and OR simultaneously so that consistency of the resulting estimator is guaranteed if any one model is correctly specified.

The rest of the paper is organized as follows. Section 2 gives the setup and a brief review of some existing methods. Section 3 presents our proposed multiply robust estimator for the ATT. Simulation studies are provided in Section 4, followed by some discussion in Section 5.

2. Setup and some existing methods

Let A denote the treatment indicator, where A = 1 for the treatment of interest and A = 0 for control. Let $Y^1$ and $Y^0$ denote the counterfactual outcomes under treatment and under control, respectively. Let Y denote the observed outcome in the data, for which we make the consistency assumption that $Y=A{Y}^1+(1-A)Y^0$. Let $\mathit{X}$ denote the pre-treatment covariates in the dataset, including potential confounders. The potential full data is $(A,Y^1,Y^0,\mathit{X})$ whereas the observed data is $(A,Y,\mathit{X})$. The causal effect we are interested in is the ATT $\tau =E(Y^1-Y^0\mid A=1)$, which is also ${\tau }_{1,1}-{\tau }_{1,0}$, where ${\tau }_{a,b}$ is the mean outcome for subjects who receive treatment a had they instead received treatment b, i.e., ${\tau }_{a,b}=E[Y^b|A=a]$.

To ensure the identifiability of τ from the observed data, we make some assumptions following the main literature (e.g., Rosenbaum & Rubin, 1983). First, we assume that all potential confounders are included in $\mathit{X}$, or in other words, the treatment assignment process and the counterfactual outcome are independent given all the covariates measured. Second, we assume that every subject has a positive probability of being assigned to either the treatment or the control group. These two assumptions are formalized as follows.

Assumption 2.1

no unmeasured confounders assumption

$A\perp (Y^0,Y^1)\mid \mathit{X}$.

Assumption 2.2

strict positivity

$0<{\sigma }_1<P(A=1\mid \mathit{X})<{\sigma }_2<1$ with probability one for some positive constants ${\sigma }_1$ and ${\sigma }_2$.

Let $\pi \left(\mathit{X}\right)=P(A=1\mid \mathit{X})$ denote the PS for treatment assignment. In the following we give a brief review of some widely used estimators of the ATT. Note that ${\tau }_{1,1}=E(Y^1\mid A=1)$ can be consistently estimated by the sample average ${\hat{\tau }}_{1,1}=\left(\sum _{i=1}^n{A}_i{Y}_i\right)/\left(\sum _{i=1}^n{A}_i\right)$ over the treatment group. Therefore, the estimation of the ATT reduces to the estimation of ${\tau }_{1,0}=E(Y^0\mid A=1)$, where the counterfactual outcome $Y^0$ is not observable for individuals in the A = 1 group.

One straightforward way to estimate ${\tau }_{1,0}$ is to impute $Y^0$ for those individuals in the A = 1 group. Suppose $d_0\left(\mathit{\gamma }\right)$ is a parametric regression model for $E(Y^0\mid \mathit{X})$ that is parametrized by $\mathit{\gamma }$. Because $E(Y^0\mid \mathit{X})=E(Y^0\mid \mathit{X},A=0)$ from Assumption 1 and $Y^0$ is fully observed in the A = 0 group, $\mathit{\gamma }$ can be estimated by $\hat{\mathit{\gamma }}$ based on individuals in the control group alone. Also from Assumption 1, we have $E(Y^0\mid \mathit{X})=E(Y^0\mid \mathit{X},A=1)$, and thus the unobserved $Y^0$ in the treatment group can be imputed by $d_0\left(\hat{\mathit{\gamma }}\right)$. Therefore, an outcome regression estimator of ${\tau }_{1,0}$ is ${\hat{\tau }}_{1,0\text{reg}}=\sum _{i=1}^n{A}_i{d}_{0i}\left(\hat{\mathit{\gamma }}\right)/\left(\sum _{i=1}^n{A}_i\right)$.

A widely used alternative method is inverse probability weighting (IPW). From Assumption 1, since $f(Y^1,Y^0,\mathit{X}\mid A=1)=\frac{\pi \left(\mathit{X}\right)f(Y^1,Y^0,\mathit{X})}{P(A=1)}\text{and}f(Y^1,Y^0,\mathit{X}\mid A=0)=\frac{\{1-\pi (\mathit{X}\left)\right\}f(Y^1,Y^0,\mathit{X})}{P(A=0)},$ where $f(\cdot )$ is a generic notation for a probability density function, we have (1) $f(Y^1,Y^0,\mathit{X}\mid A=1)=\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\frac{P(A=0)}{P(A=1)}f(Y^1,Y^0,\mathit{X}\mid A=0).$(1) Therefore, an estimator of ${\tau }_{1,0}=E(Y^0\mid A=1)$ can be constructed by properly weighting the observed $Y^0$ in the A = 0 group, and this leads to an IPW estimator (2) ${\hat{\tau }}_{1,0\text{ipw}}=\frac{1}{\sum _{s=1}^n(1-A_s)}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}\frac{\sum _{s=1}^n(1-A_s)/n}{\sum _{s=1}^n{A}_s/n}{Y}_i=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}{Y}_i,$(2) where $\hat{\pi }\left(\mathit{X}\right)$ is the estimated value of $\pi \left(\mathit{X}\right)$.

Consistency of ${\hat{\tau }}_{1,0\text{reg}}$ and ${\hat{\tau }}_{1,0\text{ipw}}$ requires correct modelling of $E(Y^0\mid \mathit{X})$ and $\pi \left(\mathit{X}\right)$, respectively. To improve the robustness against possible model misspecifications, the AIPW method combines the models for $E(Y^0\mid \mathit{X})$ and $\pi \left(\mathit{X}\right)$ so that estimation consistency is guaranteed if either model is correctly specified but not necessarily both. The AIPW estimator for ${\tau }_{1,0}$ is given as ${\hat{\tau }}_{1,0\text{aipw}}={\hat{\tau }}_{1,0\text{reg}}+\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}\{Y_i-d_{0i}(\hat{\mathit{\gamma }}\left)\right\}.$ To achieve certain covariate balancing, the EB method (Hainmueller, 2012) considers matching covariate moments between the treatment and the control groups. Specifically, let $w_i$ be a set of positive weights assigned to the control subjects $\{i:A_i=0\}$. The EB method imposes covariate balancing constraints (3) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i\mathit{h}\left({\mathit{X}}_i\right)=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=1\}}\mathit{h}\left({\mathit{X}}_i\right)$(3) on $w_i$, where $\mathit{h}\left(\mathit{X}\right)$ contains some user-specified moments of $\mathit{X}$. These constraints ensure that certain moments of $\mathit{X}$ exactly match between the two groups. The EB estimator of ${\tau }_{1,0}$ is ${\hat{\tau }}_{1,0\text{eb}}=\sum _{\{i:A_i=0\}}{\hat{w}}_{\text{eb},i}{Y}_i$ where ${\hat{w}}_{\text{eb},i}$ minimizes $\sum _{\{i:A_i=0\}}{w}_i\log w_i$ subject to the constraints in (3(3) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i\mathit{h}\left({\mathit{X}}_i\right)=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=1\}}\mathit{h}\left({\mathit{X}}_i\right)$(3) ). Although no parametric models are explicitly fitted by the EB method, Zhao and Percival (2017) showed that, implicitly, the EB method fits linear regression models for $\text{logit}\left\{\pi \right(\mathit{X}\left)\right\}$ and $E(Y^0\mid \mathit{X})$ with components of $\mathit{h}\left(\mathit{X}\right)$ as regressors. The EB estimator ${\hat{\tau }}_{1,0\text{eb}}$ is consistent if either model is correctly specified, and thus is doubly robust.

3. The proposed multiply robust estimator for ATT

Our goal is to construct an easy-to-implement estimator of the ATT that is multiply robust and achieves explicit covariate balancing as the EB estimator. Let $\mathcal{P}=\left\{{\pi }^{\left(j\right)}\right({\mathit{\alpha }}^{\left(j\right)}):j=1,\dots ,J\}$ denote a set of multiple parametric models for $\pi \left(\mathit{X}\right)$ and ${\mathcal{D}}_0=\left\{d_0^{\left(k\right)}\right({\mathit{\gamma }}^{\left(k\right)}):k=1,\dots ,K\}$ a set of multiple parametric models for $E(Y^0\mid \mathit{X})$, where ${\mathit{\alpha }}^{\left(j\right)}$ and ${\mathit{\gamma }}^{\left(k\right)}$ are the corresponding parameters. ${\mathit{\alpha }}^{\left(j\right)}$ is typically estimated by ${\hat{\mathit{\alpha }}}^{\left(j\right)}$ that maximizes the binomial likelihood (4) $\prod _{i=1}^n\left\{{\pi }_i^{\left(j\right)}\right({\mathit{\alpha }}^{\left(j\right)}\left){\}}^{A_i}\right\{1-{\pi }_i^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right){\}}^{1-A_i},$(4) and ${\mathit{\gamma }}^{\left(k\right)}$ is typically estimated by ${\hat{\mathit{\gamma }}}^{\left(k\right)}$ based on individuals in the control group because $E(Y^0\mid \mathit{X})=E(Y^0\mid \mathit{X},A=0)$ from Assumption 1 and $Y^0$ is fully observed in the control group.

From (1(1) $f(Y^1,Y^0,\mathit{X}\mid A=1)=\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\frac{P(A=0)}{P(A=1)}f(Y^1,Y^0,\mathit{X}\mid A=0).$(1) ) it is easy to see that, for any function $b\left(\mathit{X}\right)$, we have (5) $E\left\{b\right(\mathit{X})\mid A=1\}=E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}b\right(\mathit{X})\mid A=0\}\frac{P(A=0)}{P(A=1)}.$(5) In particular, taking $b\left(\mathit{X}\right)\equiv 1$ gives $1=E\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\mid A=0\}\frac{P(A=0)}{P(A=1)},$ and thus (6) $\begin{array}{rl}E\left\{b\right(\mathit{X})\mid A=1\}& =E\left\{b\right(\mathit{X})\mid A=1\}E\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\mid A=0\}\frac{P(A=0)}{P(A=1)}\\ & =E\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}E\left\{b\right(\mathit{X})\mid A=1\}\mid A=0\}\frac{P(A=0)}{P(A=1)}.\end{array}$(6) Subtracting (6(6) $\begin{array}{rl}E\left\{b\right(\mathit{X})\mid A=1\}& =E\left\{b\right(\mathit{X})\mid A=1\}E\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\mid A=0\}\frac{P(A=0)}{P(A=1)}\\ & =E\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}E\left\{b\right(\mathit{X})\mid A=1\}\mid A=0\}\frac{P(A=0)}{P(A=1)}.\end{array}$(6) ) from (5(5) $E\left\{b\right(\mathit{X})\mid A=1\}=E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}b\right(\mathit{X})\mid A=0\}\frac{P(A=0)}{P(A=1)}.$(5) ) leads to (7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) Our method starts by constructing an empirical version of (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ). Let $w_i$ be a set of positive weights assigned on the control subjects with $A_i=0$ such that $\sum _{\{i:A_i=0\}}{w}_i=1$, and then an empirical version of (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ) is $\sum _{\{i:A_i=0\}}{w}_i\left\{b\right({\mathit{X}}_i)-\frac{1}{\sum _{s=1}^n{A}_s}\sum _{s=1}^n{A}_{s}b({\mathit{X}}_s\left)\right\}=0.$ To achieve multiple robustness so that ${\tau }_{1,0}$ is consistently estimated when any one model is correctly specified, we take $b\left(\mathit{X}\right)$ to be $1/{\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$, $j=1,\dots ,J$, and $d_0^{\left(k\right)}\left({\mathit{\gamma }}^{\left(k\right)}\right)$, $k=1,\dots ,K$, and consider the constraints on the $w_i$ as (8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) where $\stackrel{\mathbf{ˆ}}{\mathit{\alpha }}=({\stackrel{\mathbf{ˆ}}{\mathit{\alpha }}}^{\left(1\right)},\dots ,{\stackrel{\mathbf{ˆ}}{\mathit{\alpha }}}^{\left(J\right)})$, $\stackrel{\mathbf{ˆ}}{\mathit{\gamma }}=({\stackrel{\mathbf{ˆ}}{\mathit{\gamma }}}^{\left(1\right)},\dots ,{\stackrel{\mathbf{ˆ}}{\mathit{\gamma }}}^{\left(J\right)})$ and $\stackrel{\mathbf{ˆ}}{\mathit{g}}(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\left(\begin{array}{c}{\displaystyle \frac{1}{{\pi }^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}}-{\displaystyle \frac{1}{\sum _{s=1}^n{A}_s}}\sum _{s=1}^n{A}_s{\displaystyle \frac{1}{{\pi }_s^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}}\\ ⋮\\ {\displaystyle \frac{1}{{\pi }^{\left(J\right)}\left({\hat{\mathit{\alpha }}}^{\left(J\right)}\right)}}-{\displaystyle \frac{1}{\sum _{s=1}^n{A}_s}}\sum _{s=1}^n{A}_s{\displaystyle \frac{1}{{\pi }_s^{\left(J\right)}\left({\hat{\mathit{\alpha }}}^{\left(J\right)}\right)}}\\ d_0^{\left(1\right)}\left({\hat{\mathit{\gamma }}}^{\left(1\right)}\right)-{\displaystyle \frac{1}{\sum _{s=1}^n{A}_s}}\sum _{s=1}^n{A}_s{d}_{0s}^{\left(1\right)}\left({\hat{\mathit{\gamma }}}^{\left(1\right)}\right)\\ ⋮\\ d_0^{\left(K\right)}\left({\hat{\mathit{\gamma }}}^{\left(K\right)}\right)-{\displaystyle \frac{1}{\sum _{s=1}^n{A}_s}}\sum _{s=1}^n{A}_s{d}_{0s}^{\left(K\right)}\left({\hat{\mathit{\gamma }}}^{\left(K\right)}\right)\\ \mathit{h}\left(\mathit{X}\right)-{\displaystyle \frac{1}{\sum _{s=1}^n{A}_s}}\sum _{s=1}^n{A}_s\mathit{h}\left({\mathit{X}}_s\right)\end{array}\right).$ Note that here $\stackrel{\mathbf{ˆ}}{\mathit{g}}(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})$ can also contain components based on taking $b\left(\mathit{X}\right)$ to be $\mathit{h}\left(\mathit{X}\right)$, a vector of user-specified moments of $\mathit{X}$. This can be used to achieve certain degree of covariate balancing between the treatment and control groups, similar to the EB method.

Subject to the constraints in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ), we consider the weights ${\hat{w}}_{\text{mr},i}$ on the subjects in the control group $\{i:A_i=0\}$ that maximize the empirical likelihood $\prod _{\{i:A_i=0\}}{w}_i$ and propose to estimate ${\tau }_{1,0}$ by ${\hat{\tau }}_{1,0\text{mr}}=\sum _{\{i:A_i=0\}}{\hat{w}}_{\text{mr},i}{Y}_i$. Following the standard empirical likelihood technique (e.g., Qin & Lawless, 1994), we have ${\hat{w}}_{\text{mr},i}=\frac{1}{\sum _{i=1}^n(1-A_i)}\frac{1}{1+{\hat{\mathit{\rho }}}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})},\text{for }i\text{ satisfying }{A}_i=0,$ where $\hat{\mathit{\rho }}$ solves (9) $\sum _{\{i:A_i=0\}}\frac{{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}{1+{\mathit{\rho }}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}=\mathbf{0}.$(9) For implementation, directly solving (9(9) $\sum _{\{i:A_i=0\}}\frac{{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}{1+{\mathit{\rho }}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}=\mathbf{0}.$(9) ) for $\hat{\mathit{\rho }}$ is not the ideal way because, as pointed out in Han (2014a), (9(9) $\sum _{\{i:A_i=0\}}\frac{{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}{1+{\mathit{\rho }}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}=\mathbf{0}.$(9) ) typically has multiple roots but only one of them makes the ${\hat{w}}_{\text{mr},i}$ positive. We recommend calculating $\hat{\mathit{\rho }}$ by minimizing $F\left(\mathit{\rho }\right)\equiv -\sum _{\{i:A_i=0\}}\log \{1+{\mathit{\rho }}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }}\left)\right\}$, which is the negative antiderivative of the left-hand side of (9(9) $\sum _{\{i:A_i=0\}}\frac{{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}{1+{\mathit{\rho }}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}=\mathbf{0}.$(9) ). As shown in Han (2014a), this is a convex minimization that always has a unique minimizer when $\mathbf{0}$ is inside the convex hull of $\left\{{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i\right(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }}):A_i=0\}$, which indeed holds at least when n is large because of the moment equality (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ) and (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) being an empirical version of (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ). The unique minimizer of $F\left(\mathit{\rho }\right)$ is $\hat{\mathit{\rho }}$ needed for calculating ${\hat{w}}_{\text{mr},i}$ and can be easily found by a Newton-Raphson-type algorithm. Refer to J. Chen et al. (2002) and Han (2014a) for such an algorithm.

In the following we provide arguments and explanations for the multiple robustness of ${\hat{\tau }}_{1,0\text{mr}}$. We do not include detailed technical conditions and mathematical proofs for these arguments so that the main ideas are easier to follow. To see the consistency of ${\hat{\tau }}_{1,0\text{mr}}$ when $\mathcal{P}$ contains a correctly specified model for $\pi \left(\mathit{X}\right)$, say ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)$ without loss of generality, let $\hat{\vartheta }=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{i=1}^n{A}_i\frac{1}{{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}$ and define $\hat{\mathit{\lambda }}$ in such a way that ${\hat{\rho }}_1=({\hat{\lambda }}_1+1)/(\hat{\vartheta }-1)$ and ${\hat{\mathit{\rho }}}_{-1}={\hat{\mathit{\lambda }}}_{-1}/(\hat{\vartheta }-1)$, where ${\hat{\rho }}_1$ and ${\hat{\lambda }}_1$ are the first component of $\hat{\mathit{\rho }}$ and $\hat{\mathit{\lambda }}$, respectively, and ${\hat{\mathit{\rho }}}_{-1}$ and ${\hat{\mathit{\lambda }}}_{-1}$ are vectors of the rest components. Then some simple algebra shows that ${\hat{w}}_{\text{mr},i}=\frac{\hat{\vartheta }-1}{\sum _{s=1}^n(1-A_s)}\frac{{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{1-{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}\frac{1}{1+{\hat{\mathit{\lambda }}}^{\mathrm{\top }}\frac{{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{1-{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})},\text{for }i\text{ satisfying }{A}_i=0$ and (9(9) $\sum _{\{i:A_i=0\}}\frac{{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}{1+{\mathit{\rho }}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}=\mathbf{0}.$(9) ) becomes an equation for $\hat{\mathit{\lambda }}$ $\sum _{\{i:A_i=0\}}\frac{\frac{{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{1-{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}{1+{\hat{\mathit{\lambda }}}^{\mathrm{\top }}\frac{{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{1-{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}=\mathbf{0}.$ From the Z-estimator theory (e.g., van der Vaart, 1998) and because of the moment equality (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ), we must have $\hat{\mathit{\lambda }}=O_p\left(n^{-1/2}\right)$, which leads to ${\hat{w}}_{\text{mr},i}=\frac{\hat{\vartheta }-1}{\sum _{s=1}^n(1-A_s)}\frac{{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}{1-{\pi }_i^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)}\{1+O_p(n^{-1/2}\left)\right\},\text{for }i\text{ satisfying }{A}_i=0.$ In other words, when $\mathcal{P}$ contains a correctly specified model for $\pi \left(\mathit{X}\right)$, this correct model is implicitly accounted for by the empirical likelihood procedure to make ${\hat{w}}_{\text{mr},i}$ the form of the IPW weight as used in (2(2) ${\hat{\tau }}_{1,0\text{ipw}}=\frac{1}{\sum _{s=1}^n(1-A_s)}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}\frac{\sum _{s=1}^n(1-A_s)/n}{\sum _{s=1}^n{A}_s/n}{Y}_i=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}{Y}_i,$(2) ), and this leads to the consistency of ${\hat{\tau }}_{1,0\text{mr}}$: $\begin{array}{rl}{\hat{\tau }}_{1,0\text{mr}}=\sum _{\{i:A_i=0\}}{\hat{w}}_{\text{mr},i}{Y}_i& ⟹pE[(E\{\frac{1}{\pi \left(\mathit{X}\right)}\mid A=1\}-1)\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}Y\mid A=0]\\ & =E[\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}{Y}^0\mid A=0]\frac{P(A=0)}{P(A=1)}=E(Y^0\mid A=1),\end{array}$ where the second last equality follows from $E[\frac{1-\pi \left(\mathit{X}\right)}{\pi \left(\mathit{X}\right)}\mid A=1]=\frac{P(A=0)}{P(A=1)}$ that can be shown by taking $b\left(\mathit{X}\right)=\{1-\pi (\mathit{X}\left)\right\}/\pi \left(\mathit{X}\right)$ in (5(5) $E\left\{b\right(\mathit{X})\mid A=1\}=E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}b\right(\mathit{X})\mid A=0\}\frac{P(A=0)}{P(A=1)}.$(5) ), and the last equality follows from (1(1) $f(Y^1,Y^0,\mathit{X}\mid A=1)=\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\frac{P(A=0)}{P(A=1)}f(Y^1,Y^0,\mathit{X}\mid A=0).$(1) ).

A difference between the constraints (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) and those used by Han and Wang (2013) and Wang (2019) for estimating ATE is that in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) it is $1/{\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$ to be balanced whereas in Han and Wang (2013) and Wang (2019) it is ${\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$. The reason for this difference is that, for consistent estimation of ATT, the calibration weight ${\hat{w}}_{\text{mr}}$ needs to implicitly be of the form $\pi \left(\mathit{X}\right)/\{1-\pi (\mathit{X}\left)\right\}$ as in (2(2) ${\hat{\tau }}_{1,0\text{ipw}}=\frac{1}{\sum _{s=1}^n(1-A_s)}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}\frac{\sum _{s=1}^n(1-A_s)/n}{\sum _{s=1}^n{A}_s/n}{Y}_i=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=0\}}\frac{\hat{\pi }\left({\mathit{X}}_i\right)}{1-\hat{\pi }\left({\mathit{X}}_i\right)}{Y}_i,$(2) ) whereas, for consistent estimation of ATE, the calibration weight needs to implicitly be of the IPW form $1/\pi \left(\mathit{X}\right)$. From the above algebra it is seen that balancing $1/{\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$ in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) ensures $\pi \left(\mathit{X}\right)/\{1-\pi (\mathit{X}\left)\right\}$ implicitly appear, whereas the algebra in Han and Wang (2013) shows that balancing ${\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$ ensures $1/\pi \left(\mathit{X}\right)$ implicitly appear. Thus, consistency of our estimator for ATT would not hold if in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) $1/{\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$ is replaced by ${\pi }^{\left(j\right)}\left({\mathit{\alpha }}^{\left(j\right)}\right)$.

When ${\mathcal{D}}_0$ contains a correctly specified model for $E(Y^0\mid \mathit{X})$, say $d_0^{\left(1\right)}\left({\mathit{\gamma }}^{\left(1\right)}\right)$ without loss of generality, the consistency of ${\hat{\tau }}_{1,0\text{mr}}$ follows from $\begin{array}{rl}\sum _{\{i:A_i=0\}}{\hat{w}}_{\text{mr},i}{Y}_i& =\sum _{\{i:A_i=0\}}{\hat{w}}_{\text{mr},i}\{Y_i-d_{0i}^{\left(1\right)}({\hat{\mathit{\gamma }}}^{\left(1\right)}\left)\right\}+\frac{1}{\sum _{s=1}^n{A}_s}\sum _{i=1}^n{A}_i{d}_{0i}^{\left(1\right)}\left({\hat{\mathit{\gamma }}}^{\left(1\right)}\right)\\ & =\frac{1}{\sum _{s=1}^n(1-A_s)}\sum _{\{i:A_i=0\}}\frac{Y_i-d_{0i}^{\left(1\right)}\left({\hat{\mathit{\gamma }}}^{\left(1\right)}\right)}{1+{\hat{\mathit{\rho }}}^{\mathrm{\top }}{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})}+E\left\{E\right(Y^0\mid \mathit{X})\mid A=1\}+o_p\left(1\right)\\ & =E\{\frac{Y^0-E(Y^0\mid \mathit{X})}{1+{\mathit{\rho }}_{\ast }^{\mathrm{\top }}{\mathit{g}}_{\ast }({\mathit{\alpha }}_{\ast },{\mathit{\gamma }}_{\ast })}\mid A=0\}+E\left\{E\right(Y^0\mid \mathit{X},A=1)\mid A=1\}+o_p\left(1\right)\\ & =E[E\{\frac{Y^0-E(Y^0\mid \mathit{X},A=0)}{1+{\mathit{\rho }}_{\ast }^{\mathrm{\top }}{\mathit{g}}_{\ast }({\mathit{\alpha }}_{\ast },{\mathit{\gamma }}_{\ast })}\mid \mathit{X},A=0\}\mid A=0]+E(Y^0\mid A=1)+o_p\left(1\right)\\ & =E(Y^0\mid A=1)+o_p\left(1\right),\end{array}$ where ${\mathit{\rho }}_{\ast }$ and ${\mathit{g}}_{\ast }({\mathit{\alpha }}_{\ast },{\mathit{\gamma }}_{\ast })$ are the probability limits of $\hat{\mathit{\rho }}$ and $\stackrel{\mathbf{ˆ}}{\mathit{g}}(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})$, respectively.

In summary, we have the following result on the multiple robustness property for the proposed estimator of the ATT.

Proposition

If $\mathcal{P}$ contains a correctly specified model for $\pi \left(\mathit{X}\right)$ or ${\mathcal{D}}_0$ contains a correctly specified model for $E\left\{Y^0|\mathit{X}\right\}$, then $\frac{1}{\sum _{i=1}^n{A}_i}\sum _{i=1}^n{A}_i{Y}_i-\sum _{\{i:A_i=0\}}{\hat{w}}_i{Y}_i⟹pE(Y^1|A=1)-E(Y^0|A=1)$ as $n\to \mathrm{\infty }$.

The proposed estimator ${\hat{\tau }}_{1,0\text{mr}}$ and the EB estimator ${\hat{\tau }}_{1,0\text{eb}}$ are both weighted averages of the observed outcomes for the control group subjects, and both are calibration-type estimators originally considered in survey sampling (Deville & Särndal, 1992) that were later extended to missing data analysis and causal inference (e.g., Chan & Yam, 2014; S. Chen & Haziza, 2017; Han & Wang, 2013; Kim, 2010; Kim & Park, 2010; Qin & Zhang, 2007; Tan, 2010; Wu & Sitter, 2001; Zhang et al., 2022). Some of the constraints in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) based on the $\mathit{h}\left(\mathit{X}\right)$ component in $\stackrel{\mathbf{ˆ}}{\mathit{g}}(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})$ are actually the constraints (3(3) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i\mathit{h}\left({\mathit{X}}_i\right)=\frac{1}{\sum _{s=1}^n{A}_s}\sum _{\{i:A_i=1\}}\mathit{h}\left({\mathit{X}}_i\right)$(3) ) used by the EB method, and thus our proposed method achieves the same degree of covariate balancing between the treatment and the control groups. The other constraints in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) based on parametric models are for multiple robustness purpose.

Note that the constraints in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) are for estimating ATT and are determined by (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ), which balance the control group with the treatment group. This is different from estimating ATE where the constraints balance the control group with the full sample or the treatment group with the full sample (e.g., Fan et al., 2023). Adding the latter type of constraints into (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) would not work, since the weight matching control to treatment is different from the weight matching control to the full sample.

Standard error of the proposed estimator is needed to make inference, such as constructing confidence intervals. Due to the presence of multiple models and the lack of knowledge of which one is correct, deriving the asymptotic distribution of the proposed estimator as an approximation to the finite sample distributions is challenging. Therefore, we recommend using bootstrap to calculate the standard error. The excellent performance of the bootstrap method for multiply robust estimators in missing data context has been demonstrated through comprehensive simulation studies (e.g., Han, 2014a, 2014b), and its effectiveness for the proposed estimator will be shown in the next section.

4. Simulation studies

In this section we conduct simulation studies to evaluate the finite sample performance of the proposed multiply robust estimator of the ATT. The simulation setting mimics that in Han and Wang (2013). The data are generated in the following way: $X\sim \text{Uniform}(-2.5,2.5)$, $Y^0\mid X\sim N\left\{d_0\right(X),4X^2+2\}$, $Y^1\mid X\sim N\left\{d_1\right(X),4X^2+2\}$ and $A\mid X\sim \text{Bernoulli}\left\{\pi \right(X\left)\right\}$, where $\pi \left(X\right)$ is either $\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$ or $1-\exp [-\exp \{0.5+0.5X-0.3\exp \left(X\right)\left\}\right]$ and $\left\{d_0\right(X),d_1(X\left)\right\}$ is either $(1+2X+3X^2,2+3X+3X^2)$ or $\{1+2X+3\exp (X),2+3X+3\exp (X\left)\right\}$. Therefore, we have in total four data generating processes. For each of them, we postulate two propensity score models ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)=\{1+\exp ({\alpha }_1^{\left(1\right)}+{\alpha }_2^{\left(1\right)}X+{\alpha }_3^{\left(1\right)}{X}^2){\}}^{-1}$ and ${\pi }^{\left(2\right)}\left({\mathit{\alpha }}^{\left(2\right)}\right)=1-\exp [-\exp \{{\alpha }_1^{\left(2\right)}+{\alpha }_2^{\left(2\right)}X+{\alpha }_3^{\left(2\right)}\exp \left(X\right)\left\}\right]$ and two regression models for $d_0\left(X\right)$, $d_0^{\left(1\right)}\left({\mathit{\gamma }}^{\left(1\right)}\right)={\gamma }_1^{\left(1\right)}+{\gamma }_2^{\left(1\right)}X+{\gamma }_3^{\left(1\right)}{X}^2$ and $d_0^{\left(2\right)}\left({\mathit{\gamma }}^{\left(2\right)}\right)={\gamma }_1^{\left(2\right)}+{\gamma }_2^{\left(2\right)}X+{\gamma }_3^{\left(2\right)}\exp \left(X\right)$.

Tables 1 and 2 contain some simulation results summarized based on 1000 replications for n = 300 and n = 800, respectively. To show the flexibility of the proposed procedure in providing a unified framework for constructing estimators, both our proposed estimators based on each possible combination of models from $\left\{{\pi }^{\left(1\right)}\right({\mathit{\alpha }}^{\left(1\right)}),{\pi }^{\left(2\right)}({\mathit{\alpha }}^{\left(2\right)}),d_0^{\left(1\right)}({\mathit{\gamma }}^{\left(1\right)}),d_0^{\left(2\right)}({\mathit{\gamma }}^{\left(2\right)}\left)\right\}$ and the corresponding OR, IPW and/or AIPW estimators, when available, are listed for comparison. Each digit of the four-digit number inside an estimator's name, from left to right, indicates if ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)$, ${\pi }^{\left(2\right)}\left({\mathit{\alpha }}^{\left(2\right)}\right)$, $d_0^{\left(1\right)}\left({\mathit{\gamma }}^{\left(1\right)}\right)$ and $d_0^{\left(2\right)}\left({\mathit{\gamma }}^{\left(2\right)}\right)$ are used, respectively. For example, MR-1101 denotes our proposed multiply robust estimator based on models ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)$, ${\pi }^{\left(2\right)}\left({\mathit{\alpha }}^{\left(2\right)}\right)$ and $d_0^{\left(2\right)}\left({\mathit{\gamma }}^{\left(2\right)}\right)$. Compared to the OR, IPW and AIPW estimators, our proposed estimators using the same parametric models have similar numerical performance. The multiple robustness property is well demonstrated by the negligible biases, especially with n = 800, of the proposed estimators MR-1100, MR-0011, MR-1110, MR-1101, MR-1011, MR-0111 and MR-1111, which are consistent under all four data generating processes. The numerical performance of MR-1100 with a small sample size may be occasionally unstable, as shown by the relatively large bias when n = 300 with data generated based on $\pi \left(X\right)=\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$, due to a high correlation between ${\pi }^{\left(1\right)}\left({\hat{\mathit{\alpha }}}^{\left(1\right)}\right)$ and ${\pi }^{\left(2\right)}\left({\hat{\mathit{\alpha }}}^{\left(2\right)}\right)$, but it improves dramatically as the sample size increases to n = 800.

Table 1. Comparison of different estimators based on 1000 replications and n = 300. Each digit of the four-digit number inside an estimator's name, from left to right, indicates if ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)$, ${\pi }^{\left(2\right)}\left({\mathit{\alpha }}^{\left(2\right)}\right)$, $d_0^{\left(1\right)}\left({\mathit{\gamma }}^{\left(1\right)}\right)$ and $d_0^{\left(2\right)}\left({\mathit{\gamma }}^{\left(2\right)}\right)$ are used, respectively. The results have been multiplied by 100.

	${\pi }^{\left(1\right)}$						${\pi }^{\left(2\right)}$
	$(d_0^{\left(1\right)},d_1^{\left(1\right)})$			$(d_0^{\left(2\right)},d_1^{\left(2\right)})$			$(d_0^{\left(1\right)},d_1^{\left(1\right)})$			$(d_0^{\left(2\right)},d_1^{\left(2\right)})$
	rbias	rmse	mae	rbias	rmse	mae	rbias	rmse	mae	rbias	rmse	mae
IPW-1000	−6	99	60	−9	79	53	11	38	26	10	38	26
IPW-0100	69	93	64	−42	68	43	−1	44	28	−1	37	24
OR-0010	−8	57	37	−134	99	68	−2	35	23	−25	47	32
OR-0001	224	131	114	−7	51	34	−19	50	34	−2	36	24
AIPW-1010	−9	65	42	−16	80	54	−2	35	23	−34	51	39
AIPW-1001	−3	85	56	−8	65	41	34	53	38	−2	35	23
AIPW-0110	−9	62	40	−69	82	53	−2	35	23	−2	40	25
AIPW-0101	63	84	56	−8	61	39	−1	42	27	−2	35	24
MR-1000	−1	71	46	−10	76	50	14	40	27	12	39	26
MR-0100	51	77	52	−59	73	48	−1	49	32	−1	39	25
MR-0010	−8	57	37	−158	108	80	−2	35	23	−30	49	35
MR-0001	229	134	116	−6	52	34	21	47	32	−2	36	24
MR-1100	−24	76	49	−56	86	52	−3	39	26	−8	41	27
MR-1010	−10	65	42	−17	73	50	−1	35	23	−27	46	32
MR-1001	−1	73	50	−9	65	42	29	51	35	−2	35	23
MR-0110	−10	64	41	−68	80	52	−2	35	24	−2	38	24
MR-0101	60	79	53	−8	61	39	−2	40	26	−2	35	24
MR-0011	−13	83	51	−12	80	43	−2	35	24	−1	38	26
MR-1110	−9	65	42	−22	72	47	−2	36	24	−1	36	25
MR-1101	4	70	46	−9	65	42	−3	38	26	−2	35	24
MR-1011	−11	73	45	−10	73	44	−2	36	24	−2	36	24
MR-0111	−11	77	47	−10	78	48	−1	36	24	−2	35	24
MR-1111	−10	75	46	−10	75	46	−2	36	24	−2	35	24
EB-1	523	282	266	352	194	177	−317	321	311	−275	279	272
EB-2	−11	65	43	−19	71	46	−2	35	23	−31	49	36
MR-1	231	166	128	197	136	107	−318	323	313	−276	280	272
MR-2	−12	91	56	57	99	64	−1	34	23	−36	52	39

${\pi }^{\left(1\right)}$: $\pi \left(X\right)=\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$. ${\pi }^{\left(2\right)}$: $\pi \left(X\right)=1-\exp [-\exp \{0.5+0.5X-0.3\exp \left(X\right)\left\}\right]$. $(d_0^{\left(1\right)},d_1^{\left(1\right)})$: $\left\{d_0\right(X),d_1(X\left)\right\}=(1+2X+3X^2,2+3X+3X^2)$. $(d_0^{\left(2\right)},d_1^{\left(2\right)})$: $\left\{d_0\right(X),d_1(X\left)\right\}=\{1+2X+3\exp (X),2+3X+3\exp (X\left)\right\}$. rbias: relative bias, bias divided by the true value. rmse: root mean square error. mae: median absolute error. IPW: inverse probability weighting. OR: outcome regression. AIPW: augmented inverse probability weighting. MR: multiply robust. EB: entropy balancing. EB-1: EB with $\mathit{h}\left(X\right)=X$. EB-2: EB with $\mathit{h}\left(X\right)=(X,X^2)$. MR-1: MR with $\mathit{h}\left(X\right)=X$ and no working models. MR-2: MR with $\mathit{h}\left(X\right)=(X,X^2)$ and no working models.

Table 2. Comparison of different estimators based on 1000 replications and n = 800. Each digit of the four-digit number inside an estimator's name, from left to right, indicates if ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)$, ${\pi }^{\left(2\right)}\left({\mathit{\alpha }}^{\left(2\right)}\right)$, $d_0^{\left(1\right)}\left({\mathit{\gamma }}^{\left(1\right)}\right)$ and $d_0^{\left(2\right)}\left({\mathit{\gamma }}^{\left(2\right)}\right)$ are used, respectively. The results have been multiplied by 100.

	${\pi }^{\left(1\right)}$						${\pi }^{\left(2\right)}$
	$(d_0^{\left(1\right)},d_1^{\left(1\right)})$			$(d_0^{\left(2\right)},d_1^{\left(2\right)})$			$(d_0^{\left(1\right)},d_1^{\left(1\right)})$			$(d_0^{\left(2\right)},d_1^{\left(2\right)})$
	rbias	rmse	mae	rbias	rmse	mae	rbias	rmse	mae	rbias	rmse	mae
IPW-1000	2	58	39	1	46	32	14	26	18	13	26	17
IPW-0100	84	66	50	−32	41	28	0	27	18	0	23	15
OR-0010	2	35	25	−127	76	62	0	22	15	−24	34	26
OR-0001	244	126	119	2	32	22	−18	34	23	0	23	15
AIPW-1010	3	40	27	0	46	30	0	22	14	−33	40	33
AIPW-1001	4	52	35	3	40	27	37	44	35	0	22	15
AIPW-0110	3	38	26	−57	52	36	0	22	15	0	25	16
AIPW-0101	75	60	43	3	37	25	0	26	17	0	22	15
MR-1000	5	44	30	1	46	32	16	28	19	15	27	18
MR-0100	63	54	39	−50	46	32	0	30	19	0	24	16
MR-0010	2	35	25	−152	86	75	0	22	14	−28	37	29
MR-0001	246	127	120	2	32	22	24	35	26	0	23	15
MR-1100	−6	46	32	−18	49	32	−1	24	16	−3	24	15
MR-1010	3	40	28	0	43	30	0	22	15	−26	34	26
MR-1001	5	46	31	3	40	28	32	41	32	0	22	15
MR-0110	3	39	28	−56	51	35	0	22	15	0	23	15
MR-0101	73	57	41	3	37	25	0	25	16	0	22	15
MR-0011	5	55	36	8	68	30	0	22	15	1	24	15
MR-1110	3	40	27	−1	42	28	0	22	15	0	22	15
MR-1101	6	44	30	3	40	28	0	24	15	0	22	15
MR-1011	4	42	29	4	42	29	0	22	15	0	22	15
MR-0111	4	46	31	7	61	31	0	22	15	0	22	15
MR-1111	4	42	29	4	42	29	0	22	15	0	22	15
EB-1	564	281	277	383	193	189	−322	316	314	−279	274	271
EB-2	2	40	28	−1	42	29	0	22	15	−30	38	29
MR-1	251	142	126	220	121	112	−323	316	314	−279	274	272
MR-2	2	73	42	92	81	52	0	21	14	−35	41	34

${\pi }^{\left(1\right)}$: $\pi \left(X\right)=\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$. ${\pi }^{\left(2\right)}$: $\pi \left(X\right)=1-\exp [-\exp \{0.5+0.5X-0.3\exp \left(X\right)\left\}\right]$. $(d_0^{\left(1\right)},d_1^{\left(1\right)})$: $\left\{d_0\right(X),d_1(X\left)\right\}=(1+2X+3X^2,2+3X+3X^2)$. $(d_0^{\left(2\right)},d_1^{\left(2\right)})$: $\left\{d_0\right(X),d_1(X\left)\right\}=\{1+2X+3\exp (X),2+3X+3\exp (X\left)\right\}$. rbias: relative bias, bias divided by the true value. rmse: root mean square error. mae: median absolute error.IPW: inverse probability weighting. OR: outcome regression. AIPW: augmented inverse probability weighting. MR: multiply robust. EB: entropy balancing. EB-1: EB with $\mathit{h}\left(X\right)=X$. EB-2: EB with $\mathit{h}\left(X\right)=(X,X^2)$. MR-1: MR with $\mathit{h}\left(X\right)=X$ and no working models. MR-2: MR with $\mathit{h}\left(X\right)=(X,X^2)$ and no working models.

The estimator EB-1 by balancing only the first moment of X has a very large bias under all four data generating processes. Replacing the exponential tilting with the empirical likelihood, the estimator MR-1 has a significantly improved performance compared to EB-1 when the data are generated based on $\pi \left(X\right)=\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$. By balancing the first two moments of X, the EB method implicitly fits linear regression models for $\text{logit}\left\{\pi \right(X\left)\right\}$ and $d_0\left(X\right)$ with regressors X and $X^2$, and this makes the estimator EB-2 consistent under the first three data generating processes and inconsistent under the last one where $\pi \left(X\right)=1-\exp [-\exp \{0.5+0.5X-0.3\exp \left(X\right)\left\}\right]$ and $\left\{d_0\right(X),d_1(X\left)\right\}=\{1+2X+3\exp (X),2+3X+3\exp (X\left)\right\}$. This theoretical conclusion is well confirmed by inspecting the bias of EB-2. The estimator MR-2 with exponential tilting replaced by empirical likelihood is inconsistent under all four data generating processes.

Table 3 contains a summary of the performance of the bootstrap method for standard error calculation. Due to similarity of the performance under different settings, we only include the results for $\pi \left(X\right)=\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$ and $\left\{d_0\right(X),d_1(X\left)\right\}=(1+2X+3X^2,2+3X+3X^2)$ with n = 300. The results are summarized based on 1000 replications with the bootstrap resampling size 200. It is seen that the average of the standard errors based on boostrap is very close to the empirical standard error. In addition, for the estimators that are consistent under this data generating process, i.e., MR-1000, MR-0010, MR-1100, MR-1010, MR-1001, MR-0110, MR-0011, MR-1110, MR-1101, MR-1011, MR-0111 and MR-1111, the coverage percentage of the $95\mathrm{\%}$ confidence intervals constructed based on bootstrap standard errors is close to the nominal level. Both observations show that the bootstrap method is reliable to calculate the standard errors for the proposed estimators.

Table 3. Performance of the bootstrap method in calculating the standard errors for the proposed estimators in the setting of $\pi \left(X\right)=\{1+\exp (0.8+0.5X-0.3X^2){\}}^{-1}$ and $\left\{d_0\right(X),d_1(X\left)\right\}=(1+2X+3X^2,2+3X+3X^2)$. Results are summarized based on 1000 replications with n = 300 and the bootstrap resampling size is 200. Each digit of the four-digit number inside an estimator's name, from left to right, indicates if ${\pi }^{\left(1\right)}\left({\mathit{\alpha }}^{\left(1\right)}\right)$, ${\pi }^{\left(2\right)}\left({\mathit{\alpha }}^{\left(2\right)}\right)$, $d_0^{\left(1\right)}\left({\mathit{\gamma }}^{\left(1\right)}\right)$ and $d_0^{\left(2\right)}\left({\mathit{\gamma }}^{\left(2\right)}\right)$ are used, respectively.

	bias	se-emp	se-bp	cp-$95\mathrm{\%}$
MR-1000	0.05	0.71	0.67	92.50
MR-0100	0.31	0.72	0.69	89.00
MR-0010	−0.01	0.56	0.58	95.30
MR-0001	1.21	0.65	0.64	52.60
MR-1100	−0.08	0.77	0.74	94.10
MR-1010	0.01	0.64	0.62	93.60
MR-1001	0.05	0.74	0.70	92.40
MR-0110	0.00	0.62	0.61	94.00
MR-0101	0.36	0.74	0.70	87.50
MR-0011	0.03	0.81	0.75	93.30
MR-1110	0.01	0.63	0.61	93.60
MR-1101	0.08	0.70	0.67	92.00
MR-1011	0.03	0.72	0.70	94.40
MR-0111	0.03	0.76	0.73	93.90
MR-1111	0.04	0.76	0.73	94.10

se-emp: empirical standard error. se-bp: averaged bootstrap standard error. cp-$95\mathrm{\%}$: coverage percentage of the $95\mathrm{\%}$ confidence interval constructed based on the bootstrap standard error. MR: multiply robust.

5. Discussion

In this paper we have proposed a multiply robust estimator for the causal effect ATT. The estimator provides more protection on estimation consistency compared to existing doubly robust estimators. It is worth pointing out that, although the proposed method can simultaneously account for multiple working models, there does not have to be multiple models to apply this method. The construction of constraints is very flexible and the method can be applied when only a single working model is available. Therefore, the proposed method provides a unified framework as an alternative to various existing methods, including the IPW and AIPW methods.

A major difference between the proposed method and the EB method is the objective function being optimized. The EB method minimizes the Shannon entropy $\sum _{\{i:A_i=0\}}{w}_i\log w_i$, or equivalently the exponential tilting function in the empirical likelihood literature (e.g., Newey & Smith, 2004), whereas we proposed to maximize the empirical likelihood $\prod _{\{i:A_i=0\}}{w}_i$. The advantage of considering an empirical likelihood is manifold. First, it improves the robustness of estimation consistency by using constraints based on parametric working models. Although the same constraints can be used by the EB method as well, a correct model does not make the EB estimator consistent, and consistency of the EB estimator is achieved only when $\text{logit}\left\{\pi \right(\mathit{X}\left)\right\}$ or $E(Y^0\mid \mathit{X})$ is a linear model with regressors $\mathit{h}\left(\mathit{X}\right)$. Second, by dividing the constraints into two parts, one based on parametric working models and one based on moments of $\mathit{X}$, the proposed method is more flexible when the dimension of $\mathit{X}$ is large. With $\pi \left(\mathit{X}\right)$ and/or $E(Y^0\mid \mathit{X})$ depending on a high dimensional $\mathit{X}$, consistency of the EB estimator requires a large number of constraints to include all the regressors for $\text{logit}\left\{\pi \right(\mathit{X}\left)\right\}$ and/or $E(Y^0\mid \mathit{X})$, and this may jeopardize the numerical performance in practice. The proposed method, on the contrary, has a separate model building step where complex working models can be built and only the ones that are deemed to be close to the truth are used in the constraints. In this case, the parametric working models help to achieve consistency, and thus the other part of constraints based on moments of $\mathit{X}$ can be chosen exclusively for the goal of achieving covariate balancing. This may considerably reduce the number of constraints compared to the EB method and thus improve the numerical performance. Third, it is well known in the empirical likelihood literature (e.g., Newey & Smith, 2004) that the empirical likelihood has smaller higher-order bias, which translates to a better numerical performance under a finite sample size, compared to exponential tilting and other alternatives. As an empirical version of the moment equality (7(7) $E\left\{\frac{\pi \left(\mathit{X}\right)}{1-\pi \left(\mathit{X}\right)}\right[b\left(\mathit{X}\right)-E\left\{b\right(\mathit{X})\mid A=1\}]\mid A=0\}=0.$(7) ), the constraints in (8(8) $w_i>0,\sum _{\{i:A_i=0\}}{w}_i=1,\sum _{\{i:A_i=0\}}{w}_i{\stackrel{\mathbf{ˆ}}{\mathit{g}}}_i(\stackrel{\mathbf{ˆ}}{\mathit{\alpha }},\stackrel{\mathbf{ˆ}}{\mathit{\gamma }})=\mathbf{0},$(8) ) are legitimate, at least when the sample size is large. Therefore, the ‘model misspecification’ problem that typically exists in the empirical likelihood literature under which the performance of empirical likelihood may be worse than that of the exponential tilting (Schennach, 2007) is not of a concern here.

Acknowledgments

We thank the Editor, Associate Editor and a referee for their valuable comments that have helped improve the quality of this paper.

Disclosure statement

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