Outlier-Resistant Estimators for Average Treatment Effect in Causal Inference

Let $(Y,T,X)$ be the observable random variables; $X$ is the outcome, $T$ is the treatment, and $X$ is the confounder. We assume that $Y$ is continuous and $T$ is binary; it is easy to extend $T$ to multiple discrete treatments. We have the observations $(Y_{i},T_{i},X_{i})_{i=1}^{n}$ drawn from the distribution of $(Y,T,X)$ in an i.i.d. manner. Denote the potential outcome under $T=t$ by $Y^{(t)}$ and let $\mu^{(t)}=\mathbb{E}[Y^{(t)}]$. $Y^{(t)}$ is uniquely defined for every treatment as a random variable, namely, well-defined. Note that i.i.d. sampling and well-definedness of the potential outcome are collectively called the stable unit treatment value assumption (SUTVA; Imbens and Rubin, 2015). The ATE is defined as $\mu^{(1)}-\mu^{(0)}$. The ATE cannot be estimated directly since we cannot observe $Y^{(1)}$ and $Y^{(0)}$ simultaneously; instead, we use the observed variables under the common assumptions (e.g. Imbens and Rubin, 2015):

1. 1.

   Conditional Unconfoundedness: $Y^{(t)}\perp\!\!\!\perp T|X$ for all $t\in\{0,1\}$,

2. 2.

   Consistency: $Y=Y^{(t)}$ if $T=t$,

3. 3.

   Positivity: $P(T=1|X)>c$ for some constant $c>0$.

The ATE can be estimated from the observed variables under these assumptions, namely, identifiable. Hereafter, we assume the triple assumption holds and focus on the estimation of $\mu^{(1)}$ for simplicity. $\mu^{(0)}$ is estimated in a similar way, and the ATE is estimated by the difference between the estimates of $\mu^{(1)}$ and $\mu^{(0)}$.

We introduce three consistent estimators of the potential mean. The IPW estimator (Rosenbaum and Rubin, 1983) is based on the propensity score (PS). Let $\pi(x;\alpha)\in(0,1)$ be the PS, which models $P(T=1|x)$. We assume the PS is correctly specified, in other words, there exists $\alpha^{*}$ such that $\pi(x;\alpha^{*})=P(T=1|x)$ for every $x$. The IPW estimator has several forms (Lunceford and Davidian, 2004), but we use the weighted average form: $\hat{\mu}^{(1)}_{IPW}=\left.\left(\sum_{i=1}^{n}T_{i}Y_{i}/\pi(X_{i};\hat{\alpha})\right)\right/\left(\sum_{i=1}^{n}T_{i}/\pi(X_{i};\hat{\alpha})\right)$, where $\hat{\alpha}$ is an estimate of $\alpha$ obtained in a consistent manner such as the maximum likelihood estimation (MLE). The IPW estimator is consistent with $\mu^{(1)}$ if the PS model is correctly specified. The IPW estimator can be seen as the root of the following estimating equation:

Outcome regression (OR) is also popular. To construct the OR estimator, we model $\mathbb{E}[Y|T=1,X]$ by some function $m_{1}(X;\beta)$. Then, the OR estimator is obtained as $n^{-1}\sum_{i=1}^{n}m_{1}(X_{i};\hat{\beta})$, where $\hat{\beta}$ is a consistent estimate of $\beta$. The IPW and OR estimators are asymptotically consistent with $\mu^{(1)}$ when the model used in each estimator is correctly specified; the consistency is not assured if the model is misspecified. The DR estimator (Scharfstein, Rotnitzky, and Robins, 1999; Bang and Robins, 2005) combines the IPW and OR estimators. Since the DR estimator is consistent with $\mu^{(1)}$ if either the PS or OR model is correctly specified, it is said to be “doubly robust.” Besides, if both models are correctly specified, the DR estimator is semiparametrically efficient (Robins and Rotnitzky, 1995; Tsiatis, 2006). Although there are many estimators equipped with double robustness, we refer the root of the following estimating equation as the DR estimator $\hat{\mu}^{(1)}_{DR}$, which is a special case of the augmented IPW estimator:

Let $\sum_{i=1}^{n}\psi(Y_{i},\theta)=0$ be an estimating equation, where $\psi$ is a known vector-valued map, and $\theta$ is the parameter of interest. An estimator $\hat{\theta}$ solving the estimating equation is called an M-estimator. M-estimator is a large class of estimators involving MLE, IPW, OR, and DR. If the estimating equation is unbiased, say $\mathbb{E}_{\theta}[\psi(Y,\theta)]=0$, the M-estimator is consistent with the truth under some regularity conditions (e.g. Chap. 5 of Van der Vaart, 2000).

