Bounding Treatment Effects by Pooling Limited Information across Observations

Bounding Treatment Effects by Pooling Limited Information across Observations

Introduction

\todo[size=tiny]{The page number starts from the title page as required; The latex option `baselinestretch' is modified to fit the main text within 45 pages.}

In many applications, causal inference hinges on strong ignorability, namely unconfoundedness and overlap imbens2015causal. The former condition is non-testable but requires that all confounders be used as covariates; the latter is a testable condition that may not be satisfied in practice.

The overlap condition has received increasing attention in the literature. In applications, it is not uncommon to have a situation where the estimated propensity scores are close to zero or one. This problem is referred to as limited overlap crump2009dealing. The existence of limited overlap may change the asymptotic behavior of the estimators Khan:Tamer:2010,Hong:et:al:2019 and may necessitate using a more robust inference method Rothe:2017,sasaki_ura_2021. DAMOUR2020 provide a cautionary tale on the overlap condition when high-dimensional covariates are adopted to make unconfoundedness more plausible.

There are several approaches in the literature to estimate treatment effects when facing limited overlap. Arguably, the most popular method is to focus on a subpopulation where the overlap condition holds crump2009dealing,Yang:Ding:18. For example, crump2009dealing recommend a simple rule of thumb to drop all observations with estimated propensity scores outside the range $[\alpha,1-\alpha]$ for some predetermined constant $\alpha$, say $\alpha = 0.1$. Alternatively, Li-et-al:2018 advocate the use of the so-called `overlap weights' to define the average treatment effect. This amounts to assigning weights equal to one minus the propensity score for the treated units and equal to the propensity score for the control units. If the treatment effects are heterogeneous, both trimming and overlap weighting change the parameter of interest from the population average treatment effect. Without changing it, nethery2019 develop a Bayesian framework by extrapolating estimates from the overlap region to the non-overlap region via a spline model. However, identification by extrapolation is subject to model misspecification.

In this paper, we start with the observation that none of the aforementioned papers would work well if the overlap condition is not satisfied at the population level and it is a priori unknown where it fails. In that case, the population average treatment effect is not point-identified and one may resort to manski1989anatomy,manski1990nonparametric's bounds, provided that the support of outcome is bounded and known. However, it may not yield tight bounds if unconfoundedness assumption is plausible, while the overlap condition being the only source of identification failure. This paper provides a systematic method to explore this possibility.

Our contributions are two-fold. First, we provide novel \todo[size=tiny]{wording change: “population bounds” to “bounds”} bounds on both average treatment effects (ATE) and average treatment effects on the treated (ATT) that are valid under an unconfoundedness assumption. Our bounds are applicable if the conditioning variables do not satisfy the overlap condition and take on a large number of different values in the observed sample. This robustness is achieved by only using limited “pooling” of information across observations. Namely, the bounds are constructed as the expectations of functions of the observed outcomes such that the contribution of each outcome only depends on the treatment status of a limited number of observations. No information pooling across observations leads to manski1989anatomy,manski1990nonparametric's bounds, which we call “first-order bounds”, while unlimited information pooling leads to standard inverse propensity score weighting. We explore the intermediate range between these two extremes by considering the setup where an applied researcher provides a reference propensity score. Our bounds are valid independent of the value of this reference propensity score, but if it happens to be close to the true propensity score, then our bounds are optimal in terms of expected width within the class of limited pooling bounds considered in this paper. The reference propensity score is therefore crucial to construct our novel treatment effect bounds uniquely, and it also allows to incorporate prior knowledge on the propensity score in a robust way.

Second, we develop estimation and inference methods for the bounds we have established under the unconfoundedness assumption. Our formal theory assumes that the observed covariates are discrete, so that multiple observations can share the same covariate value. A leading data scenario we analyze assumes that the number of distinct covariate values is large relative to the sample size, implying that for each possible covariate value only a small number of observations are available. In this scenario, it is a statistically challenging problem to provide a valid confidence interval for the treatment effects, which we tackle in this paper.

In many empirical applications, however, some or all covariates are continuous rather than discrete. To apply our method in such settings, one must first discretize the covariates, for example by binning or clustering observations with similar covariate values. This is a practical solution that we discuss in Section (ref) and employ in our Monte Carlo experiments and empirical applications. However, discretization introduces an approximation error for which we do not provide formal theory. Controlling this error would require additional smoothness assumptions on the conditional mean of the outcome variable, which we do not impose. This is a limitation of our approach in its current form, and we leave the formal analysis of discretization bias for future work.

An alternative approach to robust inference for treatment effects under unconfoundedness is provided by Armstrong:Kolesar:21. In particular, their confidence intervals are asymptotically valid under a violation of the overlap condition, as long as the researcher specifies a Lipschitz bound on the conditional mean of the outcome variable. Their approach is distinct from and complementary to ours. The approach of Armstrong:Kolesar:21 reduces to a matching estimator for the average treatment effect (e.g., abadie2006large,abadie2008failure,abadie2011bias) if the Lipschitz bound is chosen to be very large. Those matching estimators crucially require that for every observation we can find other observations with similar covariate values but opposite treatment status. This is not required in our approach. Crucially, we only pool information across observations with similar covariate values, but in contrast to Armstrong:Kolesar:21 and matching estimators, we do so completely independent of the treatment status of the observations involved. This is the key difference compared to those existing methods.

The remainder of the paper is organized as follows. In Section (ref), we describe the setup and intuition behind our approach. Section (ref) illustrates our key ideas through the simple two-unit example and introduces the main ideas, including a formal characterization of our “second-order bounds”. In Section (ref), we extend the framework to the general bounds of arbitrary order, and in Section (ref), we construct sample analogs and develop corresponding inference methods. Using those bounds we then provide asymptotically valid confidence intervals. We discuss how to cluster the covariate observations in Section (ref). The results of Monte Carlo experiments are reported in Section (ref). In Section (ref), we present two empirical applications. The first uses the well-known dataset from Connors1996's study of the efficacy of right heart catheterization (RHC), which has been extensively analyzed in the context of limited overlap (see, e.g., crump2009dealing, Rothe:2017, Li-et-al:2018). The second application uses the dataset from DehejiaWahba-JASA, which exhibits limited overlap and serves as a useful complement. These applications illustrate the practical relevance and robustness of our method. The appendices contain all proofs, technical derivations omitted from the main text, and additional results from our Monte Carlo experiments. An accompanying R package is available on the Comprehensive R Archive Network (CRAN) at \url{https://CRAN.R-project.org/package=ATbounds}.

Setup

For units $i=1,\ldots,n$, we observe treatment status $D_i \in \{0,1\}$, regressors $X_i \in {\cal X}$, where ${\cal X}$ is a discrete set, and outcome $Y_i = (1-D_i) \, Y_i(0) + D_i \, Y_i(1)$, where $Y_i(0)$ and $Y_i(1)$ are potential outcomes. While we observe the realized outcome $Y_i$, we never observe both potential outcomes for the same unit. Our main objective is to conduct inference on the average treatment effect (ATE) and average treatment effect on the treated (ATT), conditional on the covariates:

We do not assume i.i.d.\ sampling. Instead, we allow the distribution of covariates $X^{(n)} := (X_1,\ldots,X_n)$ to vary with $n$. This leads us to define estimands conditional on $X^{(n)}$, but under our assumptions, the relevant expectations depend only on $X_i$.

Assumption (ref)(i) imposes unconfoundedness, meaning that treatment assignment is effectively randomized conditional on covariates $X_i$. While we could weaken this to mean independence (i.e.\ $ \mathbb{E}\left[ Y_i(d) \, \big| \, D_i , \, X_i \right] = \mathbb{E}\left[ Y_i(d) \, \big| \, X_i \right] $) for most of our results, scenarios that justify mean independence typically also support full conditional independence. Similarly, Assumption (ref)(ii) could be relaxed to bounds on conditional expectations (i.e.\ $a_{\min} \leq \mathbb{E}\left[ Y_i(d) \, \big| \, X_i \right] \leq a_{\max}$), but in practice, known bounds usually apply directly to the potential outcomes themselves.

Assumption (ref)(iii) has two key implications for the sampling scheme. First, it ensures that $(D_i,Y_i(0),Y_i(1))$ is independently distributed across units conditional on $X^{(n)}$. Second, it requires that the conditional distribution of $(D_i,Y_i(0),Y_i(1))$ depends only on $X_i$, not on other covariates $X_j$ or the unit index $i$. This implies that both treatment effects and propensity scores depend only on individual covariates: $ \mathbb{E}\left[ Y_i(1) - Y_i(0) \, \big| \, X^{(n)} \right] = \mathbb{E}\left[ Y_i(1) - Y_i(0) \, \big| \, X_i \right] $, and $\mathbb{E}\left( D_i \, \big| \, X^{(n)} \right) = \mathbb{E}\left( D_i \, \big| \, X_i \right) =: p(x)$. If the propensity score $p(x)$ were known and satisfied $0<p(x)<1$ (overlap condition), the treatment effects would be point-identified through inverse propensity score weighting:

However, since $p(x)$ is unknown in practice and the overlap condition may fail, we can generally only obtain partial identification of ATE and ATT. This means we can construct valid large-sample confidence intervals, but these may not converge to a point as $n \rightarrow \infty$.

Let $D_{-i}$ denote the treatment statuses of all units $j \neq i$ sharing the same covariate value as unit $i$. Our approach to constructing valid confidence intervals for ATE relies on finding functions $L(Y_i, D_i, D_{-i},X_i)$ and $U(Y_i, D_i, D_{-i},X_i)$ that provide bounds on the conditional treatment effect:

Let ${\cal X}_*$ denote the set of distinct observed covariate values with cardinality $m=|{\cal X}_*|$. In asymptotic sequences where $m \rightarrow \infty$ as $n \rightarrow \infty$, equation (ref) implies:

where the probability limit is taken conditional on $X^{(n)}$. These cross-sectional averages of the bounds form the basis for constructing asymptotically valid confidence intervals for ATE.

For ATT, we similarly construct bounds using functions of the form $L(Y_i, D_i, D_{-i}, X_i)$ and $U(Y_i, D_i, D_{-i}, X_i)$ that satisfy (ref) with $\pi(X_i)$ in place of $\tau(X_i)$. This allows us to bound the numerator $\frac{1}{n} \sum_{i=1}^n \pi(X_i)$ from (ref), while the denominator $\frac{1}{n} \sum_{i=1}^n \mathbb{E}(D_i|X_i)$ can be consistently estimated by $\frac{1}{n} \sum_{i=1}^n D_i$. This approach to ATT estimation requires only that $\frac{1}{n} \sum_{i=1}^n p(X_i) > 0$, a mild condition that permits $p(X_i)=0$ for many units.

The key advantage of our approach is that the asymptotic validity of the confidence intervals for ATE and ATT relies only on Assumption (ref), with the sole additional requirement of $\frac{1}{n} \sum_{i=1}^n p(X_i) > 0$ for ATT inference. Notably, we require neither assumptions on $X^{(n)}$, nor the overlap condition ($0 < p(x) < 1$), nor knowledge or consistent estimation of $p(x)$.

For each covariate value $X_i$, the treatment assignments $(D_i,D_{-i})$ are assumed to be independent Bernoulli draws with the same mean $p(X_i)$. However, since the number of units sharing any given covariate value may be small and non-increasing asymptotically, $p(X_i)$ may not be consistently estimable under our assumptions.

Manski bounds

Manski bounds (manski1989anatomy,manski1990nonparametric) represent a simplified version of (ref) and (ref) where $L(Y_i, D_i, D_{-i}, X_i)$ and $U(Y_i, D_i, D_{-i}, X_i)$ reduce to functions $L^{(1)}(Y_i, D_i)$ and $U^{(1)}(Y_i, D_i)$ that depend only on individual outcomes and treatment status. These bounds are particularly robust as they do not require unconfoundedness (Assumption (ref)(i)) and apply even when covariate values are unique to each unit. Using the outcome bounds $a_{\min}$ and $a_{\max}$, we can establish:

where

These lead to bounds for ATE

satisfying

For ATT, defining $C^{(1)}_a(Y_i, D_i) := D_i \, (Y_i - a)$, we have:

The bounds in (ref) and (ref) are well-known, and we denote those bounds on $ \tau(X_i)$ and $\pi(X_i) $ as either Manski bounds (manski1989anatomy,manski1990nonparametric) or as “first-order bounds”.

Pooling information across observations

The Manski bounds in (ref) and (ref) are sharp when we only impose outcome boundedness (Assumption (ref)(ii)). However, under unconfoundedness (Assumption (ref)(i)) and overlap ($0<p(x)<1$), both $\tau(x)$ and $\pi(x)$ become point-identified via inverse propensity score weighting as shown in (ref).

This presents us with two extremes. The Manski bounds require minimal assumptions but pool no information across observations, leading to potentially wide bounds. In contrast, point identification through propensity score methods requires extensive pooling of information across observations to estimate $p(x)$ consistently. This latter approach demands strong data requirements and faces the curse of dimensionality as the dimension of $X_i$ increases.

This paper explores a middle ground between these extremes. Our approach pools some information across observations to tighten the bounds on average treatment effects, but requires much less pooling than needed for consistent nonparametric estimation of $p(x)$. The key idea is to use information from pairs or small groups of observations with similar or identical covariate values to construct tighter bounds.

For example, when two observations share the same covariate value ($X_i = X_j$), we can construct “second-order" bounds that use information from both observations jointly. These bounds improve upon Manski bounds by leveraging unconfoundedness and using the treatment status of both observations with the same covariate value. Similar principles can be extended to construct higher-order bounds that pool information across larger groups of observations.

Lack of overlap and curse of dimensionality

While equation (ref) shows that ATE and ATT are point-identified under Assumption (ref) and overlap ($0<p(x)<1$), estimating $p(x)= \mathbb{E}\left( D \, \big| \, X=x \right)$ in finite samples presents significant challenges, particularly due to the curse of dimensionality for multi-dimensional covariates. Consider two illustrative examples in Figure (ref), each with sample size $n=100$.

In the left example, with $X_i \sim U[0,1]$ and $p(x)=x^4$, the overlap condition is theoretically satisfied for all $x \in (0,1)$. However, we observe no treated units ($D_i=1$) for $X_i < 0.58$, making precise point estimation of ATE infeasible without strong parametric assumptions. A natural approach here would combine Manski bounds for $X_i<0.58$ with matching or inverse propensity score weighting for $X_i \geq 0.58$.

The right example, with $X_i \sim U[0,1]^2$ and constant $p(x)=0.3$, illustrates a more complex challenge. Despite perfect theoretical overlap, large regions of the covariate space (e.g., near $x=(0,0)$) contain no treated observations. From the finite sample alone, we cannot determine whether this reflects true lack of overlap or merely finite-sample variation, as these regions contain few observations overall. This challenge intensifies with higher-dimensional covariates.

In this paper, we provide asymptotically valid inference on ATE or ATT under Assumption (ref) that remains robust to any form of the unknown propensity score $p(x)$, including $p(x)=0$ and $p(x)=1$ for some $x$. While some observations share identical (or similar) covariate values and thus propensity scores, each observation typically has few such neighbors. This scenario can be modeled asymptotically by having each covariate value appear exactly $k$ times ($k=2,3,4,\ldots$) in the sample, with $m=|{\cal X}_*|=n/k$ distinct values.

Our bounds solve this inference problem where $n/m$ may stay bounded as $n \rightarrow \infty$, addressing situations where limited local sample sizes prevent reliable propensity score estimation, as in the right panel of Figure (ref). Moreover, our method can incorporate prior information about the propensity score, potentially achieving point identification when this information is correct and overlap holds, while maintaining robust confidence intervals otherwise.

Main idea for second order bounds on $\mathbb{E}[ Y_i(1) ]$

To illustrate the fundamental concept of our new bounds, we begin with the simplest non-trivial case. Consider two units $i \in \{1,2\}$ sharing identical covariate values $X_1=X_2$. For each unit, we observe an outcome $Y_i \in \mathbb{R}$ and a treatment indicator $D_i \in \{0,1\}$. As introduced in Section (ref), each unit has potential outcomes $Y_i(0)$ and $Y_i(1)$, with the observed outcome being $Y_i=Y_i(D_i)$.

All stochastic statements in this section are implicitly conditional on $(X_1,X_2)$ and $X_1=X_2$. Let $P$ denote the corresponding conditional distribution of $[(D_i,Y_i(0),Y_i(1)) \, : i \in \{1,2\}]$, and let ${\cal P}$ be the set of distributions $P$ satisfying:

These conditions are restatements of Assumptions (ref) for the case where $X_1=X_2$. We denote expectations under $P$ by $\mathbb{E}_P$.

Our goal in this section is to construct a valid upper bound function $B(Y_1,D_1,D_2) \in \mathbb{R}$ such that, for all $P \in {\cal P}$,

While equation (ref) appears to treat units asymmetrically, our final bounds symmetrize across observations through averaging: $$ \frac 1 2 \left[ B(Y_1,D_1,D_2) + B(Y_2,D_2,D_1) \right] . $$ Beyond mere validity, we seek bound functions that cannot be improved upon, as formalized in the following definition.

The proof is given in the appendix. This theorem characterizes all functions $B(Y_1,D_1,D_2)$ whose expectations provide valid and non-dominated upper bounds on $\mathbb{E}[ Y_i(1) ]$ under our assumptions, showing they can be parameterized by $p_* \in (0,1]$. For identification purposes, we could optimize over $p_*$:

The last equality follows from choosing $p_* = P(D_i=1)$ in the minimum when positive. Through such “intersection bounds” we recover the well-known identification result for $\mathbb{E}_P[ Y_i(1) ]$ under unconfoundedness.

However, replicating this identification result is not our aim. As discussed previously, a key challenge in finite-sample inference is that the true propensity score $p = P(D_i=1)$ is unknown, varies with $X_i$, and may be close to zero --- leading to potentially large variances and non-uniformity in any corresponding estimator. Therefore, we focus on bounds that can be expressed as simple sample averages, as introduced in (ref), which is the context for Theorem (ref).

An illuminating special case arises when $p_*=1$, yielding:

When symmetrized across observations, this becomes:

This special case has an intuitive interpretation: When neither unit is treated $(D_1,D_2)=(0,0)$, we can only use the worst-case bound $a_{\max}$. However, when at least one unit is treated, we can use the corresponding treatment outcome as an estimate for $\mathbb{E}_P[ Y_i(1) ]$. This requires unconfoundedness since we select between $Y_1$ and $Y_2$ based on $(D_1,D_2)$. Equation (ref) provides the simplest example of what we term a “second-order bound”, where information is pooled across two observations using unconfoundedness.

These bounds share similarities with matching estimators, where outcomes with $D_1 \neq D_2$ for units sharing the same covariate value $X_1=X_2$ are matched to obtain counterfactual outcomes. The key distinction is that we do not require $D_1 \neq D_2$, necessitating worst-case bounds when $D_1=D_2=0$. However, this relaxation allows our bounds to remain valid without requiring overlap assumptions.

While the $p_*=1$ case yields straightforward bounds, Theorem (ref) reveals a richer family of bound functions parameterized by $p_* \in (0,1]$. To better understand this family, we introduce the true propensity score $p:=P(D_i=1)$ and define the weight function $w^{(2)}: [0,1] \times (0,1] \rightarrow (-\infty,1]$ as:

For $p \neq 0$ we then have\footnote{ The condition $p \neq 0$ is only required for our discussion here, not for Theorem (ref). This is because the $1/p$ in (ref) can be canceled against the factor $p$ in $ w^{(2)}(p,p_*) = p \cdot \frac {2p_* - p} {p_*^2} $ to avoid divison by zero when $p=0$. }

and consequently

Thus, in expectation, our bounds form weighted averages between the worst-case bound $a_{\max}$ and the target parameter $\mathbb{E}_P[Y(1)]$, with weights determined by $w^{(2)} \big(p ,p_* \big)$.

Figure (ref) plots $w^{(2)}(p,p_*)$ as a function of $p$ for various values of $p_*$. When $p_*$ equals the true propensity score $p$, we have $w^{(2)} \big(p ,p_* \big)=1$ and our upper bounds are sharp. When $p \neq p_*$, we have $w^{(2)} \big(p ,p_* \big) < 1$, yielding valid but non-sharp bounds.

For comparison, the Manski bounds have $\mathbb{E}\left[ a_{\max} + D_i (Y_i - a_{\max}) \right] = [1- p(x)] \, a_{\max} + p(x) \; \mathbb{E}\left[ Y(1) \right]$, corresponding to the weight function $w^{(1)}:p \mapsto p$, also shown in Figure (ref). The figure reveals that the bounds in (ref), corresponding to $p_*(x)=1$, are the only second-order bounds that uniformly dominate the Manski bounds across all possible values of the true propensity $p$.

Since we consider only undominated second-order bounds, none of the second-order bounds dominates any other second-order bound with a different $p_*$ value uniformly across data generating processes parameterized by $p$. If we have a reliable guess (or estimate) for the propensity score $p$, it should be used for $p_*$ to ensure reasonably tight bounds. The advantage of our bounds over alternatives (like inverse propensity score weighting) is their continued validity even when our guess (or estimate) for $p$ is incorrect.

It's worth noting that for $p_*<0.5$, the weights $w^{(2)} \big(p,p_* \big)$ become negative for large values of $p$, indicating that the expected bounds can perform worse than simply reporting $a_{\max}$, but nevertheless remaining valid.

The bounds introduced in this section illustrate the middle ground discussed in Section (ref) between minimal-assumption Manski bounds and full point identification through propensity score methods. By pooling information across pairs of observations, we achieve potential improvements over Manski bounds without requiring the extensive pooling needed for consistent propensity score estimation. We believe that this approach is particularly valuable in settings with high-dimensional covariates or limited local sample sizes, where reliable propensity score estimation may be infeasible but some degree of information pooling across observations remains possible. This two-observation case serves as a building block for the general bounds developed in subsequent sections.

Generalizations to higher order, ATE and ATT

We now return to the general setting with covariates introduced in Section (ref). The basic insights gained from analyzing the case of two observations sharing the same covariate value carry over naturally to this more general context, primarily requiring adjustments to notation.

To understand how the notation from Section (ref) extends to the current setting, note that the function $B(Y_1,D_1,D_2)$ introduced there becomes $B^{(2)}_{1,a}(Y_i,D_i,p(X_i),X_i)$ in our general notation here, where the superscript $(2)$ indicates that we are dealing with second-order bounds, the subscript 1 indicates that we are bounding $Y_i(1)$, and the choice of constant $a \in \{a_{\min},a_{\max}\}$ in the second subscript depends on whether we construct upper or lower bounds. The parameter $p_*(x)$ continues to play the role of $p_*$ but can now vary with the covariate value $x$. The weight function $w^{(2)}(p,p_*)$ remains unchanged but is now applied separately for each covariate value.

Short summary of second-order bounds

The main idea of our second-order bounds was already discussed in Section (ref) for bounds on $\mathbb{E}[Y_i(1)]$. The generalization to $\mathbb{E}[Y_i(0)]$ follows naturally by symmetry, which jointly provides second-order bounds for the ATE. While the main text presents results for both ATE and ATT throughout, we defer a more detailed discussion of the second-order ATT bounds to Appendix (ref).

To implement these bounds, we require the researcher to specify a reference propensity score $p_* : {\cal X} \rightarrow (0,1)$, which can be either postulated based on prior knowledge or estimated from the data. While the resulting bounds remain valid regardless of the choice of $p_*(x)$, their sharpness depends critically on how close the true propensity score is to $p_*(x)$. For all theoretical results that follow, we assume that $p_*(x)$ is non-random, with $p_*(x)=1/2$ serving as a natural default choice in the absence of strong prior information (see also Remark (ref) for practical guidance on choosing $p_*(x)$). Using this framework, we can now formally define our second-order bounds. For every $a \in \mathbb{R}$, we define

and

Here, we introduce the second-order bound functions as dependent on the unknown propensity score \( p(X_i) \). Since \( p(X_i) \) is not observed, these bounds are initially infeasible in this form. However, because all second-order bounds are linear in \( p(X_i) \), they can be rendered feasible by substituting \( p(X_i) \) with the treatment status \( D_j \) of another unit \( j \neq i \) that shares the same covariate value $X_j = X_i$. The resulting bound function \( B^{(2)}_{1,a}(Y_i,D_i,D_j,X_i) \) then corresponds exactly to the bound derived in Theorem (ref). Expressing the bounds as functions of \( p(X_i) \) in this section is advantageous, as it eliminates the need to reference other units explicitly. The following proposition summarizes the key properties of these second-order bounds.

The proof is provided in the appendix. Once we have bounds on $ \tau(x) $ and $ \pi(x) $, then we can also construct bounds on ATE and ATT defined in (ref).

Higher-order bounds

In Section (ref), we derived second-order bounds by directly considering pairs of observations sharing the same covariate value. For higher-order bounds, it is more convenient to first formulate bounds at the “population-level” --- that is, expressing the bounds directly in terms of the unknown propensity score $p(x)$ rather than in terms of treatment indicators of other observations. This formulation can be thought of as having already taken conditional expectations over the treatment indicators of other observations, leaving us with expressions that depend on $p(X_i)$ directly.

This “population-level” formulation was already used in the second-order bounds in (ref) and (ref) above, which are linear functions of $p(X_i)$. For higher-order bounds, we generalize this idea by allowing $p(X_i)$ to enter as higher-order polynomials. The resulting bounds will be infeasible since they depend on the unknown propensity score $p(X_i)$. However, they have feasible sample analogs: when implementing the bounds in practice, each power $[p(X_i)]^r$ in our population-level expressions can be replaced by products of $r$ different treatment indicators $D_j$ from observations sharing the same covariate value. The advantage of deriving the bounds in that form is that it allows us to understand the structure of the bounds before considering their sample implementation that requires multiple observations.

Dropping the index $i$ throughout, and the arguments $(Y_i,D_i,p(X_i),X_i)$ from the bound functions, our aim is to generalize the second-order bounds in (ref) by considering, for positive integers $q$,

where the coefficients $\lambda_{r,d}(x), \lambda_r(x) \in \mathbb{R}$ still need to be determined for $q>2$. Again, the motivation for (ref) is that, we can construct unbiased estimates for $B^{(q)}_{d,a}(\lambda) $ and $C^{(q)}_a(\lambda)$ by replacing $ [p(X)]^r$ with a product of treatment indicators from $r$ different observations with the same (or similar) regressor values.

Motivated by our finding for second-order bounds we again choose a reference propensity score $p_*(x)$ to find unique solutions for the coefficients $ \lambda_{r,d}(x)$ and $ \lambda_r(x)$. Once we have chosen $p_*(x)$, then for the second-order bounds the coefficients are uniquely determined by the properties of the bounds summarized in Proposition (ref) --- namely, the bounds should be valid for all population distributions satisfying Assumption (ref), and the bounds should be binding if $p(x)=p_*(x)$. However, for $q>2$ those properties are not sufficient anymore to uniquely determine the coefficients, because we now have additional degrees of freedom in the higher-order polynomial coefficients. To make use of this additional flexibility and to obtain unique coefficients again, we therefore demand the bounds to not only have good properties when $p(x)=p_*(x)$, but also when $p(x) \in [p_*(x)-\epsilon, p_*(x) + \epsilon]$, for small $\epsilon>0$, that is, we want to have good performance in a small neighborhood around the reference propensity score $p_*(x)$.

Specifically, we choose the optimal coefficients $\lambda_{r,d}(x)$ and $\lambda_r(x)$ such that the expected widths of the bounds

are minimized not only at $p(x)=p_*(x)$, but also when considering the worst-case expected widths within an infinitesimal neighborhood of the reference propensity score $p_*(x)$, see part (iii) of Proposition (ref) below for a formalisation of this. Here, $\mathbb{E}_{p(x)}$ refers to the expectation over $(Y_i,X_i,D_i)$ with propensity score (i.e.\ distribution of $D_i|X_i=x$) specified by $p(x)$.

Once we have solved for the optimal coefficients accordingly, we obtain the following optimal $B^{(q)}_{d,a}(\lambda) $ and $ C^{(q)}_a(\lambda) $, for integers $q \geq 1$,\footnote{ Formally, for $p(x)=1$ we have $B^{(q)}_{0,a} = a + \left\{ \frac {q - p_*(x) \, \mathbbm{1} \left\{ \text{$q$ is odd} \right\} } {1-p_*(x)} \right\} (1-D) (Y-a)$ and $ C^{(q)}_a = D\, (Y-a) + \left[ \frac {q-1 + p_*(x) \, \mathbbm{1} \left\{ \text{$q$ is even} \right\}} {1-p_*(x)} \right] (1-D) (Y-a)$. For $p(x)=0$ we have $B^{(q)}_{1,a} =a + \left\{ \frac {q - [1-p_*(x)] \, \mathbbm{1} \left\{ \text{$q$ is odd} \right\} } {p_*(x)} \right\} D (Y-a)$. From the formulas in (ref) we obtain those results for $p(x)=1$ and $p(x)=0$ as limits when $p(x) \rightarrow 1$ and $p(x) \rightarrow 0$. However, the details of those special cases do not actually matter, because e.g.\ for $p(x)=1$ we also have $D=1$ with probability one, and therefore $B^{(q)}_{0,a} = a$ and $C^{(q)}_a = D\, (Y-a) $. }

where the weight functions are given by

For $q=1$ and $q=2$ the formulas in (ref) just give the same functions $B^{(q)}_{d,a}$ and $ C^{(q)}_a $ that were already discussed above. It may not be obvious from those general formulas, but $B^{(q)}_{d,a}$ and $ C^{(q)}_a $ are indeed polynomials of order $(q-1)$ in $p(X)$. For example, for $q=3$ we find

which are all second order polynomials in $p(X)$.

We now want to formally state the optimality result for these bounds. We define $L^{(q)}(Y_i,D_i,p(X_i),X_i) $ and $ U^{(q)}(Y_i,D_i,p(X_i), X_i) $ as in (ref), but with superscipt $(2)$ replaced by $(q)$, and for $\epsilon \geq 0$ and $p_*(x) \in (0,1)$, we let ${\cal B}_\epsilon(p_*(x)) := \Big\{ p(x) \in [0,1] \, \Big| \, \left| p(x) - p_*(x) \right| \leq \epsilon \Big\}$ be the $\epsilon$-ball around $p_*(x)$.

The proof is given in the appendix. To better understand the result of Proposition (ref), consider the lower bound on $\mathbb{E}\left[ Y(1)\, \big| \, X\right] $, which is given by

Thus, $ \mathbb{E}\left[ B^{(q)}_{1,a_{\min}} \, \big| \, X=x \right] $ is a weighted average between $a_{\min}$ and $\mathbb{E}\left[ Y(1)\, \big| \, X=x \right] $. The weights always satisfy $w^{(q)}(p,p_*) \leq 1$, which together with $a_{\min} \leq Y(1)$ guarantees that $ \mathbb{E}\left[ B^{(q)}_{1,a_{\min}} \, \big| \, X=x \right] \leq \mathbb{E}\left[ Y(1)\, \big| \, X=x \right] $.

Figure (ref) shows $w^{(q)}(p,p_*)$ as a function of $p$ for $p_*=0.4$ and different values of $q$. For $p=0$ we always have $w^{(q)}(p,p_*)=0$, because in that case we only have observations with $D=0$ for $X=x$, implying that we cannot learn anything about $Y(1)$ from the data. For $p = p_*$ we have $w^{(q)}(p,p_*)=1$ for $q \geq 2$, that is, the lower bound is sharp in that case. For $p$ close to $p_*$ the weights are closer to one (implying that the bounds are sharper) the larger we choose $q$. For the $k$'th derivative of $w^{(q)}(p,p_*)$ at $p=p_*$ we have $$ \frac{\partial^k w^{(q)}(p_*,p_*)} {\partial^k p} = 0 , \qquad for \; \left\{

\right. $$ which explain why for $p$ close to $p_*$ the weights are closer to one the higher we choose $q$. However, if $p$ is far away from $p_*$, then the weights $w^{(q)}(p,p_*)$ for $q \geq 2$ can be far away from one, and can even be smaller than $w^{(1)}(p)$, that is, the bounds can be worse than Manski bounds if $p$ is far away from $p_*$.

The discussion for ATT bounds is analogous. In that case we have

that is, conditional on $X$, the ATT bounds are a linear combination between their Manski bounds and the true ATT contribution for $X$. Figure (ref) also shows the weights $\widetilde w^{(q)}(p,p_*)$ as a function of $p$, for $p_*=0.4$ and various values of $q$.

Implementation of the Bounds

In this section we construct sample analogs of the bounds in Proposition (ref), and use those sample bounds to obtain asymptotically valid confidence intervals on the average treatment effects. The bounds constructed in this section are valid for both discrete and continuous covariates $X_i$. However, if the covariates are continuously distributed, then every observed value $X_i$ is typically only observed once, in which case the bounds here simply become Manski worst-case bounds.

The interesting case, for the purpose of this section, is therefore the case where the set of possible covariate values ${\cal X}$ is discrete. However, we consider an asymptotic setting where the number of covariate values grows to infinity jointly with the total sample size. This is the challenging case from the perspective of treatment effect estimation, in particular when the average number of observations available for each observed $x \in {\cal X}$ remains small.

In Section (ref) we explain how the sample bounds for discrete covariate values from this section can be generalized to continuous covariate values via clustering, that is, by approximating the continuous set ${\cal X}$ with a finite set. In that way we obtain non-trivial bounds also for the case of continuous covariates.

Sample analogs of the bounds of Section (ref)

We require some additional notation to formulate the sample bounds. Analogously to $X^{(n)} :=(X_1,\ldots,X_n)$, we also define $D^{(n)} :=(D_1,\ldots,D_n)$, the observed sample of binary treatments. Remember that ${\cal X}_* = \{X_i \, : \, i =1,\ldots,n\} \subset {\cal X}$ is the set of actually observed covariate values in the sample, and $m=| {\cal X}_* |$ is its cardinality. As already mentioned above, in our asymptotic analysis we let $m \rightarrow \infty$ as $n \rightarrow \infty$. This implies that $ {\cal X}_*$ changes with the sample size (we can allow ${\cal X}$ to change with $n$ as well), but we do not make that explicit in our notation. For $x \in {\cal X} $ we define $$ {\cal N}(x) := \left\{ i \in \{1,\ldots,n\} \, \Big| \, X_i = x \right\} , $$ the set of observations $i$ for which the observed covariate value is equal to $x$.\footnote{ ${\cal N}(x)$ is empty for $x \notin {\cal X}_*$.} Let $n(x) := \left| {\cal N}(x) \right|$ be number of observations with $X_i=x$, and let

be the number of observations with $X_i=x$, and $D_i=0$ or $D_i=1$, respectively.

To construct our sample bounds, we furthermore require the researcher to choose a “bandwidth parameter” $Q \in \{1,2,3,\ldots,\infty\}$. If $\max_{ x \in {\cal X}_*} n( x)$ remains bounded as $n \rightarrow \infty$, then we can choose $Q=\infty$, which simplifies many of the expressions in this section, and the reader may think of this case as the baseline case which makes the connection to Section (ref) most obvious.

For each covariate value $ x \in {\cal X}_*$ we need to choose the order $q( x) \in \{1,2,3,\ldots\}$ of the bounds in Proposition (ref) that we want to implement. To implement bounds of a certain order $q( x)$ we require at least that many observations for that covariate value, that is, we need to choose $q( x) \leq n( x)$. Choosing the maximal value $q( x) = n( x)$ is optimal from the perspective of expected width of the bounds, but it is not advisable in general since it can lead to upper and lower bound estimates with very large variance. In our implementation of the bounds we therefore choose

that is, we choose the maximum order that satisfies both $q( x) \leq n({ x})$ and $q( x) \leq Q$. In practice, we recommend choosing $Q$ as small as $Q=3$ or $Q=4$, but the choice $Q=\infty$ gives some theoretical optimality properties for the expected width of the bounds (but usually at the cost of higher variance).

Having chosen the order $q( x)$ for each $ x \in {\cal X}_*$, we then construct sample weights $\widehat w_0( x)$, $\widehat w_1( x)$, $\widehat v( x) $, which are functions of the chosen order $q( x)$, the chosen reference propensity score $p_*( x)$, and the values $n( x)$, $n_0( x)$, $n_1( x)$ obtained from the sample, such that

where $q( x)$ is given in (ref), and the weight functions $w^{(q)}$ and $\widetilde w^{(q)}$ on the right-hand side were defined in (ref). Here and in the following, the dependence of those sample weights on $p_*( x)$ and $q( x)$ (and thereby on $Q$) is not made explicit, and the dependence on the sample $X^{(n)} $ and $D^{(n)}$ (through $n( x)$, $n_0( x)$, $n_1( x)$) is only indicated by the “hat”. Explicit formulas for $\widehat w_0( x)$, $\widehat w_1( x)$, $\widehat v( x) $ are provided in the next subsection.

The natural sample analogs of the bounds $ B^{(q)}_{d,a} $ and $C^{(q)}_a$ in the last section are then given by

for $a \in \mathbb{R}$. Notice that the arguments $d$, $a$ were subscripts in the population analysis, but for the sample version in this section we prefer to use the unit $i$ as the subscript instead. Also, the dependence on the order $q$ is not made explicit anymore here, but we always have the choice (ref) in mind.

In view of (ref), the expressions in (ref) are direct translations of the formulas in display (ref), where the weights were replaced by sample weights, and the remaining occurrences of the unknown $1-p(x)$ and $p(x)$ were replaced by their sample analogs $n_0( x)/n( x)$ and $n_1( x)/n( x)$, respectively. In all three expressions of display (ref) the maximum function in the denominator is only included to avoid a potentially zero denominator. However, $ n_0( X_i)=0$ implies $1-D_i=0$, and $n_1( X_i)=0$ implies $D_i=0$, that is, in all cases where the maximum function is required to avoid a zero denominator, the corresponding numerator is zero anyways. In particular, we could replace $\max\{1, \ldots \}$ by $\max\{c, \ldots \}$ for any constant $0 < c \leq 1$ without changing the sample bounds in (ref) at all.

The sample analogs of the expectations over $ B^{(q)}_{d,a} $ and $ C^{(q)}_a $ are then given by

and the final upper and lower sample bounds on the ATE read

Similarly, for the ATT, the lower and upper sample bounds on $\frac 1 n \sum_{i=1}^n \pi(X_i)$ are given by $ \overline C(a_{\max}) $ and $ \overline C(a_{\min}) $, respectively. To estimate the lower- and upper bounds on the ATT itself we still need to plug-in the sample analog of the denominator $\frac 1 n \sum_{i=1}^n p(X_i)$, which gives

In Section (ref) we show that the sample bounds just constructed are unbiased and consistent estimates (as $m \rightarrow \infty$) of the corresponding population bounds from the last section, and we will also use those sample bounds to construct asymptotically valid confidence intervals for ATE and ATT.

Construction of the sample weights \texorpdfstring{$\widehat w_0(x)$, $\widehat w_1(x)$, $\widehat v(x) $}{hat w0(x), hat w1(x), hat v(x)}

A key ingredient of the sample bounds just introduced are the sample weights that satisfy (ref), and which we want to define in this section. For ease of exposition we start with the simplest case $q(x)=n(x)$, which can be even or odd, and then generalize the formulas to the case $q(x) = \min\{Q,n(x)\}$ afterwards.

Case $q(x)=n(x)$ and $n(x)$ even

Let $q(x)=n(x)$, and assume that $n(x)$ is even. We consider $\widehat w_1(x)$ first. By setting

and using that, under Assumption (ref), we have $ \mathbb{E}\left(D_i \, \big| \, X^{(n)} \right) = p(X_i)$, we find that

where we used that the set $ {\cal N}(x)$ has $n(x) = q(x)$ elements, and the definition of the population weights in (ref). Thus, $\widehat w_1(x)$ satisfies the desired result in (ref). Finally, we can rewrite equation (ref) as

which from now on will serve as our definition of $\widehat w_1(x)$ in the current case. By analogous arguments one obtains, for the current case of $q(x)=n(x)$ and $n(x)$ even, that

and one can easily verify that those expressions satisfy (ref).

Case \texorpdfstring{$q(x)=n(x)$ and $n(x)$ odd}{hat q x=n(x) and n(x) odd}

For $q(x)=n(x)$ odd we have $ w^{(q(x))}(p,p_*)= 1- \left( 1- p \right) \left( \frac{p_* - p }{ p_*} \right)^{q(x)-1} $ according to (ref), and we then need to change (ref) to

Under Assumption (ref), it is again easy to see that the approximate unbiasedness condition for $\widehat w_1(x)$ in (ref) is satisfied here. In equation (ref), the sum over $i$ only gives a contribution for the $n(x) - n_{1}(x)$ instances where $D_i=0$, in which case there still are $n_1(x)$ units $j \in {\cal N}(x) \setminus \{i\}$ with $D_j=1$. We can therefore rewrite this equation as

which from now on is our definition of $\widehat w_1(x)$ for the case $q(x)=n(x)$ odd. By analogous arguments one obtains, for the current case, that

and one can again verify that those expressions satisfy (ref).

General case \texorpdfstring{$q(x) = \min\{Q,n(x)\}$}{hat q x = min{Q,n(x)}}

For $Q=\infty$ we have $q(x) = n(x)$, in which case all the required formulas for the sample weights are already provided in Subsections (ref) and (ref) above. The generalization to finite $Q$ discussed in the following is not conceptually difficult, but it requires some combinatorial arguments. Remember that we choose the order $q(x)$ of the bounds according to (ref). For even order $q(x) =q$, we generalize the formula for $\widehat w_1(x)$ in (ref) as follows:

where the sum is over all subsets ${\cal S}_q \subset {\cal N}(x)$ with $q$ elements. For odd order $q(x) = q$, we generalize the formula for $\widehat w_1(x)$ in (ref) to

where the sum is over all subsets ${\cal S}_{q-1,i} \subset {\cal N}(x) \setminus \{i\}$ with $q-1$ elements.

Under Assumption (ref), it is again straightforward to verify that those formulas for $\widehat w_1(x)$ guarantee that $ \mathbb{E}\left[ \widehat w_1(x) \, \Big| \, X^{(n)} \right] = w^{(q(x))} \big(p(x),p_*(x) \big) $. If $q(x) < n(x)$, then alternative choices for the sample weight $\widehat w_1(x) $ exist that have the same conditional expectation -- for example, instead of averaging over ${\cal S}_q$ and ${\cal S}_{q-1,i}$, one could randomly choose one subset of $q(x)$ observations out of the set ${\cal N}(x)$ and implement the formulas in Subsections (ref) and (ref) using only that subset of observations. To avoid that ambiguity in the definition of the sample weights we have chosen the formulas in (ref) and (ref) such that the binary treatment values $D_i$ of all units $i \in {\cal N}(x)$ enter exchangeably into $\widehat w_1(x)$, that is, the sample weights remain unchanged if we swap the data of any two observations in the same cluster ${\cal N}(x)$. This requirement also guarantees that it is possible to rewrite $\widehat w_1(x)$ such that the $D_i$ only enter through their summary statistics $n_1(x) = \sum_{i \in {\cal N}(x)} D_i $ and $n(x)$. Namely, one can rewrite (ref) and (ref) as

where $ \lfloor q(x) / 2 \rfloor$ is the integer part of $q(x) / 2$, and the combinatorial coefficients $ \omega_{k,n_1(x),n(x),Q} \in [0,1]$ are implicitly determined from (ref) and (ref), and one can show that

where $n_{1}=n_{1}(x)$ and $q=q(x)$ also depend on $x$. Appendix (ref) provides a derivation of this formula for $ \omega_{k,n_1(x),n(x),Q}$. Implementing $ \widehat w_1(x)$ via (ref) and (ref) is much faster than via (ref) and (ref), and can be done quickly also for relatively large values of $n(x)$ and $n_1(x)$.

Analogously we have

where the combinatorial coefficients $\omega_{k,n_0(x),n(x),Q} \in [0,1]$ are again those in (ref), only the argument $n_1(x)$ was changed to $n_0(x)$. The equations in (ref) and (ref) provide general definitions of the sample weights that satisfy (ref).

Discussion of the sample weights

We want to briefly discuss some properties of the sample weights, again mostly focusing on $\widehat w_1(x) $ for concreteness. If we choose $Q=\infty$, then the formula for $\widehat w_1(x) $ is given in (ref) for even $n(x)$, and in (ref) for odd $n(x)$. For $p_*(x)<\frac 1 2$ we have $\left| \frac{p_*(x)-1}{p_*(x)} \right|>1$, implying that the absolute value of $\widehat w_1(x) $ grows exponentially with $n_1(x)$. Analogously, for $Q=\infty$ and $p_*(x) > \frac 1 2$ the absolute values of the weights $\widehat w_0(x) $ and $\widehat v(x)$ grow exponentially with $n_0(x)$. Only for $p_*(x) = 1/2$ are all the sample weights bounded, independent of the realization of $n_0(x)$ and $n_1(x)$.

Thus, for $Q=\infty$ the weights can take very large negative or positive values, potentially resulting in sample bounds for ATE and ATT with very large variance. This is the main reason why we introduce the bandwidth parameter $Q$, which in practice we recommend to set relative small, say $Q =3$ or $Q =4$. Once we have chosen a finite value of $Q$, then our sample weights in (ref) and (ref) are all bounded, independent of the realization of $n_0(x)$ and $n_1(x)$ --- notice that the combinatorial coefficients $\omega_{k,n_{0/1}(x),n(x),Q}$ are all bounded between zero and one.

An interesting alternative way to guarantee that the weights $\widehat w_{0}(x) $ and $\widehat w_{1}(x) $ both remain bounded is to choose $Q= \infty$, but $p_*(x)=1/2$ for all $x \in {\cal X}_*$. That is not our leading recommendation, because in many applications one might prefer values of $p_*(x)$ different from $1/2$ to obtain better bounds. If the parameter of interest is ATT, then we can choose $Q= \infty$ and $\widehat v(x)$ will remain bounded as long as $p_*(x) \leq \frac 1 2$ for all $x \in {\cal X}_*$. This could indeed be an interesting option in applications on ATT estimation. Nevertheless, the variance of the bounds will usually be smaller when a finite value of $Q$ is chosen. Furthermore, as illustrated in the following concrete examples for $\widehat w_{1}(x) $, only for finite $Q$ do the sample weights converge to the population weights as $n(x) \rightarrow \infty$.

Figure (ref) plots the weights $ \widehat w_1(x)$ for $Q=6$, $n(x) \in \{6,12\}$, and for three different values for the reference propensity score $p_*(x)$. The plot shows that as $n(x)$ becomes large the weights $ \widehat w_1(x)$ as a function of $\widehat p(x) = n_1(x)/n(x)$ converge to the population weights $w^{(q)}(p(x),p_*(x))$ as a function of $p(x)$. This, in particular, implies that $ \widehat w_1(x)$ becomes a smooth function of $n_1(x)$ for large values of $n(x)$. However, for small $n(x)=Q=6$ the weights $ \widehat w_1(x)$ heavily fluctuate as a function of $n_1(x)$. Furthermore, for $p_*(x)<0.5$ the weights $ \widehat w_1(x)$ can take on very small and very large values (notice the different scale of the plot for $p_*(x) = 0.4$), but for $p_*(x)\geq 0.5$ the weights remain within the bounded interval $[0,2]$.

Asymptotically valid confidence intervals

Remember that $m= | {\cal X}_* |$ is the number of different covariate values in our sample. Our treatment effect bounds are then based on weight functions that combine the observed treatment status $D_i$ for observations $i \in {\cal N}(x)$ of the same covariate value $x \in {\cal X}_*$ in a non-linear way. However, if we condition on realization of the covariates $X^{(n)}$, then across different covariate values the bounds are just averages of independent observations. Given that the bounds have this structure, it is useful to think of $m$ as our effective sample size, and of each $x \in {\cal X}_*$ as labelling one effective observation. It is therefore convenient to rewrite the sample bounds in (ref) not as cross-sectional averages over $i \in \{1,\ldots,n\}$, but as sample averages over $x \in {\cal X}_*$. For that purpose, for $d \in \{0,1\}$ and $a \in \mathbb{R}$, we define\footnote{ We are slightly abusing notation here, for example, $\widehat B_{x}(d,a)$ for $x=1$ (assuming $1 \in {\cal X}_*$) is not the same as $ \widehat B_i(d,a)$ for $i=1$. However, it will always be clear from the subscript letter which object is meant. }

which allows us to rewrite the sample bounds in (ref) as

Using the definitions of $\widehat B_i(d,a) $ and $\widehat C_i(a) $ in (ref) we furthermore have

where

Notice that for $n_d(x)=0$ we have $ \widehat w_d(x) =0$, and for $n_0(x)=0$ we have $ \widehat v(x) = 0$. Therefore, $\overline Y_{x}(d)$ only enters into the bounds in (ref) when $ n_d(x) >0$. In that case, $ \overline Y_{x}(d)$ is simply the average of the $n_d(x)$ observed outcomes $Y_i$ for which $X_i=x$ and $D_i=d$. However, for our theoretical discussion it is useful to also define $ \overline Y_{x}(d)$ for the case $n_d(x) =0$, because with that definition we have that, under Assumption (ref),

Equation (ref) states that $ \overline Y_{x}(d)$ is mean-independent of $D^{(n)}$ and $X^{(n)} $. The properties of $\widehat w_{0/1}(x)$ and $\widehat v(x) $ in display (ref) together with (ref) guarantee that the expected values of $ \widehat B_{x}(d,a)$ and $ \widehat C_{x}(a) $ are equal to the expectations of the population bounds $B^{(q)}_{0,a} $ and $C^{(q)}_a $ in Section (ref).

Next, we want to show consistency of those sample bounds and use them to construct confidence intervals. For that purpose, it is convenient to define

which are the four parameters of interest that we focus on in this paper after conditioning on the realization of all the covariates $X^{(n)}= (X_1,\ldots,X_n)$. For each of those parameters we have already introduced upper and lower bound estimates in (ref), (ref), (ref). For $ \theta^{(0)}$ and $\theta^{(1)}$ we now denote those bounds by

Using the above definitions we have, for $r \in \{0,1,{\rm ATE}\}$,

where

for $d \in \{0,1\}$. Our results on the “population bounds” in the last section together with (ref), (ref) and (ref) guarantee that

When comparing the last line with the definition of the actual sample bounds $ \overline L^{\rm (ATT)} $ and $\overline U^{\rm (ATT)}$ in (ref) we notice that we need to account for the randomness of the denominator term $\frac 1 n \sum_{i=1}^n D_i$ as well when constructing confidence intervals, and we therefore write those bounds as (see appendix (ref) for details)

where

Under Assumption (ref)(iii) we have that $(D_i,Y_i(0),Y_i(1))$ is independent across $i$, conditional on $X^{(n)}= (X_1,\ldots,X_n)$. This, in particular, guarantees that $ L^{(r)}_{x} $ and $U^{(r)}_{x}$, for $r \in \{0,1,{\rm ATE},\allowbreak {\rm ATT}\}$, are independent across $x \in {\cal X}_*$, conditional $X^{(n)}$. This independence is crucially used for the asymptotic convergence results stated in the following theorem. For that reason, all the stochastic statements in the theorem are conditional on $X^{(n)}$. Notice also that we have in mind a triangular array in our {asymptotic theory},\todo[size=tiny]{change from `asymptotic' to `asymptotic theory'} where the support of the regressors may change as the sample size increases.

Here, the assumptions that $Q$ is fixed and that $p_*(x)$ is bounded away from zero and one guarantee that our sample weights $\widehat w_d(x) $ and $\widehat v(x)$, and therefore also $\widehat B_{x}(d,a)$ and $ \widehat C_{x}(a)$ defined in (ref), are uniformly bounded. However, the averages $\frac 1 m \sum_{x \in {\cal X}_*} L^{(r)}_{x}$ and $\frac 1 m \sum_{x \in {\cal X}_*} U^{(r)}_{x}$ that give our bounds are over the $ L^{(d)}_{x} $ and $U^{(d)}_{x}$ defined in (ref), and those feature the additional factors $\frac{m \, n(x)} n \in [0,\infty)$. Thus, covariate values that appear often in the sample get more weight than covariate values that appear less often. Notice that $\frac{n} m = \frac 1 m \sum_{x \in {\cal X}_*} n(x)$ is the average number of observations for a given covariate value, that is, the factor $\frac{m \, n(x)} n $ simply rescales the $n(x)$ such that they average to one: $ \frac 1 m \sum_{x \in {\cal X}_*} \frac{m \, n(x)} n =1$. The assumption $ \frac 1 {m} \sum_{x \in {\cal X}_*} \left( \frac{m \, n(x)} n \right)^4 = O_P(1)$ requires that the fourth moment of $\frac{m \, n(x)} n$ remains bounded asymptotically, that is, it demands that the $n(x)$ are not distributed too heterogeneously across covariates. For example, if $X_i$ is uniformly distributed over ${\cal X}_*$, then each $n(x)$ has a Binomial distribution with parameters $n$ and $1/m$, and it is easy to verify that the assumption is satisfied. More generally, the assumption holds as long as the probabilities $P(X_i =x)$ are not too heterogeneous across $x$.

Notice that for ${M_{x}^{(r)}} \in \{ L^{(r)}_{x}, U^{(r)}_{x} \} $ we have $$ {\rm Var} \left( \frac 1 {\sqrt{m}} \sum_{x \in {\cal X}_*} M^{(r)}_{x} \, \Big| \, X^{(n)} \right) = \frac 1 m \sum_{x \in {\cal X}_*} {\rm Var} \left( M^{(r)}_{x} \, \Big| \, X^{(n)} \right), $$ that is, our assumption $ \left[\frac 1 m \sum_{x \in {\cal X}_*} {\rm Var} \left( M^{(r)}_{x} \, \Big| \, X^{(n)} \right) \right]^{-1} \, = o_P( m^{1/3} )$ simply demands that the variance of $\frac 1 {\sqrt{m}} \sum_{x \in {\cal X}_*} M^{(r)}_{x} $ is not too small. Here, $ \frac 1 {\sqrt{m}}$ is a natural rescaling, because the $M^{(r)}_{x}$ have zero mean and are independent across $x$, conditional on $X^{(n)} $. However, the $M^{(r)}_{x} $ may contribute heterogeneously to the variance because of the factors $\frac{m \, n(x)} n$ in their definition, and also because of the weights $\widehat w_d(x) $ and $\widehat v(x)$. The assumption therefore allows for the possibility that $ \frac 1 m \sum_{x \in {\cal X}_*} {\rm Var} \left( M^{(r)}_{x} \, \Big| \, X^{(n)} \right)$ converges to zero as $n,m \rightarrow \infty$, but not too fast.

Using (ref) and Theorem (ref) we obtain the following asymptotically valid confidence interval for $\theta^{(r)} $ of confidence level $(1-\alpha) \in (0,1)$,\footnote{ Here, we use the convention $[a,b]=\emptyset$ if $a>b$. }

where $\widehat \sigma^{(r)}_L := \sqrt{ {\rm SVar} \left( L^{(r)}_{x} \right) }$, $\widehat \sigma^{(r)}_U := \sqrt{ {\rm SVar} \left( U^{(r)}_{x} \right) }$. The following corollary states that $ {\rm CI}^{(r)}_{\rm basic}$ contains $\theta^{(r)} $ with probability at least $1-\alpha$ in large samples.

Thus, those confidence intervals $ {\rm CI}^{(r)}_{\rm basic}$ are asymptotically valid, but they may be conservative for three reasons: (i) the true $\theta^{(r)}$ may be an interior point of the expected bounds, implying 100% coverage in large samples; (ii) we are using an upper bound estimate for the variance of the upper and lower bounds when constructing the confidence interval, and (iii) we are using Bonferroni inequalities when dividing the statistical problem into one-sided confidence interval constructions for the upper and lower bounds --- notice the $\alpha/2$ in both the upper and lower bounds in (ref).\footnote{ One could improve on those $\alpha/2$ critical values by adapting the methods in imbens2004confidence and stoye2009more to our case. However, we want to keep the confidence interval construction simple here, and there is also the more important issue that $ {\rm CI}^{(r)}_{\rm basic}$ can be empty in our case, which we address using stoye2020. }

Here, the issues (i) and (iii) are very typical for bound estimation, and (ii) is impossible to fully overcome in our setting, unless $n_d(x)$ are sufficiently large for all $d$ and $x$. For example, if $n_d(x)=1$, then only a single outcome $Y_i$ is observed for which we have $D_i=d$ and $X_i=x$, implying that unbiased estimation of the variance of that outcome is impossible, but since $Y_i$ enters into $ \overline L^{(r)} $ and $ \overline U^{(r)} $ we can in general not expect to estimate the variances of these bounds consistently.\footnote{ Another problem is that the true propensity scores $p(x)$ are unknown, rendering the distribution of the sample weights $ \widehat w_d(x)$ also unknown. }

We therefore believe that one needs to be content with conservative confidence intervals in our setting, and that our construction so far has the advantage of being relatively simple and robust. However, a potentially more severe problem in practice is that the confidence interval $ {\rm CI}_{\rm basic}$ may be empty, that is, the lower bound may be larger than the upper bound, because nothing in our construction guarantees that $\overline L^{(r)} $ cannot be larger than $ \overline U^{(r)} $ in finite samples. While our theory guarantees that this problem cannot occur asymptotically, it is still undesirable to have a potentially empty confidence interval in applications.

We therefore use the method in stoye2020 to obtain a valid confidence interval that is never empty. The general version of that method requires knowing the correlation $\rho$ between $ \overline L^{(r)} $ and $ \overline U^{(r)} $, which we cannot estimate consistently in our setting (for the same reasons for which we can only obtain upper bounds on the variances of $ \overline L^{(r)} $ and $ \overline U^{(r)} $). We therefore apply stoye2020's method with $\rho=1$, which corresponds to the worst case: Let

and

and define the final confidence interval to be reported for $\theta$ as the union of ${\rm CI}_{\rm basic}$ and $ {\rm CI}_*$, that is, $$ {\rm CI}^{(r)}_{\theta} := {\rm CI}^{(r)}_{\rm basic} \, \cup \, {\rm CI}^{(r)}_* . $$ Then, by construction, ${\rm CI}_{\theta}$ is never empty, because ${\rm CI}_*$ is never empty, and Corollary (ref) implies that

We refer to stoye2020 for a further justification of this specific confidence interval construction. We have thus shown how to construct valid non-empty confidence intervals for all of those objects of interest.

Notice also that for the constructions of confidence intervals here we have assumed that $p_*(x)$ is non-random. If $p_*(x)$ is estimated, then the randomness of $p_*(x)$ should be accounted for when constructing those confidence intervals, either via an application of the delta method, or via a bootstrap procedure.

Clustering the covariate observations

Unconfoundedness only places restrictions on the observed data if there are at least some repeated covariate values. When each covariate vector is unique, we essentially revert to first-order (Manski) bounds, which remain valid without additional assumptions but do not leverage unconfoundedness to tighten those bounds.

A common way to exploit unconfoundedness when covariates are nearly unique is to coarsen them by binning. For example, one might group ages into years rather than days. This process discards some information but remains transparent. Alternatively, one can adopt automated methods such as clustering or nearest-neighbor matching. For clustering, we partition units into groups of similar \(X_i\) values and label each unit by its cluster identity \(\overline{X}_i\). The main steps of our approach remain unchanged, except that we substitute \(\overline{X}_i\) for \(X_i\).

Although clustering is straightforward in practice, it introduces dependence among the labeled \(\overline{X}_i\) because the clustering procedure relies on the entire sample. A fully rigorous treatment would require additional smoothness assumptions or a formal model for the underlying clustering structure. We leave these issues for future work, noting that established methods (e.g., sample-splitting or matching) can mitigate some of the complications. In summary, binning or clustering offers a practical way to address rare or unique covariates when applying our bounds, but more theoretical investigation is warranted.

The specific clustering procedure we employ in our simulations and empirical application proceeds as follows. First, we studentize each observed covariate, then use the Euclidean distance \(\| X_i - X_j \|\) to measure the closeness of observations \(i\) and \(j\). With this distance measure, we apply hierarchical, agglomerative clustering with complete linkage to the observed covariate sample \(\bigl(X_1,\ldots,X_n\bigr)\). We refer to, e.g., kaufman2009finding and everitt2011cluster for an introduction to hierarchical clustering methods, and to mullner2013fastcluster and R:cluster for software implementations. Hierarchical, agglomerative clustering begins with singleton clusters and iteratively merges pairs of clusters until all observations lie in a single cluster. A user-selected number of clusters, \(m\), can be obtained by “cutting” the resulting tree. Different forms of hierarchical clustering differ in how they measure inter-cluster distance. Complete linkage uses the maximum distance between any two points, one in each cluster, which tends to produce relatively compact clusters everitt2011cluster.

The only tuning parameter in this clustering procedure is the number of clusters \(m\in \{1,2,\dots,n\}\). This plays the same role as the number of unique covariate values in our earlier analysis. In practice, we recommend choosing

where \(L\) is a constant (e.g., \(L=10\)) and \(\lceil \cdot \rceil\) denotes the ceiling function. This ad hoc rule aims for around \(L\) observations per cluster on average, and letting \(L\) remain fixed ensures \(m\to \infty\) as \(n\to \infty\), consistent with the large-\(m\) asymptotic theory in Section (ref).

Given this partition \(\{1,\dots,n\} = \mathcal{N}_1 \cup \cdots \cup \mathcal{N}_m\), we label each cluster by its average covariate value. Concretely, for each \(g\in \{1,\dots,m\}\) and \(i\in \mathcal{N}_g\), \[ \overline{X}_i \;:=\; \frac{1}{\lvert \mathcal{N}_g\rvert} \sum_{j \in \mathcal{N}_g} X_j, \] and let \(\overline{\mathcal{X}}=\{\overline{X}_i : i=1,\dots,n\}\) be the set of these cluster averages. By construction, \(\lvert \overline{\mathcal{X}} \rvert=m\) and each \(\overline{X}_i\) uniquely identifies the cluster that observation \(i\) belongs to. For \(\overline{x}\in \overline{\mathcal{X}}\), define the corresponding cluster as \[ \mathcal{N}(\overline{x}) \;:=\; \bigl\{ i \in \{1,\dots,n\}\,\big\vert\, \overline{X}_i = \overline{x} \bigr\}, \] and let \(n(\overline{x})=\lvert \mathcal{N}(\overline{x})\rvert\) be its number of observations. Notice that if no observation is “close” to \(i\) in terms of covariates, then \(i\) may end up in a singleton cluster, i.e.\ \(n(\overline{x})=1\).

Once the partition is obtained and labeled, the construction of our sample bounds proceeds exactly as in Section (ref), except that we replace each \(X_i\) by \(\overline{X}_i\), each set \(\mathcal{X}_*\) by \(\overline{\mathcal{X}}\), and so on.

Monte Carlo Experiments

In this section, we report results from Monte Carlo experiments. The scalar covariate $X_i$ is randomly generated from $\text{Unif}[-3,3]$. The binary treatment variable $D_i$ is then obtained from the following two models:

To generate the outcome variable, define

where $V_{di} \sim N(0,1)$, $d \in \{0,1\}$, and $(V_{1i}, V_{0i})$ are independent of $(D_i, X_i)$. Finally, the observed outcome variable is generated by $$ Y_i = D_i \mathbbm{1}\{ Y_{1i}^\ast > 0 \} + (1-D_i) \mathbbm{1}\{ Y_{0i}^\ast > 0 \}. $$ To study the effect of misspecification and the lack of overlap, we take the reference propensity score $p_*(x)=0.5$. That is, under DGP A, the model is correctly specified and the overlap condition is satisfied; whereas, under DGP B, the model is misspecified and the overlap condition is not satisfied. When $X_i \leq - 2$, $p(X_i) = 1$ in DGP B. By simulation design, $a_{\min} = 0$ and $a_{\max} = 1$. In the Monte Carlo experiments, we focus on the ATT.

Recall from (ref) that in our definition, the true ATT is given by

where $p(x)=\mathbb{E}(D_i\mid X_i=x)$ and $\tau(x)=\mathbb{E}[Y_i(1)-Y_i(0)\mid X_i=x]$. To obtain the closed-form expression for $\tau(x)$ in the Monte Carlo design, let $\Phi(\cdot)$ denote the standard normal CDF. For any $x$ and $d \in \{0,1\}$,

Hence,

Then, the finite-$n$ true ATT in our setting is obtained by combining (ref) and (ref) into (ref).

Define $\widehat{p} = n^{-1} \sum_{i=1}^n D_i$. We consider the following point estimators:

Here, $\widehat{\rm ATT}_{\mathrm{Oracle}}$ is an infeasible oracle estimator of ATT, whereas $\widehat{\rm ATT}_{\mathrm{RPS}}$ is an estimator using the known (parametric) propensity score $p_*(\cdot) = 0.5$ (for both DGPs). We also consider the nearest neighbor estimator of ATT:

where $\widehat Y_{0 i}$ is the nearest neighbor estimator of $\mathbb{E}[Y \mid X=X_i, D=0]$. For the bounds, $[\mathrm{LB}(Q),\mathrm{UB}(Q)]$ denotes the $Q$th-order bounds constructed with the reference propensity score $p_*(x)=0.5$, while $[\mathrm{LBc}(Q),\mathrm{UBc}(Q)]$ denotes the alternative bounds constructed with $p_*(x)=0$.\footnote{Recall Remark (ref) regarding the use of $p_*=1$ for estimating $\mathbb{E}_P[ Y_i(1) ]$. Since our object of interest is the ATT, the relevant quantity is instead $\mathbb{E}_P[ Y_i(0) \mid D_i=1 ]$. In this case, one must consider $p_*=0$ for $\mathbb{E}_P[ Y_i(0) \mid D_i=1 ]$ if one insists on using estimators that are conditionally unbiased. See also the general weighting schemes in (ref) and (ref).} We report results for $Q=1,2,3$ throughout. The number $m$ of clusters is chosen according to (ref) with $L=10$. The sample size is $n=1{,}000$, and each design is based on $1{,}000$ Monte Carlo replications.

Table (ref) summarizes the Monte Carlo results.\footnote{In Online Appendix (ref), we report additional Monte Carlo experiments that further investigate the finite-sample performance of the proposed inference methods.} All reported statistics are computed after subtracting the true ATT from each estimator or bound. That is, for each method, the mean, median, and standard deviation summarize the distribution of the centered quantity $\widehat{\rm ATT}-{\rm ATT}$ across Monte Carlo replications. Accordingly, values close to zero indicate good finite-sample performance.

In DGP A, where the overlap condition holds and the reference propensity score is correctly specified, the oracle, RPS, NN, and higher-order bound estimators with $Q=2,3$ all have means and medians close to zero. In contrast, the first-order Manski bounds $\mathrm{LB}(1)$ and $\mathrm{UB}(1)$ are wide and centered far from zero, reflecting the fact that they do not exploit the unconfoundedness assumption.\footnote{$\mathrm{LBc}(1)$ and $\mathrm{UBc}(1)$ coincide with $\mathrm{LB}(1)$ and $\mathrm{UB}(1)$ because the reference propensity score $p_*(x)$ does not play any role when $Q=1$.} The alternative bounds $\mathrm{LBc}(Q)$ and $\mathrm{UBc}(Q)$ with $Q=2,3$ are tighter than the Manski bounds but remain relatively wide because the reference propensity score $p_*(x)=0$ is far from the true propensity.

In DGP B, the overlap condition fails and the ATT is not point identified. The NN estimator performs poorly, exhibiting substantial dispersion and a median far from zero. The RPS estimator, which relies on a misspecified propensity score, is also unreliable: its mean lies outside the range implied by the higher-order bound estimators with $Q=2,3$. In contrast, the proposed bounds with $Q=2,3$ remain informative and correctly reflect partial identification. The lower bounds become tighter as $Q$ increases, while remaining below the oracle mean.

Overall, the results from DGPs A and B illustrate that the proposed bound approach does not require the overlap condition and can substantially improve upon parametric estimators when the propensity score is misspecified. While the conservative choice $p_*(x)=0$ provides a useful worst-case benchmark, it may not be as competitive as $p_*(x)=0.5$, which is the case in the current Monte Carlo designs. These findings suggest that higher-order bound estimators constructed with a reasonable interior reference propensity offer a useful compromise between point identification under strong ignorability and worst-case Manski-type bounds.

Empirical Examples

Effects of Right Heart Catheterization

In this section, we apply our methods to Connors1996's study of the efficacy of right heart catheterization (RHC), which is a diagnostic procedure for directly measuring cardiac function in critically ill patients. This dataset has been subsequently used in the context of limited overlap by crump2009dealing, Rothe:2017, Li-et-al:2018, and Ma_Sasaki_Wang_2024 among others. The dataset is publicly available on the Vanderbilt Biostatistics website at \url{https://hbiostat.org/data/}.

In this example, the dependent variable is 1 if a patient survived after 30 days of admission, and 0 if a patient died within 30 days. The binary treatment variable is 1 if RHC was applied within 24 hours of admission, and 0 otherwise. The sample size was $n = 5735$, and 2184 patients were treated with RHC. There are a large number of covariates: hirano2001estimation constructed 72 variables from the dataset and the same number of covariates were considered in crump2009dealing and Li-et-al:2018, and Ma_Sasaki_Wang_2024, and 50 covariates were used in Rothe:2017. In our exercise, we constructed the same 72 covariates. For the purpose of illustrating our methodology, we assume that the unconfoundedness assumption holds in this example.\footnote{ bhattacharya2008treatment,BSV-2012 raise the concern that catheterized and noncatheterized patients may differ on unobserved dimensions and propose different bounds using a day of admission as an instrument for RHC.}

In this section, we focus on ATT. We first estimate ATT by the normalized inverse probability weighted estimator\footnote{See, e.g., equation (3) and discussions in busso2014new for details of the normalized inverse probability weighted ATT estimator.}:

where $W_i := \widehat{p}(X_i)/[1-\widehat{p}(X_i)]$ and $\widehat{p}(X_i)$ is the estimated propensity score for observation $i$ based on a logit model with all 72 covariates being added linearly as in the aforementioned papers. The estimator $\widehat{\text{ATT}}_{\text{PS}}$ requires that the assumed propensity score model be correctly specified and the overlap condition is satisfied. The resulting estimate is $\widehat{\text{ATT}}_{\text{PS}} = -0.0639$.\footnote{The unnormalized ATT estimate is $-0.0837$ using the same propensity scores.}

We now turn to our methods. We take the reference propensity score to be $ \widehat{p}_{\text{RPS}}(X_i) = n^{-1} \sum_{i=1}^n D_i$ for each observation $i$. That is, we assign the sample proportion of the treated to the reference propensity scores uniformly for all observations. Of course, this is likely to be misspecified; however, it has the advantage that $\widehat{p}_{\text{RPS}}(X_i)$ is never close to 0 or 1. The resulting inverse reference-propensity-score weighted ATT estimator is\footnote{When the sample proportion is used as the propensity score estimator, there is no difference between unnormalized and normalized versions of ATT estimates. In fact, it is simply the mean difference between treatment and control groups. }

None of the covariate values in the observed sample are identical among patients (that is, $n(X_i)=1$ for all observations here). We therefore implement the clustering method described in Section (ref). As recommended in Section (ref), we choose the number $m$ of clusters by (ref): $ m = \left\lceil \frac{n}{L} \right\rceil$ with $L = 5, 10, 20$. In addition, we consider $Q = 1,\ldots,4$.

Table (ref) reports the estimated ATT bounds for selected values of $L$ and $Q$ using two reference propensity scores. We first discuss Panel A, which corresponds to the reference propensity score $p_*(x)=\bar D$. When $Q=1$, our estimated bounds correspond to Manski bounds, which include zero and are wide, with interval lengths close to one for all values of $L$. Our bounds with $Q=1$ are different across $L$ because we apply hierarchical clustering before obtaining Manski bounds. With $Q=2$, the bounds shrink so that the estimated upper bound is zero for all cases of $L$; with $Q = 3$, they shrink even further so that the upper end point of the 95% confidence interval excludes zero. Among three different values of $L$, the case of $L=5$ gives the tightest confidence interval but in this case, the lower bound is larger than the upper bound, indicating that the estimates might be biased. In view of that, we take the bound estimates with $L=10$ as our preferred estimates [$-0.077, -0.039$] with the 95% confidence interval $[-0.117,-0.006]$. When $Q=4$, the lower bound estimates exceed the upper bound estimates with $L = 5, 10$. However, the estimates with $L = 20$ give an almost identical confidence interval to our preferred estimates. It seems that the pairs of $(L, Q) = (10, 3)$ and $(L, Q) = (20, 4)$ provide reasonable estimates.

Panel B reports the corresponding results obtained using the conservative reference propensity score $p_*(x)=0$. As expected, the resulting bounds are wider for all values of $L$ and $Q$, reflecting the additional conservatism of this choice.

The study of Connors1996 offered a conclusion that RHC could cause an increase in patient mortality. Based on our preferred estimates, we can exclude large beneficial effects with confidence. This conclusion is based solely on the unconfoundedness condition, but not on the overlap condition, nor on the correct specification of the logit model. Overall, our estimates seem to be consistent with the qualitative findings in Connors1996 under the maintained assumption that the unconfoundedness assumption holds.

Impacts of a Temporary Employment Program

In this section, we apply our methods to the well-known LaLonde-AER dataset, available on Rajeev Dehejia's web page at \url{http://users.nber.org/ rdehejia/nswdata2.html}. The LaLonde dataset comes from the National Supported Work Demonstration (NSW), a randomized controlled temporary employment program. The binary treatment variable indicates whether an individual is assigned to the treatment or control group. The original outcome variable (RE78) is post-experimental earnings in 1978; in our application, we define the outcome as whether an individual was employed in 1978, i.e., whether earnings in 1978 were positive (RE78 $> 0$). Because NSW is a randomized controlled trial (RCT), we first estimate the average treatment effect by computing simple mean differences, yielding a 95% confidence interval of $[0.013, 0.143]$. This interval is relatively wide but excludes zero.

Case Study Using the Dehejia-Wahba Subsample

DehejiaWahba-JASA and DehejiaWahba-RESTAT extract a subset of LaLonde's NSW experimental data that includes information on RE74 (earnings in 1974). If we assume that the Dehejia-Wahba sample preserves the initial randomization, we can impose that the reference propensity score is independent of covariates. However, this may not be the case, so our approach provides a robust method to check whether the Dehejia-Wahba sample can be viewed as a random sample from an RCT.

We define the reference propensity score to be the sample proportion of treatment in the Dehejia-Wahba sample. The covariates are: age in years, years of education, indicators for black, hispanic, married, and no degree, and earnings in 1974 as well as in 1975. If we treat this reference propensity score as a consistent estimator under preservation of randomization, the average treatment effect is again obtained by simple mean differences, producing a 95% confidence interval of $[0.026, 0.196]$, which is wide but excludes zero.

Next, we obtain our bounds on the average treatment effect (ATE). As before, it is necessary to choose $Q$ and $L$. Based on previous numerical results, we set $Q=3$ and $L=10$. Our bounds yield a 95% confidence interval of $[-0.014,0.184]$, which is similar to the interval $[0.02,0.20]$ obtained under the assumption that the Dehejia-Wahba sample is a random sample from NSW. This result suggests two points: first, there is no evidence that the random sampling assumption is violated in the Dehejia-Wahba sample; and second, our inference method does not substantially widen the confidence interval to achieve robustness, although the null effect is now included. Furthermore, our bounds remain similar if we change $Q$ to 2 or 4 or $L$ to 5 or 20, indicating that our findings are robust to the choice of tuning parameters. See details in Panel A of Table (ref).

Case Study Using the Population Survey of Income Dynamics Control Group

The Dehejia-Wahba sample can be viewed as a scenario where the propensity score is known and satisfies the overlap condition. We now turn to a different scenario where it is likely that the propensity score is unknown and may not satisfy the overlap condition. Specifically, we use one of the non-experimental comparison groups constructed by LaLonde from the Population Survey of Income Dynamics, the PSID2 controls. As in the previous subsection, we estimate the reference propensity score using the sample proportion and then obtain our bound estimates. The resulting confidence interval is $[-0.346, 0.310]$ with $Q=3$ and $L=10$, which is much larger than the interval $[-0.01,0.19]$ obtained with the Dehejia-Wahba sample. Note that the sample proportion is unlikely to be correctly specified in the NSW-treated/PSID2-control sample. Therefore, our inference method appears to produce a wider confidence interval to remain robust against possible misspecification of the propensity scores and/or a lack of overlap. As a benchmark, we also compute the Manski bounds by setting $Q=1$. The resulting confidence interval for the Manski bounds is $[-0.668, 0.770]$ with $Q=1$ and $L=10$, which is even larger. Recall that the Manski bounds do not impose the unconfoundedness assumption and do not rely on any pooling information (so the specification of the reference propensity score does not matter). See Panel B of Table (ref) for other values of $Q$ and $L$. We conclude this section by noting that our empirical findings are broadly consistent with those in Ma_Sasaki_Wang_2024. In particular, our results indicate that limited overlap is a salient concern in the NSW-treated/PSID-control sample, as evidenced by substantially wider bounds, while comparable concerns do not appear to arise in settings such as Connors1996’s study, as discussed in the previous section.