A Simple and Efficient Estimation of the Average Treatment Effect in the Presence of Unmeasured Confounders

A common approach to account for individual heterogeneity in the treatment effect literature on observational data is to assume that there exist confounders, and conditional on these confounders, there is no systematic selection into the treatment (i.e., the so-called Unconfounded Treatment Assignment condition suggested in Rosenbaum and Rubin (1983, 1984)). Under this assumption, several procedures for estimating the average treament effect (hereafter ATE) have been proposed, including the weighting procedure (Rosenbaum (1987), Hirano, Imbens, and Ridder (2003), Tan (2010), Imai and Ratkovic (2014), Chan, Yam, and Zhang (2016), Yiu and Su (2018)); the matching procedure (Rosenbaum (2002), Rosenbaum et al. (2002), Dehejia and Wahba (1999)); and the regression procedure (Heckman, Ichimura, and Todd (1997), Heckman, Ichimura, and Todd (1998), Imbens, Newey, and Ridder (2006), Chen, Hong, and Tarozzi (2008)). For example survey, see Imbens and Wooldridge (2009) and Imbens and Rubin (2015). A critical requirement in this literature is that all confounders are observed and available to researchers. In applications, however, it is often the case that some confounders are either not observed or not availale. In this case, the average treatment effect is only partially identified even with the aid of some insrumental variables (see Imbens and Angrist (1994), Angrist, Imbens, and Rubin (1996), Abadie (2003), Abadie, Angrist, and Imbens (2002), Tan (2006), Cheng, Small, Tan, and Have (2009), Ogburn, Rotnitzky, and Robins (2015) for examples).

Recently Wang and Tchetgen Tchetgen (2017) suggested a noval identification condition of ATE when some confounders are not available. Under their condition, they showed that the semiparametric efficient influence function of ATE depends on five unknown functionals. They proposed to parameterize all five functionals, estimate those functionals with appropriate parametric approaches, plug the estimated functionals into the influence function, and then estimate the ATE from the estimated influence function. They established that their estimator is consistent if certain functionals are correctly parameterized and attains the semiparametric efficiency bound if all functionals are correctly specified. In applications, it is quite possible that some or all of the five functionals are misspecified and consequently their estimator could be inefficient or worse, inconsistent. This paper proposes an alternative, intuitive and easy to compute estimation that does not require parameterization of any of the five unknown functionals. We estbalish that under some sufficient conditions the proposed estimator is consistent, asymptotically normally distributed and attains the semiparametric efficiency bound. Moreover, the proposed procedure provides a natural and convenient estimate of the asymptotic variance.

The paper is organized as follows. Section 2 describes the basic framework. Section 3 describes the proposed estimation and derives the large sample properties of the proposed estimator. Section 4 presents a consistent variance estimator. Since the proposed procedure depends on smoothing parameters, Section 5 presents a data driven method for selecting the smoothing paprameters. Section 6 reports a small scale simulation study to evaluate the finite sample performance of the proposed estimator. Some concluding remarks are in Section 7. All technical proofs are relegated to the Appendix and the supplementary material.

Let $D\in\{0,1\}$ denote the binary treatment indicator, and let $Y(1)$ and $Y(0)$ denote the potential outcomes when an individual is assigned to the treatment and control group respectively. The parameter of interest is the population average treatment effect $\tau=\mathbb{E}[Y(1)-Y(0)]$. Estimation of $\tau$ is complicated by the presence of confounders and the fact that $Y(1)$ and $Y(0)$ cannot be observed simultaneously. To distinguish observed confounders from unobserved confounders, we shall use $X$ to denote the observed confounders and use $U$ to denote the unmeasured confounders. It is well established in the literature that, when all confounders are observed, the following Unconfounded Treatment Assignment condition is sufficient to identify $\tau$:

When $U$ is unmeasured, we have the classical omitted variable problem, causing the treament indicator $D$ to be endogenous. To tackle the endogeneity problem, instrumental variable is often the preferred choice. Let $Z\in\{0,1\}$ denote the variable satisfying the following classical instrumental variable conditions:

Wang and Tchetgen Tchetgen (2017) showed that Asssumptions 2.1- 2.4 alone do not identify $\tau$, but if in addition one of the following conditions holds:

1. 1.

   there is no additive $U$-$Z$ interaction in $\mathbb{E}[D|Z,X,U]$:

   $$\mathbb{E}[D|Z=1,X,U]-\mathbb{E}[D|Z=0,X,U]=\mathbb{E}[D|Z=1,X]-\mathbb{E}[D|Z=0,X]\ .$$

2. 2.

   there is no additive $U$-$d$ interaction in $\mathbb{E}[Y(d)|X,U]$:

   $$\mathbb{E}[Y(1)-Y(0)|X,U]=\mathbb{E}[Y(1)-Y(0)|X]\ ,$$

then ATE is identified and can be expressed as

$\displaystyle\tau=\mathbb{E}[\delta(X)]=\mathbb{E}\left[\frac{\delta^{Y}(X)}{\delta^{D}(X)}\right]\ ,$

(2.1)

where

	$\displaystyle\delta^{Y}(X)=\mathbb{E}[Y|Z=1,X]-\mathbb{E}[Y|Z=0,X]\ ,$
	$\displaystyle\delta^{D}(X)=\mathbb{E}[D|Z=1,X]-\mathbb{E}[D|Z=0,X]\ ,$
	$\displaystyle\delta(X)=\delta^{Y}(X)/\delta^{D}(X)\ .$

Furthermore, Wang and Tchetgen Tchetgen (2017) derived the efficient influence function for $\tau$:

$$\varphi_{eff}(D,Z,X,Y)=\frac{2Z-1}{f_{Z|X}(Z|X)}\frac{1}{\delta^{D}(X)}\bigg{\{}Y-D\delta(X)-\mathbb{E}[Y|Z=0,X]+\mathbb{E}[D|Z=0,X]\delta(X)\bigg{\}}+\delta(X)-\tau\ ,$$

where $f_{Z|X}(Z|X)$ is the conditional probability mass function of $Z$ given $X$. Clearly, the efficient influence function depends on five unknown functionals: $\delta(X)$, $\delta^{D}(X)$, $f_{Z|X}$, $p_{0}^{Y}(X)=\mathbb{E}[Y|Z=0,X]$ and $p_{0}^{D}(X)=\mathbb{E}[D|Z=0,X]$. They proposed to parameterize all five functionals, estimate the functionals with appropriate parametric approaches, and plug the estimated functionals into the efficient influence function to estimate $\tau$. They established that their estimator of $\tau$ is consistent and asymptotically normally distributed if

• •

  either $\delta(X)$, $\delta^{D}(X)$, $p_{0}^{Y}(X)=\mathbb{E}[Y|Z=0,X]$ and $p_{0}^{D}(X)=\mathbb{E}[D|Z=0,X]$ are correctly specified

• •

  or $\delta^{D}(X)$ and $f_{Z|X}$ are correctly specified

• •

  or $\delta(X)$ and $f_{Z|X}$ are correctly specified,

and their estimator attains the semiparametric efficiency bound only when all five functionals are correctly specified. The main goal of this paper is to present an alternative, intuitive and easy approach to compute estimator that does not require parameterization of any of the functionals and is always consistent and asymptotically normal and attains the semiparametric efficiency bound.