By replacing $Y_{i}-\mu$ in (2) with an estimating function $\psi(Y_{i};\theta)$, the IPW and DR estimators can be expanded to a general M-estimator. If we are interested in the same parameter $\theta$ with respect to $Y^{(1)}$, the IPW and DR M-estimators (Tsiatis, 2006) are available:

The conditional expectation $\mathbb{E}_{\hat{q}}[\psi(Y_{i};\theta)|T=1,X_{i}]$ is calculated using the parametric OR model $q(y|T=1,x;\hat{\beta})$ via direct calculation or Monte-Carlo approximation (Hoshino, 2007). When the original M-estimating equation is unbiased, the IPW/DR estimating equations are unbiased under the proper model specification. The asymptotic properties of the IPW and DR M-estimators follow from the standard theory of M-estimators.

This section provides a brief review of the outlier-resistant estimation of the mean in a one-variable and non-causal setting. Let $\tilde{g}$ be the density function of a random variable $Z\in\mathbb{R}$. Assume that the density is contaminated as $\tilde{g}(z)=(1-\varepsilon)f_{\mu^{*}}(z)+\varepsilon\delta(z)$, where $f_{\mu^{*}}$ is the density of $Z$ without contamination equipped with the mean $\mu^{*}$, $\varepsilon$ is the contamination ratio, and $\delta$ is the density of outliers. Our goal is to estimate $\mu^{*}$ from i.i.d. observations $\{Z_{1},...,Z_{n}\}$ drawn from $\tilde{g}$. If we model the contamination in this way, the sample mean converges to $(1-\varepsilon)\mu^{*}+\varepsilon\mathbb{E}_{\delta}[Z]$; if the mean of outliers is far from $\mu^{*}$, the sample mean is asymptotically biased. To deal with contamination, many types of M-estimators are applied. The unbiasedness of the estimating equation does not usually hold under contamination because

By designing $\psi$ to eliminate or bound $\mathbb{E}_{\delta}[\psi(Z,\mu^{*})]$, the influence of outliers can be reduced. Let $\theta_{\psi}^{*}$ denote a root of $\mathbb{E}_{\tilde{g}}[\psi(Z,\theta)]=0$; then, the latent bias is defined as $\theta_{\psi}^{*}-\theta^{*}$. If $\delta$ is Dirac’s delta and $\varepsilon$ is sufficiently small, the latent bias is approximated by the IF. The IF-based discussion in Section 4 provides some insights into the outlier resistance of the estimators when the contamination ratio is small. The latent bias and M-estimators are discussed in detail elsewhere (e.g. Huber, 2004; Fujisawa, 2013; Fujisawa and Eguchi, 2008).

Next, we move to a causal setting. In this paper, we consider the vertical outliers. In other words, we assume that only the outcome $Y$ may be contaminated, and that the contamination does not affect the causal mechanism among $(Y,T,X)$. A typical example is contamination of laboratory values in medical research with foreign substances. Let $\delta_{Y|TX}$ be the conditional density of outliers given $(T,X)$, and let $\varepsilon_{t}(x)$ be the contamination ratio. Then, the contaminated conditional density given $(T,X)$ is defined as

where $g$ without tilde denotes the density without contamination; the tilde indicates that the distribution is contaminated. For simplifying the notations, we often drop the subscripts of density functions as long as there would be no confusion and write $\delta_{t}(y|x)=\delta_{Y|TX}(y|t,x)$ below. The contamination ratio and their density depend on the treatment $T$ and the confounder $X$. Since we estimate $\mu^{(t)}$ for each treatment separately, the dependence on $T$ is tractable. In contrast, the dependence on $X$ is critical in our analysis. The $X$-dependent contamination is referred to as heterogeneous contamination. We also discuss the special case in which $\varepsilon$ and $\delta$ are not dependent on $X$, called homogeneous contamination. Note that we do not assume $\varepsilon_{t}(x)$ to be small enough to be negligible, except in Section 4.

We are interested in the marginal mean of $Y^{(1)}$. Let $f_{Y^{(1)}}(y;\mu^{(1)})$ be the true marginal density of $Y^{(1)}$. It is obtained by integrating $X$ out from $g_{Y|TX}(y|T,X)$ under $T=1$:

