Model-Agnostic Covariate-Assisted Inference on Partially Identified Causal Effects

Model-Agnostic Covariate-Assisted Inference on Partially Identified Causal Effects

\begingroup \footnote{Authors are ordered alphabetically. We thank Alberto Abadie, Donald Andrews, Andres Aradillas-Lopez, P. M. Aronow, Jushan Bai, Stephen Bates, Stéphane Bonhomme, Yong Cai, Denis Chetverikov, Guilherme Duarte, Bulat Garafov, Isaac Gibbs, Patrik Guggenberger, Kevin Guo, Jinyong Hahn, Marc Henry, Keisuke Hirano, Guido Imbens, Sung Jae Jun, Vishal Kamat, Samir Khan, Yuichi Kitamura, Dean Knox, Sokbae Lee, Zhipeng Liao, Elena Manresa, Rosa Matzkin, Serena Ng, Joris Pinkse, Kirill Ponomarev, Guillaume Pouliot, Dominik Rothenh{\"a}usler, Fredrik S\"{a}vje, Vira Semenova, Azeem Shaikh, Shuyang Shen, Ruoyao Shi, Jann Spiess, Max Tabord-Meehan, Alexander Torgovitsky, Edward Vytlacil, and Martin Wainwright for helpful discussion. L.L. is grateful for the support of National Science Foundation grant DMS-2338464. A.S. was partially supported by the Two Sigma Graduate Fellowship Fund, the Citadel GQS PhD Fellowship, and a Graduate Research Fellowship from the National Science Foundation.} \addtocounter{footnote}{-1} \endgroup

Introduction

Motivation and problem statement

Many parameters of interest in econometrics and causal inference are only partially identifiable manski2003partial, tamer2010partial, molinari2020microeconometrics. Even in randomized experiments, we cannot observe the joint law of the potential outcomes $(Y_i(1), Y_i(0))$ since we observe at most one outcome per subject; thus, the law of the individual treatment effect $Y_i(1) - Y_i(0)$ is unidentifiable. However, most causal parameters of interest can be bounded using the marginal laws of $Y_i(1)$ and $Y_i(0)$. Furthermore, incorporating information from covariates $X_i \in \R^p$ can substantially reduce the width of the partially identified set.

However, partial identification bounds involving covariates can depend delicately on the relationship between the outcome and the covariates, making inference challenging. For illustration, we now give three motivating examples, although we will state a general problem formulation in Section (ref). As notation, assume that we observe $n$ i.i.d. observations $\{(X_i, W_i, Y_i)\}_{i=1}^n$ for covariates $X_i \in \mcX$, a binary treatment $W_i \in \{0,1\}$ and an outcome $Y_i \in \mcY$ with potential outcomes $Y_i(1), Y_i(0)$. This paper focuses on randomized experiments (see Assumption (ref)) where the marginal laws of $(Y_i(1), X_i)$ and $(Y_i(0), X_i)$ are identified. Thus, we say that a parameter is identified if it depends only on these marginal laws.

Given a partially identified parameter $\theta$, this paper aims to estimate sharp bounds $[\theta_L, \theta_U]$ which incorporate information from covariates. This problem is challenging because the bounds typically depend delicately on the conditional law of $Y \mid X, W$, as exemplified by Equations ((ref))-((ref)). Thus, most existing approaches to estimate $\theta_L, \theta_U$ make assumptions allowing uniformly consistent estimation of such nuisance parameters (see Section (ref) for a review). This assumption is often implausible when $X_i$ is continuous or high-dimensional, unless the researcher is willing to impose further assumptions on the conditional distributions (e.g., a parametric model, smoothness, sparsity), which may not hold in applications.

Thus, in this work, we ask the question: can we convert a working estimate of the conditional law of $Y \mid X, W$ into inferential bounds on the sharp identified set $[\theta_L, \theta_U]$ which are (i) sharp when the working estimate is consistent and (ii) conservative but valid when the working estimate is arbitrarily inaccurate?

We end this section by noting that this question is motivated by the core philosophy of partial identification. Indeed, why not simply make enough assumptions so that the parameter $\theta$ is identified? In his seminal book, manski2003partial answers this question by formulating the law of decreasing credibility:

Our objective is to enhance credibility by removing any assumptions about the accuracy of the researcher's working model of nuisance parameters without sacrificing power when the researcher's model matches the ground truth.

Contribution and overview of results

Our work introduces a framework for inference on sharp, covariate-assisted partial identification bounds on causal parameters. If $P\opt$ denotes the true joint law of $(Y_i(1), Y_i(0), X_i) \iid P\opt$ for $i \in [n]$, we consider estimands of the form

for some known function $f : \mcY^2 \times \mcX \to \R$. Many estimands can be reduced to this case, including Examples (ref)-(ref), certain conditional expectations, quantiles of treatment effects, and more (see Section (ref)). We let $\theta_L \le \theta(P\opt) \le \theta_U$ denote the sharp (population) lower and upper partial identification bounds on $\theta(P\opt)$; these quantities are defined formally in Section (ref). Our method outputs estimates $\hat\theta_L, \hat\theta_U$ of the sharp bounds $\theta_L, \theta_U$ as well as lower and upper confidence bounds $\hat\theta_\mathrm{LCB}, \hat\theta_\mathrm{UCB}$. The main idea is to leverage duality theory for optimal transport problems (reviewed in Section (ref)) to convert any estimate $\hat P_{Y \mid X, W}$ of the conditional law of the outcome into robust partial identification bounds $\hat\theta_L, \hat\theta_U$. We emphasize that this method works automatically for any estimand defined above---in our software, the analyst can specify any function $f$ and does not need to do any additional calculations to obtain the results. This eliminates the need for a closed-form representation of $\theta_L, \theta_U$.

These “dual bounds" have a few appealing properties, listed below.

1. Uniform validity. Our method allows analysts to estimate the law of $Y \mid X,W$ using any statistical or machine learning technique, e.g., quantile regression, boosting, neural networks, etc. However, in randomized experiments with known propensity scores, the resulting confidence bounds are valid even if the estimate $\hat P_{Y \mid X, W}$ is arbitrarily inaccurate relative to the ground truth $P\opt_{Y \mid X, W}$. In this sense, our method is “model-agnostic": it can leverage models for power without relying on them for validity.

Formally, $\hat\theta_L$ and $\hat\theta_U$ are always conservatively biased in the sense that $\E[\hat\theta_L] \le \theta_L$ and $\E[\hat\theta_U] \ge \theta_U$. Furthermore, the confidence bounds have uniform asymptotic coverage without any assumptions on the accuracy of $\hat P_{Y \mid X, W}$ (see Theorem (ref)). Finally our method is also doubly robust in observational studies where the propensity scores are not known (see Theorem (ref)).

2. Tightness. If one can estimate the relevant nuisance parameters at $o(n^{-1/4})$ rates, our estimators $\hat\theta_L, \hat\theta_U$ are asymptotically unbiased and $\sqrt{n}$-consistent for the sharp bounds $\theta_L, \theta_U$.

3. Easy model selection. A major question in empirical applications is (i) how to select the subset of the covariates used in the analysis and (ii) how to estimate the outcome model $Y_i \mid X_i, W_i$. Our method permits the analyst to use either nested cross-validation and/or the multiplier bootstrap chernozhukov2013multbootstrap to select the tightest bound based on different models or subsets of the covariates.

4. Computational efficiency. To compute our bounds, we propose an algorithm that is computationally efficient even when $X_i$ is high-dimensional and $Y_i$ is continuous. The python package dualbounds implements this algorithm: \url{https://dualbounds.readthedocs.io/en/latest/}.

It is noteworthy that our method achieves uniform validity and tightness simultaneously. If only the former is required, one can simply throw away all covariates and stick with covariate-independent bounds, which are by definition not tight. A common remedy is to apply a coarse stratification on a few discrete variables or to bin covariates in a data-driven fashion. However, unless the covariates are jointly discrete with a relatively small support, the former strategy could result in considerable efficiency loss and the latter is challenging to implement even with a moderate number of covariates, since one must balance the trade-off between increasing the number of bins (to improve efficiency) and ensuring there are enough observations per bin (which is necessary for inference). On the other hand, most provably tight inferential procedures crucially rely on (certain aspects of) the conditional distributions being consistently estimated and hence it is unclear if uniform validity can be achieved semenova2023classification, levis2023covariate.

Figure (ref) illustrates our contributions in a simple numerical experiment where we estimate lower Lee bounds as in Example (ref). We fit an outcome model estimate $\hat P_{Y \mid X,W}$ assuming $Y_i(k) \mid X_i$ follows a homoskedastic Gaussian linear model for $k \in \{0,1\}$. A naive estimator of $\theta_L$, which simply plugs in the estimated outcome model to Equation ((ref)), performs well when the model is well-specified. However, if the errors are made heteroskedastic, this naive “plug-in" estimator can become conservatively or anticonservatively biased (depending on the form of heteroskedasticity). In contrast, our dual bounds wrap around exactly the same estimator of the outcome model and provide provable validity under arbitrary misspecification. See Section (ref) for precise simulation details and an analogous plot showing coverage.

Related literature

Partial identification has a long history in econometrics and causal inference, and a great deal of work has been done to characterize and estimate sharp bounds in various settings manski1990treateffects, manski1997monotone, balke1997bounds, heckman1997, manski2002intervaldata, imbensmanski2004, firpo2008, molinari2008, molinari2008partial, beresteanu2008asymptotic, lee2009training, stoye2009more, fan2010, chiburis2010semiparametric, romano2010inference, beresteanu2011sharp, fan2012confidence, tetenov2012positivetreatment, andrewsshi2013, aronow2014, fan2017partial, firpo2019, kaido2019confidence, kline2021moment, russell2021sharp, jun2023persuasion, oberreynolds2023estimating,fava2024predicting, byun24a; see manski2003partial, tamer2010partial, molinari2020microeconometrics for a review. When covariates are available, the bounds can be improved by conditioning on the covariates and aggregating covariate-specific sharp bounds chernozhukov2007econometrica, chandrasekhar2012, chernozhukov2013intersection, semenova2020cate, semenova2021generalized, semenova2023classification, lee2023partial, levis2023covariate. However, unless the covariates are discrete with a few values, these methods generally either (a) make assumptions that allow the conditional distributions of the potential outcomes to be consistently estimated at semiparametric rates or (b) have to discretize covariates in a non-disciplined way at the cost of efficiency loss. In contrast, our method can handle any type of covariates without making any assumptions that enable consistent estimates of the conditional distributions.

Our key technical tool is the theory of duality in optimization. This tool is of course not new, although we use it in a novel way. In particular, many existing works use duality theory as part of an inference strategy, for example in analysis of certain linear programming problems (e.g., hsieh2022lp, andrews2023, fang2023lp) and in sensitivity analysis dornguo2021, dorn2022doublyvaliddoublysharp. Most recently, semenova2023classification independently developed a dual-based estimator for a class of intersection bounds chernozhukov2013intersection. However, they require consistent estimates of the conditional distributions at semiparametric rates uniformly over the covariate space. Moreover, to ensure a key margin condition in their proof, they can only consider intersection bounds over finite sets, which in our context requires the potential outcomes to be discrete.

Core Methodology

To aid comprehension, we mostly defer measure-theoretic details to Appendix (ref). We defer computational details to Section (ref). For brevity, we focus on the sharp lower bound $\theta_L$, but the same method can be used to estimate upper bounds by simply multiplying $\theta(P\opt)$ by negative one.

Assumptions and background on Kantorovich duality

We assume the setting of a randomized experiment, although Section (ref) relaxes the assumption that the propensity scores are known.

We also allow the analyst to optionally specify additional assumptions about $P\opt$, the joint law of $(Y(1), Y(0), X)$, via conditional moment inequalities, as defined below.

For example, $\mcP$ is the unrestricted set of all joint distributions over $\mcY^2 \times \mcX$ if $\mcW_x$ is empty for each $x \in \mcX$, in which case Assumption (ref) is always satisfied. On the other hand, in the setting of Lee bounds (Example (ref)) with compound potential outcomes $(Y_i(0), S_i(0)), (Y_i(1), S_i(1))$, the monotonicity assumption in lee2009training can be enforced by setting $\mcW_x$ to contain the single function $w((y_0, s_0), (y_1, s_0)) = \I(s_0 > s_1)$, which ensures $S_i(0) \le S_i(1)$ a.s. The conditional monotonicity assumption of semenova2021generalized is also a special case of Assumption (ref).

Given $\mcP$, the lower bound $\theta_L$ is the minimum value of $\theta(P)$ for all $P \in \mcP$ which are consistent with the true marginal distributions $P_{Y(1), X}\opt$ and $P_{Y(0),X}\opt$:

Now, we introduce the dual to this optimization problem. We refer to a collection of functions $\nu_{0,x}, \nu_{1,x} : \mcY \to \R$ indexed by $x \in \mcX$ as dual variables; we use the notation $\nu = (\nu_{0,x}, \nu_{1,x})_{x \in \mcX}$ to denote the collection of these functions. \footnote{Formally, $\nu : \mcY \times \mcX \to \R^2$ is the function defined by $\nu_{k,x}(y) = \nu_{k,x}(y)$ but to avoid confusion we mostly avoid using this notation.} Given dual variables $\nu$, the Kantorovich dual function is

Intuitively, $g(\nu)$ is an “average treatment effect" of the transformed potential outcomes $Y'(1) \defeq \nu_{1,X}(Y(1))$ and $Y'(0) \defeq - \nu_{0,X}(Y(0))$. Thus, $g(\nu)$ is easy to estimate for any fixed $\nu$.

We aim to use $g(\nu)$ as a lower bound on $\theta_L$. To ensure $g(\nu) \le \theta_L$ holds, we will enforce a collection of known constraints on $\nu$. In the simplest case where $\mcP$ is unrestricted, we require that $\nu_{0,x}(y_0) + \nu_{1,x}(y_1) \le f(y_0, y_1, x)$ for all $y_1, y_0, x \in \mcY^2 \times \mcX$. In the general case where $\mcW_x = \{w_{x,1}, \dots, w_{x,L}\}$ is nonempty, we can slightly loosen these constraints to take advantage of additional assumptions on $\mcP$. Namely, for $x \in \mcX$, we say that $\nu_{0,x}, \nu_{1,x}$ are conditionally valid at $x$ if there are a collection of nonnegative constants $\{\lambda_{x,\ell}\}_{\ell=1}^{L}$ such that the following holds:

and we let $\mcV_x \subset \{\mcY \to \R^2\}$ denote the set of all pairs of functions satisfying this condition. Finally, we say that the full set of dual variables $\nu = (\nu_{0,x}, \nu_{1,x})_{x \in \mcX}$ are fully valid or “dual-feasible" if $\nu_{0,x}, \nu_{1,x} \in \mcV_x$ are conditionally valid for every $x \in \mcX$, and we let $\mcV \subset \{\mcY \times \mcX \to \R^2\}$ denote the set of all valid dual variables.

Computational issues aside (see Section (ref)), we emphasize that $\mcV$ is a known set which does not depend on $P\opt$. Satisfying this known constraint ensures that weak duality holds, i.e., $g(\nu) \le \theta_L$. The theorem below states this formally; it also states a strong duality result and gives a useful characterization of the optimal dual variables $\nu\opt \in \argmax_{\nu \in \mcV} g(\nu)$. \footnote{As notation, $\nu\opt \in \argmax_{\nu \in \mcV} g(\nu)$ denotes any “optimal" dual variables; when the $\argmax$ is not unique, $\nu\opt$ represents an arbitrary choice of maximizer.} To ease readability, we defer technical regularity conditions regarding measurability and the proof of Theorem (ref) to Appendix (ref).

This theorem has two statistical implications. First, to estimate $\theta_L$, we need only (i) estimate $\nu\opt$ and (ii) estimate $g(\nu\opt)$. Second, $\nu\opt$ is only a functional of $\PCopt$ and does not depend on $P_X\opt$. We now use these insights to estimate $\theta_L$.

Inference via dual bounds

To motivate our method, recall that the dual function $g(\nu)$ is easy to estimate for any fixed choice of $\nu \in \mcV$ using an inverse probability weighting (IPW) estimator. Thus, if only we knew the value of $\nu\opt$, we could easily estimate $\theta_L = g(\nu\opt)$. The main idea is to use the first split of the data to estimate $\hat \nu \approx \nu\opt$, and the second split of the data to estimate $g(\hat \nu)$. Crucially, even if our first-stage estimate $\hat \nu$ is poor, our inference will be conservative but valid, since weak duality ensures $g(\hat \nu) \le \theta_L$. And as we will see in Section (ref), if $\hat \nu$ is close to $\nu\opt$, then our confidence interval will be tight.

We will see in Section (ref) that this procedure is uniformly valid (in randomized experiments) and that it provides an asymptotically exact and sharp lower confidence bound if we can estimate $\hatPC$ at semiparametric rates. The main drawbacks of this procedure are that it requires splitting the data and that Eq. ((ref)) assumes the propensity scores are known. In Section (ref), we overcome these drawbacks by employing cross-fitting and by plugging in estimates $\hat \pi$ of the propensity scores in observational data. Before presenting these additional results, however, we first give a few guidelines and examples of how to apply this procedure.

Model selection via the multiplier bootstrap

To compute dual bounds, one must estimate the dual variables $\nu\opt$---in practice, we recommend first estimating $\hatPC$ and then computing $\hat \nu$ as per Eq. ((ref)). However, there are many ways to estimate $\hatPC$. E.g., analysts may prefer to use only a subset of the covariates to predict $Y$, but even after observing $\mcD_1$, it is not clear which subset of the covariates to choose. And even after making this decision, as discussed in Section (ref), there are still countless existing methods to estimate $\hatPC$. This raises the question: in practice, how should analysts choose between $K$ candidate estimates $\hat \nu^{(1)}, \dots, \hat \nu^{(K)} \in \mcV$ of the dual variables $\nu\opt$? Or more colloquially, how should we perform model selection?

One solution is to perform cross-validation within the first fold ($\mcD_1$) and pick the best-performing model. This approach is clearly valid since the final estimated dual variables $\hat \nu$ still depend only on $\mcD_1$, satisfying Definition (ref). In Section (ref), we recommend this approach for observational studies, where the validity of the final bounds may depend on the accuracy of the outcome model. However, in randomized experiments, we can improve upon this method.

In particular, let $\tilde{\theta}_L^{(k)} = g(\hat \nu^{(k)})$ denote the dual lower bound on $\theta_L$ implied by the estimate $\hat \nu^{(k)}$, for $k=1, \dots, K$. We will estimate $\max_{k\in [K]} \tilde{\theta}_L^{(k)}$, the tightest possible lower bound on $\theta_L$ based on $\{\hat \nu^{(k)}\}_{k=1}^K$, using the Gaussian multiplier bootstrap chernozhukov2013multbootstrap, as defined below.

The multiplier bootstrap is well-suited to this problem for two reasons. First, our bounds are valid no matter which model we select, i.e., $\tilde{\theta}^{(k)}_L \le \theta_L$ always holds. This may not be true in other problems---for example, when estimating regression coefficients, selecting different subsets of covariates may lead to anticonservative bias, but in our setting, any bias from misspecification is conservative. Second, after estimating $\{\hat \nu^{(k)}\}_{k=1}^K$, the dual bounds $\{\tilde{\theta}_L^{(k)}\}_{k=1}^K$ can be expressed as marginal moments and estimating them does not require (e.g.) any complicated M-estimation. As a result, in Section (ref), we conclude that the multiplier bootstrap quantile $\hat q_{1-\alpha}$ is consistent even if $K$ grows exponentially with a power of $n$ chernozhukov2018multibootstrap.

Guidelines on estimating the conditional distributions $\hatPC$

The first step in computing a dual bound $\hat\theta_\mathrm{LCB}$ is to estimate $\hatPC$, or equivalently, to estimate the conditional law of $Y_i \mid X_i, W_i$. An immense literature exists on this modeling problem koenker1978regression, chernozhukov2010quantile, chernozhukov2013inference, friedman2020contrast, and any choice will yield valid inferences. However, we make a few recommendations here.

To start, note that it is usually insufficient to model the conditional mean $\E_{P\opt}[Y_i \mid X_i, W_i]$, since the sharp lower bound $\theta_L$ may depend on the whole conditional law (e.g. Examples (ref)- (ref)). Instead, we can apply distributional regression is devoted to the task of estimating the law $Y_i \mid X_i, W_i$ (see kneib2023distreg for a review). One way to do this is to fit many quantile regressions. Another simple method is to assume a Gaussian linear model, i.e.,

where $\phi(X_i, W_i) \in \R^d$ is some feature transformation of $X_i, W_i$ and $\epsilon_i \iid \mcN(0, \sigma^2)$. To fit this model, one can (i) adaptively fit the feature representation $\phi$ using the first fold $\mcD_1$, (ii) fit a regularized estimate $\hat \beta$ of $\beta$ using (e.g.) a cross-validated lasso on $\mcD_1$, and (iii) estimate $\sigma^2$ using the usual OLS estimator of the residual variance. Of course, the Gaussian assumption may not always be realistic. Instead, our default implementation in dualbounds fits the same coefficients $\hat \beta$ and uses the empirical residuals $\hat \epsilon_i \defeq Y_i - \phi(X_i, W_i)^T \hat \beta$ to nonparametrically estimate the law of $\epsilon_i$. Similarly, in the presence of heteroskedasticity, we can estimate $\var(Y_i \mid X_i, W_i)$ using a nonparametric estimator like a random forest; clearly, the possibilities are endless. The main point is that misspecification of these models will not affect the validity of $\hat\theta_\mathrm{LCB}$, although better models will yield tighter estimates and confidence intervals.

Examples

In this section, we give a few examples of estimands that fit into the framework from Section (ref).

\begingroup

\addtocounter{example}{-1} \endgroup

\begingroup

\addtocounter{example}{-1} \endgroup

We now return to the case of Lee bounds (Example (ref)) from Section (ref).

\begingroup

\addtocounter{example}{-1} \endgroup

The ideas in Example (ref) apply to any quasilinear function of $P$. Two examples are given below.

Theory

Uniform validity

For expositional convenience, we assume $|\mcD_2| \ge cn$ for some constant $c > 0$ throughout the section. Throughout, $\E[\cdot \mid \mcD_1]$ denotes an expectation conditional on the first fold of data. All proofs will be presented in Appendix (ref).

Our first main theoretical result is that in randomized experiments, $\hat\theta_\mathrm{LCB}$ is a valid $1-\alpha$ lower confidence bound on $\theta_L$ under arbitrary model misspecification. Note that the following result allows for the analyst to use any method to estimate the optimal dual variables $\hat \nu$ as long as $\hat \nu \in \mcV$ are dual-feasible. It also places no restrictions on the relationship between the potential outcomes and $X_i$, although we do require the following moment condition on $\hat \nu$.

Assumption (ref) is weak, since in practice one could always “clip" $\hat \nu$ below some large value to ensure its moments exist without violating dual feasibility. It can also be substantially relaxed at the cost of a more technical statement (see Appendix (ref), Remark (ref)). All we need is for the moments of $S_i$ to be sufficiently regular such that we can apply a univariate central limit theorem (CLT) to $\{S_i\}_{i \in \mcD_2}$ conditional on $\mcD_1$.

Of course, by multiplying $\theta(P)$ by negative one, these theorems prove that we can get a $1-\alpha$ upper confidence bound $\hat \theta_{\mathrm{UCB}}$ on the sharp upper bound $\theta_U$. These bounds can be combined to cover either the partially identified set or the parameter $\theta(P\opt)$ imbensmanski2004, stoye2009more (see Section (ref)).

Theorem (ref) has two key ingredients---(i) weak duality plus (ii) the fact that $\tilde{\theta}_L$ has a representation as a marginal moment, which allows us to apply the CLT. As discussed in Section (ref), these properties also allow us to use the multiplier bootstrap to select a “good" choice of $\hat \nu$. In particular, the multiplier bootstrap is asymptotically valid as long as the central moments of the IPW summands $S_i^{(k)}$ do not grow too quickly with $n$ and $K$, as stated formally below.

This assumption is weak and is standard in the literature. When $\hat{\nu}_1^{(k)}(Y_i, X_i), \hat{\nu}_0^{(k)}(Y_i, X_i)$ are uniformly bounded, it is satisfied if $\log K = O(n^{1/7 - \epsilon})$ for some $\epsilon > 0$, meaning that we can select from many different models without sacrificing validity.

Tightness

Each of the previous results relies on the weak duality result that $\tilde{\theta}_L \le \theta_L$ holds deterministically. Although this ensures that $\hat\theta_\mathrm{LCB}$ and $\hat\theta_\mathrm{LCB}\mb$ are valid lower confidence bounds, one might worry that it will make inference too conservative. We investigate this question in this subsection.

General analysis

We now give high-level conditions under which $\hat\theta_\mathrm{LCB}$ converges to $\theta_L$ at oracle rates. The main intuition follows from the decomposition

The univariate CLT suggests that the second term is asymptotically exact. Thus, the main question is how large the first-stage bias is.

The following theorem tells us that the first stage bias is bounded by the product of the errors in estimating $(\PZC\opt,\POC\opt)$ and $\nu\opt$. Thus, if the product of the errors decays at an $o(n^{-1/2})$ rate, the first stage bias will be negligible compared to the variance from the univariate CLT. As notation, let $p_0\opt(y_0 \mid x), p_1\opt(y_1 \mid x)$ denote the conditional densities of $Y(0) \mid X$ and $Y(1) \mid X$ with respect to some base measure $\psi$ on $\mcY$ \footnote{We choose $\psi$ to be the Lebesgue measure for continuous outcomes and the counting measure for discrete outcomes.}; similarly, let $\hat p_1(y_1 \mid x), \hat p_0(y_0 \mid x)$ denote the estimated densities under $\hatPC$.

For each $x \in \mcX$, we define $\mathrm{error}_P(x)$ to be the $\ell_2$ distance between $(p_0\opt(\cdot \mid x), p_1\opt(\cdot \mid x))$ and $(\hat p_0(\cdot \mid x), \hat p_1(\cdot \mid x))$:

Similarly, we define $\mathrm{error}_\nu(x)$ to be the corresponding $\ell_2$ distance between $\hat \nu$ and $\nu\opt$:

Overall, Theorem (ref) gives intuition that if strong duality holds, the first stage bias should decay at a faster rate than $\E_{X \sim P_X\opt}[\mathrm{error}_P(X) \mid \mcD_1]$, which represents the error in estimating the outcome model. Intuitively, this is because if $\hatPC$ is close to $\PCopt$, then $\hat \nu$ should be close to $\nu\opt$, since $\hat \nu$ maximizes the empirical dual based on $\hatPC$, and $\nu\opt$ solves the population dual based on $\PCopt$.

We emphasize that as long as strong duality holds, Theorem (ref) makes no assumptions whatsoever about the form of $\theta(P)$, the dimension of the covariates $X$, or the model class $\mcP$---furthermore, it is a finite-sample result with no “hidden" constants.

Refined theory for discrete potential outcomes

Previously, we used Theorem (ref) to argue that the first-stage bias of dual bounds decays faster than the estimation error of the outcome model $\mathrm{error}_P(X)$. Now, we formalize this intuition in the case where $Y$ has finite support and $\hat{\nu}$ are chosen as the dual variables corresponding to $(\hat{P}_{Y(0)\mid X}, \hat{P}_{Y(1)\mid X})$. In particular, we use a technical tool called Hoffman constants hoffman1952, which measure the stability of linear programs. We provide a detailed discussion of Hoffman constants in Appendix (ref). Lemma (ref) now shows that for each $x \in \mcX$, the error in estimating the dual variables decays linearly in $\mathrm{error}_P(x)$.

Note that Lemma (ref) allows for settings where the optimal dual variables $\nu\opt$ are not unique. Nonetheless, there always exists some choice of $\nu\opt \in \argmax_{\nu} g(\nu)$ such that Lemma (ref) holds.

Combining Theorem (ref) and Lemma (ref) establishes that if the error in estimating the outcome model, $\mathrm{error}_P(x)$, decays at $o(n^{-1/4})$ rates, then the effective estimand $\tilde{\theta}_L = g(\hat{\nu})$ is statistically indistinguishable from the sharp lower bound $\theta_L$. To state this result, we denote $Z_n = o_{L_k}(a_n)$ for a sequence of random variables $Z_n$ and fixed numbers $a_n$ if $(\E[Z_n^k])^{1/k} = o(a_n)$.

Theorem (ref) shows that as long as one can estimate the conditional laws of the potential outcomes at semiparametric rates, then $\hat\theta_L$ is asymptotically unbiased. Furthermore, the proof of Theorem (ref) shows that $\hat \theta_L$ is asymptotically equivalent to the “oracle" estimator which has perfect knowledge of the outcome model and uses the optimal dual variables $\nu\opt$ in place of $\hat \nu$. Please see Appendix (ref) for further details.

Cross fitting

This subsection shows that employing cross-fitting can recover the factor of two lost by sample splitting, without sacrificing validity under most forms of outcome model misspecification or tightness when the outcome model can be estimated at $o(n^{-1/4})$ rates. As notation, let $\hat\theta_L^\mathrm{swap}$ denote the same estimator as $\hat \theta_L$ but with the roles of $\mcD_1$ and $\mcD_2$ swapped. The cross-fit estimator is then

A cross-fit lower confidence bound can be computed as follows. Let $\hat \nu$ and $\hat \nu^\mathrm{swap}$ denote the estimated dual variables from $\mcD_1$ and $\mcD_2$, respectively. For ease of exposition, we assume $n$ is even and $|\mcD_1| = |\mcD_2| = n/2$. \footnote{The results in this section can be easily extended to $M$-fold cross-fitting for $M > 2$.} Let $S_i = \frac{\hat\nu^\mathrm{swap}_{1,X_i}(Y_i) W_i}{\pi(X_i)} + \frac{\hat\nu^\mathrm{swap}_{0,X_i}(Y_i) (1-W_i)}{1 - \pi(X_i)}$ if $i \in \mcD_1$. If $i \in \mcD_2$, let $S_i$ be defined analogously but with $\hat \nu^\mathrm{swap}$ replaced with $\hat \nu$. Then if $\hat \sigma_s^\mathrm{crossfit}$ is the empirical standard deviation of $\{S_i\}_{i=1}^n$, the cross-fit lower confidence bound is

We first establish validity when $\hat{\nu}$ is potentially inconsistent. Due to the dependence introduced by cross-fitting, we need more regularity conditions to show an analogue of Theorem (ref). Interestingly, we show that $\hat\theta_\mathrm{LCB}^\mathrm{crossfit}$ is valid under two separate and non-nested conditions.

The first condition of Theorem (ref) shows that if the estimated dual functions $\hat \nu_0, \hat \nu_1$ are asymptotically deterministic, though the limits may differ from ($\nu_0\opt, \nu_1\opt)$, $\hat\theta_\mathrm{LCB}^\mathrm{crossfit}$ is a valid lower confidence bound. Similar conditions on estimated nuisance parameters have been studied in other contexts chernozhukov2020adversarial, arkhangelsky2021double. A strength of this result is that it allows $\hat \nu_0, \hat \nu_1$ to converge at arbitrarily slow rates. Indeed, the proof technique for this result is based on a novel argument leveraging weak duality; it is not necessarily true that under Condition 1, $\hat\theta_\mathrm{LCB}^\mathrm{crossfit}$ is equivalent to an “oracle" confidence bound of any form. The second condition suggests that even if this is not true and the fluctuations of $\hat \nu_0, \hat \nu_1$ do not vanish asymptotically, cross-fitting can be valid if the first-stage bias is sufficiently large, making $\hat\theta_\mathrm{LCB}^\mathrm{crossfit}$ conservative but valid. Except in pathological examples, we expect the second condition to hold whenever the first condition does not. Intuitively, if $\hat\nu$ has non-vanishing fluctuations, this suggests that $\hat\nu$ is not consistently estimating $\nu\opt$, in which case we should expect a substantial conservative bias, satisfying Condition 2. Thus, in practice, we recommend using cross-fitting.

Now we turn to tightness of cross-fitting when the outcome model can be estimated at $o(n^{-1/4})$ rates. Analogous to $\tilde{\theta}_L$, we define the effective estimand of $\hat\theta_\mathrm{LCB}^\mathrm{crossfit}$ as

We now prove that under the same conditions as in Theorem (ref), namely discrete potential outcomes and semiparametric convergence rate of $\mathrm{error}_P(X)$, $\tilde{\theta}_L^\mathrm{crossfit}$ is statistically indistinguishable from the sharp lower bound $\theta_L$.

Dual bounds for observational studies

So far, our theory has assumed that the propensity scores are known. However, when $\pi(X_i)$ is unknown, we can replace the IPW estimator with an augmented IPW (AIPW) estimator to increase robustness. In particular, define the conditional mean of the estimated dual variables $\hat\nu$ as

\[c_0(x) \defeq \E_{P\opt_{Y(0)\mid X = x}}[\hat \nu_{0,x}(Y(0))] \text{ and } c_1(x) \defeq \E_{P\opt_{Y(1)\mid X = x}}[\hat \nu_{1,x}(Y(1))]\]

so $c_k(X_i)$ is the conditional mean of $\hat \nu_{k,X_i}(Y(k))$ given $X_i$ and $\mcD_1$. Also, let $\hat c_0(x), \hat c_1(x)$ denote estimators of $c_0(x), c_1(x)$ fit on $\mcD_1$; for example, one can automatically compute $\hat c_0(x), \hat c_1(x)$ by plugging in $\hatPC$. Lastly, for any $i \in \mcD_2$, define the AIPW summand

where $\hat \pi$ are propensity scores estimated on $\mcD_1$. Then, if $\hat \sigma_s^\mathrm{aug}$ is the sample standard deviation of $\{S_i\}_{i \in \mcD_2}$ on $\mcD_2$, the “augmented" version of $\hat\theta_\mathrm{LCB}$ is

We can now prove a validity result for $\hat\theta_\mathrm{LCB}^\mathrm{aug}$. There are two cases. In the first case, we assume that the product of estimation errors for the outcome model and propensity scores decays faster than $o(1/n)$, in which case $\hat\theta_\mathrm{LCB}^\mathrm{aug}$ will be a valid lower confidence bound for $\tilde{\theta}_L$ based on standard results for the AIPW estimator robins1994aipw. However, even outside this standard regime, $\hat\theta_\mathrm{LCB}^\mathrm{aug}$ may still be valid. In the second case, we assume the outcome model is sufficiently misspecified such that the first stage bias $\theta_L - \tilde{\theta}_L$ dominates either the error in estimating $\pi$ or the error in estimating $c$. In this situation, the fluctuations of $\hat\theta_\mathrm{LCB}^\mathrm{aug}$ around $\tilde{\theta}_L$ are of smaller order than the first-stage bias.

See Appendix (ref) for a proof.

Computation

General strategy and ensuring validity

In this section, we discuss how to compute the dual bounds in Definition (ref). Computation is straightforward except for two questions:

We now outline a general strategy to answer these questions based on two key observations. Note that for simplicity, in this section, we assume the response $Y$ is real-valued.

Observation 1: the problem separates in $\mcX$. Theorem (ref) makes clear that to compute $\hat\nu \approx \argmax_{\nu \in \mcV} g(\nu)$, it suffices to repeatedly solve the problem conditional on $x$. The solutions to these problems are independent in the sense that the value of $\hat\nu_{0,x},\hat\nu_{1,x}$ does not affect the value of $\hat\nu_{0,x'},\hat\nu_{1,x'}$ for some $x' \ne x$. Similarly, by definition we have that $\hat \nu \in \mcV$ is dual-feasible if and only if $\hat\nu_{0,x}, \hat\nu_{1,x} \in \mcV_x$ for all $x \in \mcX$.

Observation 2: Only compute what we need. To apply dual bounds, we need to only compute $\{\hat\nu_{0,x}, \hat\nu_{1,x}\}_{x \in \{X_i : i \in \mcD_2\}}$ to compute the IPW estimator $\hat\theta_L$ and lower confidence bound $\hat\theta_\mathrm{LCB}$---i.e., we do not need to solve Eq. ((ref)) for all $x \in \mcX$.

These observations have two implications.

Implication 1: ensuring validity. Given any initial estimate $\hat\nu\init$ which may or may not be dual-feasible, we can convert $\hat\nu\init$ into dual-feasible estimators as follows:

For each $x \in \mcX$, $c_x$ can be computed using a two-dimensional grid search---crucially, because this grid search is low-dimensional, we can accurately compute $c_x$. Furthermore, Observation 2 implies that we only need to compute $c_x$ for $\{X_i : i \in \mcD_2\}$. As a result, the steps above represent a generic algorithm to convert any initial estimates $\hat\nu\init$ into valid dual estimates $\hat\nu\in\mcV$ via $|\mcD_2|$ grid searches.

Implication 2: a generic strategy for computing optimal dual variables. Similarly, to compute $\hat\nu \in \argmax_{\nu \in \mcV} g(\nu)$, we have the following general strategy:

In other words, we need to only solve the conditional problem $|\mcD_2|$ times to compute the dual bounds. We discuss how to do this in the next section.

Finding conditionally optimal dual variables

We suggest a discretization-based method to approximately solve this conditional problem (ref) and obtain initial estimates $\hat\nu_{0,x}\init, \hat\nu_{1,x}\init$.\footnote{Our software implements this method by default, although Appendix (ref) discusses an alternative approach based on series estimators.} The idea is to approximate $\hat P_{Y(k) \mid X = x}$ as a discrete distribution with support $\{y_{k,1,x}, \dots, y_{k,\nvals,x}\}$ and probability mass function (PMF) $p_{k,1,x}, \dots, p_{k,\nvals,x} \in (0,1)$ so that

where $\delta_z$ denotes the point mass on $z \in \R$. In particular, we suggest taking $y_{k,j,x}$ as the $\frac{j}{\nvals+1}$th quantile of $\hat P_{Y(k) \mid X = x}$ and setting $p_{k,j,x} = \frac{1}{\nvals}$ for $k \in \{0,1\}, j\in \{1, \ldots, \nvals\}$. The conditional optimization problem then becomes a discrete linear program with $2\nvals+L$ variables and $\nvals^2+L$ constraints:

where the optimization variables are $\{\nu_{0,x}(y_{0,j,x})\}_{j=1}^{\nvals}, \{\nu_{1,x}(y_{1,i,x})\}_{i=1}^{\nvals}$ and $\lambda_{x,1}, \dots, \lambda_{x,L}$. This problem can be solved efficiently using off-the-shelf LP solvers if, e.g., $\nvals \le 100$. Furthermore, when $\mcW_x = \emptyset$, this is the dual to a standard optimal transport problem, so it can be solved even more efficiently using specialized solvers such as the network simplex algorithm flamary2021pot. After solving this problem, we obtain initial values $\{\hat\nu_{0,x}\init(y_{0,j,x})\}_{j=1}^{\nvals}, \{\hat\nu\init_{1,x}(y_{1,i,x})\}_{i=1}^{\nvals}$ and we define the full functions $\hat\nu_{0,x}\init, \hat\nu_{1,x}\init : \R \to \R$ via linear interpolation. Then, as described in Section (ref), we can use a two-dimensional grid search to obtain valid dual variables $\hat\nu_{0,x}, \hat\nu_{1,x}$. As discussed in Remark (ref), this gridsearch ensures that the final confidence bounds are valid even if the discretization yields an inaccurate initial solution $\hat\nu_{0,x}\init, \hat\nu_{1,x}\init$.

Empirical applications

We now illustrate our method in applications to two randomized experiments and one observational study. Code and data are publicly available at \url{https://github.com/amspector100/dual_bounds_paper}.

Persuasion effects of political news

We first analyze data from gerber2009, who in 2005 randomly assigned a set of individuals in Prince William County, Virginia, to receive a free subscription offer for the Washington Post.\footnote{The original experiment had a third treatment condition, namely to receive a free subscription offer for the Washington Times. For simplicity, we follow jun2023persuasion and only analyze subjects in the Washington Post or control treatment groups.} Using administrative data, they also determined whether each subject voted in the November 2006 elections. Thus, for $n=2400$ individuals, $W_i \in \{0,1\}$ denotes whether individual $i$ received a free subscription to the Washington Post, and $Y_i \in \{0,1\}$ denotes whether individual $i$ voted in the 2006 elections.

In this context, jun2023persuasion (henceforth JL) studied the “persuasion effect" of the treatment, defined as the probability that the treatment causes an individual who would not otherwise have voted:

This estimand is also known as the Probability of Sufficiency pearl1999probabilities. As noted by JL, without covariates, the sharp bounds on $\theta(P\opt)$ are rescaled \fh\, bounds:

However, gerber2009 also collected a rich set of covariate information, including demographic information, political preferences, and previous voter turnout data. Furthermore, $\theta(P\opt)$ takes the form of an unidentifiable expectation divided by an identifiable expectation (since $P\opt(Y(0) = 0)$ is identified in Eq. ((ref))). Thus, we can use our methodology to form covariate-assisted estimates of the numerator and apply the bivariate delta method to perform inference on $\theta(P\opt)$, as described in Appendix (ref).

To form the dual bounds, we estimate the conditional laws of $Y(1) \mid X$ and $Y(0) \mid X$ using three outcome models: a cross-validated logistic ridge regression, a random forest, and a k-nearest neighbors (KNN) classifier, where the covariates are the $43$ baseline covariates from gerber2009 plus interaction terms with the treatment. For each outcome model, we form dual bounds following the methodology from Sections (ref) and (ref) using $10$-fold cross-fitting. We also compute non-robust plug-in bounds, which plug in the estimated conditional distributions and the empirical law of $X$ into Eq. ((ref)); unlike dual bounds, these bounds can be anti-conservatively biased. We also aggregate the results across all dual bounds using the multiplier bootstrap-like procedure detailed in Appendix (ref).

Table (ref) shows the results, from which we report three main findings. First, the covariate-assisted dual bounds are more than twice as narrow as the covariate-free bounds. Second, the dual bounds appear to be more reliable than the covariate-assisted plug-in bounds. For example, the KNN and random forest outcome models produce plug-in lower bounds larger than $15\%$. This is implausible because the ATE point estimate is $0.029$ and not significant; indeed, we do not even have power to reject the sharp null that $Y(1) = Y(0)$ with probability one. In contrast, dual bounds can leverage each outcome model to provide provably valid confidence bounds without assuming that the outcome model is accurate. Third, the multiplier bootstrap method successfully selects the tightest lower and upper bounds while providing rigorous uncertainty quantification.

Estimating intensive margins

carranza2022 conducted a randomized experiment in South Africa where treated individuals received assessment results that they could share with potential employers. They found that treated individuals had higher employment rates and higher earnings, suggesting that the tests provided useful information about workers' skills. However, we might wonder: is the treatment effect driven by increases in employment (extensive margin), or does the treatment increase hours worked for individuals who would have been employed with or without the treatment (intensive margin)?

To estimate the intensive margin, chenroth2023 (henceforth CR) analyzed the following quantities:

where above, the outcome $Y$ measures the average hours worked per week post-treatment, and the logs in the latter estimand ensure that it is scale-invariant and can roughly be interpreted as a “percentage" effect. CR bounded these quantities using the methodology from lee2009training, which assumes that $Y(1) > 0$ holds whenever $Y(0) > 0$, i.e., the treatment does not cause any individual to be unemployed. To defend this assumption, CR noted that individuals with poor test results likely did not share them with their employers, and we agree that this assumption seems plausible in this setting.

However, the dataset from carranza2022 contains a rich set of pre-treatment covariates, including baseline earnings, demographic information, and educational history. Thus, we produce covariate-assisted variants of the bounds from CR. To fit the outcome model, we use the default settings in the dualbounds package, which employs a linear model with interactions:

To estimate $\beta$ and $\gamma$, we use a cross-validated ridge regression. We estimate the law of $\epsilon_i \mid X_i, W_i$ as the empirical law of the estimated residuals $\{Y_i - X_i^T\hat\beta - W_i X_i^T \hat\gamma : i \in [n], W_i = w\}$.\footnote{This estimate severely restricts the heteroskedasticity pattern, since it asserts that the residuals are independent of the covariates given the treatment, i.e., $\epsilon_i \Perp X_i \mid W_i$. That said, we emphasize that the final dual bounds are valid even if the model for the law of $\epsilon_i \mid X_i, W_i$ is completely inaccurate.} We then convert this outcome model into a cross-fit dual bound using the methodology from Sections (ref) and (ref).

Table (ref) shows the results: for both the logged and non-logged outcome, the covariate-assisted bounds are only $\approx 60\%$ as wide as the covariate-free bounds. Although the bounds are still quite wide, this analysis nonetheless shows that covariate adjustment can substantially sharpen partial identification bounds without requiring additional assumptions.

401k eligibility

We now study how 401(k) eligibility impacts wealth. An extensive literature argues that 401(k) eligibility is essentially exogenous conditional on covariates poterba1995, poterba1998, poterba2000, chernhansen2004, since workers likely choose employers based on job characteristics besides 401(k) eligibility, e.g., income. We adopt this assumption; thus, the outcome $Y \in \R$ measures total household wealth and the treatment $W \in \{0,1\}$ indicates 401(k) eligibility.

We obtain data from chernozhukov2018dml, who estimated average treatment effects (ATEs) using a sample of households from the 4th wave of the 1990 Survey of Income and Program Participation. Yet the literature emphasizes that treatment effects may be highly heterogeneous. For instance, 401(k) eligibility may reduce wealth for households who otherwise would participate in a different retirement plan or whose 401(k) contributions adversely reduce their liquidity. To study whether negative effects contribute to the overall ATE, we now bound the positive treatment effect:

Our analysis uses the same covariates as chernozhukov2018dml, including income, demographics, and financial indicators such as homeownership status. Following chernozhukov2018dml, we use the raw covariates except when fitting regularized GLMs, where we include polynomial transformations and pairwise interactions. We estimate cross-fit propensity scores using a cross-validated logistic elastic net.

The outcome model takes the form $Y_i = \E[Y_i \mid X_i, W_i] + \epsilon_i$. To estimate $\E[Y_i \mid X_i, W_i]$, we use (i) a cross-validated elastic net, (ii) a KNN regressor, (iii) an sklearn histogram gradient boosting (HGBoost) regressor sklearn, lightgbm2017 with the constraint that $\E[Y_i \mid X_i, W_i]$ is increasing in $W_i$,\footnote{Note that while the “HGBoost monotone" model asserts that the conditional average treatment effect $\E[Y(1) - Y(0) \mid X]$ is nonnegative, it nonetheless allows the positive treatment effect $\theta(P\opt) = \E[\max(Y(1) - Y(0), 0)]$ to differ from the ATE, for example, because there may be unobserved covariates $U$ such that $\E[Y(1) - Y(0) \mid X, U]$ may be negative.} and (iv) an intercept-only model as a baseline. We estimate the law of $\epsilon_i \mid X_i, W_i$ as in Section (ref). For each model, we report cross-fit AIPW dual bounds (as described in Sections (ref)-(ref)) as well as non-robust “plug-in" bounds, which plug $\hat P_{Y \mid X, W}$ and the empirical law of the covariates into Eq. ((ref)).

Table (ref) shows the out-of-sample $R^2$, cross-fit AIPW ATE estimates, dual bounds, and plug-in bounds for each outcome model. We report two main conclusions.

1. Incorporating covariates improves robustness. It is known that in observational studies, accurate outcome models can reduce the bias of ATE estimates. E.g., in Table (ref), more accurate outcome models yield smaller ATE estimates, ranging from $\approx \$11$K to $\approx \$6$K. Table (ref) suggests that the same logic applies to partially identified estimands (see Theorem (ref)), since the dual lower bounds decrease with the ATE estimates.

2. Plug-in bounds can be anticonservative. The KNN plug-in lower bound is $\approx \$10$K. This value seems implausible, since it is twice as large as the corresponding dual bound and $50\%$ larger than the ATE estimate from the best-performing model. Indeed, covariate-assisted plug-in bounds rely entirely on the accuracy of the outcome model, whereas dual bounds are doubly robust as per Theorem (ref). That said, it is reassuring that the best-performing model (HGBoost) yields similar plug-in and dual bounds.

A Monte-Carlo simulation

In this section, we run simulations to demonstrate the power, validity, and computational efficiency of dual bounds. Throughout, we consider randomized experiments where the propensity scores $\pi(x) = \frac{1}{2}$ are known. Replication code is available at \url{https://github.com/amspector100/dual_bounds_paper/}.

We perform simulations where we estimate lower Lee bounds (Example (ref)). We sample covariates $X_{i} \iid \mcN(0,I_p)$ for $p=20$ covariates and draw $Y_i \mid X_i$ from a homoskedastic Gaussian linear model:

for variance $\sigma^2 = 1$, coefficients $\beta \in \R^p$ chosen such that $\var(Y_i(1)) = \var(Y_i(0)) = 10$, and average treatment effect $\tau = 2$. We sample the selection events $S_i \mid X_i$ from a logistic regression model:

with $\|\beta_S\|_2 = 1$ and $\tau_{S,0} = 0, \tau_{S,1} = 1$. Following general practice in the literature, our simulations enforce the monotonicity condition $S(1) \ge S(0)$ a.s., and we assume that the practitioner knows this a-priori. We compare three methods for estimating the sharp bound $\theta_L$ in this problem:

We also consider the performance of each method in two misspecified settings where $Y_i \mid X_i$ is actually heteroskedastic:

for the functions $\sigma_1^2(X) = \sigma_1^2 \|X\|_2^2, \sigma_0^2(X) = \sigma_0^2 \|X\|_2^2$ for constants $\sigma_1, \sigma_0 \ge 0$. In this case, for both the naive plug-in and dual crossfit methods, the estimated outcome model is misspecified, since it incorrectly assumes homoskedasticity. In the first setting (labelled as “Heteroskedasticity (I)"), we set $\sigma_1 / \sigma_0 = 3$; in the other setting (“Heteroskedasticity (II)"), we set $\sigma_0 / \sigma_1 = 0.3$.

Figures (ref) and (ref) show the results with $n \in \{100, 200, 600, 1000\}$. Figure (ref) shows the average value of the estimate $\hat\theta_L$ and the lower confidence bound $\hat\theta_\mathrm{LCB}$; it shows that the naive plug-in estimator is biased when $n$ is small (due to the effect of regularization) and when the model is misspecified. In contrast, the cross-fit dual bounds are (i) guaranteed to be conservatively biased at worst and (ii) less sensitive to errors in estimating the outcome model, yielding valid and reasonably sharp inference in all three settings. Figure (ref) confirms that dual bounds provide $\ge 95\%$ coverage in all settings, whereas the naive plug-in method can be quite conservative or anticonservative, depending on the form of heteroskedasticity. Overall, in this setting, cross-fit dual bounds perform well even in small samples.

Discussion

This paper introduces a dual bound method to estimate and perform inference on a class of partially identified causal parameters. The method can leverage any statistical and machine learning techniques to learn the conditional distribution of the outcome $Y$ given the covariates $X$. In randomized experiments, the resulting bounds are always valid regardless of whether the estimates are consistent and asymptotically sharp when the conditional distributions are estimated at semiparametric rates. In addition, one can apply the multiplier bootstrap to perform model selection. For observational studies, the method can be easily extended to be doubly robust. In all settings, the dual bounds can be computed efficiently.

Our analysis leaves open many questions. For example, a few of our theoretical results require $Y$ to be discrete, and it would be interesting to investigate if the same results hold when $Y$ is continuous. Perhaps the most pressing question is whether the techniques developed in this paper can be applied more generally. In particular, we use a duality argument to guarantee the robustness of our method. Does this same argument apply to settings beyond causal inference? In the next two sections, we begin to address this question. Then, Section (ref) discusses an alternative computational strategy, and Section (ref) discuss two-sided intervals.

Extensions beyond causal effects

In this section, we discuss whether our method can be extended to settings that cannot be reduced to estimating an expectation of the form $\E[f(Y(0), Y(1), X)]$. Indeed, the ideas in Section (ref) apply to many estimands in economics which can be written as the optimal value of an optimization problem. We describe two classes of problems below.

Cost of robustness when the propensity scores are unknown

However, dual bounds as defined in Section (ref) are not appropriate for every problem. For example, one might hope that a simple modification of Definition (ref) can produce always valid bounds on the average treatment effect when the propensity scores $\pi(X_i)$ are not known. Unfortunately, our method yields valid yet trivial lower and upper bounds due to the lack of strong duality. This is consistent with aronow2021nonparametric, who prove that no uniformly consistent estimator of ATE exists under strong ignorability and strict overlap without further assumptions if one of the covariates is continuous.

Suppose the vectors $(X_i, W_i, Y_i(0), Y_i(1))$ are sampled i.i.d. from some population distribution $P\opt \in \mcP$, where $\mcP$ is the set of distributions on $\mcX \times \{0,1\} \times \mcY^2$ satisfying unconfoundedness, i.e., $\{Y_i(1), Y_i(0)\} \Perp W_i \mid X_i$, and strict overlap, i.e., $0 < \pi_P(x) < 1$ for all $x\in \mcX$ and $P \in \mcP$, where $\pi_P(X) \defeq \E_P[W \mid X]$.

Given i.i.d. observations $(X_i, W_i, Y_i)$, we seek to form a lower bound on the average treatment effect $\theta(P\opt) \defeq \E_{P\opt}[Y_i(1) - Y_i(0)]$ which is valid even under arbitrary misspecification of $\pi_P(X_i)$ and the outcome model. Although $\theta$ is identifiable, it can still be written as the solution to the optimization problem

where $P_{X,W,Y}$ is the law of $(X, W, Y)$ under $P$ and $P\opt_{X,W,Y}$ is the true law of $(X, W, Y)$. Note that the optimization variable is $P$, which is a joint law over $(X,W,Y(0),Y(1))$, and $P_{X,W,Y}$ is a functional of $P$. For any $h : \mcX \times \{0,1\} \times \mcY \to \R$, the Lagrange dual to this problem is

where $\kappa(h) \defeq \inf_{P \in \mcP} \E_P[Y(1) - Y(0) - h(X, W, Y)]$ is a known constant depending on $h$. For any $h$, $g(h)$ is a valid lower bound on $\theta(P\opt)$ by weak duality, but unfortunately, strong duality does not hold. In particular, for any $h$, we have that

See Appendix (ref) for a proof.

This result tells us that any dual bound (in the sense of Def. (ref)) on the ATE which is valid under arbitrary misspecification must also be trivial, since it must impute $Y_i(1)$ to have the minimum possible value whenever it is not observed, and it must impute $Y_i(0)$ to have the maximum possible value when it is not observed. Thus, applied in this way, our method reduces to the nonparametric bounds from manski1989, even though we have made the extra assumptions of unconfoundedness and strict overlap, which are not made in manski1989.

An alternative computational strategy

To compute our recommended estimator $\hat \nu$ of the optimal dual variables, one must first model the conditional laws $\PCopt$. This procedure may not be feasible when conditional distributions are hard to model. For example, when $X$ includes unstructured data such as images (e.g., profile pictures as in athey2022smiles), texts (e.g., resumes as in vafa2022career), and embeddings vafa2024estimating, existing machine learning algorithms may not be able to provide distribution estimates. In Appendix (ref), we present Deep Dual Bounds, an alternative approach that directly learns the optimal dual variables by solving the dual problem. This approach parametrizes $\nu_{0,X}(Y(0)), \nu_{1,X}(Y(1))$ by neural networks, which, unlike the two-stage approach, can exploit the smoothness of the dual functions in $X$.

Although this end-to-end formulation is conceptually clear, standard gradient-based algorithms cannot be directly applied since potential outcomes cannot be observed simultaneously. We resolve this issue by matching treated and control units and optimizing an approximate objective function. We emphasize that the approximation error does not affect the validity of Dual Bounds, as only dual feasibility is required. The details and experimental results of the algorithm are discussed in Appendix (ref) and Table (ref).

Two-sided confidence intervals

In previous sections we focused on one-sided confidence bounds on the sharp population bounds $\theta_L, \theta_U$. To cover the full identified set, we can simply construct $(1-\alpha/2)$ lower/upper confidence bounds on the lower/upper bounds horowitz2000nonparametric. However, in many applications, it suffices to cover the true parameter. It is well-known imbensmanski2004, stoye2009more that tighter uniform confidence intervals can be constructed by estimating the gap between upper and lower bounds and/or the correlation between two estimators.

With data splitting, the upper and lower dual bounds are both empirical moments: \[\hat{\theta}_L = \frac{1}{|\mathcal{D}_2|}\sum_{i\in \mathcal{D}_2}S_i^{L}, \quad \hat{\theta}_U = \frac{1}{|\mathcal{D}_2|}\sum_{i\in \mathcal{D}_2}S_i^{U}.\] Under mild regularity assumptions on the marginal moments of $S_i^L$ and $S_i^U$ discussed in Section (ref), $(\hat{\theta}_L, \hat{\theta}_U)$ is asymptotically bivariate Gaussian and the empirical covariance matrix of $(S_i^L, S_i^U)_{i\in \mathcal{D}_2}$ is a consistent estimate of the true asymptotic covariance matrix. While the superefficiency assumption in imbensmanski2004 does not necessarily hold in our case, we can apply the construction studied in Proposition 3 of stoye2009more to guarantee the uniform coverage of the true parameter.