To motivate our estimation procedure, we rewrite the treatment effect coefficient. Applying the tower law of conditional expectation, we obtain:

	$\displaystyle\tau=$	$\displaystyle\mathbb{E}\left[\frac{\delta^{Y}(X)}{\delta^{D}(X)}\right]=\mathbb{E}\left[\frac{\mathbb{E}[Y|Z=1,X]}{\delta^{D}(X)}-\frac{\mathbb{E}[Y|Z=0,X]}{\delta^{D}(X)}\right]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left[\frac{Z}{f_{Z|X}(1|X)}\cdot\frac{\mathbb{E}[Y|Z=1,X]}{\delta^{D}(X)}-\frac{1-Z}{f_{Z|X}(0|X)}\cdot\frac{\mathbb{E}[Y|Z=0,X]}{\delta^{D}(X)}\right]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left[\frac{Z}{f_{Z|X}(1|X)}\cdot\frac{\mathbb{E}[Y|Z,X]}{\delta^{D}(X)}-\frac{1-Z}{f_{Z|X}(0|X)}\cdot\frac{\mathbb{E}[Y|Z,X]}{\delta^{D}(X)}\right]$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left[\left\{\frac{2Z-1}{f_{Z|X}(Z|X)}\right\}\frac{Y}{\delta^{D}(X)}\right]\ .$		(3.1)

The above expression suggests a natural and intuitive plugin estimation, with $f_{Z|X}(Z|X)$ and $\delta^{D}(X)$ replaced by some consistent estimates. There are many approaches to estimate these functionals including parametric and nonparametric approaches, but as noted by Hirano, Imbens, and Ridder (2003), not all estimates can lead to efficient estimation of $\tau$. In this paper, we present an intuitive and easy way to compute estimates of functionals that ensure efficiency of the plugin estimation of $\tau$. To illustrate our procedure, we notice that the following conditions hold for any integrable functions $u_{1}(X)$ and $u_{2}(X)$:

	$\displaystyle\mathbb{E}\left[\frac{Z}{f_{Z|X}(1|X)}u_{1}(X)\right]=\mathbb{E}[u_{1}(X)]\ =\mathbb{E}\left[\frac{1-Z}{f_{Z|X}(0|X)}u_{1}(X)\right],$		(3.2)
	$\displaystyle\mathbb{E}\left[D\left\{\frac{2Z-1}{f_{Z|X}(Z|X)}\right\}u_{2}(X)\right]=\mathbb{E}\left[\delta^{D}(X)u_{2}(X)\right]\ ,$		(3.3)

and (3.2) and (3.3) uniquely determine $f_{Z|X}(Z|X)$ and $\delta^{D}(X)$. These conditions impose restrictions on the unknown functionals and they must be taken into account when estimating those functionals. One difficulty with these conditions is that they must be imposed in an infinite dimmensional functional space. To overcome this difficulty, we propose to impose the conditions on a smaller sieve space. Specifically, let $u_{K}(X)=(u_{K,1}(X),\ldots,u_{K,K}(X))^{\top}$ denote a known basis functions that can approximate any suitable function $u(X)$ arbitrarily well (see Chen (2007) or Appendix A.1 for further dicussion). Conditions (3.2) and (3.3) imply for any integers $K_{1}$ and $K_{2}$:

$$\mathbb{E}\left[\frac{Z}{f_{Z|X}(1|X)}u_{K_{1}}(X)\right]=\mathbb{E}[u_{K_{1}}(X)]=\mathbb{E}\left[\frac{1-Z}{f_{Z|X}(0|X)}u_{K_{1}}(X)\right]\$$

(3.4)

and

$$\mathbb{E}\left[D\left\{\frac{Z}{f_{Z|X}(1|X)}-\frac{1-Z}{f_{Z|X}(0|X)}\right\}u_{K_{2}}(X)\right]=\mathbb{E}[\delta^{D}(X)u_{K_{2}}(X)].\$$

(3.5)

We shall construct estimates of the functionals by imposing the above conditions. To ensure consistency, we shall allow $K_{1}$ and $K_{2}$ to increase with sample size at appropriate rates.

To conduct the statistical inference on $\tau$, we need a consistent estimator of the asymptotic variance of $\widehat{\tau}$. Note that the asymptotic varance of $\widehat{\tau}$,