The second equality holds from the triple assumption in Section 2.1. We often write $f_{Y^{(1)}}(y;\mu^{(1)})$ as $f_{1}(y)$ for simplicity.

Under contamination, the IPW estimating equation is severely biased even if the true PS is obtained as $\pi(X|\alpha^{*})=P(T=1|X)$:

It is found that the remaining term contains the expectation of $Y$ with respect to $\delta$, which implies the estimating equation is severely affected by outliers. The DR estimating equation is also biased. To estimate $\mu^{(1)}$ accurately, we have to remove the influence of contamination.

First, we introduce an extension of the IPW estimator, called the density power inverse probability weighting (DP-IPW) estimator. The DP-IPW estimator is defined as a root of the following estimating equation:

Under no contamination, the DP-IPW estimating equation is unbiased.

Now we consider the contaminated case. The bias of the DP-IPW estimating equation takes a different form from (8).

Next, we introduce the density power doubly robust (DP-DR) estimator. The DP-DR estimator is a straight application of the DR M-estimator and defined as a root of the following estimating equation:

As we have discussed in Section 2.1, $\mathbb{E}_{\hat{q}}\left[h(Y;\mu)^{\gamma}(Y-\mu)|T=1,X\right]$ is obtained by direct calculation or Monte Carlo approximation based on the parametric OR model $\hat{q}:=q(y|T=1,X;\hat{\beta})$. In the Appendix, we present the explicit forms of $\mathbb{E}_{\hat{q}}\left[h(Y;\mu)^{\gamma}(Y-\mu)|T=1,X\right]$ when $h$ and $q$ are assumed to be Gaussian. The parameter $\beta$ is usually estimated in an outlier-resistant manner: Huber regression (Huber, 2004, Chap.7), MM estimator (Yohai, 1987), density power regression (Basu et al., 1998; Kanamori and Fujisawa, 2015), and $\gamma$-regression (Fujisawa and Eguchi, 2008; Kawashima and Fujisawa, 2017), for example.

The DP-DR estimator is doubly robust under no contamination as with the general DR M-estimator.

Now, we evaluate the bias of the DP-DR estimating equation under contamination.

In the OR-correct case, the reason why DP-DR is biased under homogeneous contamination is as follows. Under Assumption 1, the expectation of the DP-DR estimating function becomes

where we denote the density power estimating function by $\psi$. In the last formula, it is found that the terms in the curly brackets do not cancel because the first term is reduced by $1-\varepsilon_{1}$. Based on this consideration, we propose a bias-corrected version of DP-DR, called the $\varepsilon$DP-DR estimator. $\varepsilon$DP-DR is designed to cancel the dominant bias under homogeneous contamination. The $\varepsilon$DP-DR estimator is a root of the following estimating equation:

where $\hat{\varepsilon}_{1}$ is a consistent estimator of the expected contamination ratio $\overline{\varepsilon}_{1}=\int\varepsilon_{1}(x)g(x)dx$. $\hat{\varepsilon}_{1}$ can be obtained simultaneously with the parametric OR model by the unnormalized modeling with the density power score (Kanamori and Fujisawa, 2015), for example. While DP-DR is a special case of the DR M-estimator, $\varepsilon$DP-DR goes beyond this framework by the bias correction. Under no contamination, the $\varepsilon$DP-DR estimating equation is asymptotically identical to the DP-DR estimating equation. The $\varepsilon$DP-DR estimating equation is also biased under heterogeneous contamination; however, the bias takes a different form.

Similar to (19), the second term of (21) is approximately zero if we assume that $\pi(\cdot;\alpha)$ is bounded away from 0 and 1.

We have proposed three types of outlier resistant semiparametric estimators: DP-IPW, DP-DR, and $\varepsilon$DP-DR. Table 1 shows the bias of the estimating equations under the conditions discussed above. Under heterogeneous contamination, all estimators are biased, but the biases are hardly influenced by the absolute value of outliers. Furthermore, as discussed in Section 4, outliers have negligible influence if the contamination ratio is sufficiently small. $\varepsilon$DP-DR improves DP-DR in the OR-correct case under homogeneous contamination, but we continue to discuss DP-DR for three reasons: the contamination ratio is sometimes hard to estimate, the bias (19) is not serious if $\pi(X;\alpha)$ is close to $P(T=1|X)$, and the simulation results presented in Section 6 indicate that DP-DR remains better than the existing methods even in the OR-correct case.

