Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
85,444 characters · 12 sections · 28 citation commands
Partial identification via conditional linear programs: estimation and policy learning
\if00 {
\thispagestyle{empty} }\fi
\if10 \thispagestyle{empty} \fi
\pagenumbering{gobble}
\pagenumbering{arabic} \onehalfspacing
Important estimands are oftentimes only partially identified from the observable data. Even an infinite amount of data cannot uniquely determine the estimand of interest, and instead can only limit it to be within a set of plausible values. For example, in causal inference settings it may not be possible to point-identify average treatment effects if treatment assignment is confounded; even with auxiliary information, such as an instrumental variable, the average treatment effect is not necessarily point identified. There is also often interest in complex counterfactual estimands that go beyond average treatment effects, for instance the proportion of individuals who benefit from treatment or characteristics of the spread of treatment effects across the population.
This paper develops a unified framework for covariate-assisted estimation, inference, and decision making in a broad class of partial identification problems where the parameter of interest satisfies a series of linear constraints. In such settings, bounds for the parameter of interest are solutions to conditional linear programs (LPs) that optimize a linear function of unknown parameters subject to linear constraints, where both the objective function and the constraints may depend on additional covariates and need to be estimated. This encompasses a wide range of settings; as illustrative examples we consider estimands related to the joint distribution of potential outcomes, policy learning with inequality-aware value functions, and instrumental variable settings (Section (ref)).
We develop two de-biased estimators of bounds on the parameter of interest that are robust to errors in estimating the nuisance functions that define the objective and constraints (Section (ref)). The first approach directly solves the optimization problems using plugin estimates of the conditional constraints and objectives and leverages properties of optimal solutions of LPs to de-bias these estimates using the output of standard LP solvers. In contrast to prior work, this estimator does not require derivation of analytical solutions or enumeration of all possible solutions, which can be computationally intractable for large problems. The second approach leverages recent advances in the theory of regularized linear programs to create smooth approximations to the conditional LPs using entropic regularization and estimates the smoothed bounds in a de-biased manner. This approach trades off a small amount of approximation error for gains in estimation and computational efficiency. In contrast to alternative smoothing schemes, entropic regularization introduces approximation error that rapidly decays with the choice of hyper-parameter, and, crucially, does not require explicit enumeration of all feasible solutions.
For both estimators, we establish the conditions for asymptotic normality and construct standard, Wald-type approximate confidence intervals for the partially identified parameter. These de-biased estimators are robust to first-order errors in estimating the covariate-conditional constraints and objective vectors, allowing for valid inference even when these nuisance functions are estimated using flexible machine learning methods. We also extend these results to decision-making problems where the goal is to estimate an optimal decision policy that optimizes the lower or upper bound of a partially identified quantity (Section (ref)). We demonstrate these estimators by applying them to a study of the effects of Medicaid enrollment on health outcomes (Section (ref)).
\paragraph{Related literature.}
This paper builds on three key strands of the literature on identification and estimation of partially identified parameters. One line of work focuses on derivation of analytic or symbolic bounds balke_bounds_1997, sachs_general_2023, gabriel_sharp_2024. This work typically characterizes sharp bounds on partially identified parameters as solutions to linear programs, then enumerate all feasible solutions (i.e. all vertices of the feasible set) to obtain the bounds. Because the number of solutions can grow combinatorial with the size of the problem, this approach can be computationally intractable for larger problems. The approaches we develop do not require vertex enumeration and so can computationally scale to larger problems. duarte_automated_2024 extend these ideas to polynomial programs and use numerical approaches to bound the parameter of interest. Our approach similarly arrives at the bounds computationally, rather than analytically, but for a more limited class of problems for which we are able to provide de-biased estimators that allow us to construct straightforward Wald-type asymptotic confidence intervals.
The second line of work focuses on de-biased, covariate-assisted estimation and inference for partially identified parameters. Some of this work develops de-biased estimators of aggregations of conditional versions of analytic bounds on the parameter of interest Kallus2022_harm,levis_assisted_2023, semenova_aggregated_2024, benmichael_ai_2024, or targets smooth approximations to these bounds levis_assisted_2023, lanners_data_2025. Again, the methods we develop here differ in that they do not require analytical bounds, nor do they require enumeration of all possible solutions. ji_model-agnostic_2024 consider estimating covariate-assisted bounds on estimands involving the joint distribution of potential outcomes jiang_semiparametric_2025. Unlike the approaches we develop here, this work allows for continuous potential outcomes; however, we consider a general setup that allows for multi-level treatments, instrumental variables, and other forms of partial identification problems.
A third line of work focuses on learning optimal decision rules in settings where the value of a decision policy is only partially identified. This work builds on ideas from decision theory Manski2005_partial, and considers settings where there are unmeasured confounding variables, treatment-control overlap violations, or the value function involves counterfactual quantities Kallus2021,Pu2021,Han2021_partial, Cui2021_partial,benmichael2021_safe,Zhang2022_safe,benmichael_asymmetric_2024,DAdamo2023. The approach to policy learning here builds directly on this work and generalizes it to a broad class of problems.
We begin by describing the generic estimand and associated upper and lower bounds defined by conditional linear programs. We then discuss several specific examples of such estimands related to partially identified parameters in causal inference and policy learning settings.
Consider a set of covariates $X \in \mathcal{X}$. Our goal will be to estimate a linear functional of some true non-negative function $\bm{p}^\ast:\mathcal{X} \to \ensuremath{\mathbb{R}}^K_+$ defined by a function $\bm{c}:\mathcal{X} \to \ensuremath{\mathbb{R}}^K$, i.e.
Our focus will be on settings where the objective vector $\bm{c}(x)$ is not necessarily known but is point-identifiable from observable data while $\bm{p}^\ast(x)$ is not necessarily point-identified for all $x \in \mathcal{X}$. Rather, we can only observe $J \leq K$ linear functions of $\bm{p}^\ast(x)$, i.e. $A \bm{p}^\ast(x) = \bm{b}(x)$ for some known matrix $A \in \ensuremath{\mathbb{R}}^{J \times K}$ and some unknown function $\bm{b}:\mathcal{X} \to \ensuremath{\mathbb{R}}^J$ that is identifiable from observable data. For each value of the covariates $x \in \mathcal{X}$, these linear constraints along with non-negativity imply that the true function $\bm{p}^\ast(x)$ lies in a convex set $\mathcal{P}(x) = \{\bm{p} \in \ensuremath{\mathbb{R}}^K_+ \mid A \bm{p} = \bm{b}(x)\}$, which we will refer to as the feasible set. Note that this framework also allows for inequality constraints via standard LP transformations.
In many cases, the feasible set $\mathcal{P}(x)$ contains more than one element (e.g. because there are more equations than unknowns), and $\bm{p}^\ast(x)$ is not point-identified. This means that the parameter of interest $\theta$ may also not be point identified and so it may not be possible to estimate it precisely. We will instead focus on finding lower and upper bounds $\theta_L$ and $\theta_U$ such that $\theta \in [\theta_L, \theta_U]$. To construct these, we define $X$-conditional linear programs that minimize and maximize the linear function of $\bm{p}^\ast(x)$ over the feasible set $\mathcal{P}(x)$:
The overall bounds are then the expected solutions to the conditional linear programs, i.e. $\theta_L = \ensuremath{\mathbb{E}}[\theta_L(X)]$ and $\theta_U = \ensuremath{\mathbb{E}}[\theta_U(X)]$. In order for these bounds to be meaningful, we will assume that the the feasible set $\mathcal{P}(x)$ is both not empty and bounded for all $x \in \mathcal{X}$ and that $\|\bm{b}\|_\infty < \infty$, where for a vector-valued function $\bm{f}(x)$, we use denote $\|\bm{f}\|_\infty = \max_i \sup_{x} |f_i(x)|$. To simplify the exposition, we will also assume that the constraint matrix $A$ is sufficiently pre-processed to remove redundant constraints and so has full row rank equal to $J$. Finally, note that we have not ruled out that $\theta$ is point identified even though the function $\bm{p}^\ast(x)$ is not point identified. In such cases the lower and upper bounds match, i.e. $\theta_L = \theta_U = \theta$. As an additional bit of notation, for a scalar-valued function $f: \mathcal{X} \to \ensuremath{\mathbb{R}}$ we will denote $\|f\|_2^2 = \ensuremath{\mathbb{E}}[f(X)^2]$, and for a vector-valued function $\bm{f}:\mathcal{X} \to \ensuremath{\mathbb{R}}^p$ we will denote $\|\bm{f}\|_2^2 = \sum_{i=1}^p \ensuremath{\mathbb{E}}[f_i(X)^2] = \sum_{i=1}^p \|f_i\|_2^2$.
Next, we describe several specific examples of estimands that can be expressed in this way. \paragraph{Functions of the joint distributions of potential outcomes.} Consider a setting where we are interested in understanding the causal relationship between a discrete treatment variable $D \in \{0,\ldots, M-1\}$ and a discrete outcome variable $Y \in \mathcal{Y}$ where $|\mathcal{Y}| = L$. For example, $D$ may represent a multi-level treatment or different doses. We posit the existence of $M$ potential outcomes $(Y(0),\ldots,Y(M-1))$, one for each treatment level. In this setting, our goal is to estimate some moment of the joint distribution of potential outcomes, i.e. $\theta = \ensuremath{\mathbb{E}}[f(Y(0), \ldots, Y(M-1); X)]$ for some function $f:\mathcal{Y}^M \times \mathcal{X} \to \ensuremath{\mathbb{R}}$. We observe iid data $(X, D, \{\ensuremath{\mathbbm{1}}\{D = d\} Y(d) \}_{d = 0}^{M-1})$ from some distribution $P$.
Using the tower property, we can write the estimand as
In this case the unknown function $\bm{p}^\ast(x)$ corresponds to the joint distribution of the potential outcomes given that the covariates $X = x$, and the function $\bm{c}(x)$ is a vector collecting the function values $f(y_0, \ldots, y_{M-1}; x)$. For example, $f(y_0, \ldots, y_{M-1}; x) = y_d - y_d'$ would correspond to the average treatment effect of $d$ versus $d'$---which is identified---and $f(y_0, \ldots, y_{M-1}; x) = \max_{d} y_d$ would correspond to average outcome under an oracle treatment assignment rule that maximizes the potential outcome---this is not necessarily identified. Another unidentified parameter is the proportion of individuals who are not assigned an optimal treatment, i.e. $P(Y \neq \max_{d} Y(d))$. This corresponds to choosing $f(y_0, \ldots, y_{M-1}; x) = \sum_{d'}P(D = d'\mid X = x)\ensuremath{\mathbbm{1}}\{y_{d'} \neq \max_{d} y_d\}$. Note that if the propensity score $P(D = d' \mid X = x)$ is unknown, then the objective vector $\bm{c}(x)$ is also unknown but can be estimated.
Under the standard strong ignorability assumption that there is treatment overlap and the potential outcomes are independent of the treatment given the covariates Rosenbaum1983, the marginal distribution of each potential outcome given the covariates is identified as $P(Y(d) = y_\ell \mid X = x) = P(Y = y_\ell \mid D = d, X = x)$ for all $d \in \{0,\ldots,M-1\}$ and $\ell \in \mathcal{Y}$. This gives a set of marginal constraints on the conditional joint distribution of potential outcomes: \[ \sum_{\substack{y_0,\ldots,y_{d - 1}, \\ y_{d + 1}, \ldots, y_{M-1} \in \mathcal{Y}}} P(Y(0) = y_0, \ldots, Y(d) = y_\ell,\ldots Y(M - 1) = y_{M-1} \mid X = x) = P(Y = y_\ell \mid D = d, X = x). \] Since the probabilities sum to one, each margin gives $L - 1$ constraints. There is also an overall sum-to-one constraint on the joint distribution. This leads to a total of $J = (L-1)M + 1$ constraints, on $K = L^M$ variables (the number elements in the joint distribution) which can be collected into an integral constraint matrix $A \in \{0,1\}^{J \times K}$. The vector corresponding to the right hand side of the constraints, $\bm{b}(x)$, consists of the marginal probability distributions for each potential outcome, as well as 1 for the sum-to-one constraint.
\paragraph{Policy learning with collective utility functions.} A special case of the above is when we have a binary treatment $D \in \{0,1\}$ and we want to learn a decision policy $\pi:\mathcal{X} \to [0,1]$ that assigns $D=1$ with probability $\pi(X)$ and achieves some notion of optimality based on a function $u:\mathcal{Y} \to \ensuremath{\mathbb{R}}^+$ that assigns a utility to an outcome $y$ of $u(y)$. We can write the expected utility for an individual under decision policy $\pi$ as $\pi(X) \times u(Y(1)) + (1 - \pi(X)) \times u(Y(0))$.
To choose a policy $\pi$, we can use power-law collective utility functions that aggregate the individual utilities into a single number moulin_axioms_1988:
Here the value of the power-law parameter $\lambda$ determines the level of emphasis placed on individuals with lower or higher utilities under $\pi$. When $\lambda = 1$, the value function is simply the expected individual utility under $\pi$: $V^1(\pi) = \ensuremath{\mathbb{E}}[U(Y(0)) + \pi(X) (U(Y(1)) - U(Y(0)))]$. Values of $\lambda < 1$ will lead to inequality-averse value functions. When $\lambda = 0$, the value function is the logarithm of the geometric mean of the individual utilities: $V^0(\pi) = \ensuremath{\mathbb{E}}[\log(U(Y(0)) + \pi(X) (U(Y(1)) - U(Y(0)))]$. As $\lambda \to -\infty$, $V^\lambda(\pi)$ goes towards prioritizing only the individuals with the smallest utilities under $\pi$; choosing $\pi$ to maximize this results in the egalitarian rule that maximizes the minimum utility. Finally, values of $\lambda > 1$ will lead to inequality-seeking value functions. Note that when $\lambda \neq 1$, the power-law value function is not point-identifiable from the observed data.
To create a notion of optimality, we can compare the value under the policy $\pi$ to the value under an infeasible oracle policy that has access to the potential outcomes and assigns the treatment that maximizes the individual utility, i.e. $\ensuremath{\mathbbm{1}}\{u(Y(1)) \geq u(Y(0))\}$. The regret of policy $\pi$ relative to the oracle is thus: \[ R^\lambda(\pi) = \frac{1}{\lambda}\ensuremath{\mathbb{E}}\left[\sum_{y_0, y_1 \in \mathcal{Y}} \left(\max\{u(y_0), u(y_1)\} ^\lambda - \left\{u(y_0) + \pi(X) (u(y_1) - u(y_0))\right\}^\lambda\right) P(Y(0) = y_0, Y(1) = y_1 \mid X)\right]. \] This is a special case of Equation (ref) with $M = 2$ and $f(y_0, y_1; x) = \frac{1}{\lambda}\max\{u(y_0), u(y_1)\} ^\lambda - \frac{1}{\lambda}\left(u(y_0) + \pi(x) (u(y_1) - u(y_0))\right)^\lambda$. We will discuss finding optimal policies $\pi$ that minimize the maximum regret or maximize the minimum value in Section (ref).
\paragraph{Instrumental variables.} Now consider the case of a binary treatment $D$ and a binary instrumental variable $Z \in \{0,1\}$ that affects takeup of the treatment, where the outcome space $\mathcal{Y}$ still contains $L$ distinct levels (this discussion generalizes to non-binary treatments and instruments). The instrument leads to two potential treatments $D(0)$ and $D(1)$. Under an exclusion restriction that the instrument does not have a direct effect on the outcome variable, we then have two potential outcomes $Y(0)$ and $Y(1)$. In this case, we are interested in estimating a property of the joint distribution of potential treatments and potential outcomes $\theta = \ensuremath{\mathbb{E}}\left[f(Y(0), Y(1), D(0), D(1))\right]$ for some function $f:\mathcal{Y}^2\times \{0,1\}^2 \to \ensuremath{\mathbb{R}}$, The observable data consists of iid samples of $(X, Z, (1-Z) D(0), Z D(1), (1-D(Z)) Y(0), D(Z) Y(1))$.
We can write the estimand as \[
\] Now the unknown function $\bm{p}^\ast(x)$ corresponds to the joint distribution of the potential outcomes and potential treatments given the covariates $X = x$, and the function $\bm{c}(x)$ is a vector that collects the values of the function $ f(y_0, y_1, d_0, d_1; x) $. For example, $f(y_0, y_1, d_0, d_1; x) = y_1 - y_0$ would again correspond to the average treatment effect, while $f(y_0, y_1, d_0, d_1; x) = \ensuremath{\mathbbm{1}}\{d_0 < d_1\}$ would correspond to the proportion of compliers. Another example is the proportion of individuals who are not assigned an optimal treatment under the treatment distribution induced by a particular level of the instrument: $P(Y(D(z)) \neq \max_{d} Y(d))$. This corresponds to choosing $f(y_0, y_1, d_0, d_1; x) = \sum_{d'} P(D = d' \mid Z = z, X = x)\ensuremath{\mathbbm{1}}\{y_{d'} \neq \max_{d} y_d\}$.
Under ignorable assignment of the instrument, the joint distribution of the observed treatment and outcome satisfies \[ P(Y = y, D = d \mid Z = z, X) = \sum_{y_{1-d} \in \mathcal{Y}}\sum_{d_{1-z}=0}^1 P(Y(d) = y, Y(1-d) = y_{1-d}, D(z) = d, D(1-z) = d_{1-z} \mid X). \] Accounting for the sum-to-one constraints gives $J = 4(L-1) + 1$ constraints with $K=4L^2$ decision variables. Here the constraint matrix $A \in \{0,1\}^{J \times K}$ encodes the marginal constraints and the vector $\bm{b}(x)$ encodes the conditional joint distribution of the observed outcome and treatment variables along with the sum to one constraint.
We now turn to estimating the lower and upper bounds $\theta_L$ and $\theta_U$. Throughout, we will assume that we have iid samples of random variables $O$ that include the covariates $X$ as well as variables that give information about the constraint and objective vectors, respectively. Our principle assumption is that we can estimate the constraint vector $\bm{b}(x)$ and the objective vector $\bm{c}(x)$ in a de-biased manner using the information contained in the observed data $O$.
This structure allows us to estimate the lower and upper bounds in a de-biased manner that will be robust to errors in estimating the constraint and objective vectors. We will assume that we have estimates of the constraint and objective vectors $\hat{\bm{b}}(\cdot)$ and $\hat{\bm{c}}(\cdot)$, along with estimates of the de-biasing functions $\hat{\bm{\varphi}}^{(b)}(\cdot)$ and $\hat{\bm{\varphi}}^{(c)}(\cdot)$, that are fit on a separate sample. When it is clear from context we will use $\varphi_j^{(b)}(O) \equiv \varphi_j^{(b)}(O, b_j)$ and $\varphi_k^{(c)}(O) \equiv \varphi_k^{(c)}(O, c_k)$ to denote the true de-biasing function using the true constraint or objective functions, and $\hat{\varphi}_j^{(b)}(O) \equiv \hat{\varphi}_j^{(b)}(O, \hat{b}_j)$ and $\hat{\varphi}_k^{(c)}(O) \equiv \hat{\varphi}_k^{(c)}(O, \hat{c}_k)$ to denote the estimated de-biasing function using the estimated constraint or objective function. Our analysis can also be extended to cross-fit estimates, which we use in practice to avoid losing sample size.
We assume that the estimates of the constraint and objective vectors and their de-biasing functions are consistent, and that the rate of convergence for de-biased estimates of linear functions of the constraint and objective vectors is controlled by a rate $r_n$.
In Appendix (ref) we give the forms of these de-biasing functions for the examples in Section (ref). Broadly, these de-biasing functions can be written in terms of weighted residuals of the form $\varphi_{j}^{(b)}(O; b_j) = g_{j}^{(b)}(O; b_j)(Y_j - b_j(X))$ and $\varphi_{k}^{(c)}(O; c_k) = g_{k}^{(c)}(O; c_k)(C_k - c_k(X))$ where $Y_j$ and $C_k$ are elements of $O$ and $g_{j}^{(b)}(\cdot)$, $g_{k}^{(c)}(\cdot)$ are corresponding Riesz representers (e.g. inverse propensity weights). To estimate the de-biasing functions, we can estimate the Riesz representers, giving a double robust or double machine learning-style estimate Robins1994, chernozhukov_doubledebiased_2018. In this case the rates will depend on the product of the errors in the estimates of the constraints $\bm{b}(\cdot)$ and objectives $\bm{c}(\cdot)$ with the errors in their respective Riesz representers---e.g. $r_n = \max\{\|\hat{\bm{b}} - \bm{b}\|_2\|\hat{\bm{g}}^{(b)} - \bm{g}^{(b)}\|_2, \|\hat{\bm{c}} - \bm{c}\|_2\|\hat{\bm{g}}^{(b)} - \bm{g}^{(b)}\|_2\}$ where $\hat{\bm{g}}^{(b)}$ and $\hat{\bm{g}}^{(c)}$ are estimates of the Riesz representers. Generally, such de-biasing functions can be constructed by deriving efficient influence functions for the expected constraint and objective vectors.
To begin, we will re-frame the conditional linear programs in Equation (ref) as discrete optimization problems. Because the matrix of marginal constraints $A \in \ensuremath{\mathbb{R}}^{J \times K}$ has rank equal to $J$, and the feasible set $\mathcal{P}(x)$ is non-empty by assumption, any solution to this LP must lie on the boundary of the feasible set $\mathcal{P}(x)$. In particular there exists a solution that is a basic feasible solution (BFS) with $J$ active, non-negative variables and $K - J$ inactive variables on the boundary at 0 klatt_limit_2022.
For a set of $J$ unique indices $B = \{i_1,\ldots, i_{J}\}$, let $A_B$ denote the sub-matrix of $A$ that includes all rows and the $J$ columns defined by $B$, and similarly for a vector $v$ let $v_B$ denote the sub-vector of $v$ that includes the $J$ elements defined by $B$. Then, a basic solution to the linear program is given by $\bm{p}_B = A_B^{-1} \bm{b}(x)$, and $p_i = 0$ for all $i \not \in B$. If $\bm{p}_B \geq 0$, then it is also a basic feasible solution. Abusing notation, we will use $A_B^{-1}$ to denote the re-ordered inverse matrix that places the rows of the inverse into $i_1,\ldots,i_J$, and sets the other rows equal to zero, so that $p = A_B^{-1}\bm{b}(x)$ is the BFS corresponding to the basis $B$.
We can then write the LPs in Equation (ref) as optimization problems over $\mathcal{B}$, the set of basic feasible solutions. Letting $\mathcal{B}^\ast_L(x) = \text{argmin}_{B \in \mathcal{B}} \langle\bm{c}(x), A_B^{-1} \bm{b}(x)\rangle$ and $\mathcal{B}^\ast_U(x) = \text{argmax}_{B \in \mathcal{B}} \langle\bm{c}(x), A_B^{-1} \bm{b}(x)\rangle$ denote the sets of optimal bases for the conditional minimization and maximization problems respectively, which may not be unique, we can write the bounds in Equation (ref) as $\theta_L = \ensuremath{\mathbb{E}}\left[\langle\bm{c}(x) ,A^{-1}_{B^\ast_L(X)} \bm{b}(x) \rangle\right]$ and $\theta_U = \ensuremath{\mathbb{E}}\left[\langle\bm{c}(x), A^{-1}_{B^\ast_U(X)} \bm{b}(x)\rangle\right]$, where $B^\ast_L(x) \in \mathcal{B}_L^\ast(x)$ and $B^\ast_U(x) \in \mathcal{B}_U^\ast(x)$ are any particular optimal bases.
If we knew these true optimal bases $B_L^\ast(x)$ and $B_U^\ast(x)$ for each value of the covariates $x$, then we could estimate the upper and lower bounds $\theta_L$ and $\theta_U$ using de-biased estimators, since they are expectations of bi-linear forms of the nuisance functions $\bm{b}(x)$ and $\bm{c}(x)$. However, this estimator is infeasible because we do not know the sets of optimal bases. Instead, we will use a plug-in approach, where we estimate the optimal bases using the estimates of the conditional constraints $\hat{\bm{b}}(x)$ and the conditional objective $\hat{\bm{c}}(x)$:
The corresponding plugin estimates of the basic feasible solutions are $\hat{\bm{p}}_L(x) \equiv A_{\widehat{B}_L(x)}^{-1} \hat{\bm{b}}(x)$ and $\hat{\bm{p}}_U(x) \equiv A_{\widehat{B}_U(x)}^{-1} \hat{\bm{b}}(x)$. We will proceed with estimating the lower and upper bounds $\theta_L$ and $\theta_U$ using de-biased BFS estimators as if the optimal bases were known:
Below, we will focus on characterizing the properties of the lower bound estimator; the upper bound is analogous.
We note that written this way, this estimator is a form of the covariate-assisted intersection bounds estimator proposed by semenova_aggregated_2024. However, even though there are potentially $\binom{K}{J}$ possible bases, we need not perform an exhaustive search over all possibilities. Indeed, we need not even enumerate nor analytically characterize the set of bases $\mathcal{B}$. Instead, we can use the simplex algorithm dantzig_maximization_1951 to determine an optimal basic feasible solution, which has polynomial run-time in typical cases and is efficient in practice, though it has exponential run-time in the worst-case spielman_smoothed_2004. This provides a substantial computational speedup over analytic vertex enumeration approaches that require enumerating all possible bases, which is particularly helpful because we must find the optimal bases for each unique value of the covariates $x$ in the sample, potentially up to $n$ times.
The quality of the plug-in estimates of the optimal bases will depend on the difficulty of picking out an optimal basis from the set of all possible bases. To quantify this, we will define the sub-optimality gap $\Delta_L(x)$ as the minimum difference between an optimal value of the conditional LP and the value of the LP at any other non-optimal vertex of the constraint set $\mathcal{P}(x)$. This can be written in terms of the non-optimal bases $\mathcal{B}\setminus \mathcal{B}_L^\ast$:
It may be that all feasible solutions are optimal and so $\mathcal{B}\setminus\mathcal{B}_L^\ast = \emptyset$, implying that the parameter of interest is point-identified at that value of the covariates. For such cases, any choice of basis will suffice, and we will define the sub-optimality gaps to be $\Delta_L(x) = 0$. We then make the following margin condition assumption about the distribution of the sub-optimality gaps across the covariate space $\mathcal{X}$.
The margin parameter $\alpha$ determines the difficulty of the problem. If $\alpha$ is small then there is a high probability across the covariate space $\mathcal{X}$ that there are suboptimal bases that are close to the optimal ones; if $\alpha$ is large then this is less likely. Margin conditions such as Assumption (ref) have been used to analyze the performance of plugin estimators in classification settings Audibert2007, estimators for classification in causal inference settings Luedtke2016,kennedy_sharp_2020, and estimating bounds on partially identified parameters as we do here levis_assisted_2023,DAdamo2023, semenova_aggregated_2024. Note that the definitions of the sub-optimality gap and the margin condition do not require that the optimal bases are unique. They also allow for the possibility that the parameter of interest is point-identified for some or all values of the covariates and so $0 < P(\Delta_L(X) = 0)$.
If the optimal bases were known, then the estimates $\hat{\theta}_L$ and $\hat{\theta}_U$ would be standard, de-biased estimators of the lower and upper bounds. However, since we do not know the optimal bases, we will need to account for the estimation error in their plug-in estimates. The following theorem gives the rate of convergence of the lower bound estimator $\hat{\theta}_L$; the properties of the upper bound are analogous.
Theorem (ref) shows that the estimation error $\hat{\theta}_L - \theta_L$ is equal to a mean-zero term $\tilde{\theta}_L - \ensuremath{\mathbb{E}}[\tilde{\theta}_L]$ plus additional bias terms that depend on the estimation error of the nuisance functions $\hat{\bm{b}}$, $\hat{\bm{c}}$ and the de-biased rate $r_n$. The first bias term is due to the sub-optimality of using plugin estimates of an optimal basis $\widehat{B}_L(x)$ relative to a true optimal basis $B_L^\ast(x)$. This decreases at a faster rate than the worst-case error of the nuisance functions by a factor of the margin parameter $\alpha$. Roughly, this is because in cases where the optimal bases are mis-classified, the sub-optimality gap must be smaller than the estimation error in the nuisance functions (note this can be extended to error in other norms as well). The next term comes from de-biasing, and the final term involves cross terms between the estimation errors of the constraint and objective vectors that arise from the bi-linear form of the estimator.
If the estimates of the nuisance functions converge quickly enough---though potentially slower than parametric rates---then the estimation error $\sqrt{n}(\hat{\theta}_L - \theta_L)$ is asymptotically equivalent to a centered version of the de-biased estimator using the plugin estimates of the optimal basis $\widehat{B}_L(x)$, $\sqrt{n}(\tilde{\theta}_L - \ensuremath{\mathbb{E}}[\tilde{\theta}_L])$. Therefore, if this is asymptotically normally distributed, then so is the estimator $\hat{\theta}_L$, allowing for straightforward asymptotic inference.
This result involves the variance of the plugin estimates of the optimal basis and requires that the estimates are sufficiently well-behaved to allow for $\tilde{\theta}_L$ to be asymptotically normal, which may be difficult to verify. Ideally, we would like to find an asymptotic expansion that relies only on population quantities, such as the set of true optimal bases $\mathcal{B}^\ast_L(x)$. The challenge is that there are potentially several such true optimal bases, and even though we can control the probability that the plugin estimates are sub-optimal, we cannot guarantee that the estimated basis corresponds to any particular true optimal basis. Indeed, it may be that the plugin estimate shifts between different optimal bases depending on the particular error structure of the nuisance functions.
If for (almost) all $x \in \mathcal{X}$ there is only one unique solution (i.e. $|\mathcal{B}^\ast_L(x)| = 1$), then the asymptotic expansion of $\hat{\theta}_L$ involves the true optimal basis. However, even if there are multiple optimal solutions to the conditional LP, if at least one solution is non-degenerate with basis values that are strictly greater than zero (i.e. $A_B^{-1}\bm{b}(x) > 0$), then the asymptotic expansion of $\hat{\theta}_L$ will be invariant to the choice of optimal basis if the objective vector $\bm{c}(x)$ is known. This is because a non-degenerate basis implies that the there is a unique solution to the Lagrangian dual of the conditional LP, and so the dual solution corresponding to an optimal basis, $A_{B_L^\ast(x)}^{-1\prime}\bm{c}(x)$, is the same regardless of which optimal basis is chosen.
Corollary (ref) slightly generalizes the above discussion to allow for cases where the objective vector is not known, but the de-biasing function is a linear function of the objective vector. This is the case, for instance, when the objective function is known up to a scaling factor such as the propensity score. Corollary (ref) also allows for cases where there is neither a unique nor a non-degenerate solution, but the variance of the de-biasing function is zero for such cases, a stringent condition that may be satisfied in some cases. Finally, note that Corollary (ref) explicitly allows for point-identified cases when the objective vector $\bm{c}(x)$ is known up to a scaling factor, as long as there is a non-degenerate solution. In the examples in Section (ref), non-degeneracy is a mild condition that will be satisfied if there is a non-zero probability of observing each outcome or treatment level at each value of the covariates $x$. It can also be possible to simplify the program to avoid degenerate solutions to remove such zero probabilities.
To construct confidence intervals, we estimate the variances with
for the lower bound, and analogously for the upper bound. Following imbens_confidence_2004, we can construct a Wald-type approximate $(1 - \alpha)$ confidence interval for the partially identified parameter $\theta$ by combining one-sided confidence intervals for the lower and upper bounds as $\left[\hat{\theta}_L - z_{1-\alpha} \sqrt{\frac{\hat{V}_L}{n}}, \hat{\theta}_U + z_{1-\alpha} \sqrt{\frac{\hat{V}_U}{n}}\right]$, where $z_{1-\alpha}$ is the $1-\alpha$ quantile of a standard normal distribution.
Due to the non-smoothness of the conditional linear programs, small errors in either the constraint or objective vectors can lead to large changes in the optimal solution of the conditional LP. As an alternative, we propose to target a smoothed version of the estimand that adds an entropy penalty to the conditional linear program. As we will see, this will induce some approximation error, but we will typically be able to tune the level of regularization so that the approximation error will be small relative to the standard error of the estimator.
For a vector $\bm{p} \in \ensuremath{\mathbb{R}}_+^K$ define the solutions to the entropic conditional linear programs as
We then define the entropic regularized lower and upper bounds as the expected value using the entropic solutions: $ \theta_L^\eta \equiv \ensuremath{\mathbb{E}}\left[\left\langle\bm{c}(X), \bm{p}_L^{\eta}(X)\right\rangle\right]$ and $ \theta_U^\eta \equiv \ensuremath{\mathbb{E}}\left[\left\langle\bm{c}(X), \bm{p}_U^{\eta}(X)\right\rangle\right]$.
The hyperparameter $\eta > 0$ controls the level of regularization. As $\eta \to 0$, the entropy penalty in Equation (ref) dominates the objective function, and so the entropic solutions $\bm{p}_L^\eta(x)$ and $\bm{p}_U^\eta(x)$ will converge to the maximum entropy distribution that satisfies the constraints $A\bm{p} = \bm{b}(x)$. In this case, the lower and upper bounds will be equal, $\theta_L^0 = \theta_U^0$, and correspond to point-identifying the estimand $\theta$ by making a maximum entropy assumption. In the potential outcomes setting from Section (ref), sending $\eta \to 0$ corresponds to identifying $\theta$ by assuming that the potential outcomes are mutually independent given the covariates. On the other hand, as $\eta \to \infty$, the entropic solutions will converge to the solutions of the unregularized conditional linear programs. In this way, $\eta$ can act as a sensitivity parameter, moving smoothly between the bounds under no additional assumptions and the point-identified case under the maximum entropy assumption.
There are several benefits of using the entropic regularized bounds in place of solving the linear program as in Section (ref). The computational complexity required to solve the conditional LPs can scale unfavorably with aspects of the problem, such as the number of outcome or treatment levels. Because the conditional LPs need to be solved for each data point, this can lead to a large computational burden if the number of outcome levels or especially the number of treatment levels is large in potential outcome settings. In certain problems such as cases with a binary treatment, entropic regularization can lead to significant computational speedups by using the sinkhorn_diagonal_1967 algorithm cuturi_sinkhorn_2013.
Furthermore, the regularizer is strongly convex, so there exists a unique solution to both the lower and upper entropic bounds for all values of the covariates $x \in \mathcal{X}$, and the bounds are smooth functions of the nuisance parameters $\bm{b}(x)$ and $\bm{c}(x)$. In particular, by taking the Lagrangian dual of the conditional entropic linear programs (ref), for any constraint and objective vector $\bm{b} \in \ensuremath{\mathbb{R}}^J$ and $\bm{c} \in \ensuremath{\mathbb{R}}^K$, we can write the solutions as $\bm{p}^\eta_L(\bm{b}, \bm{c}) = \exp\left(-A'\bm{\lambda}_L^\eta(\bm{b}, \bm{c}) + \eta \bm{c}\right)$ and $\bm{p}^\eta_U(\bm{b}, \bm{c}) = \exp\left(A'\bm{\lambda}_U^\eta(\bm{b}, \bm{c}) + \eta \bm{c}\right)$, where the dual variables $\bm{\lambda}_L^\eta(\bm{b}, \bm{c})$ and $\bm{\lambda}_U^\eta(\bm{b}, \bm{c})$ are the solutions to the entropic regularized dual problems:
where $A_{\cdot k}$ denotes the $k$\textsuperscript{th} column of the matrix $A$.
This removes the ambiguity of multiple solutions, and allows us to construct a de-biased estimator using standard semi-parametric theory because the dual variables are an implicitly defined function of the nuisance parameters $\bm{b}(x)$ and $\bm{c}(x)$, sidestepping the issue of selecting an optimal basis that arises in the unregularized linear programs. To do so, we can compute the Jacobian of the entropic regularized solutions with respect to the vectors $\bm{b}(x)$ and $\bm{c}(x)$ and use a Taylor expansion of the true solution around the plugin estimates (i.e. $\hat{\bm{p}}_L^\eta(x) \equiv \bm{p}_L^\eta(\hat{\bm{b}}(x), \hat{\bm{c}}(x))$ and $\hat{\bm{p}}_U^\eta(x) \equiv \bm{p}_U^\eta(\hat{\bm{b}}(x), \hat{\bm{c}}(x))$). Then, we can approximately de-bias the plugin estimates as
where $\nabla_{\bm{b}} \hat{\bm{p}}^\eta_L(x) \equiv \nabla_{\bm{b}}\bm{p}^\eta_L(\hat{\bm{b}}(x), \hat{\bm{c}}(x))$ is shorthand for the Jacobian of the entropic regularized lower bound with respect to the constraint vector $\bm{b}$ evaluated at the plugin estimates of the nuisance functions $\hat{\bm{b}}(x)$ and $\hat{\bm{c}}(x)$; $\nabla_{\bm{c}} \hat{\bm{p}}^\eta_L(x), \nabla_{\bm{b}} \hat{\bm{p}}^\eta_U(x), \nabla_{\bm{c}} \hat{\bm{p}}^\eta_U(x) $ are similarly defined. Lemma (ref) in the Appendix gives the explicit form of these Jacobians. These estimators are similar to the de-biased BFS estimators in Section (ref), but now rather than relying on plugin estimates of the optimal basic feasible solutions, we can use the differentiability of the entropic solutions to make the estimator more robust to errors in the nuisance function estimates. Again, we focus on the properties of the lower bound estimator for brevity, but the upper bound is analogous.
In Appendix Theorem (ref), we show that the asymptotic bias of $\hat{\theta}_L^\eta$ relative to the true entropic regularized bound $\theta_L^\eta$ behaves as we would expect, decreasing with the de-biased rate $r_{n}$ and second order errors in the nuisance functions $\hat{\bm{b}}$ and $\hat{\bm{c}}$, and so there is no first-order dependence on the nuisance functions. Of course, changing the estimand by including the entropy penalty will introduce approximation error relative to the original bounds because $\theta_L^\eta \geq \theta_L$. However, the approximation error induced via entropic regularization has been shown to decrease exponentially with the hyperparameter $\eta$ when $\eta$ is large enough relative to the sub-optimality gap weed_explicit_2018. The following result incorporates the approximation error in a setting where we choose the hyperparameter $\eta$ to grow with the sample size.
Theorem (ref) includes the approximation error induced by regularization. The key requirement is that the hyperparameter $\eta$ is large---and so the level of regularization is small---relative to the ratio of the size of the feasible set $\mathcal{P}(X_i)$ and the sub-optimality gap $\Delta_L(X_i)$ for each observed covariate $X_i$ in the sample for which the sub-optimality gap is not zero (if it is zero, then any feasible solution is optimal). This ensures that the overall level of regularization $\frac{1}{\eta}$ is small enough so that the approximation error decreases exponentially with $\eta$ for each observed $X_i$ in the sample weed_explicit_2018. Because these aspects of the conditional linear programs are independent of the sample size, we can expect that to find an $\eta$ that is larger than this minimal value for a finite sample.
On the other hand, the first and second derivatives of the entropic solutions with respect to the nuisance functions $\bm{b}(x)$ and $\bm{c}(x)$ scale with $\eta$: larger values of $\eta$ can increase the magnitude of the derivatives and increase the scale of the bias quadratically with $\eta$. Theorem (ref) summarizes sufficient conditions for this increase in magnitude to be outweighed by the decrease in approximation error so that the overall bias converges faster than $n^{-1/2}$---and so is asymptotically smaller than the standard error and will not affect our confidence intervals. The faster $\eta$ grows, the more quickly the nuisance function errors need to converge (the maximum being a parametric convergence rate when $\eta = O(n^{1/4})$), but it is sufficient for $\eta$ to scale logarithmically with the sample size, $\eta = O(\log n)$, a relatively slow growth.
Theorem (ref) relies on being able to choose a large enough $\eta$ to ensure that the approximation error decays exponentially. However, depending on the tail behavior of the sub-optimality gaps $\Delta_L(X_i)$, it may be possible that under some data generating processes the minimal size of $\eta$ grows too quickly and the magnitude of the derivatives outweighs the decrease in estimation error. While it is not clear under which data generating processes this will occur, in Appendix Corollary (ref) we consider a more pessimistic case where we have no guarantee that $\eta$ is large enough to have a small bias. In this case, the approximation error will not decay exponentially, but rather polynomially with $\eta$: if $\eta$ is too small, then the approximation error may only decay with the inverse of $\eta$ weed_explicit_2018, but under the margin condition in Assumption (ref), this occurs with probability $O(\eta^{-\alpha})$, where $\alpha$ is the margin parameter, leading to an approximation error that decays with $\eta^{-(1+\alpha)}$. This affects the minimum rate of growth for $\eta$ and the required convergence rates of the nuisance functions. For example, if we choose $\eta = O(n^\beta)$ for some $0 < \beta < \frac{1}{4}$, then we need $\beta > \frac{1}{2(1+\alpha)}$. If $\alpha$ is just larger than 1, then this would require that the nuisance functions converge at a parametric rate, but for larger values of $\alpha$ this would decrease. In Appendix (ref) we inspect the impact of different growth rates for $\eta$ via simulation and find that the de-biased entropic estimator is not particularly sensitive to the choice of hyper-parameter. As a baseline heuristic, one can first compute the de-biased BFS estimate, then choose $\eta$ to be large enough so that the difference between the de-biased BFS and entropic estimates is not substantively meaningful (e.g. an order of magnitude or two smaller than the estimates).
Finally, we can construct confidence intervals for the entropic bounds by first estimating the variances with \[
\] then constructing $1-\alpha$ level confidence intervals for the entropic bounds as $\hat{\theta}_L^\eta - z_{1-\alpha} \sqrt{\frac{\hat{V}_L^\eta}{n}}$ and $\hat{\theta}_U^\eta + z_{1-\alpha} \sqrt{\frac{\hat{V}_U^\eta}{n}}$.
We now turn to applying the de-biased BFS and entropic estimators to the problem of decision making, where the goal is to find a decision rule that optimizes a parameter that is partially identified via conditional linear programs. In particular, consider a decision rule $\pi: \mathcal{X} \to [0,1]$ that maps covariates to a continuous action between 0 and 1 (e.g. a treatment probability). We will focus on settings where the covariate-conditional objective function is given by $\bm{c}(x, \pi(x))$, so that the conditional lower bound is $\theta_L(x, \pi) = \min_{\bm{p} \in \mathcal{P}(x)} \langle\bm{c}(x, \pi(x)), \bm{p}\rangle,$ and the overall expected lower bound on the objective value for a given decision rule $\pi$ is given by $\theta_L(\pi) = \ensuremath{\mathbb{E}}[\theta_L(X, \pi)]$.
Our goal is to find a decision rule that maximizes this lower bound across all decision rules $\pi$ in a policy class $\Pi$, i.e. $\pi^\ast \in \text{argmax}_{\pi \in \Pi} \ \theta_L(\pi)$. This setup accommodates a wide variety of decision making problems under uncertainty Manski2011 including maximin value and minimax regret rules (i.e. by maximizing the minimum negative regret). In this section we will consider finding such policies by maximizing an estimate of the lower bound $\hat{\theta}_L(\pi)$, using either the the de-biased BFS or entropic estimators from the previous section.
Throughout, we will focus on settings where the objective function $\bm{c}(x, \pi(x))$ is a known function of the covariates $x$ and the decision rule $\pi(x)$ and we will make several regularity assumptions.
These assumptions ensure that the objective is bounded and Lipschitz continuous in the decision rule $\pi(x)$. As a final piece of setup, we will characterize the complexity of the policy class $\Pi$ via the Rademacher complexity $\mathcal{R}_n(\Pi) = \ensuremath{\mathbb{E}}\left[\sup_{\pi \in \Pi} \left|\frac{1}{n}\sum_{i=1}^n \pi(X_i) \varepsilon_i\right|\right],$ where $\varepsilon_i \in \{-1,1\}$ are independent Rademacher random variables.
For a policy $\pi$, let $\mathcal{B}_L(x, \pi(x))$ denote the set of optimal bases given the objective vector $\bm{c}(x, \pi(x))$ and let $\widehat{B}_L(x, \pi(x))$ be the plugin estimate of an optimal basis. Then, we estimate a policy $\hat{\pi}$ by maximizing the de-biased BFS estimator of the lower bound on the expected objective value. Denoting $\hat{\bm{p}}(x, \pi(x)) \equiv A^{-1}_{\widehat{B}_L(x, \pi(x))} \hat{\bm{b}}(x)$, the minimax optimal policy $\hat{\pi}$ is given by
To evaluate the quality of the estimated policy $\hat{\pi}$, we compare it to the best-in-class policy that maximizes the true lower bound on the expected objective value, i.e. $\pi^\ast \in \text{argmax}_{\pi \in \Pi} \theta_L(\pi)$, using the excess regret, $\theta_L(\pi^\ast) - \theta_L(\hat{\pi})$. Controlling this excess regret involves (i) linking the excess regret to the bias in the de-biased BFS estimator, and (ii) applying standard results on empirical risk minimization.
To do so, we first make a stronger version of the margin condition in Assumption (ref).
Assumption (ref) is primarily a requirement on the form of the objective function $\bm{c}(x, \pi(x))$. It restricts the number of cases where $\bm{c}(x, \cdot)$ is such that it is possible to choose a policy value $\pi(x)$ with a small sub-optimality gap $\Delta_L(x, \pi(x))$. If the optimal BFS does not depend on the policy value $\pi(x)$ (e.g. if the objective function is linear in $\pi(x)$), then this “strong” margin condition is is not stronger than the margin condition in Assumption (ref).
Theorem (ref) shows that the excess regret of the estimated policy $\hat{\pi}$ inherits most of the properties of the de-biased BFS estimator in Theorem (ref). The error in the estimated constraint vector $\hat{\bm{b}}(x)$ again enters, with the impact of the error lessened by the margin parameter $\alpha$, and the de-biased rate $r_{n}$ enters as well. In this case the objective vector is known, and so there are no estimation errors to consider. There is an additional term in the excess regret that depends on the Rademacher complexity of the policy class $\Pi$, which is standard for empirical risk minimization problems. In settings where it would be possible to estimate the lower bound $\theta_L(\pi)$ using standard, asymptoic normality-based inference as in Corollary (ref), the excess regret would be primarily impacted by the Rademacher complexity of the policy class $\Pi$, a quantity that is controllable by the analyst by choosing a simpler or more expressive policy class.
Just as with the de-biased BFS estimator, we can use the entropic estimator to estimate the lower bound on the expected objective value for a policy $\pi$, and then find the policy that maximizes this estimate. Fixing a regularization hyperparameter $\eta$, and denoting $\hat{\bm{p}}_L^\eta(x, \pi(x)) \equiv \bm{p}_L^\eta(\hat{\bm{b}}(x), \bm{c}(x, \pi(x)))$ as the solution to the entropic conditional linear program with policy value $\pi(x)$, we can find a policy $\hat{\pi}^\eta$ by solving
It is more computationally straightforward to solve this optimization problem than for the de-biased BFS estimator in (ref). If the optimal BFS depends on the policy value $\pi(x)$, optimizing the de-biased BFS estimator in (ref) would require using gradient-free optimization methods. In contrast, the entropic estimator of the lower bound $\hat{\theta}_L^\eta(\pi)$ is differentiable with respect to the policy value $\pi(x)$ and so it can be optimized using standard gradient-based methods---though the overall objective is non-convex. Note, however, that each objective value and gradient computation requires solving the dual to the entropic conditional linear program for each unique value of the covariates $x$ in the sample, which can be computationally expensive for iterative optimization methods.
We will compare this estimated decision policy $\hat{\pi}^\eta$ to two notions of an optimal policy: (i) the best-in-class policy that maximizes the true lower bound on the expected objective value, $\pi^\ast \in \text{argmax}_{\pi \in \Pi} \ \theta_L(\pi)$, and (ii) the best-in-class policy that maximizes the entropic lower bound on the expected objective value, $\pi^\ast_\eta \in \text{argmax}_{\pi \in \Pi} \ \theta_L^\eta(\pi)$. In each case we will measure the excess regret relative to the corresponding bound.
Theorem (ref) shows that relative to the optimal decision policy for the true entropic lower bound, $\pi^\ast_\eta$, the excess regret of the estimated policy $\hat{\pi}^\eta$ is impacted by estimation errors in the constraint vector $\hat{\bm{b}}(x)$ only via the de-biased rate $r_{n}$ and squared errors $\|\hat{\bm{b}} - \bm{b}\|_2^2$. This means that the estimated policy is relatively robust to such errors and we would typically expect the complexity of the policy class $\Pi$ to drive the excess regret. The same is true for the excess regret relative to the optimal policy for the true lower bound, $\pi^\ast$, along with the additional scaling with the regularization hyperparameter $\eta$ as with Theorem (ref). The excess regret bound also requires that $\eta$ be chosen sufficiently large across all potential policies $\pi$ in the policy class $\Pi$. As in Section (ref), we can also extend these results to make no assumptions on $\eta$ but require the strong margin condition in Assumption (ref) to hold, in which case the approximation error would decay like $\eta^{-1(1 + \alpha)}$.
In this section we apply the general procedure outlined in the previous sections to measure the impact of Medicaid access and enrollment on emergency department visits finkelstein_oregon_2012, taubman_medicaid_2014. In Appendix (ref) we also conduct a simulation study to understand the finite sample behavior of the estimators.
This natural experiment occurred in 2008 when the state of Oregon offered a Medicaid program to individuals by randomly selecting names off a waiting list. Here, being chosen off the waitlist is an instrument $Z$, while enrolling in Medicaid is the treatment $D$. We will focus on the number of emergency department (ED) visits as the outcome $Y$, it is discrete with $L=23$ levels corresponding to the number of ED visits in the past 6 months, ranging from 0 to 22 visits (the maximum censored value). Our dataset consists of the 13,019 individuals who lived in areas served by one of 12 hospitals for which the authors received ED records, were sampled to be part of an in-person follow-up survey, and had no other members in the household. For this group, the probability of being selected off the waitlist was uniform across individuals (with 48% selected). We include basic demographic variables, pre-period measures of health and health-care utilization, and pre-period enrollment in public assistance programs in our covariate vector $X$. Here, strong ignorability of $Z$ is satisfied by design; we also assume the exclusion restriction that being selected from the waitlist does not affect the outcome except through Medicaid enrollment.
\paragraph{Impacts of Medicaid access.} First, we focus on the impacts of access to Medicaid by being chosen off the waitlist ($Z = 1$) versus not ($Z = 0$). We compare the observed randomization procedure to an oracle policy that could observe each individual's potential outcomes (i.e. $Y(D(0))$ and $Y(D(1))$) and can select individuals from the waitlist in order to minimize the number of ED visits (i.e. when $Y(D(1)) \leq Y(D(0))$). We do this by estimating bounds on (i) the number of additional ED visits under the observed distribution relative to the counterfactually optimal policy, $E[Y] - \ensuremath{\mathbb{E}}[\min_z Y(D(z))]$, (ii) the proportion of individuals who experience any additional ED visits under the observed distribution relative to the oracle optimal policy---$P(Y \neq \min_{z} Y(D(z)))$, and (iii) the regret under the powerlaw value function from Section (ref), ${V^\lambda(\ensuremath{\mathbbm{1}}\{Y(D(1)) \leq Y(D(0))\}) - V^\lambda(P(Z = 1 \mid X))}$. Because we would like to minimize the number of ED visits, we set the utility of $Y$ ED visits to be the maximum number of ED visits minus $Y$, $U(Y) = 22 - Y$. We estimate the nuisance functions using multinomial gradient boosted decision trees and use 3-fold cross-fitting. To make the hyperparameter $\eta$ comparable across estimands, we scale the objective of each linear program to have a maximum value of 1 i.e. $\|\bm{c}\|_\infty = 1$.
Figure (ref) shows the estimated lower and upper bounds on the first two estimands using the de-biased entropic estimator with varying levels of entropic regularization $\eta$. Setting $\eta$ to be small assumes that the potential outcomes are nearly independent, and we estimate bounds narrowly centered around an additional 0.75 ED visits per person under the randomization policy vs the ED-minimizing waitlist policy (relative to an observed average of $1.04 \pm 0.04$ visits per person) and 25% of individuals experiencing any more ED visits relative to the oracle. Decreasing the level of regularization, the estimates stabilize at between $[0.02, 1]$ additional ED visits per person (95% CI $[0, 1.04]$) and between [0%, 34%] of individuals experiencing any more ED visits (95% CI $[0\%, 34.4\%]$).
Figure (ref) shows the estimated regret of the randomization policy versus the ED-minimizing waitlist policy using the power-law collective utility function as the power parameter $\lambda$ changes. We see that the upper bound on the regret is larger the closer the utility function is to the expected utility ($\lambda \to 1$), and it decreases as $\lambda$ decreases to $-1$ and is more inequality-averse. This indicates that randomizing incurs less potential regret for more inequality-averse the utility functions. These bounds do not rule out that randomly selecting individuals from the waitlist led to higher than the minimal number of ED visits. However, they do indicate that there is potentially substantial room for improvement on these policies at moderate to low levels of inequality aversion.
These ranges also mask considerable variation in the level of uncertainty. Figure (ref) shows the estimated bounds using the least amount of regularization, stratified by education and whether the individual had any pre-randomization ED visits. We see that for individuals that had no ED visits in the pre-period, the upper bound on the expected number of additional ED visits is much smaller than for individuals that had at least one visit; a similar trend holds for the proportion of individuals with any additional visits. There is limited heterogeneity by level of education. In Appendix (ref) we consider estimating bounds conditional on continuous measures using a pseudo-outcome regression approach.
Taken together, these results imply that there is potential for an optimized waitlist rule. To evaluate this, we consider a stylized policy learning problem where the goal is to reduce ED visits (this is not the only goal of Medicaid). We learn a minimax regret waitlist allocation rule relative to the oracle policy under power-law collective utility functions using the de-biased entropic estimator of the regret relative to the oracle with $\eta = 100$. We parameterize the rule as a logistic function of the number of pre-period ED visits, and find the minimax optimal parameters by solving Equation (ref). Figure (ref) shows the estimated minimax allocation rules using a range of power-law collective utility functions. With $\lambda = 1$, the regret is minimized by minimizing the conditional average effect of $Z$ on $Y$; this selects nobody off the waitlist because the average effect is negative across the range of pre-period ED visits. As $\lambda$ decreases, the minimax regret rule is no longer point-identifiable, and the estimated minimax regret waitlist policies assign a higher selection probability to individuals with a higher number of pre-period ED visits, with the slope flattening for more inequality-averse utility functions.
\paragraph{Impacts of Medicaid enrollment.}
Next, we consider the same question but for enrollment in Medicaid $D$ rather than selection from the waitlist $Z$. In particular, if there were no capacity constraints, and we could enroll individuals into Medicaid directly based off their potential outcomes ($Y(0)$ and $Y(1)$), how would this compare to simply taking everyone off the waitlist (i.e. $Y(D(1))$)? We make this comparison in two ways: (i) the number of additional ED visits $E[Y(D(1))] - \ensuremath{\mathbb{E}}[\min_d Y(d)]$, and the proportion of individuals with additional ED visits under the no-waitlist policy, $P(Y(D(1)) \neq \min_{d} Y(d))$. Here we use the de-biased BFS estimator with nuisance functions again fit using multinomial gradient boosted decision trees and 3-fold cross-fitting. For $14\%$ of units in the sample, the LP was not feasible, implying the exclusion restriction is not possible with those estimated values. We exclude these units from the analysis.
Overall, we estimate that if everyone were to be taken off the waitlist there would be an additional $[0.06, 1.23]$ ED visits per person relative to the oracle policy that enrolls individuals into Medicaid based on their potential outcomes (95% CI $[0.04, 1.28]$), and between $[0\%, 48.3\%]$ of individuals would experience any more ED visits relative to the oracle policy (95% CI $[0\%, 49\%]$). Again there is substantial heterogeneity in these upper bounds. Figure (ref) shows the estimated bounds stratified by education and whether the individual had an ED visits in the 14 months pre-randomization. As with selection off the waitlist, we see that the bounds are much smaller for individuals with no ED visits.
In this paper we proposed a general framework for estimating bounds on partially-identified parameters that satisfy a series of linear constraints, conditional on auxiliary covariates, leading to a series of linear programs that are conditional on the covariates. We developed two estimators for these bounds. The de-biased BFS estimator directly solves the conditional linear programs for each unique covariate value in the data, and then extracts readily available information from standard LP solvers to construct de-biased estimates of the bounds. This estimator is related to existing methods for covariate-assisted bounds, but eschews the need to analytically derive the bound or fully enumerate a combinatorial number of potential solutions.
A key issue with solving the linear programs directly is that the solutions are not smooth in the constraint and objective vectors, so small errors in these nuisance parameters can lead to large changes in the estimated bounds. The entropic regularized estimator instead targets a smoothed version of the bounds, including an entropy penalty to the conditional linear programs that ensures that the solutions are smooth in the nuisance parameters and so we can estimate them in a de-biased manner. Although this induces approximation error, we can likely expect the approximation error to be small relative to the standard error of the estimates if we allow the level of regularization to shrink appropriately. Finally, we extend these results to decision making problems where the value of a decision policy is only partially identified by a set of conditional linear constraints. We show that taking an empirical risk minimization approach and optimizing bounds estimated via the de-biased BFS or entropic estimators results in a policy with excess regret that has an analogous structure to the bias of the estimators.
An important limitation of these approaches is that they require a finite number of decision variables and so for instance cannot accommodate continuous outcomes. While it may be possible to approximate some continuous outcomes with a finite number of discrete levels, this may not be feasible in all cases. One potential path forward is to leverage the duality of the linear programs to adapt the de-biased BFS and entropic estimators to return valid bounds even if an approximation is used, similar in spirit to the dual estimator proposed by ji_model-agnostic_2024. In addition, while we have used an entropy penalty to regularize the LPs, other forms of regularization may have different benefits. For example, it may be possible to use barrier-based interior point methods to solve the LPs, and directly account for regularization bias via a double de-biased estimator following recent results in liu_beyond_2025. Such an approach could lessen the need to select a regularization hyperparameter. Finally, another limitation is that we have sidestepped the question whether and when conditioning on covariates information leads to narrower bounds.