$$\mathbb{E}\left[\left(\frac{2Z-1}{f_{Z|X}(Z|X)}\frac{1}{\delta^{D}(X)}\bigg{\{}Y-D\delta(X)-\mathbb{E}[Y|Z=0,X]+\mathbb{E}[D|Z=0,X]\delta(X)\bigg{\}}+\delta(X)-\tau\right)^{2}\right]\ ,$$

depends on five unknown functionals. Direct estimation of the variance requires replacing the five unknown functionals with consistent estimates. In this section, we present an alternative estimation that does not require estimation of those functionals.

To illustrate the idea, we denote:

	$\displaystyle g_{1}(Z,{X};\lambda)\triangleq Z\rho^{\prime}\left(\lambda^{\top}u_{K_{1}}({X})\right)u_{K_{1}}({X})-u_{K_{1}}({X})\ ,$
	$\displaystyle g_{2}(Z,{X};\beta)\triangleq(1-Z)\rho^{\prime}\left(\beta^{\top}u_{K_{1}}({X})\right)u_{K_{1}}({X})-u_{K_{1}}({X})\ ,$
	$\displaystyle g_{3}(Z,D,{X};\lambda,\beta,\gamma)\triangleq D\left\{Z\cdot\rho^{\prime}\left(\lambda^{\top}u_{K_{1}}(X)\right)-(1-Z)\cdot\rho^{\prime}\left(\beta^{\top}u_{K_{1}}(X)\right)\right\}u_{K_{2}}({X})-f^{\prime}\left(\gamma^{\top}u_{K_{2}}({X})\right)u_{K_{2}}({X})\ ,$
	$\displaystyle g_{4}(Z,D,X,Y;\lambda,\beta,\gamma,\tau)\triangleq\left\{Z\cdot\rho^{\prime}\left(\lambda^{\top}u_{K_{1}}(X)\right)-(1-Z)\cdot\rho^{\prime}\left(\beta^{\top}u_{K_{1}}(X)\right)\right\}Y/f^{\prime}\left(\gamma^{\top}u_{K_{2}}(X)\right)-\tau\ ,$

and

$$g(Z,D,X,Y;\theta)\triangleq\left(\begin{array}[]{c}g_{1}(Z,{X};\lambda)\\ g_{2}(Z,{X};\beta)\\ g_{3}(Z,D,{X};\lambda,\beta,\gamma)\\ g_{4}(Z,D,X,Y;\lambda,\beta,\gamma,\tau)\end{array}\right)$$

with $\theta\triangleq(\lambda,\beta,\gamma,\tau)^{\top}$. Let $\hat{\theta}\triangleq(\hat{\lambda}_{K_{1}},\hat{\beta}_{K_{1}},\hat{\gamma}_{K_{2}},\hat{\tau})^{\top}$ and ${\theta}^{\ast}\triangleq({\lambda}_{K_{1}}^{\ast},{\beta}_{K_{1}}^{\ast},{\gamma}_{K_{2}}^{\ast},{\tau})^{\top}$. Then $\hat{\theta}$ is the moment estimator solving the following moment condition:

$$\frac{1}{N}\sum_{i=1}^{N}g(Z_{i},D_{i},X_{i},Y_{i};\hat{\theta})=0.$$

(4.1)

Applying Mean Value Theorem, we obtain

$$0=\frac{1}{N}\sum_{i=1}^{N}g(Z_{i},D_{i},X_{i},Y_{i};\theta^{\ast})+\frac{1}{N}\sum_{i=1}^{N}\frac{\partial g(Z_{i},D_{i},X_{i},Y_{i};\tilde{\theta})}{\partial\theta}(\hat{\theta}-\theta^{\ast})$$

(4.2)

where $\tilde{\theta}=(\tilde{\lambda}_{K_{1}},\tilde{\beta}_{K_{1}},\tilde{\gamma}_{K_{2}},\tilde{\tau})^{\top}$ lies on the line joining $\hat{\theta}$ and $\theta^{\ast}$. We show in the supplemental material that