Since the proposed estimating equations cannot be solved explicitly, we develop an iterative algorithm. Various algorithms are available, but we propose a standard algorithm for M-estimators (Huber, 2004; Hampel et al., 2011). Detailed algorithm is available in the Appendix. Hereafter, we suppose $h$ and $q$ are Gaussian. We also provide explicit updating formulae in this case. Note that some additional parameters of $h$ should be estimated in a roughly unbiased and outlier-resistant way.

We simulated random observations based on a simple causal setting. The confounders $(X_{1},X_{2})$ were independently drawn from a Gaussian or uniform distribution with mean zero and unit variance. The treatment $T$ was assigned along with the conditional probability $P(T=1|X_{1},X_{2})$ that was defined as a sigmoid function of $0.8X_{1}+0.2X_{2}$. The potential outcomes $(Y^{(1)},Y^{(0)})$ were generated according to a linear function of $(X_{1},X_{2})$ with Gaussian error: $Y^{(1)}=\mu^{(1)}+1.2X_{1}+0.3X_{2}+e$ and $Y^{(0)}=\mu^{(0)}+1.2X_{1}+0.3X_{2}+e$. $\mu^{(1)}$ and $\mu^{(0)}$ were set to 3 and 0, respectively. The standard deviation (SD) of $e$ was set to $\sqrt{0.72}$; then, $\mathrm{SD}[Y^{(1)}]=\mathrm{SD}[Y^{(0)}]=1.5$. When the confounders were not Gaussian, the potential outcomes were not Gaussian. The observed outcome $Y$ was defined as $Y=TY^{(1)}+(1-T)Y^{(1)}$ under no contamination. Outliers were drawn from $\mathcal{N}(\mu^{(t)}+10\sigma^{(t)},1)$, with $\sigma^{(t)}=\mathrm{SD}[Y^{(t)}]=1.5$. For the homogeneous contamination settings, the contamination ratio was set to be a constant $\varepsilon_{t}$. For the heterogeneous contamination settings, the contamination ratio was set to be $1.5\varepsilon_{t}$ if $X_{1}+X_{2}\leq 0$ and $0.5\varepsilon_{t}$ if $X_{1}+X_{2}>0$. The average contamination ratio is set to $\varepsilon_{t}\in\{0,0.05,0.1,0.2\}$. Then, the observations of $Y$ were randomly replaced with outliers according to the contamination ratio. The sample size was fixed to $n=100$ throughout the Monte Carlo simulations. Furthermore, we generated datasets in which the outcome follows a symmetric and heavy-tailed distribution. We drew the error term of $Y^{(t)}$ from the standard Cauchy distribution instead of inserting outliers.

First, we performed a comparative study. The potential mean $\mu^{(1)}$ was estimated using the proposed and comparative methods. In this experiment, we used all settings illustrated in the previous section. The propensity score was estimated by logistic regression. The parametric OR was conducted in two ways: Gaussian MLE with non-outliers or unnormalized Gaussian modeling (the tuning parameter was set to 0.5) (Kanamori and Fujisawa, 2015). For the DR estimators, we investigated three patterns of model misspecification: PS-correct/OR-correct, PS-correct/OR-incorrect, and PS-incorrect/OR-correct. For the model-correct case, we included an intercept and $(X_{1},X_{2})$ as covariates. For the model-incorrect case, we included only an intercept and $X_{2}$. We performed 10,000 simulations for every setting and method. Tables 2 and 3 show the results of the comparative study when the covariates were Gaussian and the OR for the DR-type estimators was the Gaussian MLE with non-outliers. The estimation error was measured by the root mean square error (RMSE). The mean and SD of all estimates, the mean computation time, and the results for the other settings are provided in the Appendix. In Table 2, the naive IPW estimator had a significantly larger RMSE under contamination. Both the median-based methods and DP-IPW dramatically reduced the RMSE. As the contamination ratio increased, the RMSE increased. The RMSE tended to be larger for heterogeneous contamination than for homogeneous contamination. When the optimal $\gamma$ was chosen, the proposed method outperformed the comparative methods and had the smallest RMSE for all settings. Looking at Table 3, the results for the DR-type estimators were similar to those for the IPW estimators. The proposed method with a proper $\gamma$ outperformed the comparative methods and had the smallest RMSE in all settings. DP-DR and $\varepsilon$DP-DR performed similarly, although $\varepsilon$DP-DR was slightly superior in many settings. Among the median-based methods, TMLE performed better, but it took much more time than the other methods, including the proposed methods, and occasionally ($<1\%$) failed to converge.