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.
109,493 characters · 34 sections · 82 citation commands
Assessing Sensitivity to Unconfoundedness: Estimation and Inference
JEL classification: C14; C18; C21; C51
Keywords: Treatment Effects, Conditional Independence, Unconfoundedness, Selection on Observables, Sensitivity Analysis, Nonparametric Identification, Partial Identification
\onehalfspacing
A core goal of causal inference is to identify and estimate effects of a treatment variable on an outcome variable. A common assumption used to identify such effects is unconfoundedness, which says that potential outcomes are independent of treatment conditional on covariates. This assumption is also known as conditional independence, selection on observables, ignorability, or exogenous selection; see Imbens2004 for a survey. This assumption is not refutable, meaning that the data alone cannot tell us whether it is true. Nonetheless, empirical researchers may wonder: How important is this assumption in their analyses? Put differently: How sensitive are their results to failures of the unconfoundedness assumption?
A large literature on sensitivity analysis has developed to answer this question. Moreover, researchers widely acknowledge that answering this question is an important step in empirical research. For example, in their figure 1, CaliendoKopeinig2008 describe the workflow of a standard analysis using selection on observables. Their fifth and final step in this workflow is to perform sensitivity analysis to the unconfoundedness assumption. ImbensWooldridge2009, ImbensRubin2015, and AtheyImbens2017 all also recommend that researchers conduct sensitivity analyses to assess the importance of non-refutable identifying assumptions. In particular, AtheyImbens2017 describe these methods as “a systematic way of doing the sensitivity analyses that are routinely done in empirical work, but often in an unsystematic way.”
Most of the existing approaches to assessing unconfoundedness rely on strong auxiliary assumptions, however. For example, they often assume treatment effects are homogeneous and that all unobserved confounding arises due to a single unobserved variable whose distribution is parametrically specified, like a binary or normal distribution. They also often assume a parametric functional form for potential outcomes, like a logit model for binary potential outcomes or a linear model for continuous potential outcomes. These assumptions---which are not needed for identification of the baseline model when unconfoundedness holds---raise a new question: Are the findings of these sensitivity analyses themselves sensitive to these extra auxiliary assumptions?
In this paper, we provide a set of tools for assessing the sensitivity of the unconfoundedness assumption which do not rely on strong auxiliary assumptions that are not used for the baseline analysis. We do this by studying nonparametric relaxations of the unconfoundedness assumption. Specifically, we apply the identification results of MastenPoirier2018, who consider a class of assumptions called conditional $c$-dependence. This class measures relaxations of conditional independence by a single scalar parameter $c \in [0,1]$. This parameter $c$ is the largest difference between the propensity score and the probability of treatment conditional on covariates and an unobserved potential outcome. Hence it has a straightforward interpretation as a deviation from conditional independence, as measured in probability units. For any positive $c$, conditional independence only partially holds, and so we cannot learn the exact value of our treatment effect parameters, like the average treatment effect (ATE) or the average effect of treatment on the treated (ATT). Instead, we only get bounds. MastenPoirier2018 derive closed-form expressions for these bounds as a function of $c$. Setting $c= 0$ yields the baseline model where unconfoundedness holds. Setting $c= 1$ yields the other extreme where no assumptions on selection are made, and hence gives the no assumption bounds as in Manski1990. The bounds are monotonic in $c$, so that small values of $c$ give narrow bounds while larger values of $c$ give wider bounds. Just how wide these bounds are---and hence how sensitive one's results are---depends on the data.
While MastenPoirier2018 studied identification of treatment effects under nonparametric relaxations of unconfoundedness, they did not study estimation or inference. We do that in this paper. First we propose sample analog estimators of the bounds on the conditional quantile treatment effect (CQTE), the conditional average treatment effect (CATE), the ATE, and the ATT. We do this using flexible parametric first step estimators of the propensity score and the conditional quantile function of the observed outcomes given treatment and covariates. Although such parametric restrictions are not required for our identification theory, the analysis of inference is complicated and non-standard even with these parametric first step estimators. Doing inference based on fully nonparametric first step estimators will likely require deriving and applying more general asymptotic theory for non-Hadmard differentiable functionals than currently exists. Hence we leave that to future work. Moreover, note that our approach of using nonparametric identification results paired with flexible parametric estimators is analogous to what is commonly done in the baseline model which imposes unconfoundedness: Identification is shown nonparametrically but many commonly used estimators are based on flexible parametric first step estimators. For example, see chapter 13 in ImbensRubin2015.
We derive the asymptotic distributions of our bound estimators using the delta method for Hadamard directionally differentiable functionals from FangSantos2014. We then show consistency of a non-standard bootstrap based on estimating the analytical Hadamard directional derivatives of our bound functionals. This step again involves using the recent results of FangSantos2014. We show how to construct confidence bands for the bound functions which are uniform over all values of $c \in [0,1]$. We also provide a sufficient condition on the propensity score and the distribution of the covariates under which we can do inference using the standard nonparametric bootstrap. Finally, we show how to implement our analysis in an empirical illustration to the National Supported Work Demonstration program (MDRC MDRC1983). Using the techniques developed in this paper, and implemented in an accompanying Stata module, researchers can quantify the robustness of treatment effects estimated using the unconfoundedness assumption.
The rest of this paper is organized as follows. In the rest of this section we briefly discuss the related literature. In section (ref) we summarize the identification results from MastenPoirier2018. We also discuss how to use and interpret these results in practice. Section (ref) describes the definition of our bound estimators. Section (ref) provides the corresponding asymptotic estimation and inference theory for these estimators. Section (ref) describes how to use these inference results to conduct bootstrap based inference. In section (ref) we give sufficient conditions under which standard bootstrap approaches are valid. Section (ref) shows how to use our methods in an empirical illustration. Appendix (ref) contains theoretical results and proofs for our first step estimators. Appendices (ref), (ref), and (ref) have proofs for our main results. Appendix (ref) gives the full expressions for various analytical Hadamard directional derivatives used in our analysis. Appendix (ref) provides several additional results.
We conclude this section with a brief literature review. As mentioned earlier, there is a large existing literature that studies how to relax unconfoundedness. This includes RosenbaumRubin1983sensitivity, Mauro1990, Rosenbaum1995, Rosenbaum2002, RobinsRotnitzkyScharfstein2000, Imbens2003, AltonjiElderTaber2005, AltonjiElderTaber2008, IchinoMealliNannicini2008, HosmanHansenHolland2010, Krauth2016, KallusMaoZhou2019, Oster2019, and CinelliHazlett2020, among others. Here we discuss the most closely related work and several recent papers. For further details about the related literature, see section 1 in MastenPoirier2018 for identification and Appendix D in MastenPoirier2020 for estimation and inference.
A key feature of our results is that they are based on the fully nonparametric analysis of MastenPoirier2018. There are only a few other alternative nonparametric analyses available in the literature. The first is IchinoMealliNannicini2008, who require that all variables are discretely distributed. In contrast, we allow for continuous outcomes, covariates, and unobservables. Their approach requires picking a vector of sensitivity parameters that determines the joint distribution of the discrete observable and unobservable variables. In contrast, our approach uses a scalar sensitivity parameter. Finally, unlike us, they do not provide any formal results for doing estimation or inference. The second is Rosenbaum1995,Rosenbaum2002, who proposed a sensitivity analysis for unconfoundedness within the context of doing randomization inference based on the sharp null hypothesis of no unit level treatment effects for all units in the data set. Like our approach, he only uses a scalar sensitivity parameter and also does not rely on a parametric model for outcomes or treatment assignment probabilities. His approach, however, is based on finite sample randomization inference (for more discussion, see chapter 5 of ImbensRubin2015). This approach to inference is conceptually distinct from the approach we use based on repeated sampling from a large population. For this reason, we view these different approaches to inference in sensitivity analyses as complementary. Finally, KallusMaoZhou2019 study bounds on CATE under the same nonparametric relaxations defined by Rosenbaum1995,Rosenbaum2002. Unlike him, however, they take a large population view. They propose sample analog kernel estimators based on an implicit characterization of the identified set using extrema. They show consistency of these estimators, but they do not provide any inference results. As we discuss later, this is a key distinction because inference in this setting is non-standard.
A few recent papers provide methods for assesssing unconfoundedness in parametric linear models. This includes Oster2019 and CinelliHazlett2020. These results rely on the assumption that outcomes are linear functions of treatment and covariates, among other parametric assumptions. In contrast, we build on the selection on observables literature that has emphasized nonparametric identification. That literature emphasizes that identification by functional form is often implausible. Sensitivity analyses that rely on functional form assumptions are subject to the same criticism: Findings that one's results are robust to violations of unconfoundedness can be driven primarily from the parametric functional form restrictions. To address this, our estimation and inference results are based on nonparametric sensitivity analyses that do not require parametric assumptions.
Finally, we discuss the relationship with our own previous work. As noted earlier, our paper provides estimation and inference results for population bounds derived in MastenPoirier2018. That paper did not provide any estimation or inference theory. MastenPoirier2020 builds on those results in several ways: First, they extend the identification analysis to identification of distributional treatment effect parameters, with a focus on assessing the importance of the rank invariance assumption. Second, they provide some asymptotic distributional results for sample analog estimators of the average treatment effect (ATE), the conditional average treatment effect (CATE), and the conditional quantile treatment effect (CQTE), among other results. Those results are limited in a variety of ways, which we discuss next.
Specifically, our paper differs from the results in MastenPoirier2020 in several important ways: (1) Our paper allows for both discrete and continuous covariates, whereas that paper focused on the case where all covariates are discrete. In particular, to allow for continuous covariates we develop a different estimator of the bound functions. This is important since many empirical applications, like ours in section (ref), use continuous covariates. (2) Our results allow for all possible values of $c \in [0,1]$, whereas that paper restricted attention to small values of the sensitivity parameter $c$ (see their assumption A2.1). This is also important for practice and requires a substantial amount of new theoretical work. (3) Our results use the FangSantos2014 bootstrap based on estimators of analytical Hadamard directional derivatives to do inference. That paper instead used the numerical delta method bootstrap of HongLi2015. Our approach allows us to avoid choosing the step size tuning parameter required for the numerical delta method bootstrap, although our estimators of the analytical Hadamard directional derivatives also have tuning parameters. (4) Unlike that paper, we also discuss inference on the average effect of treatment on the treated (ATT). (5) In this paper we provide a new companion Stata module implementing our results.
In this section we describe the model and review standard results on point identification of treatment effects under unconfoundedness. We then describe how we relax unconfoundedness. Finally, we review the bounds on treatment effects derived by MastenPoirier2018 when unconfoundedness is relaxed.
We use the standard potential outcomes model. Let $X \in \{0, 1 \}$ be an observed binary treatment. Let $Y_1$ and $Y_0$ denote the unobserved potential outcomes. The observed outcome is
Let $W \in \ensuremath{\mathbb{R}}^{d_W}$ denote a vector of observed covariates, which may be discrete, continuous, or mixed. Let $\mathcal{W} = \operatorname*{supp}(W)$ denote the support of $W$. Let \[ p_{x \mid w} = \ensuremath{\mathbb{P}}(X=x \mid W=w) \] denote the observed generalized propensity score.
It is well known that the conditional distributions of potential outcomes $Y_1 \mid W$ and $Y_0 \mid W$ are point identified under the following two assumptions:
Consequently, any functional of the distributions of $Y_1 \mid W$ and $Y_0 \mid W$ is also point identified. We focus on two leading examples: The average treatment effect, $\text{ATE} = \ensuremath{\mathbb{E}}(Y_1 - Y_0)$ and the average treatment effect for the treated, $\text{ATT} = \ensuremath{\mathbb{E}}(Y_1 - Y_0 \mid X=1)$. We also consider the conditional quantile treatment effects $\text{CQTE}(\tau \mid w) = Q_{Y_1 \mid W}(\tau \mid w) - Q_{Y_0 \mid W}(\tau \mid w)$ and the conditional average treatment effect $\text{CATE}(w) = \ensuremath{\mathbb{E}}(Y_1 - Y_0 \mid W=w)$.
As discussed in section (ref), the overlap assumption is refutable and hence can be directly verified from the data. The unconfoundedness assumption, however, is not refutable. Consequently, like much of the literature reviewed in section (ref), we perform a sensitivity analysis. This entails replacing unconfoundedness with a weaker assumption and investigating how this changes the conclusions we can draw about our parameter of interest. Specifically, we define the following class of assumptions, which we call conditional $c$-dependence (MastenPoirier2018):
When $c = 0$, conditional $c$-dependence is equivalent to $X \mathbin{ \mathpalette{\@indep}{} } Y_x \mid W$. For $c > 0$, however, we allow for violations of unconfoundedness by allowing the unobserved conditional probability \[ \ensuremath{\mathbb{P}}(X=1 \mid Y_x=y_x, W=w) \] to differ from the observed propensity score \[ \ensuremath{\mathbb{P}}(X=1 \mid W=w) \] by at most $c$. Thus we actually allow for some selection on unobservables, since treatment assignment may depend on $Y_x$, but in a constrained manner. For sufficiently large $c$, however, conditional $c$-dependence imposes no constraints on the relationship between $Y_x$ and $X$. This happens when $c \geq \overline{C}$ where $\overline{C} =\sup_{w \in \mathcal{W}} \max \{ p_{1 \mid w}, p_{0 \mid w} \}$. When $c \in (0,\overline{C})$, conditional $c$-dependence imposes some constraints on treatment assignment, but it does not require conditional independence to hold exactly. For this reason, we call it a conditional partial independence assumption. Thus our sensitivity analysis replaces unconfoundedness with
By relaxing conditional independence our main parameters of interest---ATE and ATT---are no longer point identified. Instead they are partially identified: We can bound them from above and from below. As $c$ gets close to zero, however, these bounds collapse to a point. Hence for small $c$ these bounds can be quite narrow. The goal of a sensitivity analysis is to understand how the shape and width of these bounds changes as $c$ varies from 0 to 1.
These bounds were derived in MastenPoirier2018, which we summarize here. Although that paper studied both continuous and binary outcomes, here we only summarize the results for continuous $Y_x$. All of our parameters of interest can be written in terms of bounds on the quantile regressions $Q_{Y_x \mid W}(\tau \mid w)$. Under the conditional partial independence assumption stated above and some regularity conditions, MastenPoirier2018 showed that $[\underline{Q}^c_{Y_x \mid W}(\tau \mid w), \overline{Q}^c_{Y_x \mid W}(\tau \mid w)]$ are are sharp bounds on this quantile regression, uniformly in $\tau$, $x$, and $w$, where
and
Taking differences of these bounds for $x=1$ and $x=0$ yields sharp bounds on the conditional quantile treatment effect $\text{CQTE}(\tau \mid w)$, uniformly in $\tau$ and $w$:
Integrating these bounds over $\tau$ yields sharp bounds on $\text{CATE}(w)$, uniformly in $w$: \[ \left[\underline{\text{CATE}}^c(w),\overline{\text{CATE}}^c(w) \right] \equiv \left[\int_0^1 \underline{\text{CQTE}}^c(\tau \mid w) \; d\tau, \int_0^1 \overline{\text{CQTE}}^c(\tau \mid w) \; d\tau \right]. \] Further integrating over the marginal distribution of $W$ yields sharp bounds on ATE: \[ \left[\underline{\text{ATE}}^c,\overline{\text{ATE}}^c \right] \equiv \Big[ \ensuremath{\mathbb{E}} \big( \underline{\text{CATE}}^c(W) \big), \, \ensuremath{\mathbb{E}} \big( \overline{\text{CATE}}^c(W) \big) \Big] \] To obtain bounds on ATT, let \[ \underline{E}_x^c(w) = \int_0^1 \underline{Q}_{Y_x}^c(\tau \mid w) \; d\tau \qquad \text{and} \qquad \overline{E}_x^c(w) = \int_0^1 \overline{Q}_{Y_x}^c(\tau \mid w) \; d\tau \] denote bounds on $\ensuremath{\mathbb{E}}(Y_x \mid W=w)$. Averaging these over the marginal distribution of $W$ yields bounds on $\ensuremath{\mathbb{E}}(Y_x)$, denoted by \[ \underline{E}_x^c = \ensuremath{\mathbb{E}} \big( \underline{E}_x^c(W) \big) \qquad \text{and} \qquad \overline{E}_x^c = \ensuremath{\mathbb{E}} \big( \overline{E}_x^c(W) \big). \] This yields the following bounds on ATT:
where $p_x = \ensuremath{\mathbb{P}}(X=x)$ for $x \in \{0,1\}$. Finally, note that all of these bounds are sharp.
So far we've discussed sharp bounds on various parameters of interest as a function of the sensitivity parameter $c$. In addition to the bounds themselves, it is common to analyze breakdown points for various conclusions of interest. For example, suppose that under the baseline model ($c=0$) we find that $\text{ATE} > 0$. We then ask: How much can we relax unconfoundedness while still being able to conclude that the ATE is nonnegative? To answer this question, define the breakdown point for the conclusion that the ATE is nonnegative as
This number is a quantitative measure of the robustness of the conclusion that ATE is positive to relaxations of the key identifying assumption of unconfoundedness. Breakdown points can be defined for other parameters and conclusions as well. See MastenPoirier2020 for more discussion and additional references.
We conclude this section by giving some suggestions for how to interpret conditional $c$-dependence in practice. In particular, what values of $c$ are large? What values are small? Here we summarize and extend the discussion on page 321 of MastenPoirier2018. We illustrate these interpretations in our empirical analysis in section (ref).
Let $W_k$ denote a component of $W$. Denote the propensity score by \[ p_{1 \mid W}(w_{-k},w_k) = \ensuremath{\mathbb{P}}(X=1 \mid W=(w_{-k},w_k) ). \] Let \[ p_{1 \mid W_{-k}}(w_{-k}) = \ensuremath{\mathbb{P}}(X=1 \mid W_{-k}=w_{-k}) \] denote the leave-out-variable-$k$ propensity score. This is just the proportion of the population who are treated, conditional on only $W_{-k}$. Consider the random variable \[ \Delta_k = | p_{1 \mid W}(W_{-k}, W_k) - p_{1 \mid W_{-k}}(W_{-k}) |. \] This difference is a measure of the impact on the observed propensity score of adding $W_k$, given that we already included $W_{-k}$. Conditional $c$-dependence is defined by a similar difference, except there we add the unobservable $Y_x$ given that we already included $W$. Hence we suggest using the distribution of $\Delta_k$ to calibrate values of $c$. For example, you could examine the 50th, 75th, and 90th quantiles of $\Delta_k$, along with the upper bound on the support, $\bar{c}_k = \max \operatorname*{supp}(\Delta_k)$. You may also find it useful to plot an estimate of the density of $\Delta_k$. All of these reference values can be compared to the breakdown point $c_\textsc{bp}$ for a specific conclusion of interest. Specifically, if $c_\textsc{bp}$ is larger than the chosen reference value, then the conclusion of interest could be considered robust. In contrast, if $c_\textsc{bp}$ is smaller than the chosen reference value, then the conclusion of interest could be considered sensitive. You may also want to see where $c_\textsc{bp}$ lies relative to the distribution of $\Delta_k$. This can be done by computing $F_{\Delta_k}(c_\textsc{bp})$.
While you could do this for all covariates $k$, it may be helpful to restrict attention to covariates $k$ that have a sufficiently large impact on the baseline point estimates. For example, suppose we are interested in the ATE. Let $\text{ATE}_{-k}$ denote the ATE estimand obtained in the baseline selection on observables model using only the covariates $W_{-k}$. Let $\text{ATE}$ denote the ATE estimand obtained in the baseline model using all the covariates. Then \[ \left| \frac{\text{ATE} - \text{ATE}_{-k} }{\text{ATE}} \right| \] denotes the effect of omitting covariate $k$ on the ATE point estimand, as a percentage of the baseline estimand that uses all covariates in $W$. You may want to restrict attention to covariates $k$ for which this ratio is relatively large. We illustrate this approach in our empirical analysis in section (ref).
In the previous section we assumed the entire population distribution of $(Y,X,W)$ was known. In practice we only have a finite sample $\{ (Y_i,X_i,W_i) \}_{i=1}^n$ from this distribution. In this section we explain how to use this finite sample data to estimate the population bounds of section (ref). We give the corresponding asymptotic theory in section (ref) where we obtain the joint limiting distribution of treatment effects bounds. We describe how to perform bootstrap based inference on these bounds in section (ref).
As shown in section (ref), all of our bounds can be constructed from the marginal distribution of $W$ and the bounds on $Q_{Y_x \mid W}$ given in equations (ref) and (ref). These bounds on $Q_{Y_x \mid W}$, in turn, depend on just two features of the data:
In both cases, we can use parametric, semiparametric, or nonparametric estimation methods. In this paper we focus on flexible parametric approaches. Even in this case the asymptotic distribution theory is non-standard and quite complicated. We discuss this point further in the conclusion, section (ref).
In section (ref) we describe our first step estimators of these two functions. Given these estimators, we then construct sample analog estimates of our bound functions in a second step. We describe these estimators in section (ref).
We estimate $Q_{Y \mid X,W}$ by a linear quantile regression of $Y$ on flexible functions of $(X,W)$ that we denote by $q(X,W) \in \ensuremath{\mathbb{R}}^{d_q}$. For example, $q(x,w)$ could be $(1,x,w)$, $(1,x, w, x \cdot w)$, or could contain additional interactions between the treatment indicator $X$ and functions of the covariates $W$. For $\tau \in (0,1)$, let
be the estimated coefficients from a linear quantile regression of $Y$ on $q(X,W)$ at the quantile $\tau$. Here $\rho_\tau(s) = s (\tau - \ensuremath{\mathbbm{1}}(s<0))$ is the check function. Let $\widehat{Q}_{Y \mid X,W}(\tau \mid x,w) = q(x,w)'\widehat{\gamma}(\tau)$ denote this estimator.
We estimate the propensity score by maximum likelihood. In particular, specify the parametric model \[ \ensuremath{\mathbb{P}}(X=1 \mid W=w) = F( r(w)'\beta_0) \] where $F$ is a known cdf, $r(w)$ is a known vector function, and $\beta_0$ is an unknown constant vector. The functions $r(w)$ could simply be $(1,w)$ or may contain functions of $w$, like squared or interaction terms. For notational simplicity, we will assume throughout the paper that $r(w) = w$. Given this assumption, the dimension of $\beta_0$ is $d_W$, the length of $W$. Suppose $\beta_0$ lies in the parameter space $\mathcal{B} \subseteq \ensuremath{\mathbb{R}}^{d_W}$.
This specification for the propensity score includes the probit and logit estimators as special cases. Those estimators are commonly used in the literature; for example, see chapter 13 of ImbensRubin2015. Let $\widehat{\beta}$ denote the maximum likelihood estimate of $\beta$: \[ \widehat{\beta} = \operatorname*{argmax}_{\beta \in \mathcal{B}} \sum_{i=1}^n \log L(X_i,W_i'\beta) \] where \[ L(x,w'\beta) = F(w'\beta)^x(1-F(w'\beta))^{1-x}. \] For each $x\in \{0,1\}$, let $\widehat{p}_{x|w} = L(x,w'\widehat\beta)$ denote our propensity score estimator.
Given the first step estimators from section (ref), we obtain the following sample analog estimators of the CQTE bound functions defined in equations (ref) and (ref):
and
As discussed in section (ref), averaging these over $\tau \in (0,1)$ yields sample analog estimates of bounds on $\text{CATE}(w)$, which we can then use to get bounds on $\text{ATE}$. This approach requires estimation of extremal quantiles, however---estimation for $\tau$'s close to 0 or 1. This is well known to be a delicate problem (see ChernozhukovFernandezValKaji2017 for details). So in this paper we use a common solution: Fixed trimming of the extremal quantiles. We do this by modifying the quantile bound estimators to ensure that the quantile index lies in $[\varepsilon,1-\varepsilon]$ for some fixed and known $\varepsilon \in (0,0.5)$. Specifically, this yields the trimmed estimators of the quantile bounds
and
We use these estimators for the rest of the paper. Common choices of $\varepsilon$ are $0.05$ or $0.01$. In our asymptotic analysis we hold $\varepsilon$ fixed with sample size. In principle we could generalize the results to allow $\varepsilon \rightarrow 0$ as $n \rightarrow \infty$, but this would complicate the analysis of inference, which is already non-standard for other reasons. Since we fix $\varepsilon$ throughout, we omit $\varepsilon$ from the notation for brevity, except when necessary.
We next estimate the CQTE bounds by taking differences of the quantile bound estimators:
Since our CATE bounds are simply the integral of the CQTE bounds over all the quantiles $\tau$, we can estimate them by \[ \left[\widehat{\underline{\text{CATE}}}^c(w),\widehat{\overline{\text{CATE}}}^c(w) \right] \equiv \left[\int_0^1 \widehat{\underline{\text{CQTE}}}^c(\tau \mid w) \; d\tau, \int_0^1 \widehat{\overline{\text{CQTE}}}^c(\tau \mid w) \; d\tau \right]. \] A second integration over $w$ with respect to the marginal distribution of $W$ yields bounds on $\text{ATE}$. Like much of the literature, we use the empirical distribution of $W$ to estimate the marginal distribution of $W$. This yields the following estimator of our $\text{ATE}$ bounds: \[ \left[ \widehat{\underline{\text{ATE}}}^c, \widehat{\overline{\text{ATE}}}^c \right] = \left[ \frac{1}{n}\sum_{i=1}^n \widehat{\underline{\text{CATE}}}^c(W_i), \, \frac{1}{n}\sum_{i=1}^n \widehat{\overline{\text{CATE}}}^c(W_i) \right]. \] Next consider the estimation of the $\text{ATT}$ bounds. Let \[ \widehat{\underline{E}}_0^c = \frac{1}{n}\sum_{i=1}^n \int_0^1 \widehat{\underline{Q}}_{Y_0}^c(\tau \mid W_i) \; d\tau \qquad \text{and} \qquad \widehat{\overline{E}}_0^c = \frac{1}{n}\sum_{i=1}^n \int_0^1 \widehat{\overline{Q}}_{Y_0}^c(\tau \mid W_i) \; d\tau. \] For $x \in \{0,1\}$ let \[ \widehat{\ensuremath{\mathbb{E}}}(Y \mid X=x) = \frac{\sum_{i=1}^n Y_i \ensuremath{\mathbbm{1}}(X_i=x)}{\sum_{i=1}^n \ensuremath{\mathbbm{1}}(X_i=x)} \qquad \text{and} \qquad \widehat{p}_x = \frac{1}{n}\sum_{i=1}^n \ensuremath{\mathbbm{1}}(X_i = x). \] We can then estimate the ATT bounds by replacing the population quantities in (ref) with their estimators that we just defined.
For $c = 0$, our estimated upper and lower bounds are equal and give point estimates of the various parameters of interest. For $c > 0$, our bounds have positive width. To use these bounds in a sensitivity analysis, we recommend producing the following plot: Pick a grid $\{ c_1,\ldots,c_K \} \subseteq [0,1]$ of values for $c$. Compute our bound estimates on this grid and plot them against these values of $c$. Then compute and plot confidence bands for these bound estimates against $c$ as well; we describe how to compute these bands in section (ref). We illustrate all of these steps in our empirical analysis of section (ref).
In this section we provide formal results on the consistency and limiting distributions of the estimators we described in section (ref). In section (ref) we show how to use these results to do inference based on a non-standard bootstrap. In section (ref) we provide sufficient conditions under which standard bootstrap inference is valid.
Throughout this paper we assume that we observe a random sample.
Our first step estimators are standard in the literature. Hence we only briefly review the main assumptions and results for these estimators. For completeness, we provide a formal analysis in appendix (ref).
We assume that both the propensity score and quantile regression functions are correctly specified: \[ \ensuremath{\mathbb{P}}(X=x\mid W=w) = L(x,w'\beta_0) \] and \[ Q_{Y \mid X,W}(\tau \mid x,w) = q(x,w)'\gamma_0(\tau) \] for all $\tau \in [\varepsilon,1-\varepsilon]$.
Since the first step estimators consist of linear quantile regression and maximum likelihood estimation, their $\sqrt{n}$-convergence to Gaussian elements can be shown under standard assumptions and arguments. For example, see NeweyMcFadden1994. Moreover, the convergence of $\widehat{\gamma}(\tau)$ to $\gamma_0(\tau)$ is uniform over $\tau \in [\varepsilon,1-\varepsilon]$. Formally, as we show in appendix (ref) lemma (ref), \[ \sqrt{n}
\rightsquigarrow \mathbf{Z}_1(\tau), \] where $\mathbf{Z}_1(\cdot)$ is a mean-zero Gaussian process in $\ensuremath{\mathbb{R}}^{d_W} \times \ell^\infty([\varepsilon,1-\varepsilon],\ensuremath{\mathbb{R}}^{d_q})$ with continuous paths. The covariance kernel of this process is defined in appendix (ref), equation (ref). Also see appendix (ref) for the formal assumptions under which this result holds.
Next we consider the limiting distribution of our various second step estimators.
We start with equations (ref) and (ref), our estimators of the conditional quantile bounds. The population conditional quantile bounds, equations (ref) and (ref), are known functions of $\theta_0 = (\beta_0,\gamma_0)$. Define
Throughout the paper, we let $\Gamma_j = (\overline{\Gamma}_j, \underline{\Gamma}_j)$ for $j\geq 1$. Evaluating these at $\theta_0$ gives the trimmed population conditional quantile bounds. Evaluating these at $\widehat{\theta}$ gives their sample analog estimators. Define \[ \overline{\Gamma}_2(x,w,\theta) = \int_0^1 \overline{\Gamma}_1(x,w,\tau,\theta) \;d\tau \qquad \text{ and } \qquad \underline{\Gamma}_2(x,w,\theta) = \int_0^1 \underline{\Gamma}_1(x,w,\tau,\theta) \;d\tau. \] Then \[ \left[\underline{\text{CATE}}_\varepsilon^c(w),\overline{\text{CATE}}_\varepsilon^c(w) \right] = \Big[ \underline{\Gamma}_2(1,w,\theta_0)- \overline{\Gamma}_2(0,w,\theta_0), \, \overline{\Gamma}_2(1,w,\theta_0) - \underline{\Gamma}_2(0,w,\theta_0) \Big] \] are the trimmed population CATE bounds. We estimate them by \[ \left[\widehat{\underline{\text{CATE}}}^c(w),\widehat{\overline{\text{CATE}}}^c(w) \right] \equiv \left[\underline{\Gamma}_2(1,w,\widehat\theta)- \overline{\Gamma}_2(0,w,\widehat\theta), \, \overline{\Gamma}_2(1,w,\widehat\theta) - \underline{\Gamma}_2(0,w,\widehat\theta)\right]. \] If these mappings were Hadamard differentiable in $\theta$ at $\theta_0$, we could use the functional delta method to show that the above estimators have limiting Gaussian distributions and converge at $\sqrt{n}$ rates. Because they depend on the $\min$ and $\max$ functions these mappings are not Hadamard differentiable. They are, however, Hadamard directionally differentiable (HDD); see definition (ref) in appendix (ref). It turns out that this weaker version of differentiability is sufficient to establish their (non-Gaussian) limiting distribution.
To formally derive the limiting distribution of the CATE estimators, we show that the mapping \[ \Gamma_2(x,w,\cdot): \ensuremath{\mathbb{R}}^{d_W} \times \ell^\infty([\varepsilon,1-\varepsilon],\ensuremath{\mathbb{R}}^{d_q}) \to \ensuremath{\mathbb{R}}^2 \] is Hadamard directionally differentiable at $\theta_0$ tangentially to $\ensuremath{\mathbb{R}}^{d_W} \times \mathscr{C}([\varepsilon,1-\varepsilon],\ensuremath{\mathbb{R}}^{d_q})$. Here $\mathscr{C}(A,B)$ is the set of continuous functions from $A$ to $B$.
As a technical assumption, we restrict the complexity of the space that the quantile regression coefficient $\gamma_0(\cdot)$ lives in. Specifically, we assume that it is in a H\"older ball. To define this parameter space precisely, let $\mathscr{C}_m(\mathcal{D})$ denote the set of $m$-times continuously differentiable functions $f : \mathcal{D} \rightarrow \ensuremath{\mathbb{R}}$, where $m$ is an integer and $\mathcal{D}$ be an open subset of $\ensuremath{\mathbb{R}}^{d_q}$. Denote the differential operator by \[ \nabla^\lambda = \frac{\partial^{|\lambda |}}{\partial x_1^{\lambda_1} \cdots \partial x_{d_q}^{\lambda_{d_q}}} \] where $\lambda = (\lambda_1,\ldots,\lambda_{d_q})$ is a $d_q$-tuple of nonnegative integers and $| \lambda | = \lambda_1 + \cdots + \lambda_{d_q}$. Let $\nu \in (0,1]$. Let $\| \cdot \|$ without any subscripts denote the $\ensuremath{\mathbb{R}}^{d_q}$-Euclidean norm. Define the H\"older norm of $f: \mathcal{D} \rightarrow \ensuremath{\mathbb{R}}$ by \[ \| f \|_{m,\infty,\nu} = \max_{| \lambda | \leq m} \sup_{x \in \text{int}(\mathcal{D})} | \nabla^{\lambda} f(x) | + \max_{| \lambda | = m} \sup_{x,y \in \text{int} (\mathcal{D}), x \neq y} \frac{ | \nabla^\lambda f(x) - \nabla^\lambda f(y) |}{\| x - y \|^\nu}. \] For any $B > 0$, let $\mathscr{C}_{m,\nu}^B(\mathcal{D}) = \{ f \in \mathscr{C}_m(\mathcal{D}) : \| f \|_{m, \infty, \nu} \leq B \}$ denote a H\"older ball.
In this assumption we assume $m \geq 3$ to obtain bounded third derivatives of $\gamma_0$. In appendix (ref) we state several additional standard regularity conditions that we use to obtain asymptotic normality of the first step estimators; see assumptions A(ref)--A(ref) starting on page (ref). We continue to maintain these assumptions here. As a first preliminary result, we use these assumptions to derive the limiting distribution of the CQTE bound estimators; see proposition (ref) in appendix (ref). Using that result, we can then derive the limiting distribution of the CATE bound estimators.
In the statement of this result we deferred the full characterization of $\mathbf{Z}_\text{CATE}(w)$ to the proof. To get a brief idea of what it looks like, however, consider the first component. It is \[ \textbf{Z}_{\text{CATE}}^{(1)}(w) = \overline{\Gamma}_{2,\theta_0}'(1,w,\mathbf{Z}_1) - \underline{\Gamma}_{2,\theta_0}'(0,w,\mathbf{Z}_1) \] where $\overline{\Gamma}_{2,\theta_0}'(x,w,\mathbf{Z}_1)$ is the Hadamard directional derivative of $\overline{\Gamma}_{2}$ evaluated at $\mathbf{Z}_1$, the limiting distribution of the first step estimators. See page (ref) for the expression for $\overline{\Gamma}_{2,\theta_0}'$. Likewise, $\underline{\Gamma}_{2,\theta_0}'(x,w,\mathbf{Z}_1)$ is the Hadamard directional derivative of $\underline{\Gamma}_2$ evaluated at $\mathbf{Z}_1$. Although $\mathbf{Z}_1$ is Gaussian, the HDDs are continuous but generally nonlinear functionals. Hence the distribution of $\mathbf{Z}_\text{CATE}(w)$ is non-Gaussian. In section (ref) we show how to use a non-standard bootstrap to approximate its distribution.
Next we derive the limiting distribution of our ATE bound estimators. Let \[ \overline\Gamma_3(x,\theta) = \int_\mathcal{W} \overline{\Gamma}_2(x,w,\theta) \; dF_W(w) \quad \text{ and } \qquad \underline\Gamma_3(x,\theta) = \int_\mathcal{W} \underline{\Gamma}_2(x,w,\theta) \; dF_W(w). \] Then \[ [\underline{\text{ATE}}_\varepsilon^c, \overline{\text{ATE}}_\varepsilon^c ] = \left[ \underline{\Gamma}_3(1,\theta_0) - \overline{\Gamma}_3(0,\theta_0), \,\overline{\Gamma}_3(1,\theta_0) - \underline{\Gamma}_3(0,\theta_0) \right] \] are the trimmed population ATE bounds. We estimate them by \[ \left[\widehat{\underline{\text{ATE}}}^c,\widehat{\overline{\text{ATE}}}^c \right] \equiv \left[\frac{1}{n}\sum_{i=1}^n\left(\underline{\Gamma}_2(1,W_i,\widehat\theta)- \overline{\Gamma}_2(0,W_i,\widehat\theta)\right), \, \frac{1}{n}\sum_{i=1}^n\left(\overline{\Gamma}_2(1,W_i,\widehat\theta) - \underline{\Gamma}_2(0,W_i,\widehat\theta)\right)\right]. \] Unlike $\Gamma_2$, the $\Gamma_3$ mapping depends on $F_W$, which is unknown. Here we estimate it by the empirical distribution of $W$.
Next, let $\delta > 0$ and define \[ \mathcal{B}_\delta = \{ \beta \in \mathcal{B} : \|\beta - \beta_0\| \leq \delta\} \qquad \text{and} \qquad L_\beta(x,w'\beta) = \frac{\partial}{\partial\beta}L(x,w'\beta). \] The following assumption bounds the inverse ratio of the squared propensity score to its derivative with respect to the parameter $\beta$. This assumption holds under common parametric specifications for the propensity score, like logit or probit. It also holds if strong overlap holds. Moreover, under our other assumptions, note that strong overlap holds when $W$ has finite support.
Under these assumptions, we show the following result.
Like the CATE bound estimators, the limiting distribution of the ATE bound estimators is non-Gaussian. To understand this limiting distribution, first recall that we denote our bounds on the means $\ensuremath{\mathbb{E}}(Y_x)$ by \[ \overline{E}_{x,\varepsilon}^c = \overline{\Gamma}_3(x,\theta_0) = \ensuremath{\mathbb{E}}[ \overline{\Gamma}_2(x,W,\theta_0)] \qquad \text{and} \qquad \underline{E}_{x,\varepsilon}^c = \underline{\Gamma}_3(x,\theta_0) = \ensuremath{\mathbb{E}}[ \underline{\Gamma}_2(x,W,\theta_0)]. \] We estimate them by \[ \widehat{\overline{E}}_x^c = \frac{1}{n}\sum_{i=1}^n \overline{\Gamma}_2(x,W_i,\widehat{\theta}) \quad \text{ and } \quad \widehat{\underline{E}}_x^c = \frac{1}{n}\sum_{i=1}^n \underline{\Gamma}_2(x,W_i,\widehat{\theta}). \] In the proof of theorem (ref), we show that the following asymptotic expansion holds:
where $\Gamma_{3,\theta_0}'$ is the Hadamard directional derivative of $\Gamma_3$, which we define in the proof of theorem (ref). The first term in this expansion comes from the sample variation in the first step estimators: the propensity score $\widehat{p}_{x|w} = L(x,w'\widehat{\beta})$ and quantile function $\widehat{Q}(\tau \mid x,w) = p(x,w)'\widehat{\gamma}(\tau)$. The functional $\Gamma_{3,\theta_0}'(x,\cdot)$ is nonlinear in $\widehat{\beta}$. Therefore, since $\sqrt{n}(\widehat{\beta}-\beta_0)$ converges in distribution to a Gaussian limiting process, the limiting distribution of this functional is non-Gaussian. If $\beta_0$ was known---and hence the propensity score was known---then this component would follow a Gaussian distribution since the remaining component $\widehat{\gamma}$ is asymptotically Gaussian and enters $\Gamma_{3,\theta_0}'(x,\cdot)$ linearly.
The second term in this expansion comes from the variation of the CATE bounds over the values of the covariates $W$. It follows a limiting Gaussian distribution by the central limit theorem. This term is asymptotically independent of the sampling variation in the first step estimators $\widehat{\theta}$ since the influence function of $\widehat{\theta}$ is mean independent of $W$. Overall, we see that the limiting distribution of the ATE bounds is the sum of two independent random vectors, one Gaussian and one non-Gaussian. We approximate the distribution of these two random vectors using two separate bootstraps in section (ref).
Finally we study the limiting properties of our ATT bound estimators. Our trimmed population ATT bounds are \[ [\underline{\text{ATT}}_\varepsilon^c, \overline{\text{ATT}}_\varepsilon^c] = \left[\ensuremath{\mathbb{E}}(Y \mid X=1) - \frac{\overline{E}_{0,\varepsilon}^c - p_0 \ensuremath{\mathbb{E}}(Y \mid X=0)}{p_1}, \, \ensuremath{\mathbb{E}}(Y \mid X=1) - \frac{\underline{E}_{0,\varepsilon}^c - p_0\ensuremath{\mathbb{E}}(Y \mid X=0)}{p_1}\right]. \] We estimate them by \[ [\widehat{\underline{\text{ATT}}}^c, \widehat{\overline{\text{ATT}}}^c] = \left[\widehat\ensuremath{\mathbb{E}} (Y \mid X=1) - \frac{\widehat{\overline{E}}_0^c - \widehat{p}_0 \widehat\ensuremath{\mathbb{E}} (Y \mid X=0)}{\widehat{p}_1}, \, \widehat\ensuremath{\mathbb{E}} (Y \mid X=1) - \frac{\widehat{\underline{E}}_0^c - \widehat{p}_0 \widehat\ensuremath{\mathbb{E}} (Y \mid X=0)}{\widehat{p}_1} \right], \] where $\widehat{\ensuremath{\mathbb{E}}}(Y \mid X=0)$ and $\widehat{p}_0$ are defined in section (ref).
$\widehat{\ensuremath{\mathbb{E}}}(Y \mid X=0)$ and $\widehat{p}_0$ are asymptotically Gaussian. Like our analysis of the ATE bounds, however, $\widehat{\overline{E}}_0^c$ and $\widehat{\underline{E}}_0^c$ have non-Gaussian asymptotic distributions. Overall, $\mathbf{Z}_\text{ATT}$, the asymptotic distribution of our ATT bound estimators, is a linear combination of Gaussian and non-Gaussian random variables.
We now show how to conduct inference on our bounds for CATE, ATE, and ATT. Earlier we noted that these bounds are generally not ordinary Hadamard differentiable mappings of the underlying parameters $\theta_0$. By corollary 3.1 in FangSantos2014, this implies that standard bootstrap approaches cannot be used for these bounds. We instead use the non-standard bootstrap approach developed by FangSantos2014. For brevity we focus on ATE and ATT in this section. We provide analogous results for CQTE and CATE in lemmas (ref) and (ref) in appendix (ref).
The bounds for ATE and ATT can be written in terms of bounds on $\ensuremath{\mathbb{E}}(Y_x)$. In this section we describe how to do inference on bounds for these means. We'll then use these results to do inference on our ATE and ATT bounds in the next subsection. Recall that our bounds on $\ensuremath{\mathbb{E}}(Y_x)$ can be written as a functional of $\theta$. This functional is Hadamard directionally differentiable in $\theta$, but it is generally not ordinary Hadamard differentiable. Theorem 3.1 of FangSantos2014 shows how to do bootstrap inference by consistently estimating the Hadamard directional derivative (HDD). This can be done by using analytical estimators or by using a numerical derivative as described in HongLi2015. Here we use analytical estimates of the HDD. This approach explicitly uses the functional form of the HDD to estimate it. It allows us to avoid picking the numerical derivative step size, although other tuning parameters are used to estimate the HDDs analytically.
Next we define some general notation. Let $Z_i = (Y_i,X_i,W_i)$ and $Z^n = \{Z_1,\ldots,Z_n\}$. Let $\vartheta_0$ denote some parameter of interest and let $\widehat{\vartheta}$ be an estimator of $\vartheta_0$ based on the data $Z^n$. Let $\mathbf{A}_n^*$ denote $\sqrt{n}(\widehat{\vartheta}^* - \widehat{\vartheta})$ where $\widehat{\vartheta}^*$ is a draw from the nonparametric bootstrap distribution of $\widehat{\vartheta}$. Suppose $\ensuremath{\mathbf{A}}$ is the tight limiting process of $\sqrt{n} ( \widehat{\vartheta} - \vartheta_0)$. Denote bootstrap consistency by $\mathbf{A}_n^* \overset{P}{\rightsquigarrow} \mathbf{A}$ where $\overset{P}{\rightsquigarrow}$ denotes weak convergence in probability, conditional on the data $Z^n$. Weak convergence in probability conditional on $Z^n$ is defined as \[ \sup_{h\in \text{BL}_1} \left| \ensuremath{\mathbb{E}} [h(\mathbf{A}_n^*) \mid Z^n] - \ensuremath{\mathbb{E}} [h(\mathbf{A})] \right| = o_p(1) \] where $\text{BL}_1$ denotes the set of Lipschitz functions into $\ensuremath{\mathbb{R}}$ with Lipschitz constant no greater than 1. We leave the domain of these functions and its associated norm implicit.
We focus on the choices $\vartheta_0 = \theta_0$ and $\widehat{\vartheta} = \widehat{\theta}$. For these choices, let $\mathbf{Z}_n^* = \sqrt{n}(\widehat{\theta}^* - \widehat{\theta})$. Let $\mathbf{Z}_1$ denote the limiting distribution of $\sqrt{n}(\widehat{\theta} - \theta_0)$; see lemma (ref) in appendix (ref). Theorem 3.6.1 of VaartWellner1996 implies that $\mathbf{Z}_n^* \overset{P}{\rightsquigarrow} \mathbf{Z}_1$. Our parameters of interest are all functionals $\Gamma$ of $\theta_0$. In particular, in section (ref) we showed that \[ \sqrt{n}(\Gamma(\widehat{\theta}) - \Gamma(\theta_0)) \rightsquigarrow \Gamma'_{\theta_0}(\mathbf{Z}_1) \] for a variety of functionals $\Gamma$. To do inference on these functionals, we therefore want to estimate the distribution of $\Gamma_{\theta_0}'(\mathbf{Z}_1)$. FangSantos2014 show that \[ \widehat{\Gamma}_{\theta_0}'(\mathbf{Z}_n^*) \overset{P}{\rightsquigarrow} \Gamma'_{\theta_0}(\mathbf{Z}_1) \] where $\widehat{\Gamma}_{\theta_0}'$ is a suitable estimator of the Hadamard directional derivative $\Gamma'_{\theta_0}$. In this section we construct the estimators $\widehat{\Gamma}_{\theta_0}'$ and show that they can be used in this bootstrap.
Next, recall the asymptotic expansion in equation (ref) on page (ref). As we will show, the second term in this expansion can be approximated using standard bootstrap approaches and replacing $\theta_0$ by $\widehat{\theta}$. The first term requires estimating the HDDs $\overline{\Gamma}_{3,\theta_0}'$ and $\underline{\Gamma}_{3,\theta_0}'$. The formulas for our estimators of these HDDs are long, and so we describe them in appendix (ref). Denote these estimators by $\widehat{\overline{\Gamma}}_{3,\theta_0}'$ and $\widehat{\underline{\Gamma}}_{3,\theta_0}'$. They require choosing two scalar tuning parameters, $\kappa_n$ and $\eta_n$. $\kappa_n$ is a slackness parameter and $\eta_n$ is a step size parameter used to compute numerical derivatives of $\widehat{\gamma}(\cdot)$. Although not used in our proof, the asymptotic independence of the two components implies that approximating their respective marginal distributions is sufficient to obtain their joint distribution.
As we just mentioned, we'll use the standard nonparametric bootstrap to approximate the second term of equation (ref). To formalize this, let $\mathbb{G}_n^*$ denote the nonparametric bootstrap empirical process: \[\mathbb{G}_n^* = \frac{1}{\sqrt{n}}\sum_{i=1}^n (M_{n,i} - 1) \delta_{Z_i}\] where $(M_{n,1},\ldots,M_{n,n})$ are multinomially distributed with parameters $(1/n,\ldots,1/n)$ independently of $Z^n$, and where $\delta_{Z_i}$ is a distribution which assigns probability one to the value ${Z_i} \in \ensuremath{\mathbb{R}}^{2+d_W}$. Then for any function $g$, \[ \mathbb{G}_n^* g(Z) = \sqrt{n}\left(\frac{1}{n} \sum_{i=1}^n g(Z_{i}^*) - \overline{g(Z)}\right) \] where $Z_{i}^*$, $i=1,\ldots,n$ are drawn independently with replacement from $\{Z_1,\ldots,Z_n\}$ and $\overline{g(Z)} = \frac{1}{n}\sum_{i=1}^n g(Z_i)$. In particular, we'll study the asymptotic distribution of \[ \mathbb{G}_n^* \Gamma_2(x,W,\widehat{\theta}) = \sqrt{n}\left(\frac{1}{n} \sum_{i=1}^n \Gamma_2(x,W_i^*, \widehat{\theta}) - \frac{1}{n} \sum_{i=1}^n \Gamma_2(x,W_i,\widehat{\theta}) \right). \]
The following proposition is our main bootstrap consistency result.
This result shows how to use the bootstrap to approximate the joint limiting distribution of upper and lower bounds of $\ensuremath{\mathbb{E}}(Y_x)$, $x\in\{0,1\}$. As we show in section (ref) below, these approximations can be used to conduct pointwise or uniform-in-$c$ inference on the ATE bounds. As part of the proof, we show that $\widehat{\Gamma}_{3,\theta_0}'(x, \sqrt{n}(\widehat{\theta}^* - \widehat{\theta}))$ weakly converges in probability conditional on the data to $\Gamma_{3,\theta_0}'(x,\ensuremath{\mathbf{Z}}_1)$, a non-Gaussian vector which reflects the sample variation in the first step estimators. We also show weak convergence in probability conditional on the data of $\mathbb{G}_n^* \Gamma_2(x,W,\widehat{\theta})$ to $\mathbb{G} \Gamma_2(x,W,\theta_0) \sim \ensuremath{\mathcal{N}}(0, \operatorname*{var}(\Gamma_2(x,W,\theta_0))$, a bivariate Gaussian vector which reflects the variation of the CATE bounds over $W$. This variation can be approximated using the standard nonparametric bootstrap. Hence the bounds' limiting distribution is approximated by a combination of standard and non-standard bootstraps. Note that the two bootstrap distributions can be computed from a unique sequence of draws $Z_i^*$ and so has the same computational burden as a single bootstrap.
Next we show how to use proposition (ref) to do inference on our ATE bounds $[\underline{\text{ATE}}_\varepsilon^c, \overline{\text{ATE}}_\varepsilon^c]$. We first consider inference pointwise in $c$. We then construct confidence bands that are uniform over $c$.
An immediate corollary of proposition (ref) is
Thus we can also use this specific bootstrap to approximate the asymptotic distribution of our ATE bounds estimators. Given this result, we can construct a $100(1-\alpha)$% confidence set for the ATE identified set under $c$-dependence as follows. Let \[ \text{CI}_\text{ATE}^c(1-\alpha)= \left[\widehat{\underline{\text{ATE}}}^c - \frac{\widehat{d}_\alpha}{\sqrt{n}}, \, \widehat{\overline{\text{ATE}}}^c + \frac{\widehat{d}_\alpha}{\sqrt{n}}\right] \] where \[ \widehat{d}_\alpha = \inf\left\{z\in\ensuremath{\mathbb{R}}: \ensuremath{\mathbb{P}} \left(\sqrt{n} (\widehat{\underline{\text{ATE}}}^{c,*} - \widehat{\underline{\text{ATE}}}^c ) \leq -z \; \text{ and } \; \sqrt{n} (\widehat{\overline{\text{ATE}}}^{c,*} - \widehat{\overline{\text{ATE}}}^c ) \geq z \mid Z^n \right) \geq 1-\alpha \right\}. \] The probability in this expression can be approximated by taking a large number of bootstrap draws according to equation (ref). Proposition (ref) then implies that \[ \liminf_{n \rightarrow \infty} \ensuremath{\mathbb{P}} \Big( \text{CI}_\text{ATE}^c(1-\alpha) \supseteq [\underline{\text{ATE}}_\varepsilon^c,\overline{\text{ATE}}_\varepsilon^c] \Big) \geq 1-\alpha. \]
Let \[d_\alpha = \inf\left\{z\in\ensuremath{\mathbb{R}}: \ensuremath{\mathbb{P}} \left(\textbf{Z}_\text{ATE}^{(2)} \leq -z \; \text{ and } \; \textbf{Z}_\text{ATE}^{(1)} \geq z\right) \geq 1-\alpha\right\}.\] If $\ensuremath{\mathbb{P}} (\textbf{Z}_\text{ATE}^{(2)} \leq -z \; \text{ and } \; \textbf{Z}_\text{ATE}^{(1)} \geq z)$ is continuous and strictly increasing in a neighborhood of $d_\alpha$, corollary 3.2 in FangSantos2015workingPaper yields $\widehat{d}_\alpha = d_\alpha + o_p(1)$ and hence \[ \lim_{n \rightarrow \infty} \ensuremath{\mathbb{P}} \Big( \text{CI}_\text{ATE}^c(1-\alpha) \supseteq [\underline{\text{ATE}}_\varepsilon^c,\overline{\text{ATE}}_\varepsilon^c] \Big) = 1-\alpha. \]
We just described how to use proposition (ref) to do inference on the ATE bounds for any fixed $c$. Those results can be immediately extended to do inference the ATE bounds for any finite grid of $c$'s. In this section we show how to construct confidence bands that are uniform over all $c \in [0,1]$. We do this by using monotonicity of the ATE bound functions in $c$. This lets us extrapolate bands that are uniform on a finite grid in such a way that they have uniform coverage. A related procedure is described in corollary 1 of MastenPoirier2020.
Although $\overline{\text{ATE}}_\varepsilon^c$ is nondecreasing in $c$, its estimate $\widehat{\overline{\text{ATE}}}^c$ may be nonmonotonic in $c$ because of the quantile crossing problem with linear quantile regression. In that case, we could monotonize the estimated ATE bound function by using the rearrangement procedure of ChernozhukovFernandez-ValGalichon2010, for example. As they show, the rearrangement operator is Hadamard directionally differentiable, and thus can be accommodated in our inferential results. Likewise, $\underline{\text{ATE}}_\varepsilon^c$ is nonincreasing in $c$ and $\widehat{\underline{\text{ATE}}}^c$ is also nonincreasing after applying a suitable rearrangement. From here on we assume our bound estimators have been monotonized. Note that this does not affect the asymptotic distribution of the estimators, by corollary 1 of ChernozhukovFernandez-ValGalichon2010.
Next consider a grid of values $\mathcal{C} = \{c_1,\ldots,c_K\}$ such that $0 = c_1 <\cdots < c_K = 1$. All of our analysis works with $0 < c_1$ and $c_K < 1$, but typically researchers will want to include the endpoints of $[0,1]$ so we do that from here on. Using methods similar to those in section (ref), for all $c \in \mathcal{C}$ let \[ \text{CI}_\text{ATE}^c(1-\alpha) = \left[\widehat{\underline{\text{ATE}}}^c - \frac{\widehat{d}_\alpha(c)}{\sqrt{n}},\widehat{\overline{\text{ATE}}}^c + \frac{\widehat{d}_\alpha(c)}{\sqrt{n}}\right] \] be $100(1-\alpha)$% confidence sets for $[\underline{\text{ATE}}_\varepsilon^c,\overline{\text{ATE}}_\varepsilon^c]$ where the critical values $\widehat{d}_\alpha(c)$ are chosen such that these sets are uniform over $c$ in the finite grid $\mathcal{C}$. That is, \[ \ensuremath{\mathbb{P}}\left( \text{CI}_\text{ATE}^c(1-\alpha) \supseteq [\underline{\text{ATE}}_\varepsilon^c,\overline{\text{ATE}}_\varepsilon^c] \text{ for all $c \in \mathcal{C}$} \right) \rightarrow 1-\alpha \] as $n \rightarrow \infty$. Finally, let $\min(c) = \inf \{ c_k \in \mathcal{C} : c \leq c_k \}$ denote the smallest element in the grid $\mathcal{C}$ that is still larger than $c$. For $c \in [0,1]$ define \[ \widehat{\text{UB}}(c) = \widehat{\overline{\text{ATE}}}^{\min(c)} + \frac{\widehat{d}_\alpha(\min(c))}{\sqrt{n}} \qquad \text{and} \qquad \widehat{\text{LB}}(c) = \widehat{\underline{\text{ATE}}}^{\min(c)} - \frac{\widehat{d}_\alpha(\min(c))}{\sqrt{n}}. \] $\widehat{\text{UB}}(c)$ is the greatest monotonic interpolation of the upper bounds of the confidence intervals on the grid $\mathcal{C}$. $\widehat{\text{LB}}(c)$ is the least monotonic interpolation of the lower bounds of the confidence intervals on the grid $\mathcal{C}$. By the definition of these interpolated bands and by monotonicity of the population ATE bounds,
as $n \rightarrow \infty$.
In this subsection we've shown that, although we cannot obtain the limiting distribution of the ATE bounds uniformly over $c \in [0,1]$, the fact that these bounds are monotonic lets us nonetheless do inference on them uniformly over $[0,1]$. This monotonicity comes from the nested nature of $c$-dependence: $c_1$-dependence implies $c_2$-dependence when $c_1 \leq c_2$. This kind of monotonicity is common in many other approaches to sensitivity analysis and hence greatest and least monotonic interpolations can likely be used more broadly to construct uniform confidence bands.
Bootstrap inference on the ATT bounds is quite similar to that on the ATE bounds. By examining the ATT bounds' limiting distribution (see the proof of proposition (ref)) we see that it depends on two types of terms:
Combining all the terms gives \[
\overset{P}{\rightsquigarrow} Z_ATT. \] We can use this result to construct pointwise confidence sets for the ATT bounds for a fixed $c$, or to construct confidence bands that are uniform on a finite grid $\mathcal{C}$. Like the ATE bounds, the ATT bounds are monotonic in $c$. Thus a similar interpolation can be used to construct confidence bands for the ATT bounds that are uniform over $c \in [0,1]$.
We conclude this section be showing how to use the confidence bands we just described to do inference on breakdown points. For brevity we focus on the breakdown point for the conclusion that the ATE is nonnegative, which we defined earlier in equation (ref). Inference on other breakdown points for other conclusions can be done similarly.
Let $\text{CI}_\text{ATE}^c(1-\alpha)$ be a pointwise-in-$c$ confidence band for the ATE bounds, as described in section (ref). Define \[ c_L = \sup \{ c \in [0,1] : \text{CI}_\text{ATE}^c(1-\alpha) \subseteq [0,\infty) \}. \] This is simply the value at which the confidence band first intersects the horizontal line at zero. By proposition S.2 in Appendix D of MastenPoirier2020, \[ \lim_{n \rightarrow \infty} \ensuremath{\mathbb{P}}( c_L \leq c_\textsc{bp} ) \geq 1 - \alpha. \] Thus $[c_L, 1]$ is a valid one-sided lower confidence interval for the breakdown point $c_\textsc{bp}$.
In the previous section we showed how to use a non-standard bootstrap method to conduct inference on the CATE, ATE, and ATT bounds. The key technical problem was that these bounds are not necessarily Hadamard differentiable functionals of the first step estimators; they are only Hadamard directionally differentiable. In this section, we provide simple sufficient conditions on the propensity score under which the CATE, ATE, and ATT bounds are in fact Hadamard differentiable. Under this condition, the methods in section (ref) are still valid, but so is the standard nonparametric bootstrap. After stating the formal result, we discuss when this sufficient condition holds and when it does not.
First consider the average treatment effect. Recall that the ATE bounds depend on the functional $\Gamma_3(x,\theta)$. We will show that its Hadamard directional derivative $\Gamma_{3,\theta_0}'(x, h)$ is linear in $h$ under a condition on the value of $c$, the propensity score $p_{1 | w}$, and the distribution of $W$. By proposition 2.1 in FangSantos2014, this linearity is equivalent to Hadamard differentiability. By theorem 3.9.11 in VaartWellner1996, this linearity also implies that the bootstrap process $\sqrt{n}(\Gamma_{3,\theta_0}(x, \widehat{\theta}^*) - \Gamma_{3,\theta_0}(x, \widehat{\theta}))$ converges weakly in probability conditional on the data to $\Gamma_{3,\theta_0}'(x, \mathbf{Z}_1)$, a Gaussian vector. In other words, we can conduct inference on the ATE bounds using the standard nonparametric bootstrap: Take $n$ independent draws from the data with replacement, compute the bound estimates in this bootstrap sample, and then use the distribution of these bound estimates across many such bootstrap samples to approximate the sampling distribution of the bound estimators. The following theorem provides the explicit sufficient condition for validity of this bootstrap. In this result, let $p_{1 \mid W} = \ensuremath{\mathbb{P}}(X=1 \mid W)$ denote the random variable obtained by evaluating the propensity score at the random vector $W$.
In the proof of this result we show that when the propensity score does not contain a point mass on either $c$ or $1-c$, the mapping $\Gamma_3(x,\theta_0)$ is Hadamard differentiable for $x \in \{0,1\}$. Hence the nonparametric bootstrap is valid. Although it is not formally stated in the theorem, we conjecture that our sufficient condition for validity of the standard bootstrap is also a necessary condition. That is, we expect the standard bootstrap to be invalid when $c$ or $1-c$ are point masses of the propensity score's distribution. From the proof of theorem (ref), we see when this condition on the propensity score fails, $\Gamma'_{1,\theta_0}(x,w,\tau,h)$ is nonlinear in $h$ on a set of $(\tau,w)$ values of positive measure. Since the HDD of $\Gamma_3(x,\cdot)$ is the integral over $(\tau,w)$ of the HDD of $\Gamma_1(x,w,\tau,\cdot)$, we expect that $\Gamma_{3,\theta_0}'(x,h)$ will also be nonlinear in $h$, a failure of Hadamard differentiability.
By further examining the proof of theorem (ref), we can also show that the CATE bounds are Hadamard differentiable at covariate values $w$ and sensitivity parameter values $c$ such that $p_{1|w} \notin \{c,1-c\}$. Finally, the following proposition gives a similar result for the ATT, using slightly weaker assumptions.
The ATT bounds only depend on our bounds for $\ensuremath{\mathbb{E}}(Y_0)$, and not our bounds for $\ensuremath{\mathbb{E}}(Y_1)$. Hence we only need to examine $\Gamma_3(x,\theta_0)$ for $x=0$. So the proof of this proposition proceeds by showing that $\Gamma_3(0,\theta_0)$ is Hadamard differentiable when the propensity score does not have a point mass at $c$.
The sufficient conditions in theorem (ref) and proposition (ref) depend on the support of the propensity score $p_{1 \mid W}$. If the propensity score's distribution is absolutely continuous, it contains no point masses and therefore these condition holds. This holds when one covariate $W_k$ has nonzero coefficient $\beta_{0,k}$ and has a continuous distribution conditional on the other covariates $W_{-k}$. However, if all covariates are discrete or mixed, the support of the propensity score will generally contain point masses. The nonparametric bootstrap may not be valid whenever $c$ coincides with these points. Even when $p_{1|W}$ has point masses, the nonparametric bootstrap is valid for $c$ outside of these points. To use this bootstrap, one could in principle estimate the support of $p_{1|W}$ to determine at which values of $c$ inference might be invalid, and select sensitivity parameters outside of this support.
The nonparametric bootstrap has the advantage of being computationally simple and does not require the choice of tuning parameters. While more involved, the bootstrap technique detailed in section (ref) is valid regardless of the support of the propensity score. For example, in our empirical analysis in section (ref), all of the covariates are either mixed or discrete. Given our analysis above, we therefore use the non-standard bootstrap in our empirical analysis since the standard bootstrap may fail at some values of $c$.
In this section we illustrate our methods using data on the National Supported Work (NSW) demonstration project studied by LaLonde1986. Since this is a highly studied and well-known program, we only briefly summarize it here. See, for example, HLS1999 for further details. We use LaLonde's data as reconstructed by DehejiaWahba1999.
The NSW demonstration project randomly assigned participants to either receive a guaranteed job for 9 to 18 months along with frequent counselor meetings or to be left in the labor market by themselves. We use the DehejiaWahba1999 sample, which are all males in LaLonde's NSW dataset and where earnings are observed in 1974, 1975, and 1978. This dataset has 445 people: 185 in the treatment group and 260 in the control group. Like Imbens2003, we use this experimental sample primarily as an illustration; in experiments where treatment was truly randomized it is not necessary to assess sensitivity to unconfoundedness. Our results may be useful for assessing the impact of randomization failure in experiments, but that is not our focus here.
In addition to this experimental sample, we construct a sample using observational data. This sample combines the 185 people in the NSW treatment group with 2490 people in a control group constructed from the Panel Study of Income Dynamics (PSID). This control group, called PSID-1 by LaLonde, consists of all male household heads observed in all years between 1975 and 1978 who were less than 55 years old and who did not classify themselves as retired. We further drop observations with earnings above \$5,000 in 1974, 1975, or both. This leaves 148 treated units (out of 185) and 242 untreated units (out of 2490). This observational sample was also considered by Imbens2003.
The outcome of interest is earnings in 1978. There are also nine covariates: Earnings in 1974, earnings in 1975, years of education, age, indicators for race (Black, Hispanic, other), an indicator for marriage, an indicator for having a high school degree, and an indicator for treatment. All earnings variables are measured in 1982 dollars. Table (ref) shows the summary statistics, as reported in table 1 of Imbens2003.
Table (ref) shows the baseline point estimates of both ATE and ATT under the unconfoundedness assumption in the two samples we consider. These estimates are all computed by inverse probability weighting (IPW) using a parametric logit propensity score estimator. We do not consider other estimators, since our goal is to illustrate sensitivity to identifying assumptions, rather than finite sample sensitivity to the choice of estimator.
Figure (ref) shows our main results. These are estimated treatment effect bounds under $c$-dependence, along with corresponding pointwise confidence bands, as described in sections (ref)--(ref). The top plot shows bounds on ATE while the bottom plot shows bounds on ATT. The solid lines are bounds for the observational dataset while the dashed lines are bounds for the experimental dataset. The light dotted lines are confidence bands for the observational dataset while the light dashed-dotted lines are confidence bands for the experimental dataset. These bands are constructed to have nominal 95% coverage probability pointwise in $c$ based on our non-standard bootstrap results in section (ref). For the tuning parameters we use $\varepsilon = 0.05$, $\eta_n = 0.05 n^{-1/4}$, and $\kappa_n = n^{-1/3}$. Note that the sufficient conditions for validity of the standard bootstrap that we gave in section (ref) do not apply here, since the distribution of the propensity score variable $p_{1 \mid W}$ has point masses. This occurs because seven of the nine covariates are discrete, while the other two mixed discrete-continuous. The mixed variables are earnings in 1974 and earnings in 1975, which have point masses at zero since many people in the sample did not work in those years.
For both datasets, at $c=0$ the bounds collapse to the baseline point estimate. When $c>0$, we allow for some selection on unobservables. Comparing the shape of the bounds for both datasets we see that the experimental data are substantially more robust to relaxations from the baseline assumptions than the observational data. Specifically, for most values of $c$ the bounds for the experimental data are substantially tighter than the bounds for the observational data. Even the no assumptions bounds ($c = 1$) are tighter for the experimental data than for the observational data.
A second way to measure robustness uses breakdown points. MastenPoirier2020 discuss these in detail and give additional references. In the current context, the breakdown point is simply the largest value of $c$ such that we can no longer draw a specific conclusion about some parameter. Specifically, in the next two subsections we consider two conclusions: The conclusion that ATE is nonnegative, and the conclusion that ATT is less than the per participant program cost.
First consider the conclusion that ATE is nonnegative. Our point estimates support this conclusion, but does it still hold if the baseline unconfoundedness assumption fails? In the experimental dataset, the estimated breakdown point is $0.082$. This is simply the value of $c$ such that the lower bound function in figure (ref) intersects the horizontal axis. For all $c \leq 0.082$, the estimated identified sets for ATE only contains nonnegative values. For $c > 0.082$, the estimated identified sets contain both positive and negative values. Hence, for such relaxations of unconfoundedness, we cannot be sure that the average treatment effect is positive.
For the observational dataset, the estimated breakdown point for the conclusion that ATE is nonnegative is $0.037$. This is more than twice as small as the breakdown point for the experimental dataset. Hence again we see that conclusions about ATE from the experimental dataset are substantially more robust than the observational dataset. The same conclusion holds for ATT: The point estimates in both datasets suggest that it is positive. But how robust is that conclusion? The estimated breakdown point for the conclusion that ATT is nonnegative in the experimental data is $0.123$ while it is $0.049$ for the observational dataset. By this measure, the conclusion that ATT is positive is more than twice as robust using the experimental data compared to the observational data.
Thus far we have compared the robustness of results obtained from the experimental data with results obtained from the observational data. Next we discuss whether either of these results are robust in an absolute sense. To do this, we use the leave-out-variable $k$ propensity score analysis discussed in section (ref).
First consider table (ref), which uses data from the experimental sample. For each variable $k$, listed in the rows of this table, we compute four summary statistics from the estimated distribution of \[ \Delta_k = | p_{1 \mid W}(W_{-k},W_k) - p_{1 \mid W_{-k}}(W_{-k}) |. \] Specifically, we estimate the 50th, 75th, and 90th percentiles of $\Delta_k$, along with the maximum observed value, denoted $\bar{c}_k$. As discussed in section (ref), these quantities tell us about the marginal impact of covariate $k$ on treatment assignment. $c$-dependence constrains the maximum value of the marginal impact of the unobserved potential outcome on treatment assignment, above and beyond the observed covariates. Thus the values in table (ref) can help us calibrate $c$. Specifically, we will compare the breakdown point to the values in this table. These values could be interpreted as upper bounds on the magnitude of selection on unobservables that we might think is present. Thus, for a given reference value from this table, if the breakdown point is larger than the reference value, we could consider the conclusion of interest to be robust to failure of unconfoundedness. In contrast, if the breakdown point is smaller than the reference value, we could consider the conclusion of interest to be sensitive to failure of unconfoundedness.
Recall that the estimated breakdown point for the conclusion that ATE is nonnegative is $0.082$. This is larger than three of the $\bar{c}_k$ values and on the same order of magnitude as four more. If we look at a less stringent comparison, the 90th percentile, we see that the estimated breakdown point is now larger than all but one of the rows, corresponding to the indicator for Hispanic. Let's examine this variable more closely. Figure (ref) plots the density $\Delta_k$ for $k=$ Hispanic indicator. Here we see that there is a small proportion of mass who have values larger than $0.082$, but most people have values well below the breakdown point. Next suppose we weaken the criterion even more by considering the 75th percentile column in table (ref). The breakdown point is larger than all values in this column.
The leave-out-variable $k$ propensity score analysis focuses on the relationship between observed covariates and treatment assignment. It does not use data on outcomes. A less conservative analysis is to only worry about covariates $k$ which have large values in table (ref) and which also affect our outcomes in some way. Specifically, we next consider leave-out-variable $k$ IPW estimates of ATE under the baseline unconfoundedness assumption. Table (ref) shows the effect of leaving out a single variable on the ATE point estimates for both datasets. Continue to consider just the experimental dataset. Here we first see that omitting any covariate at most changes the point estimate by 5.4%. Moreover, recall the main variable we were concerned about before: the indicator for Hispanic. Omitting this variable only changes the ATE point estimate by 1.5%.
Overall, the leave-out-variable $k$ analysis suggests that, on an absolute scale, the conclusion that ATE is nonnegative using the experimental data is quite robust. A similar analysis applies to conclusions about ATT.
Next consider the observational data. Table (ref) shows the leave-out-variable $k$ propensity score analysis. Recall that the estimated breakdown point for the conclusion that ATE is nonnegative in the observational dataset is 0.037. By any of these measures the conclusion that ATE is nonnegative is not robust. Suppose we only consider variables which also substantially change the point estimates, as shown in table (ref). Even then we still find that the results are sensitive. For example, the indicator for Black changes the ATE point estimate by 14% and also has substantial marginal impact on the propensity score, with its 50th percentile in table (ref) about 1.5 times as large as the estimated ATE breakdown point. Thus, using these as absolute measures of robustness, we find that the conclusion that ATE is positive using the observational data is not robust.
This conclusion that findings based on the observational dataset are not robust contrasts with the sensitivity analysis of Imbens2003, who finds that the same observational dataset yields relatively robust results. Imbens' analysis relied importantly on fully parametric assumptions about the joint distribution of the observables and unobservables. In particular, he assumed outcomes were normally distributed, that the treatment effect is homogeneous, and that any selection on unobservables arises due to an omitted binary variable. Our identification analysis does not require any of these assumptions. As discussed in section (ref), we do impose some parametric assumptions to simplify estimation, but even these assumptions are substantially weaker than those used by Imbens. Given that we are making weaker auxiliary assumptions, it is not surprising that our analysis shows the findings to be more sensitive than the analysis in Imbens2003. Nonetheless, even with these weaker assumptions, we continue to find that conclusion from the experimental dataset remain robust.
Finally, note that all of our discussion thus far has focused on the point estimates of the breakdown points. In section (ref) we showed that the value at which the pointwise confidence band intersects the horizontal axis is a valid one-sided lower confidence interval for the breakdown point. For the ATE with experimental data, this gives a confidence set of $[0.0156,1]$, with a point estimate of $0.082$. For the ATE with observational data, this gives a confidence set of $[0.009,1]$, with a point estimate of $0.037$. Thus the lower bound of the confidence interval for the experimental data is almost twice as large as the lower bound for the observational data. So the relative comparison of the two datasets continues to hold once we account for sampling uncertainty. Unfortunately, the lower bound of $0.0156$ for the experimental dataset is quite small, if we compare it to the variation in the leave-out-variable-$k$ propensity scores. This is not surprising though, given that there is a substantial amount of sampling uncertainty---even the lower bound of the confidence intervals for the baseline estimates are quite close to zero.
In the previous subsection we studied the sensitivity of the conclusion that the ATE is nonnegative. In practice, however, this is not necessarily the most policy relevant conclusion. For example, HeckmanSmith1998 give a model where the socially optimal decision whether to continue a small scale program or to shut it down can be computed by comparing the ATT with the program's per participant cost. In this subsection we show how our sensitivity analysis can be used in these kinds of cost-benefit analyses. Specifically, we consider the conclusion that the ATT is less than the per participant program cost. Under the model in HeckmanSmith1998, the program should be shut down when this conclusion holds.
Chapter 8 of MDRC MDRC1983 reports NSW per participant program costs. For males, total costs ranges between \$4,637 and \$5,218 (tables 8-2, 8-3, and 8-4) in 1976 dollars. Our treatment effect estimates are in 1982 dollars. Adjusting these reported costs to 1982 dollars (CPI-U series) gives a range of \$7,865 to \$8,850.
First consider the experimental dataset. The conclusion of interest holds for the baseline estimate: The ATT of \$1,738 is far less than the per participant cost. Suppose, however, that a supporter of the program claims that this baseline estimate is implausible due to selection on unobservables. How strong does selection on unobservables need to be to allow for the possibility that the program is cost effective? Formally, what is the smallest $c$ such that the identified set for ATT includes values that are larger than the per participant cost? If we look at the bounds' point estimates, there are no values of $c$ under which the program is cost effective. Accounting for sampling uncertainty by examining the confidence bands, we need $c$ to be at least about 0.5 before it is possible that the program is cost effective. As we argued earlier, these are very large values, so it is unlikely that selection on unobservables is this strong.
Next consider the observational dataset. Here again the conclusion of interest holds for the baseline estimate: The ATT of \$4,001 is smaller than the per participant cost. Next consider the breakdown point for this conclusion. Since the bounds for the observational dataset are larger than those for the experimental dataset, we need less selection on unobservables to allow for possibly large values of the ATT. Despite this, there are still no values of $c$ under which the program is cost effective, based on the bounds' point estimates. This is largely because the uncertainty due to the impact of selection on unobservables is asymmetric in this example: The lower bound grows much faster in $c$ than the upper bound does. Hence conclusions about the largest possible value of the ATT are more robust to relaxations of unconfoundedness than conclusions about the smallest possible value of the ATT. If we account for sampling uncertainty by examining the confidence bands, then we need $c$ to be at least about 0.08 before the confidence intervals contain ATT values larger than the per participant costs. This is a relatively large value, although it is smaller than a decent number of the leave-out-variable $k$ propensity score values in table (ref). That, however, likely just reflects the large amount of sampling uncertainty in this data.
Overall, our analysis suggests that the program does not pass a cost-benefit analysis, even if we allow for a large amount of selection on unobservables. Hence the conclusion that the ATT is less than the per participant cost, and hence that the program should be shut down, is quite robust to failures of unconfoundedness.
Finally, note that our analysis here is primarily illustrative. A more comprehensive cost-benefit analysis would require examining many other program outcomes besides just short run post-program earnings. For example, see the analysis in chapter 8 of MDRC MDRC1983 and section 10 of HLS1999. Note, however, that given data on these additional outcomes, our methods could then be used to analyze the sensitivity of total program impacts to failures of unconfoundedness.
Identification, estimation, and inference on treatment effects under unconfoundedness has been widely studied and applied. This approach uses two assumptions: Unconfoundedness and Overlap. The overlap assumption is refutable, and many tools have been developed for checking this assumption in practice. For example, Stata's built-in package teffects has commands for checking overlap. In this paper, we provide a complementary suite of tools for assessing the unconfoundedness assumption. There are two key distinctions between our results and the previous literature. First, we begin from fully nonparametric bounds. In contrast, most of the previous literature relies on parametric assumptions for their identification analysis. Second, we provide tools for inference. This is important because, just like baseline estimators, sensitivity analyses are also subject to sampling uncertainty.
We conclude by discussing several extensions and directions for future work. As we just mentioned, a key distinguishing feature of our sensitivity analysis is that we begin from fully nonparametric bounds. We then estimated these bounds using flexible parametric estimators of the propensity score and the quantile regression of outcomes on treatment and covariates. These estimators can include quadratic terms, cubic terms, and interactions, for example, but they are not fully nonparametric. We restricted attention to parametric estimators for one reason: Even in this case, the asymptotic distribution theory is non-standard, complicated, and at the frontier of current research. This difficulty comes from the fact that our estimands are not Hadamard differentiable. Extending our analysis to first step nonparametric estimators is an important next step, but doing so will likely require both deriving and applying more general asymptotic theory for non-Hadamard differentiable functionals than currently exists. Hence we leave that analysis to future work.
A second extension is to consider additional parameters of interest. In this paper we focus on estimation and inference on the ATE and ATT bounds. We also developed analogous results for the CQTE and CATE. The conditional average treatment effect for the treated, $\text{CATT}(w) = \ensuremath{\mathbb{E}}(Y_1 - Y_0 \mid X=1, W=w)$, can be studied with the same tools we use in section (ref). We omit that analysis for brevity. MastenPoirier2018 also derive sharp bounds on unconditional quantile treatment effects (QTEs). Estimation and inference on the QTE bounds is more complicated than the ATE and ATT bounds. The reason is identical to the explanation Vaart2000 gives when discussing inference on unconditional sample quantiles: “to derive the asymptotic normality of even a single quantile estimator $\widehat{F}_n^{-1}(p)$, we need to know that the estimators $\widehat{F}_n$ are asymptotically normal as a process, in a neighborhood of $F^{-1}(p)$.” In our case, performing inference on the QTE bounds requires showing convergence of corresponding bounds on the unconditional potential outcome cdfs as a process in a neighborhood of the quantile of interest.\footnote{MastenPoirier2020 prove some results along these lines; see their lemma 1. Those results are only valid for sufficiently small values of $c$ and with discrete $W$, which substantially simplifies the analysis.} For this reason, we leave estimation and inference on the QTE bounds to a separate paper.
\singlespacing
\allowdisplaybreaks