$$\frac{1}{N}\sum_{i=1}^{N}\frac{\partial g(Z_{i},D_{i},X_{i},Y_{i};\tilde{\theta})}{\partial\theta}=\mathbb{E}\left[\frac{\partial g(Z,D,X,Y;\theta^{\ast})}{\partial\theta}\right]+o_{p}(1)$$

(4.3)

Note that

$\displaystyle\hat{\tau}-\tau=\mathbf{e}_{2K_{1}+K_{2}+1}^{\top}(\hat{\theta}-\theta^{\ast})\ ,$

(4.4)

where $\mathbf{e}_{2K_{1}+K_{2}+1}$ is a $(2K_{1}+K_{2}+1)$-dimensional column vector whose last element is $1$ and other components are all of $0$’s.

Combining (4.2), (4.3) and (4.4), we obtain

$$\sqrt{N}(\hat{\tau}-\tau)=-\mathbf{e}_{2K_{1}+K_{2}+1}^{\top}\left\{\mathbb{E}\left[\frac{\partial g(Z,D,X,Y;\theta^{\ast})}{\partial\theta}\right]+o_{p}(1)\right\}^{-1}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}g(Z_{i},D_{i},X_{i},Y_{i};\theta^{\ast})\ ,$$

which in turn implies

$$V_{eff}=\lim_{N\rightarrow\infty}Var(\sqrt{N}(\hat{\tau}-\tau))=\lim_{N\rightarrow\infty}\mathbf{e}_{2K_{1}+K_{2}+1}^{\top}\left\{L\cdot\Omega\cdot(L^{-1})^{\top}\right\}\mathbf{e}_{2K_{1}+K_{2}+1}\ .$$

where

	$\displaystyle L=\mathbb{E}\left[\frac{\partial g(Z,D,X,Y;\theta^{\ast})}{\partial\theta}\right]\ ,$
	$\displaystyle\Omega=\mathbb{E}\left[g(Z,D,X,Y;\theta^{\ast})g(Z,D,X,Y;\theta^{\ast})^{\top}\right]\ .$

Therefore, we can define the sandwich estimator for the efficient variance $V_{eff}$ by

$$\hat{V}=\mathbf{e}_{2K_{1}+K_{2}+1}^{\top}\left\{\hat{L}^{-1}\cdot\hat{\Omega}\cdot(\hat{L}^{-1})^{\top}\right\}\mathbf{e}_{2K_{1}+K_{2}+1}\ ,$$

where

	$\displaystyle\hat{L}=\frac{1}{N}\sum_{i=1}^{N}\frac{\partial g(Z_{i},D_{i},X_{i},Y_{i};\hat{\theta})}{\partial\theta};$
	$\displaystyle\hat{\Omega}=\frac{1}{N}\sum_{i=1}^{N}g(Z_{i},D_{i},X_{i},Y_{i};\hat{\theta})g(Z_{i},D_{i},X_{i},Y_{i};\hat{\theta})^{\top}.$

The large sample properties of the proposed estimator permit a wide range of values of $K_{1}$ and $K_{2}$. This presents a dilemma for applied researchers who have only one finite sample and would like to have some guidance on the selection of smoothing parameters. In this section, we present a data-driven approach to select $K_{1}$ and $K_{2}$. Notice that $f_{Z|X}(1|X)^{-1}$, $f_{Z|X}(0|X)^{-1}$ and $\delta^{D}(X)$ satisfy the following regression equations:

	$\displaystyle\mathbb{E}\left[Zf_{Z|X}(1|X)^{-1}\bigg{|}X\right]=1\ ,$
	$\displaystyle\mathbb{E}\left[(1-Z)f_{Z|X}(0|X)^{-1}\bigg{|}X\right]=1\ ,$
	$\displaystyle\mathbb{E}\left[D\left\{Zf_{Z|X}(1|X)^{-1}-(1-Z)f_{Z|X}(0|X)^{-1}\right\}\bigg{|}X\right]=\delta^{D}(X)\ .$

