Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
22,379 characters · 7 sections · 11 citation commands
Semiparametric estimation of heterogeneous treatment effects under the nonignorable assignment condition
In causal inference, treatment effects such as average treatment effect (ATE) or average treatment effect on the treated (ATT) have been of primary interest in the literature rubin_estimating_1974. These parameters are, as the names stand for, an averaged treatment effect over a population of interest. However, a treatment effect may differ among units depending on the covariates and outcomes. Such heterogeneity of treatment effects has been intensively studied in recent years. For example, wager_estimation_2018 proposed a method for finding subgroups in which the treatment effect is similar using the random forest. Understanding heterogeneity of treatment effects aids not only to more detailed analysis of a population of interest but also to a more efficient policy-making where an intervention is costly.
While most studies have concerned heterogeneity over the covariates which are fully observed, takahata_identification_2018 (henceforth TH) studied the heterogeneity over the untreated potential outcome, which is defined as
where $y_{1}$ and $y_{0}$ are the outcome when receiving the treatment and control condition respectively. Following TH, we call ((ref)) as the heterogeneous treatment effect (HTE), which is also of interest in this paper. To estimate the HTE, we have to deal with non-ignorable missingness because we need to estimate $p(y_{1}|y_{0})$ but $y_{1}$ and $y_{0}$ are never observed simultaneously. It is known that identification is not trivial in non-ignorable missing models (e.g., miao_identifiability_2016). TH provided the sufficient condition for the identification of the HTE using the information of the marginal distribution of $y_{0}$, $p(y_{0})$. Although the identification condition is rather general, estimation of the HTE is difficult in that it is necessary to solve some integral equation; the integral equation that we need to consider is a Fredholm integral equation of the first kind, which is known to be a ill-posed problem. In general, appropriate regularization methods are needed to obtain a stable solution to such equation. In TH, this problem was avoided using a parametric Bayesian modeling, but for wide applicability, a more flexible approaches is desired.
In this paper, we propose a semiparametric two-stage least square (2SLS) estimator for the HTE. Our approach relies on the quadratic programming method proposed by newey_instrumental_2003, in which they concerned the estimation of a nonparametric instrumental variable model. The function of interest is approximated by series of a finite number and then the integral equation reduces to a constrained least square problem with the regressors replaced by its expectation. To overcome instability of the solution due to the ill-posedness, certain bounds are placed on the coefficients of the series to make the inverse mapping of integral to be continuous. The numerical experiments show that the proposed method correctly estimate the HTE.
We follow the same setup as TH. The HTE (eq. ((ref))) is rewritten as
where $x\in \mathbb{R}^{d}$ is $d$-dimensional covariate. From this formula we observe that, for the identification of the HTE, it is sufficient to identify $p(x|y_{0})$ and $p(y_{1}|y_{0},x)$. Let $z\in\left\{0,1\right\}$ be the binary indicator, which is equal to 1 when assigned to the treatment condition. If $z=1$, then only $y_{1}$ is observed and $y_{0}$ is missing, and vice versa. TH showed that the following two assumptions play a primary role for the identification:
Assumption (A1) is called weak ignorability. Intuitively, this assumption implies that we can identify the HTE by observing the difference of two groups in which the assignment probability depends on $y_{0}$. Therefore, in a situation where strong ignorability, $p(z|y_{1},y_{0},x)=p(z|x)$, is satisfied, our approach is not applicable. Assumption (A2) is needed for identify $p(z|y_{0},x)$ hirano_combining_2001. In addition to these conditions, several constraints on the functional form and parameter space of $p(z|y_{0},x)$ is needed; for more detailed discussion on the identification condition, see TH. In what follows, we suppose that all the conditions mentioned in Theorem 2 in TH are satisfied. Note that, in this setup, the identification of ATE is trivial because
Consider the integral equation
where the second equality holds by weak ignorability. Under the identification condition, it can be proved that the solution to eq. ((ref)) for $\mathbb{E}_{}\mathopen{}\left[y_{1}|y_{0},x\mathopen{}\right]$ is unique, that is, $\mathbb{E}_{}\mathopen{}\left[y_{1}|y_{0},x\mathopen{}\right]$ is identified. Then our goal is to obtain an actual solution to eq. ((ref)). However, if we employ a nonparametric method for estimating $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0},x\mathopen{}\right]$, a solution suffers from the instability due to the discontinuity of the inverse mapping of integral. Then we need to take an appropriate regularization method to overcome the ill-posedness of eq. ((ref)). We address this problem in the next section.
In this section we propose a two-stage least square estimator (2SLSE) for $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0},x\mathopen{}\right]$ based on newey_instrumental_2003's method. The strategy is that (i) we approximate $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0},x\mathopen{}\right]$ by a finite number of orthogonal basis functions, (ii) take expectation of them with respect to $p(y_{0}|x,z=1)$, and (iii) do least square estimation under the constraint to make the inverse mapping of the integral to be continuous and eliminate the ill-posedness of eq. ((ref)). For simplicity we suppose $x\in\mathbb{R}$ in the rest of the paper.
We consider approximating $\phi(y_{0},x)\equiv \mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0},x\mathopen{}\right]$ with a finite number of orthogonal basis functions, $\{q_{j}\}_{j=0}^{J}$, as
We specify $p(z=1|y_{0},x)$ by the logistic regression such that it satisfies the identification condition from TH:
where the additivity holds between $y_{0}$ and $x$. Expand $p(y_{0}|x,z=1)$ as
where $c(x) = \exp(k_{0} + k_{x}(x))p(z=0)/p(x,z=1)$. Plugging eq. ((ref)) and eq. ((ref)) into eq. ((ref)) yields
where $s_{j_{1}}(x) = \int q_{j_{1}}(y_{0})\exp(k_{y_{0}}(y_{0}))p(y_{0},x|z=0)dy_{0}$. The estimation of the missing mechanism will be discussed later.
Here, we consider estimating $p(y_{0},x|z=0)$ by the kernel density estimator,
where $h_{y_{0}}$ and $h_{x}$ are the bandwidths and $N_{0}$ is the sample size of the control group. In this case, $s_{j_{1}}(x)$ can be estimated by
where $\hat{t}_{j_{1}}(y_{i0}) = \int q_{j_{1}}(y_{0})\exp(\hat{k}_{y_{0}}(y_{0}))K\left(\frac{y_{0}-y_{i0}}{\hat{h}_{y_{0}}}\right)dy_{0}$. Similarly, we obtain an estimator for $c(x)$, $\hat{c}(x)$, by the kernel density estimator,
By inserting eq. ((ref)) and eq. ((ref)) into eq. ((ref)), we obtain
where $N_{1}$ is the sample size of the treatment group. Therefore, the least square estimator for $\gamma_{j_{1} j_{2}}$ is obtained by the following quadratic problem:
where $\gamma = (\gamma_{11},\gamma_{12},\dots,\gamma_{j_{1}j_{2}},\dots,\gamma_{JJ})'$, $B_{\gamma}$ is a positive constant. $\Lambda_{J}$ is the Sobolev norm of the basis functions, $\left\{q_{j_{1}}(\cdot)q_{j_{2}}(\cdot)\right\}_{(j_{1},j_{2})}\, (j_{1}=0,\dots,J, \, j_{2}=0,\dots,J)$, which imposes compactness on both the true and estimated functions of $\phi(y_{0},x)$. This compactness eliminates the ill-posedness of the inverse problem (eq. ((ref))) because the inverse operator of the integral becomes continuous mapping newey_instrumental_2003.
Finally, by integrating out $x$ in $\hat{\phi}$, we obtain $\hat{E}[y_{1}|y_{0}]$:
Note that we can calculate $\hat{p}(x|y_{0})$ by plugging corresponding estimators into the following formula:
We estimate the missing mechanism (eq. ((ref))) referring to nevo_using_2003, who proposed a generalized method of moments (GMM) estimator for a nonignorable missing model. Suppose that the auxiliary moment condition, $\mathbb{E}_{}\mathopen{}\left[m(x,y_{0})\mathopen{}\right]=0$, is available. From assumption (A2), we can calculate any moments of $y_{0}$ up to infinite dimension, but for simplicity, the dimension of the moment condition is set to be equal to the sum of the dimension of the parameters of $k_{y_{0}}(\cdot)$ and $k_{x}(\cdot)$. For example, if $k_{x}(x)=\beta_{0}x$ and $k_{y_{0}}(y_{0})=\beta_{1}y_{0}+\beta_{2}y_{0}^{2}$, then we may set the moment function as
Note that
Therefore, the solution of the following system of equations,
is an unbiased estimator for $(k_{0},\beta_{0},\beta_{1},\beta_{2})$. $N$ is the total sample size. The second equation implies the normalization of the weights, which is due to $1/N\cdot \sum_{i:z_{i}=0}1/p(z_{i}=0|y_{i0},x_{i})\stackrel{p}{\longrightarrow}1$. For detailed identification conditions, see nevo_using_2003.
We conduct simulations to examine the performance of the estimator shown in the previous section. Each data set is generated from
where $\sigma_{0} = 1/5$, $\sigma_{1} = 1/2$, $\rho = 1/2$, $\mu_{0}(x) = -3x/5 - 1/10$, $\mu_{1}(x) = -(x-1)^{2}/10 + 1$ and the sample size is $N=3000$ (Figure (ref), (ref)). In this study, we use Legendre polynomials up to the third-order ($J=3$) as the basis functions in eq. ((ref)):
Then, we make an appropriate linear transformation so that each variable is included in $[-1,1]$ because Legendre polynomials are orthogonal on this interval (hereafter, the variables denote the transformed values. After estimation, the inverse transformation is made to calculate ATE). We specify the functions in the missing mechanism (eq. ((ref))) as $k_{x}(x) = \beta_{0} x$ and $k_{y_{0}}(y_{0}) = \beta_{1} y_{0} + \beta_{2}y_{0}^{2}$ and $z_{i}$ is drawn from $\mathrm{B}(1,p_{i})$ where $p_{i} = g(k_{0}+k_{y_{0}}(y_{i0}) + k_{x}(x_{i}))$. We set $k_{0} = -3/2$, $\beta_{0} = -2$, $\beta_{1} = -2$, $\beta_{2} = 1$ and the mean of the probability of being assigned to the treatment group becomes about $30\%$. For the units where $z_{i} = 1$, only $y_{i1}$ and $x_{i}$ are used for estimation, and conversely, $y_{i0}$ and $x_{i}$ are used where $z_{i} = 0$. Bandwidths for the kernel density estimators are chosen by Scott's normal reference rule of thumb. For estimation, we first estimate the missing mechanism, and then $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0},x\mathopen{}\right]$ given the former estimator.
Figure (ref) shows the result of the estimation of $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0}\mathopen{}\right]$. The horizontal axis is $y_{0}$ and the vertical axis is $y_{1}$. The dashed line shows the theoretical value of $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0}\mathopen{}\right]$. The solid line and the gray region show the mean of the estimator and the 90% confidence interval from 1000 replications respectively. As the figure shows, the performance of the estimator is substantially influenced by $B_{\gamma}$; if we set a small value as $B_{\gamma}$, the variance of the estimator also becomes small, whereas the confidence interval may not include the true curve and the expectation of the estimator may be apart from it. This problem is particularly serious on the edge, where only a small number of samples is observed. However, it is notable that we can estimate $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0}\mathopen{}\right]$ to some extent although none of the pair $(y_{i1},y_{i0})$ is observed and there are large overlaps in the distributions (see Figure (ref)).
Table (ref) shows the result of estimating ATE by integrating out $y_{0}$ in $\hat{E}[y_{1}|y_{0}]$. The theoretical value of the ATE in our setting is 0.900. As the estimation of $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0}\mathopen{}\right]$, the variance becomes larger as $B_{\gamma}$ becomes larger. On the other hand, the mean of the estimator for the ATE gets closer to the true value when $B_{\gamma}$ is small. Although we recognize its importance, how to determine $B_{\gamma}$ is beyond the scope of this paper.
We proposed a semiparametric two-stage least square estimator for the HTE and examine its properties through a simple simulation study, showing the availability of estimating $\mathbb{E}_{}\mathopen{}\left[y_{1} \middle| y_{0}\mathopen{}\right]$ even though none of the pair $(y_{1},y_{0})$ is observed. As mentioned in Section 4, the performance of the estimator shown in this paper is influenced substantially by the constraint parameter, which we have to tune. In addition, although we use Legendre polynomials up to the third-order in the simulation, the order actually needs to be determined reflecting the characteristics of the target population. Although there is literature on this issue (e.g., horowitz_adaptive_2014), a decisive method has not been developed. More importantly, our approach would not work for a multivariate case because it uses the kernel density estimator, that is, our method suffers from the curse of dimensionality. Similarly, an approximation using the orthogonal basis functions as eq. ((ref)) would be a problem because the number of parameters grows with $J^{d}$ (“approximation order" to the power “the number of dimensions of covariates") rate. Considering this issue, a nonparametric Bayesian approach may be helpful. We are planning to address this in future work.