Since $N\hat{p}(X)$, $N\hat{q}(X)$ and $\hat{\delta}^{D}(X)$ are consistent estimators of $f_{Z|X}(1|X)^{-1}$, $f_{Z|X}(0|X)^{-1}$ and $\delta^{D}(X)$ respectively, the mean-squared-error (MSE) of the nuisance parameters $(\hat{\lambda}_{K_{1}},\hat{\beta}_{K_{1}})$ and $\hat{\gamma}_{K_{2}}$ are defined by

	$\displaystyle MSE_{1}(K_{1})=$	$\displaystyle\sum_{i=1}^{N}\left\{Z_{i}N\hat{p}(X_{i})-1\right\}^{2}+\sum_{i=1}^{N}\left\{(1-Z_{i})N\hat{q}(X_{i})-1\right\}^{2}\ ,$
	$\displaystyle MSE_{2}(K_{1},K_{2})=$	$\displaystyle\sum_{i=1}^{N}\left\{D_{i}\left\{Z_{i}N\hat{p}(X_{i})-(1-Z_{i})N\hat{q}(X_{i})\right\}-\hat{\delta}^{D}(X_{i})\right\}^{2}\ .$

The smoothing parameters $K_{1}$ and $K_{2}$ shall be chosen to minimize $MSE_{1}$ and $MSE_{2}$. Specifically, denote the upper bounds of $K_{1}$ and $K_{2}$ by $\bar{K}_{1}$ and $\bar{K}_{2}$ (e.g. $\bar{K}_{1}=\bar{K}_{2}=5$ in our simulation studies). The data-driven $K_{1}$ and $K_{2}$ are given by

	$\displaystyle\hat{K}_{1}=\arg\min_{K_{1}\in\{1,...,\bar{K}_{1}\}}MSE_{1}(K_{1})\ ,$
	$\displaystyle\hat{K}_{2}=\arg\min_{K_{2}\in\{1,...,\bar{K}_{2}\}}MSE_{2}(\hat{K}_{1},K_{2})\ .$

In this section, we conduct a small scale simulation study to evaluate the finite sample performance of the proposed estimator. To evaluate the performance of our estimator against the existing alternatives, particularly the estimators proposed by Wang and Tchetgen Tchetgen (2017), we adopt the exact same design (i.e., the same data generating processes (DGP)). In each Monte Carlo run, we generate sample of data from DGP for two sizes: $N=500$ and $N=1000$ respectively, and from each sample we compute our estimator and other existing estimators. We then repeat the Monte Carlo runs for $500$ times.

The observed baseline covariates are $X=(1,X_{2})$, where $X$ include an intercept term and a continuous random variable $X_{2}$ uniformly distributed on the interval $(-1,-0.5)\cup(0.5,1)$. The unmeasured confounder $U$ is a Bernoulli random variable with mean 0.5. The instrumental variable $Z$, treatment variable $D$ and outcomes variable $Y\in\{0,1\}$ are generated according to the simulation design of Wang and Tchetgen Tchetgen (2017). The true value of the average treatment effect is $\tau=0.087$.

We compute the proposed estimator (cbe), the naive estimator, the multiply robust estimator (mr) and the bounded multiply robust estimator (b-mr) proposed by Wang and Tchetgen Tchetgen (2017). Details of calculations are given below.

1. 1.

   the proposed estimator (cbe) is computed with $\rho(v)=\log(1+v)$;

2. 2.

   the naive estimator is computed by the difference of group means between treatment and control groups;

3. 3.

   the multiply robust estimator (mr) and the bounded multiply robust estimator (b-mr) are computed by the procedures proposed by Wang and Tchetgen Tchetgen (2017).