The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
90,645 characters
Robust Design and Evaluation of Predictive Algorithms under Unobserved Confounding
{\singlespacing
\maketitle
\thispagestyle{empty}
\setcounter{page}{0}
\begin{abstract}
Predictive algorithms inform consequential decisions in settings with selective labels: outcomes are observed only for units selected by past decision makers.
This creates an identification problem under unobserved confounding --- when selected and unselected units differ in unobserved ways that affect outcomes.
We propose a framework for robust design and evaluation of predictive algorithms that bounds how much outcomes may differ between selected and unselected units with the same observed characteristics.
These bounds formalize common empirical strategies including proxy outcomes and instrumental variables.
Our estimators work across bounding strategies and performance measures such as conditional likelihoods, mean square error, and true/false positive rates.
Using administrative data from a large Australian financial institution, we show that varying confounding assumptions substantially affects credit risk predictions and fairness evaluations across income groups.
\medskip
\noindent \textbf{Keywords}: predictive algorithms, missing data, confounding, partial identification. \\
\noindent \textbf{JEL Codes}: C14, C52, G20
\end{abstract}
}
\newpage \clearpage
\section{Introduction}
Predictive algorithms inform high-stakes decisions in pretrial release, consumer lending, health care, and many other domains.
Designing and evaluating such tools is often complicated by the ``selective labels'' problem: outcomes are observed only for units selected by past decision makers.
For example, pretrial risk assessments predict the likelihood defendants would commit pretrial misconduct if they were released before trial; but we only observe whether a defendant committed misconduct if a judge decided to release them \citep[e.g.,][]{KLLLM(18), Rambachan(21), AngelovaEtAl(22)}.
Consumer credit scores predict the likelihood of default if an applicant were granted a loan; but we only observe whether an applicant defaulted if the financial institution approved them and they accepted the loan terms \citep[e.g.,][]{BlattnerNelson(21), CostonEtAl(21)-Rashomon}.\footnote{The selective labels problem arises when analyzing predictive algorithms in many other empirical settings as well, such as medical diagnosis \citep[e.g.,][]{MullainathanObermeyer(21)}, child protective services \citep[e.g.,][]{ChouldechovaEtAl(18)}, hiring \citep[e.g.,][]{LiRaymondBergman(20)}, education \citep[][]{BergmanEtAl(21)}, and tax audits \citep[e.g.,][]{BlackEtAl(22)-TaxAudits, BattagliniEtAl(22), ElzaynEtAl(23)}.}
The selective labels problem creates an identification challenge under unobserved confounding--when selected and unselected units differ in ways not captured by observed covariates.
Judges may release defendants based on courtroom interactions not in administrative data; loan applicants may accept offers based on competing credit options unavailable to researchers.
Ignoring unobserved confounding leads to inaccurate predictions and misleading evaluations of predictive algorithms.\footnote{While there exists methods for designing and evaluating predictive algorithms in computer science from selectively observed data \citep[e.g.,][]{Coston(20)-Evaluation, CostonEtAl(21)-Rashomon, MishlerEtAl(21), mishler2021fade, guerdan2023counterfactual}, they tackle this challenge assuming that the historical selection decision was unconfounded conditional on observed covariates.}
Existing solutions are limited.
Applied researchers rely on informal heuristics---imputing missing outcomes using proxy variables \citep[][]{ChouldechovaEtAl(18), BlattnerNelson(21), MullainathanObermeyer(21)} or applying ad-hoc adjustments to observed outcomes among selected units.
Recent work uses instrumental variables but typically focuses on specific estimands \citep[][]{lakkaraju2017selective, ArnoldDobbieHull(20)-discrimination_alg} or predictive algorithms that are linear in covariates \citep[][]{ArnoldEtAl(24)}.
We develop a framework for robust design and evaluation of predictive algorithms under unobserved confounding.
We bound how much outcomes may differ between selected and unselected units with the same observed characteristics.
These bounds formalize common empirical strategies --- including proxy outcomes and instrumental variables --- and partially identify key predictive estimands such as conditional likelihoods, mean square error, calibration, and true/false positive rates.
Our estimation approach works across multiple bounding strategies and performance measures, flexibly incorporating machine learning methods.
We consider data $(X_i, D_i, Y_i)$ for $i = 1, \hdots, n$ drawn i.i.d. from a joint distribution $\mathbb{P}(\cdot)$, where $X_i$ is a vector of covariates and $D_i \in \{0, 1\}$ is a selection decision by past decision makers (e.g., judges or loan applicants).
The observed outcome satisfies $Y_i = D_i \times Y_i^*$ for some true outcome $Y_i^* \in \{0, 1\}$; for example, $Y_i^*$ is whether a defendant would commit pretrial misconduct if released before trial or whether an applicant would default if granted a loan.
Researchers have two objectives.
First, design \textit{predictive algorithms} $s(X_i)$ that predict the probability $Y_i^* = 1$ given the covariates $X_i$ and construct data-driven decision rules.
Second, evaluate existing predictive algorithms through performance measures, such as \textit{overall performance} given by $\mbox{perf}(s; \beta) = \mathbb{E}[\beta_0(X_i; s) + \beta_1(X_i; s) Y_i^*]$,
and \textit{class-specific performance} given by $\mbox{perf}_{+}(s; \beta) = \mathbb{E}[\beta_0(X_i; s) \mid Y_i^* = 1]$ and $\mbox{perf}_{-}(s; \beta) = \mathbb{E}[\beta_0(X_i; s) \mid Y_i^* = 0]$, where $\beta_0(X_i; s), \beta_1(X_i; s)$ are researcher-specified functions.
Alternative choices recover popular diagnostics such as mean square error, true positive rate, and false positive rate.
Comparing these measures across sensitive attributes evaluates fairness properties \citep[e.g.,][]{MitchellEtAl(19), barocas-hardt-narayanan}.
Rather than assuming unconfounded selection, we bound the confounding function $\delta(x) := \mathbb{P}(Y_i^* = 1 \mid D_i = 0, X_i = x) - \mathbb{P}(Y_i^* = 1 \mid D_i = 1, X_i = x)$ using researcher-specified bounding functions that may depend on covariates and identified nuisance parameters.
The confounding function measures how much unselected units differ from selected units with the same characteristics; bounding it formalizes assumptions about unobserved confounding severity without requiring a full selection model.
Under these bounds, the conditional likelihood and predictive performance measures are partially identified.
Our framework encompasses popular strategies for addressing selective labels as special cases.
When researchers observe proxy outcomes correlated with $Y_i^*$ --- such as default on other credit products \citep{BlattnerNelson(21)} or long-term health outcomes \citep{MullainathanObermeyer(21), ChanGentzkowYu(21)} --- this implies specific bounds.
When researchers have instrumental variables like randomly assigned decision makers \citep{lakkaraju2017selective, KLLLM(18)}, this implies alternative bounds.
Our framework accommodates other plausible identifying assumptions.
We apply this framework to two objectives.
First, we develop estimators for the bounds on the conditional likelihood $P(Y_i^* = 1 \mid X_i = x)$ (e.g., the pretrial misconduct rate or default rate given covariates), enabling design of predictive algorithms robust to unobserved confounding.
Building on advances in causal machine learning \citep[][]{NieWager(20), kennedy2020optimal}, we propose pseudo-outcome regression estimators with sample splitting: estimate nuisance parameters on one fold, construct influence-function-based pseudo-outcomes, then regress them on covariates using machine learning on a second fold.
A key innovation is an oracle inequality extending \citet{kennedy2020optimal}'s pointwise analysis to integrated mean square error: our feasible estimators approximate an infeasible oracle up to a smoothed, doubly-robust product of first-stage errors.
This applies to any regression method satisfying a mild stability condition, including random forests, kernel methods, and series estimators.
We leverage this to characterize integrated mean square error convergence rates and establish regret bounds for plug-in decision rules, connecting estimation guarantees to decision-theoretic performance.
Second, we develop estimators for bounds on the performance of a given predictive algorithm $s(\cdot)$ under unobserved confounding.
For overall performance measures $\mbox{perf}(s; \beta)$ --- which nest mean square error, calibration, and precision --- we construct debiased estimators using crossfitting and debiased machine learning \citep[e.g.,][]{RobinsEtAl(08), ZhengvanderLaan(11), ChernozhukovEtAl(18)-DML}.
Under mild regularity conditions, these estimators are $\sqrt{n}$-consistent and jointly asymptotically normal.
For class-specific performance $\mbox{perf}_{+}(s; \beta)$ and $\mbox{perf}_{-}(s; \beta)$ --- such as true and false positive rates --- we show sharp bounds correspond to optimal values of linear fractional programs with nuisances in both the objective and constraints.
We estimate using two folds: estimate nuisances and bounding functions on the first fold, then solve the sample optimization program on the second fold.
We establish partial double robustness: estimation error decomposes into sampling error, a doubly-robust product for objective nuisances, and constraint estimation errors entering through root mean square error.
Proving this exploits the optimization structure of linear fractional programs.\footnote{Since linear fractional programs arise in several related sensitivity analysis frameworks in causal inference \citep[e.g.,][]{AronowLee(2013), MiratrixWagerZubizaretta(18), KallusMaoZhou(18)-IntervalEstimation, ZhaoSmallBhattacharya(19), KallusZhou(21)}, our results may be broadly applicable beyond the selective labels setting.}
Our estimation strategies accommodate multiple bounding strategies, flexibly incorporating domain-specific knowledge while maintaining statistical guarantees.
As an empirical illustration, we apply our framework to administrative data on personal loan applications from a large Australian bank.
First, we construct confounding-robust credit risk scores under varying confounding assumptions.
The resulting algorithms produce substantially different risk rankings that could meaningfully affect credit access.
Second, inspired by recent work in consumer finance \citep[][]{BlattnerNelson(21), FusterEtAl(22)predictably, DiMaggioEtAl(22)-Fintech}, we evaluate a benchmark credit score's performance across income groups.
An analysis ignoring unobserved confounding suggests the score is more accurate (in mean square error) for higher-income applicants.
However, we show this disparity disappears at plausible levels of unobserved confounding, illustrating how researchers can assess the robustness of fairness conclusions.
Our framework prioritizes breadth and applicability. Rather than bespoke methods for each combination of bounding strategy and performance measure, we provide a unified approach for a variety of assumptions and performance measures.
We deliver practical estimation procedures that enable researchers to robustly design and evaluate predictive algorithms under unobserved confounding.
\vspace{-1em}
\paragraph{Related Work:}
Our analysis builds on foundational work on partial identification with selective labels \citep[e.g.,][]{Manski(89), Manski(90), Manski(94)-SampleSelectionBounds}.
We generalize \citet[][]{Manski(03)}'s ``approximate mean independence'' assumption by incorporating covariates and nuisance parameters into the bounding functions, enabling flexible formalization of empirical heuristics using proxy outcomes and other strategies. See also \cite{diaz2013sensitivity, luedtke2015statistics, diaz2018sensitivity} for similar bounding assumptions without covariates.
We characterize identified sets for novel estimands arising in predictive algorithm evaluation.
Instrumental variable strategies are widely used to address unobserved confounding in this setting. \citet{lakkaraju2017selective} and \citet{KLLLM(18)} evaluate specific performance measures by returning a single point in the identified set.
\cite{ArnoldDobbieHull(20)-discrimination_alg, ArnoldEtAl(24)} extrapolate across instrument values using identification-at-infinity arguments. We instead derive covariate-dependent bounds following \cite{Manski(94)-SampleSelectionBounds}, accommodating IVs as a special case while also encompassing other identifying strategies.\footnote{In related work, \cite{LevisEtAl(23)-CovariateAssistedIV} analyze covariate-assisted Balke-Pearl bounds on the average treatment effect.
\cite{semenova2023generalized} generalizes the \cite[][]{Lee(09)} bounds on the average treatment effect on the always takers, which requires a binary instrument that monotonically affects selection, to incorporate covariate information. Monotonicity is particularly unlikely to hold when the instrument arises from the random assignment of decision makers. See, for example, \citet[][]{FrandsenLefgrenLeslie(19), ArnoldDobbieHull(20)-discrimination_bail, ChanGentzkowYu(21)}.}
Our work relates to sensitivity analysis in causal inference, where researchers bound treatment effects under relaxations of unconfoundedness.
Influential frameworks include the marginal sensitivity model \citep[e.g.,][]{Tan(06), KallusMaoZhou(18)-IntervalEstimation, ZhaoSmallBhattacharya(19), DornGuoKallus(21), KallusZhou(21), DornGuo(21)}, Rosenbaum's sensitivity model \citep[e.g.,][]{Rosenbaum(02), yadlowsky2018bounds}, and partial $c$-dependence \citep[e.g.,][]{MastenPoirier(18), MastenPoirier(20)}.
These models bound how unobservables affect selection propensities—assumptions that may be difficult to formulate when decisions involve human judgment. We instead bound outcome differences, which may be more intuitive in applied settings. Nonetheless, these models imply specific bounding functions in our framework, as we discuss in Remark \ref{remark:causal}.
\citet{chernozhukov2022long} develop omitted variable bounds for linear functionals; this covers our overall performance measures but not class-specific performance.
\section{Setup and Identification Framework}\label{section: model set-up}
We consider data $(X_i, D_i, Y_i)$ for $i = 1, \hdots, n$ drawn i.i.d. from a joint distribution $\mathbb{P}(\cdot)$, where $X_i \in \mathcal{X} \subseteq \mathbb{R}^{d}$ are covariates and $D_i \in \{0, 1\}$ is the decision.
There is a ``selective labels'' problem: the observed outcome satisfies $Y_i = D_i \, Y_i^*$ for some true outcome $Y_i^* \in \{0, 1\}$ \citep[e.g.,][]{KLLLM(18)}.
We assume $0 < \mathbb{P}(Y_i^* = 1) < 1$ so that we face a non-trivial prediction problem and that the decision satisfies strict overlap (i.e., there exists $\epsilon > 0$ such that $\mathbb{P}(D_i = 1 \mid X_i = x) \geq \epsilon$ with probability one).
\vspace{-1em}
\paragraph{Example: lending decisions} We observe data on loan applications, where $X_i$ contains information reported on the application such as income and prior credit history, $D_i$ is whether the application was granted a loan, and $Y_i = D_i \, Y_i^*$ is whether the applicant defaulted if granted the loan.
We study credit scores that predict whether a new applicant would default $Y_i^* = 1$ \citep[e.g.,][]{BlattnerNelson(21), DiMaggioEtAl(22)-Fintech, FusterEtAl(22)predictably}. $\blacktriangle$
\vspace{-1em}
\paragraph{Example: pretrial release} We observe data on defendants, where $X_i$ contains information about the defendant such as their current charge and prior conviction history, $D_i$ is whether the defendant was released before trial, and $Y_i = D_i \, Y_i^*$ is whether the defendant committed pretrial misconduct if released.
We study pretrial risk tools that predict whether a new defendant would commit pretrial misconduct $Y_i^* = 1$ \citep[e.g.,][]{KLLLM(18), ArnoldDobbieHull(20)-discrimination_alg, AngelovaEtAl(22)}. $\blacktriangle$
\vspace{-1em}
\paragraph{Example: medical testing} We observe data on patients, where $X_i$ contains patient information such as vital signs at admission and prior medical history, $D_i$ is whether a doctor decided to conduct a diagnostic test such as a cardiac stress test or lab test, and $Y_i = D_i \, Y_i^*$ is whether the patient was diagnosed with a heart attack or bacterial infection.
We study medical risk tools that predict whether a new patient suffered a heart attack or bacterial infection $Y_i^* = 1$ \citep[e.g.,][]{MullainathanObermeyer(21), HuangEtAl(22)}. $\blacktriangle$
\medskip
Let $\pi_d(x) := \mathbb{P}(D_i = d \mid X_i = x)$ for $d \in \{0, 1\}$ be the propensity score and $\mu_{1}(x) := \mathbb{P}(Y_i^* = 1 \mid D_i = 1, X_i = x)$.
Averages of a random variable $V_i$ are $\mathbb{E}_n[V_i] := n^{-1} \sum_{i=1}^{n} V_i$.
We let $\| \cdot \|$ denote the appropriate $L_2$-norm by context, where $\|f \| = \left( \int f(v)^2 d\mathbb{P}(v) \right)^{1/2}$ for a measurable function $f(\cdot)$ and $\|v\| = \left( \sum_{j=1}^{k} v_j^2 \right)^{1/2}$ for a vector $v \in \mathbb{R}^k$.
\subsection{Predictive Algorithms and Performance Measures}\label{section: estimands}
A \textit{predictive algorithm} $s(\cdot) \colon \mathcal{X} \rightarrow [0, 1]$ predicts the probability $Y_i^* = 1$ given covariates $X_i$, therefore estimating the \textit{conditional likelihood} $\mu^*(x) := \mathbb{P}(Y_i^* = 1 \mid X_i = x)$.
Given a predictive algorithm, researchers evaluate its performance in two ways.
First, given a predictive algorithm, researchers construct decision rules $d(\cdot) \colon \mathcal{X} \rightarrow \{0, 1\}$ by applying a threshold rule; that is, defining $d(X_i) = 1\{s(X_i) \geq \tau\}$ for some threshold $\tau \in [0, 1]$.
Alternative decision rules are evaluated based on their \textit{expected payoff}
\begin{equation}\label{eqn: definiton of expected welfare}
U(d) := \mathbb{E}[ (-u_{1,1}(X_i) Y_i^* + u_{1,0}(X_i) (1 - Y_i^*)) d(X_i) + (-u_{0,0}(X_i) (1-Y_i^*) + u_{0,1}(X_i) Y_i^*) (1 - d(X_i)],
\end{equation}
where each $u_{d,y^*}(\cdot) \geq 0$ for $d,y^* \in \{0,1\}^2$ specifies the known payoff associated with a combination of decision $D_i$ and outcome $Y_i^*$ at features $X_i$.
Payoffs are normalized so that $\sum_{d,y^* \in \{0,1\}^2} u_{d,y^*}(x) = 1$ for all $x \in \mathcal{X}$.
In consumer lending, the profitability of approving customers that would not default $u_{1,0}(\cdot)$ may vary based on features like the requested loan size.
In pretrial release, the societal benefits of releasing defendants that would not commit misconduct $u_{1,0}(\cdot)$ may vary based on the defendant's age, charge severity, or prior history of pretrial misconduct \citep[e.g., ][]{Baughman(17)}.
Second, given a predictive algorithm $s(\cdot)$, researchers evaluate its performance over the full population or specific groups.
The target parameters are its \textit{overall performance}
\begin{equation}\label{eqn: overall performance}
\mbox{perf}(s; \beta) = \mathbb{E}[\beta_0(X_i; s) + \beta_1(X_i; s) Y_i^*]
\end{equation}
and its \textit{positive class} and \textit{negative class performance}
\begin{equation}\label{eqn: class performance}
\mbox{perf}_{+}(s; \beta) = \mathbb{E}[\beta_0(X_i; s) \mid Y_i^* = 1] \mbox{ and } \mbox{perf}_{-}(s; \beta) = \mathbb{E}[\beta_0(X_i; s) \mid Y_i^* = 0],
\end{equation}
where $\beta_0(X_i; s), \beta_1(X_i; s)$ are researcher-specified functions that may depend on covariates and the predictive algorithm.
Alternative choices of $\beta_0(X_i; s), \beta_1(X_i ; s)$ recover popular evaluations used in empirical research.
We discuss several leading examples.
\vspace{-1em}
\paragraph{Example: mean square error and calibration} \textit{Mean square error} is $\mathbb{E}[\left( s(X_i) - Y_i^* \right)^2] = \mbox{perf}(s; \beta)$ for $\beta_0(X_i; s) = s^2(X_i)$, $\beta_1(X_i; s) = 1 - 2s(X_i)$.
\textit{Calibration} at prediction bin $[r_1, r_2] \subseteq [0,1]$ is $\mathbb{E}[Y_i^* \mid r_1 \leq s(X_i) \leq r_2] = \mbox{perf}(s; \beta)$ for $\beta_0(X_i; s) := 0$, $\beta_1(X_i; s) := \frac{1\{r_1 \leq s(X_i) \leq r_2\}}{\mathbb{P}\left(r_1 \leq s(X_i) \leq r_2 \right)}$ assuming $\mathbb{P}(r_1 \leq s(X_i) \leq r_2) > 0$. $\blacktriangle$
\vspace{-1em}
\paragraph{Example: generalized true positive and false positive rates}
Generalized \textit{true positive rate} is $\mathbb{E}[s(X_i) \mid Y_i^* = 1] = \mbox{perf}_{+}(s; \beta)$ for $\beta_0(X_i; s) = s(X_i)$, and generalized \textit{false positive rate} is $\mathbb{E}[s(X_i) \mid Y_i^* = 0] = \mbox{perf}_{-}(s; \beta)$ for $\beta_0(X_i; s) = s(X_i)$. $\blacktriangle$
\vspace{-1em}
\paragraph{Example: ROC curve}
The \textit{true positive rate} or \textit{recall} at threshold $\tau$ is $\mathbb{P}\left( s(X_i) \geq \tau \mid Y_i^* = 1\right) = \mbox{perf}_{+}(s; \beta)$ for $\beta_{0}(X_i; s) = 1\{s(X_i) \geq \tau\}$.
The \textit{false positive rate} at threshold $\tau$ is $\mathbb{P}\left(s(X_i) \geq \tau \mid Y_i^* = 0\right) = \mbox{perf}_{-}(s; \beta)$.
The \textit{Receiver Operating Characteristic (ROC) curve} $\{( \mathbb{P}(s(X_i) \geq \tau \mid Y_i^* = 0 ), \mathbb{P}( s(X_i) \geq \tau \mid Y_i^* = 1 ) ) \colon \tau \in [0, 1] \}$ summarizes the predictive algorithm's ability to separate the positive class $Y_i^* = 1$ from the negative class $Y_i^* = 0$ as $\tau$ varies. $\blacktriangle$
\vspace{-1em}
\paragraph{Example: precision-recall curve}
The \textit{precision} at threshold $\tau$ is $\mathbb{P}\left( Y_i^* = 1 \mid s(X_i) \geq \tau \right) := \mbox{perf}(s; \beta)$ for $\beta_0(X_i; s) = 0$ and $\beta_1(X_i; s) = \frac{1\{ s(X_i) \geq \tau \}}{\mathbb{P}\left( s(X_i) \geq \tau \right)}$ assuming $\mathbb{P}(s(X_i) \geq \tau) > 0$.
The \textit{precision-recall curve} $\{ \left( \mathbb{P}\left( s(X_i) \geq \tau \mid Y_i^* = 1\right), \mathbb{P}\left( Y_i^* = 1 \mid s(X_i) \geq \tau \right) \right) \colon \tau \in [0, 1] \}$ summarizes the predictive algorithm's ability to accurately classify $Y_i^* = 1$ as $\tau$ varies. $\blacktriangle$
\medskip
\noindent Returning to the earlier examples, \cite{BlattnerNelson(21)}, \cite{FusterEtAl(22)predictably}, and \cite{DiMaggioEtAl(22)-Fintech} compare mean square error, ROC curves, and precision-recall curves of alternative predictive algorithms for default likelihood in consumer finance.
\cite{KLLLM(18)} and \cite{MullainathanObermeyer(21)} analyze predictive algorithms in pretrial release and medical diagnosis by reporting calibration at various prediction bins and ROC curves.\footnote{In analyzing pretrial risk assessments, \cite{lakkaraju2017selective} and \cite{AngelovaEtAl(22)} consider the ``failure rate'' or ``counterfactual misconduct rate'' $\mathbb{P}(Y_i^* = 1, s(X_i) \leq \tau)$ at threshold $\tau \in [0,1]$. This is recovered by the overall performance $\mbox{perf}(s; \beta)$ for $\beta_0(X_i; s) = 0, \beta_1(X_i;s) = 1\{s(X_i) \leq \tau\}$.}
\subsection{Identification}\label{section: identification of estimands}
Selective labels create an identification problem: since $Y_i^*$ is only observed when $D_i = 1$, the conditional likelihood $\mathbb{P}(Y_i^* = 1 \mid D_i = 0, X_i)$ is not point identified from the data alone.
Consequently, neither are the conditional likelihood, expected payoffs for any decision rule, nor predictive performance measures defined in Section \ref{section: estimands}.
A natural approach would be to assume unconfounded selection --- that is, $Y_i^* \perp \!\!\!\perp D_i \mid X_i$ \citep[e.g.,][]{Coston(20)-Evaluation, mishler2021fade}.
Under unconfoundedness, the conditional likelihood among unselected units equals that among selected units with the same covariates: $P(Y_i^* = 1 \mid D_i = 0, X_i) = P(Y_i^* = 1 \mid D_i = 1, X_i)$.
However, this assumption is often implausible in practice.
Selected and unselected units may differ systematically in ways not captured by observed covariates $X_i$.
In consumer lending, applicants' decisions to accept loan terms may depend on competing offers from other lenders or their private assessment of their ability to repay --- information unavailable to the researcher.
In pretrial release, judges may rely on courtroom interactions or details from case files that are not recorded in administrative data.
Consequently, unobserved confounding --- unobservables correlated with both the outcome $Y_i^*$ and the decision $D_i$ --- is a central challenge.
Rather than imposing unconfoundedness, we partially identify the conditional likelihood, expected payoffs, and predictive performance measures by placing restrictions on how much outcomes can differ between unselected and selected units with the same covariates.
These restrictions formalize researcher beliefs about the severity of unobserved confounding in a transparent way and nest popular empirical strategies as special cases.
The key to our approach is bounding the \textit{confounding function} $\delta(x) := \mathbb{P}(Y_i^* = 1 \mid D_i = 0, X_i = x) - \mathbb{P}(Y_i^* = 1 \mid D_i = 1, X_i = x)$, which summarizes the difference in outcomes between unselected and selected units conditional on covariates.
When $\delta(x) = 0$ for all $x \in \mathcal{X}$, the selection decision is unconfounded.
When $\delta(x) \neq 0$, unobserved confounding is present: units with the same observed covariates have systematically different outcomes depending on whether they were selected.
\begin{assumption}\label{asm: bounding assumption}
For researcher-specified bounding functions $\underline{\delta}(\cdot; \eta), \overline{\delta}(\cdot; \eta)$ with (possibly infinite-dimensional) nuisance parameters $\eta$, the confounding function satisfies
\begin{equation}\label{eqn: bounding assumption}
\underline{\delta}(x; \eta) \leq \delta(x) \leq \overline{\delta}(x; \eta) \mbox{ for all } x \in \mathcal{X}.
\end{equation}
Let $\Delta$ denote the set of all confounding functions satisfying Equation \eqref{eqn: bounding assumption}.
\end{assumption}
\noindent Assumption \ref{asm: bounding assumption} restricts how much outcomes can differ between unselected and selected units with the same observed covariates.
In consumer lending, this bounds the difference in default rates between unfunded and funded applicants conditional on application characteristics.
In pretrial release, this bounds the difference in misconduct rates between detained and released defendants conditional on charge and criminal history.
Assumption \ref{asm: bounding assumption} allows researchers to encode beliefs about unobserved confounding without committing to a specific selection model.
It nests two important special cases: unconfoundedness with $\underline{\delta}(x; \eta) = \overline{\delta}(x; \eta) = 0$ for all $x \in \mathcal{X}$, and the assumption-free, worst-case bounds with $\underline{\delta}(x; \eta) = -\mathbb{P}(Y_i^* = 1 \mid D_i = 1, X_i = x)$, $\overline{\delta}(x; \eta) = 1 - \mathbb{P}(Y_i^* = 1 \mid D_i= 1, X_i = x)$ \citep[e.g.,][]{Manski(89), Manski(94)-SampleSelectionBounds}.
Thus Assumption \ref{asm: bounding assumption} interpolates between strong identifying assumptions and possibly uninformative bounds.
More broadly, Assumption \ref{asm: bounding assumption} follows the structure of widely-used partial identification frameworks in econometrics, most notably \citet[][]{Manski(03)}'s "approximate mean independence" assumption.
It extends this by conditioning on covariates and allowing bounds to depend on identified nuisance parameters.
In this sense, Assumption \ref{asm: bounding assumption} formalizes restrictions in terms that align with domain knowledge: researchers in these settings naturally reason about outcome differences rather than, say, how unobservables affect selection propensities.
As shown in Section \ref{section: choice of bounding function}, the bounding functions can be chosen through direct calibration based on domain expertise or derived indirectly from additional data features such as proxy outcomes or instrumental variables.
Each strategy is instantiated by specifying both the functional form of the bounding functions and the collection of nuisance parameters $\eta$ (possibly infinite-dimensional) to be estimated.
Under Assumption \ref{asm: bounding assumption}, the conditional likelihood and evaluation estimands are partially identified.
We assume throughout that the chosen bounding functions yield valid conditional probabilities. Specifically, for all $\delta \in \Delta$ and $x \in \mathcal{X}$, $\mu_1(x) + \pi_0 \delta(x; \eta) \in [0, 1]$ and $\delta(x; \eta) \in [-1, 1]$.
This is a mild restriction in practice; as we discuss below, for common bounding functions, these constraints can be enforced through the researcher's specification of sensitivity parameters or are satisfied automatically, as we will discuss next.
If these constraints are violated and the researcher proceeds with our identification framework, the resulting bounds will be conservative (non-sharp) but may still be informative.
Since $\mu^*(x) = \mu_1(x) + \pi_0(x) \delta(x)$,
the \textit{identified set} is the set of all values consistent with the bounds on $\delta(x)$,
\begin{equation*}
\mathcal{H}(\mu^{*}(x)) = \left\{ \mu_{1}(x) + \delta(x) \pi_0(x) \mbox{ such that } \delta(x) \in \Delta \right\}.
\end{equation*}
The identified sets for the performance measures, $\mathcal{H}(\mbox{perf}(s; \beta))$ and $\mathcal{H}(\mbox{perf}_{+}(s; \beta))$, and for the expected payoff of a decision rule $\mathcal{H}(U(d))$ are defined analogously.
As shorthand, let $\beta_{0,i} := \beta_0(X_i; s), \beta_{1,i} = \beta_{1}(X_i; s)$ and $\underline{\delta}_i := \underline{\delta}(X_i; \eta)$, $\overline{\delta}_i := \overline{\delta}_i(X_i; \eta)$.
\begin{lemma}\label{lem: bounds under MOSM}
For all $x \in \mathcal{X}$, $\mathcal{H}(\mu^*(x)) = \left[\underline{\mu}^{*}(x), \overline{\mu}^{*}(x) \right],$
where $\overline{\mu}^{*}(x) = \mu_{1}(x) + \pi_0(x) \overline{\delta}(x; \eta) $, $\underline{\mu}^{*}(x) = \mu_1(x) + \pi_0(x) \underline{\delta}(x; \eta)$.
Furthermore,
\begin{align*}
& \mathcal{H}(\mbox{perf}(s; \beta)) = \left[ \underline{\mbox{perf}}(s; \beta), \overline{\mbox{perf}}(s; \beta) \right] \mbox{ and }\mathcal{H}(\mbox{perf}_{+}(s; \beta)) = \left[ \underline{\mbox{perf}}_{+}(s; \beta), \overline{\mbox{perf}}_{+}(s; \beta) \right],
\end{align*}
where $\overline{\mbox{perf}}(s; \beta) = \mathbb{E}[\beta_{0,i} + \beta_{1,i} \mu_{1}(X_i) + \beta_{1,i} \pi_0(X_i) \left( 1\{\beta_{1,i} > 0\} \overline{\delta}_{i} + 1\{\beta_{1,i} \leq 0\} \underline{\delta}_{i} \right) ]$, $\underline{\mbox{perf}}(s; \beta) = \mathbb{E}[\beta_{0,i} + \beta_{1,i} \mu_{1}(X_i) + \beta_{1,i} \pi_0(X_i) \left( 1\{\beta_{1,i} \leq 0\} \overline{\delta}_{i} + 1\{\beta_{1,i} > 0\} \underline{\delta}_{i} \right) ]$, and
\begin{align*}
& \overline{\mbox{perf}}_{+}(s; \beta) = \sup_{\delta \in \Delta} \mathbb{E}[\mu_{1}(X_i) + \pi_0(X_i) \delta(X_i)]^{-1} \mathbb{E}[ \beta_{0,i} \mu_{1}(X_i) + \beta_{0,i} \pi_0(X_i) \delta(X_i) ], \\
& \underline{\mbox{perf}}_{+}(s; \beta) = \inf_{\delta \in \Delta} \mathbb{E}[\mu_{1}(X_i) + \pi_0(X_i) \delta(X_i)]^{-1} \mathbb{E}[ \beta_{0,i} \mu_{1}(X_i) + \beta_{0,i} \pi_0(X_i) \delta(X_i) ].
\end{align*}
\end{lemma}
\begin{lemma}\label{lem: sharp bounds on ECU under MOSM}
For $d(\cdot) \colon \mathcal{X} \rightarrow \{0, 1\}$ and $u_{d, y^*,i} = u_{d,y^*}(X_i)$, $\mathcal{H}(U(d)) = [\underline{U}(d), \overline{U}(d)]$ for $\underline{U}(d), \overline{U}(d)$ defined in the proof.
\end{lemma}
The remainder of the paper applies this partial identification framework to the two objectives for predictive algorithms commonly seen in applied work.
Section \ref{section: DR Learner, main text} addresses algorithm design: we estimate bounds on the conditional likelihood, thereby constructing predictive algorithms robust to unobserved confounding.
Section \ref{section: robust evaluation, main text} addresses algorithm evaluation: we estimate bounds on performance measures for a given predictive algorithm.
Before proceeding to these estimation results, we next discuss how researchers might choose the bounding functions in Assumption \ref{asm: bounding assumption}.
\begin{remark}[Connection to Algorithmic Fairness]\label{remark: fairness}
Differences in predictive performance across sensitive attributes (e.g., race, gender) formalize violations of popular fairness criteria studied in the algorithmic fairness literature \citep[see][for reviews]{MitchellEtAl(19), barocas-hardt-narayanan}.
Specifically, differences in positive and negative class performance across groups characterize violations of ``equality of opportunity'' \citep{hardt2016equality}, ``equalized odds'' \citep{chouldechova2017fair, KleinbergMullainathanRaghavan(17)}, and ``group-level discrimination'' \citep{ArnoldDobbieHull(20)-discrimination_alg, ArnoldDobbieHull(20)-discrimination_bail}. Differences in overall performance across groups characterize violations of ``bounded group loss'' \citep{AgarwalEtAl(19)-FairRegression}. In Online Appendix \ref{section: identification and estimation of predictive disparities}, we partially identify these performance disparities under Assumption \ref{asm: bounding assumption}, enabling robust fairness audits.
\end{remark}
\subsection{Choice of Bounding Functions}\label{section: choice of bounding function}
How should the bounding functions $\underline{\delta}(\cdot; \eta)$ and $\overline{\delta}(\cdot; \eta)$ be chosen?
The choice encodes assumptions about how much unselected units might differ from selected units with the same observed characteristics. In consumer lending: How much more likely are unfunded applicants to default than funded applicants with the same credit profile? In pretrial release: How much more likely are detained defendants to commit misconduct than released defendants with the same criminal history?
We next show how alternative bounding functions formalize three common strategies for addressing selective labels in applied research.
\subsubsection{Observed Outcome Bounds}\label{section: observed outcome bounds}
A natural starting point is to assume that the unobserved outcome rate among unselected units is similar to the observed outcome rate among selected units with the same characteristics.
For instance, in pretrial release, detained defendants with a given criminal history might be 1.5 times as likely to commit misconduct as released defendants with the same history.
In consumer lending, unfunded applicants might be twice as likely to default as funded applicants with the same credit profile.
\textit{Observed outcome bounds} formalize this intuition by specifying, for researcher-chosen constants $\underline{\Gamma}, \overline{\Gamma} > 0$ and all $x \in \mathcal{X}$,
\begin{equation}
\underline{\delta}(x; \eta) := (\underline{\Gamma} - 1) \mu_1(x) \mbox{ and } \overline{\delta}(x; \eta) := (\overline{\Gamma} - 1) \mu_1(x).
\end{equation}
Under this choice, the unobserved conditional probability among unselected units satisfies $\underline{\Gamma} \mu_1(x) \leq P(Y_i^* = 1 \mid D_i = 0, X_i = x) \leq \overline{\Gamma} \mu_1(x)$ --- that is, it is bounded by constant multiples of the observed conditional probability $\mu_1(x)$.
Setting $\underline{\Gamma} = \overline{\Gamma} = 1$ recovers unconfoundedness.
For observed outcome bounds, the probability constraints are satisfied for all choices of $\underline{\Gamma}, \overline{\Gamma}$ in a neighborhood of one under strict overlap on the decision and observed outcome probabilities being bounded away from zero and one.
This approach formalizes strategies commonly used in applied work. \citet[][]{KLLLM(18)} evaluate a constructed pretrial risk tool under varying assumptions on how much misconduct rates among detained defendants differ from released defendants (their Table 5).
In consumer finance, researchers impute default rates for rejected applicants by applying multiplicative adjustments to observed default rates among accepted applicants \citep[e.g.,][]{HandHenley(93), ZengZhao(14)}. Observed outcome bounds flexibly encode such adjustments while allowing the imputed values to vary richly with covariates.
\begin{remark}[Connection to Sensitivity Analysis Models in Causal Inference]\label{remark:causal}
Several influential sensitivity analysis frameworks --- including the marginal sensitivity model \citep{Tan(06), ZhaoSmallBhattacharya(19), DornGuo(21)}, Rosenbaum's $\Gamma$-sensitivity model \citep{Rosenbaum(02), Yadlowsky(21)evaluating}, and partial $c$-dependence \citep{MastenPoirier(18), MastenPoirier(20)} --- place bounds on how unobservables affect selection propensities.
While conceptually distinct from bounding outcome differences, these models imply observed outcome bounds with specific choices of $\underline{\Gamma}, \overline{\Gamma}$ (see Online Appendix \ref{section: connections to existing sensitivity analysis models}).
The connection between outcome-based bounds and selection-based sensitivity models arises through Bayes' rule; though the observed outcome bounds we derive may not exhaust all implications of the original sensitivity model.
\end{remark}
\subsubsection{Proxy Outcome Bounds}\label{section: proxy variable bounds}
Sometimes researchers observe a \textit{proxy outcome} $\widetilde{Y}_i \in \{0, 1\}$ for both selected and unselected units that is statistically related to the true outcome $Y_i^*$.
This proxy can help discipline assumptions about unobserved confounding.
For example, in consumer lending, whether an applicant defaulted on other credit products (e.g., credit cards) serves as a proxy for whether they would default on the loan of interest \citep[][]{BlattnerNelson(21)}.
In medical testing, long-term health outcomes or total medical spending serves as a proxy for whether a patient suffered an acute condition \citep[][]{ObermeyerEtAl(19), ChanGentzkowYu(21), MullainathanObermeyer(21)}.
The key assumption is that the proxy outcome is equally informative about the true outcome for both selected and unselected units.
Formally, suppose the proxy outcome satisfies, for all $x \in \mathcal{X}$,
\begin{equation}\label{eqn: proxy assumption}
P(Y^*_i = \tilde{Y}_i \mid D_i = 0, X_i = x) = P(Y^*_i = \tilde{Y}_i \mid D_i = 1, X_i = x).
\end{equation}
Conditional on the covariates, the probability that the true and proxy outcome agree is the same regardless of selection status.
This assumption, combined with the observed distribution of $\widetilde{Y}_i$ among unselected units, pins down specific bounding functions on $\delta(x)$.
Under mild restrictions satisfied by the proxy outcomes used in \cite{BlattnerNelson(21)} and \cite{MullainathanObermeyer(21)}, the bounding functions simplify to
\begin{align*}
& \underline{\delta}(x; \eta) = 1 - \gamma_1(x) - \widetilde{\mu}_0(x) - \mu_1(x), \\
& \overline{\delta}(x; \eta) = 1 - \gamma_1(x) + \widetilde{\mu}_0(x) - \mu_1(x),
\end{align*}
where $\gamma_1(x) := \mathbb{P}(Y_i^* = \tilde{Y}_i \mid D_i = 1, X_i = x)$ is the probability the true and proxy outcome agree among selected units and $\widetilde{\mu}_0(x) := \mathbb{P}(\tilde{Y}_i = 1 \mid D_i = 0, X_i = x)$ is the observed proxy outcome rate among unselected units.
These bounding functions satisfy the probability constraints under the maintained proxy assumptions under strict overlap on the decision and observed outcome probabilities being bounded away from zero and one.
Online Appendix \ref{section: proxy outcome bounds} provides the general characterization without these simplifying restrictions.
To our knowledge, proxy outcome assumptions of this form have not been explored in the causal inference sensitivity analysis literature.
\subsubsection{Instrumental Variable Bounds}\label{section: instrumental variable bounds}
Sometimes researchers observe an \textit{instrumental variable}
$Z_i \in \mathcal{Z}$ with finite support that generates quasi-random variation in selection decisions but does not directly affect outcomes (i.e., $Y_i^*$ independent of $Z_i$ given $X_i$).
For example, in pretrial release, judges randomly assigned to cases vary in their propensity to release defendants \citep[][]{KLLLM(18), ArnoldDobbieHull(20)-discrimination_bail, Rambachan(21), AngelovaEtAl(22)}.
The day-of-week affects whether certain tests are ordered \citep[e.g.,][]{MullainathanObermeyer(21)}, and state-level variation in banking regulations affects loan approval rates \citep[e.g.,][]{BlattnerNelson(21)}.
An instrument $Z_i$ implies specific bounding functions through a classic result of \citet[][]{Manski(94)-SampleSelectionBounds}.
For any fixed instrument value $z \in \mathcal{Z}$, the confounding function satisfies, for all $x \in \mathcal{X}$,
$$
\underline{\delta}_{z}(x; \eta) \leq \delta(x) \leq \overline{\delta}_{z}(x; \eta),
$$
where $\underline{\delta}_{z}(x) = \left( \mathbb{E}[Y_i D_i \mid X_i = x, Z_i = z] - \mu_1(x) \right)/\pi_0(x)$ and $\overline{\delta}_z(x) := ( \pi_0(x, z) + \mathbb{E}[Y_i D_i \mid X_i = x, Z_i = z] - \mu_1(x) )/\pi_0(x)$.
Our identification and estimation results apply to these bounds for any fixed $z \in \mathcal{Z}$.
Researchers can obtain sharper bounds by aggregating across all instrument values: $\max_{z \in \mathcal{Z}} \ \underline{\delta}_{z}(x; \eta) \leq \delta(x) \leq \min_{z \in \mathcal{Z}} \ \overline{\delta}_{z}(x; \eta)$.
These intersection bounds \citep[][]{Chernozhukov-IntersectionBounds} are tighter but involve maximum and minimum operators, creating non-smoothness that complicates estimation.
Following \citet{LevisEtAl(23)-CovariateAssistedIV}, we address this challenge using smooth approximations to the intersection bounds.
Online Appendix \ref{section: log-sum-exp bounds for IV} discusses both approaches: bounds based on a fixed instrument value and smooth approximations to intersection bounds.\footnote{\citet{LevisEtAl(23)-CovariateAssistedIV} and \citet{semenova2023adaptive} consider an alternative strategy for dealing with pointwise intersection bounds that invokes a ``margin condition'' on the separation of the maximal value of the intersection bounds.
Under such a margin condition, existing results would apply to bounds on overall performance.
Extending this approach to estimators for the conditional likelihood bounds or positive class performance bounds is not immediate.
We leave the full analysis of margin conditions for future work.}
\section{Robust Design of Predictive Algorithms}\label{section: DR Learner, main text}
In this section, we propose estimators for the bounds on the conditional likelihood $\mu^*(\cdot)$ under Assumption \ref{asm: bounding assumption}.
We build on recent work studying the estimation of heterogeneous treatment effects in causal inference \citep[e.g.,][]{KunzelEtAl(19), NieWager(20), Kennedy(22)-towardsoptimal}, and our estimators use cross-fitting and doubly-robust bias corrections.
Since the bounds on the conditional likelihood are an infinite-dimensional parameter, we establish finite-sample bounds on our estimator's integrated mean square error (MSE) convergence using a novel oracle inequality for pseudo-outcome regression procedures.
These convergence guarantees enable us to construct plug-in decision rules with provable regret bounds.
To simplify exposition, we develop our estimator assuming the researcher specifies observed outcome bounds (Section \ref{section: observed outcome bounds}) and illustrate it with two folds.
The analysis for multiple folds and cross-fitting is straightforward.
We extend to proxy outcome bounds in Online Appendix \ref{section: proxy outcome bounds} and instrumental variable bounds in Online Appendix \ref{section: log-sum-exp bounds for IV}.
\vspace{-1em}
\paragraph{Estimation Procedure for Observed Outcome Bounds:}
Under Assumption \ref{asm: bounding assumption}, the unobserved conditional probability among unselected units can be decomposed as $P(Y^*_i = 1 | D_i = 0, X_i = x) = \mu_1(x) + \delta(x)$, where $\mu_1(x)$ is the observed rate among selected units and $\delta(x)$ measures confounding.
Our estimation strategy constructs pseudo-outcomes based on influence functions that combine observed outcomes among selected units with extrapolation to unselected units.
The pseudo-outcomes correct for both observable selection (via inverse propensity weighting) and unobservable confounding (via the sensitivity parameters $\underline{\Gamma}, \overline{\Gamma}$).
We estimate the bounds on the conditional likelihood under observed outcome bounds with the following steps.
Split the data into two folds and estimate the nuisance functions $\hat \eta = \left( \hat{\mu}_1(\cdot), \hat{\pi}_1(\cdot) \right)$ on the first fold.
On the second fold, construct the pseudo-outcomes $\phi_{\mu,i}(\hat \eta) + \left( \overline{\Gamma} - 1 \right) \phi_{\pi \mu, i}(\hat \eta)$ for the upper bound and $\phi_{\mu,i}(\hat \eta) + \left( \underline{\Gamma} - 1 \right) \phi_{\pi \mu, i}(\hat \eta)$ for the lower bound, where $\phi_{\mu,i}(\eta) := \mu_1(X_i) + \frac{D_i}{\pi_1(X_i)} (Y_i - \mu_1(X_i))$ and $\phi_{\pi \mu, i}(\eta) := \left((1-D_i) - \pi_0(X_i)\right) \mu_1(X_i) + \frac{D_i}{\pi_1(X_i)} (Y_i - \mu_1(X_i)) \pi_0(X_i) + \pi_0(X_i) \mu_1(X_i)$ are the influence functions for $\mathbb{E}\left[ \mu_1(X_i)\right]$ and $\mathbb{E}\left[ \pi_0(X_i) \mu_1(X_i) \right]$ respectively.
These can be derived using standard arguments \citep[e.g.,][]{Hines_2022}.
We regress the estimated pseudo-outcomes on the covariates $X_i$ using a researcher-specified nonparametric regression procedure in the second fold.
This yields our estimators $\widehat{\underline{\mu}}(\cdot)$, $\widehat{\overline{\mu}}(\cdot)$.
\subsection{Bound on Integrated Mean Square Error Convergence}\label{section: main text, L2 inequality}
We analyze the integrated MSE convergence rate of our proposed estimators by comparing them against an infeasible oracle nonparametric regression.
This oracle has access to the true nuisance functions $\eta$ and estimates the conditional likelihood bounds by regressing the true pseudo-outcomes $\phi_{\mu,i}(\eta) + \left( \overline{\Gamma} - 1 \right) \phi_{\pi \mu, i}(\eta)$ for the upper bound and $\phi_{\mu,i}(\eta) + \left( \underline{\Gamma} - 1 \right) \phi_{\pi \mu, i}(\eta)$ for the lower bound on $X_i$ in the second fold using the same nonparametric regression procedure specified by the researcher.
We denote these oracle estimators by $\widehat{\underline{\mu}}_{oracle}(\cdot), \widehat{\overline{\mu}}_{oracle}(\cdot)$.
This oracle comparison is useful because it decomposes the integrated MSE of our proposed estimators into two sources.
The first is approximation error from using a finite-sample nonparametric regression procedure to estimate the infinite-dimensional conditional likelihood bounds --- this is captured by the oracle's integrated MSE.
The second is estimation error from first-stage nuisance function estimation.
Our next result shows that the latter only affects our proposed estimator through a smoothed, doubly-robust remainder term.
\begin{proposition}\label{prop: oracle result for learning, outcome regression bounds}
Let $\widehat{\mathbb{E}}_n[\cdot \mid X_i = x]$ denote the second-stage pseudo-outcome regression estimator.
Suppose $\widehat{\mathbb{E}}_n[\cdot \mid X_i = x]$ satisfies the $L_2(\mathbb{P})$-stability condition (Assumption \ref{asm: L2 stability condition}), and $\mathbb{P}(\epsilon \leq \hat{\pi}_{1}(X_i) \leq 1 - \epsilon) = 1$ for some $\epsilon > 0$.
Define $\tilde{R}(x) = \widehat{\mathbb{E}}_n[ (\pi_1(X_i) - \hat{\pi}_1(X_i)) (\mu_{1}(X_i) - \hat{\mu}_1(X_i)) \mid X_i = x ]$, and $R_{oracle}^2 = \mathbb{E}[ \| \widehat{\overline{\mu}}_{oracle}(\cdot) - \overline{\mu}^*(\cdot) \|^2]$.
Then,
\begin{align*}
& \| \widehat{\overline{\mu}}(\cdot) - \overline{\mu}^*(\cdot) \| \leq \|\widehat{\overline{\mu}}_{oracle}(\cdot) - \overline{\mu}^*(\cdot) \| + \epsilon^{-1} \sqrt{2} ( \overline{\Gamma} - 1 ) \| \tilde{R}(\cdot) \| + o_\mathbb{P}(R_{oracle})
\end{align*}
The analogous result holds for $\widehat{\underline{\mu}}(x)$.
\end{proposition}
\noindent Proposition \ref{prop: oracle result for learning, outcome regression bounds} delivers two insights.
First, our analysis is modular: researchers can pair choices of nuisance function estimators in the first step with any choice of nonparametric regression method (satisfying Assumption \ref{asm: L2 stability condition}) in the second step, and the integrated MSE is characterized by the oracle rate plus a doubly-robust remainder term.
Second, Proposition \ref{prop: oracle result for learning, outcome regression bounds}
shows that the cost of estimating the nuisance functions in the first step only enters through the product of nuisance function errors.
To prove this result, we establish an oracle inequality on the $L_2(\mathbb{P})$-error of pseudo-outcome regression procedures (Lemma \ref{lemma: L2 oracle inequality}), extending \cite{Kennedy(22)-towardsoptimal}'s pointwise analysis.
This result may be of independent interest.
The $L_2(\mathbb{P})$-stability condition (Assumption \ref{asm: L2 stability condition}) on $\widehat{\mathbb{E}}_n[\cdot \mid X_i = x]$ is quite mild in practice.
It is satisfied by a variety of generic linear smoothers such as linear regression, series regression, nearest neighbor matching, random forest models, and several others (Proposition \ref{prop: linear smoothers satisfy L2 stability}).
Once again, in this sense, Proposition \ref{prop: oracle result for learning, outcome regression bounds} is agnostic to both the researcher's choice of ${\widehat{\mathbb{E}}_n[\cdot \mid X_i = x]}$ and nuisance function estimators.
Furthermore, Proposition \ref{prop: oracle result for learning, outcome regression bounds} can be applied in settings where nuisance functions satisfy additional smoothness or sparsity conditions, and for particular choices of the second-stage regression estimator and nuisance function estimators.
Known results on mean-squared error convergence rates of nonparametric regression procedures can then be used to characterize $\|\widehat{\overline{\mu}}_{oracle}(\cdot) - \overline{\mu}^*(\cdot) \|$ and $\| \tilde{R}(\cdot) \|$, yielding specific rates of convergence in terms of primitives such as underlying smoothness assumptions and sample size.
This follows in the spirit of recent work on optimal estimation of heterogeneous treatment effects \citep[][]{Kennedy(22)-towardsoptimal, KennedyEtAl(23)-Minimax}.
\subsection{Regret Bound for Plug-In Decision Rules}
Researchers often apply a threshold rule to an estimated predictive algorithm to make decisions.
For example, \cite{KLLLM(18), Rambachan(21), AngelovaEtAl(22)} compare outcomes under alternative decisions rules that threshold a predictive algorithm in pretrial release.
We show that plug-in decision rules based on our estimators of the conditional likelihood bounds achieve low regret, ensuring that the statistical uncertainty from estimation does not lead to substantially worse decisions than if the true conditional probabilities were known.
Since the expected payoff $U(d)$ of any decision rule $d(\cdot) \colon \mathcal{X} \rightarrow \{0, 1\}$ is partially identified (Lemma \ref{lem: sharp bounds on ECU under MOSM}), we evaluate decision rules by comparing their worst-case expected welfare \citep[e.g.,][]{Manski(07)}.
An immediate consequence of Lemma \ref{lem: sharp bounds on ECU under MOSM} is that the optimal max-min decision rule $d^*(\cdot) \in \arg \max_{d(\cdot) \colon \mathcal{X} \rightarrow [0,1]} \underline{U}(d)$ is a threshold rule:
$$
d^*(X_i) = 1\{ \widetilde{\mu}^*(X_i) \leq u_{1,0,i} + u_{0,0,i}\},
$$
where $\widetilde{\mu}^*(x) = (u_{1,1,i} + u_{1,0,i}) \overline{\mu}^*(x) + (u_{0,0,i} + u_{0,1,i}) \underline{\mu}^*(x)$ is a welfare-weighted average of the conditional probability bounds.
We construct a feasible plug-in decision rule $\widehat{d}(X_i)$ by substituting in our nonparametric estimators $\widehat{\overline{\mu}}(\cdot)$, $\widehat{\underline{\mu}}(\cdot)$.
The regret of the plug-in decision rule is $R(\hat{d}) = \underline{U}(d^*) - \underline{U}(\widehat{d})$, which measures the welfare loss from using estimated rather than true conditional probabilities.
Our next result bounds squared regret by the integrated MSE of the oracle estimators plus the smoothed, doubly-robust remainder term.
\begin{proposition}\label{prop: regret bound for feasible plug-in decision rule}
Under the same conditions as Proposition \ref{prop: oracle result for learning, outcome regression bounds},
$$
R(\widehat{d})^2 \leq 2 \| \widehat{\overline{\mu}}_{oracle}(\cdot) - \overline{\mu}^*(\cdot) \| + 2 \| \widehat{\underline{\mu}}_{oracle}(\cdot) - \underline{\mu}^*(\cdot) \| + 4 \epsilon^{-1} \sqrt{2} \| \tilde{R}(x) \| + o_\mathbb{P}( R_{oracle} ).
$$
\end{proposition}
\noindent The proof proceeds in two steps. First, we show that the regret of the plug-in decision rule is bounded by the integrated MSE of our estimators for the conditional probability bounds --- this follows because the optimal decision rule is a threshold rule, and threshold rules are sensitive to estimation errors measured in $L_2(\mathbb{P})$ norm.
Second, we apply Proposition \ref{prop: oracle result for learning, outcome regression bounds} to decompose this integrated MSE into the oracle rate plus the doubly-robust remainder term.
Consequently, worst-case regret converges to zero whenever the oracle's integrated MSE converges to zero.
Plug-in decision rules based on our nonparametric estimators are therefore robust: they achieve nearly optimal welfare despite estimation uncertainty.
\section{Robust Evaluation of Predictive Algorithms}\label{section: robust evaluation, main text}
In this section, we construct estimators for the bounds on the performance of a given predictive algorithm $s(\cdot)$ under Assumption \ref{asm: bounding assumption}.
This setting arises empirically in two scenarios: (i) when researchers evaluate an externally provided algorithm (e.g., a proprietary credit score or a commercial risk assessment tool), and (ii) when researchers conduct hold-out evaluation of an algorithm trained on separate data.
In both cases, the algorithm $s(\cdot)$ is taken as fixed, and the goal is to assess its performance under varying assumptions about unobserved confounding.
Our estimators are consistent for the bounds on the identified set as sample size grows large, and we characterize their convergence rates in terms of errors in the first-step estimation of nuisance parameters.
\subsection{Bounds on Overall Performance}\label{section: estimating overall performance bounds}
We first construct estimators for the bounds on overall performance $\mbox{perf}(s; \beta)$.
To simplify exposition, we develop estimators for observed outcome bounds (Section \ref{section: observed outcome bounds}).
We extend to proxy outcome bounds and instrumental variable bounds in Online Appendices \ref{section: proxy outcome bounds}-\ref{section: log-sum-exp bounds for IV} respectively.
Under observed outcome bounds, Lemma \ref{lem: bounds under MOSM} implies that the upper bound can be written as
$$
\overline{\mbox{perf}}(s; \beta) := \mathbb{E}\left[ \beta_{0,i} + \beta_{1,i} \mu_1(X_i) + \beta_{1,i} \left( 1\{ \beta_{1,i} > 0 \} \left( \overline{\Gamma} - 1 \right) + 1\{ \beta_{1,i} \leq 0 \} (\underline{\Gamma} - 1) \right) \pi_0(X_i) \mu_1(X_i) \right].
$$
The lower bound $\underline{\mbox{perf}}(s; \beta)$ has an analogous expression.
Both are linear functionals of known functions of the data and identified nuisance parameters $\eta = (\mu_1, \pi_0)$.
This linearity enables us to construct debiased estimators using standard arguments based on influence functions and cross-fitting \citep[e.g.,][]{ChernozhukovEtAl(18)-DML, ChernozhukovEtAl(22)-LocallRobust, Kennedy(22)-IFreview}.
\vspace{-1em}
\paragraph{Estimation Procedure for Observed Outcome Bounds:} For completeness, we sketch the construction of our estimators based on $K$-fold cross-fitting.
We randomly split the data into $K$ disjoint folds.
For each fold $k$, we estimate the nuisance functions $\widehat{\eta}_{-k}$ using all observations not in the $k$-th fold and construct
$$
\overline{\mbox{perf}}_{i}(\hat \eta_{-k}) := \beta_{0,i} + \beta_{1,i} \phi_{\mu,i}(\hat \eta_{-k}) + \beta_{1,i} \left( 1\{ \beta_{1,i} > 0 \} \left( \overline{\Gamma} - 1 \right) + 1\{ \beta_{1,i} \leq 0 \} (\underline{\Gamma} - 1) \right) \phi_{\pi \mu, i}(\hat \eta_{-k})
$$
for each observation $i$ in the $k$-th fold, where $\phi_{\mu,i}(\eta)$ and $\phi_{\pi \mu, i}(\eta)$ are defined as before.
We then average over all observations $\widehat{\overline{\mbox{perf}}}(s; \beta) := \mathbb{E}_n\left[ \overline{\mbox{perf}}_i(\hat \eta_{-K_i}) \right],$
where $K_i$ is the fold containing observation $i$.
Our estimator for the lower bound $\widehat{\underline{\mbox{perf}}}(s;\beta)$ is defined analogously.
Our next result characterizes the convergence rate and asymptotic distribution of our estimators for the bounds on overall performance under observed outcome bounds.
\begin{proposition}\label{prop: overall performance estimator, nonparametric outcome bounds}
Define $R_{1,n}^{k} := \|\hat{\mu}_{1,-k}(\cdot) - \mu_{1}(\cdot)\| \|\hat{\pi}_{1,-k}(\cdot) - \pi_1(\cdot)\|$ for each fold $k$.
Assume (i) there exists $M < \infty$ such that $\| \beta_1(\cdot) \| \leq M$; (ii) $\mathbb{P}(\pi_1(X_i) \geq \delta) = 1$ for some $\delta > 0$, (iii) there exists $\epsilon > 0$ such that $\mathbb{P}(\hat{\pi}_{1,-k}(X_i) \geq \epsilon) = 1$ for each fold $k$, and (iv) $\|\hat{\mu}_{1, -k}(\cdot) - \mu_{1}(\cdot)\| = o_P(1)$ and $\|\hat{\pi}_{1, -k}(\cdot) - \pi_1(\cdot)\| = o_P(1)$ for each fold $k$.
Then,
$$\left| \widehat{\overline{\mbox{perf}}}(s; \beta) - \overline{\mbox{perf}}(s; \beta) \right| = O_\mathbb{P}\left(1 / \sqrt{n} + \sum_{k=1}^{K} R_{1,n}^{k}\right),$$
and the analogous result holds for $\widehat{\underline{\mbox{perf}}}(s; \beta)$.
If further $R_{1,n}^{k} = o_\mathbb{P}(1/\sqrt{n})$ for each fold $k$, then
\begin{equation*}
\sqrt{n} \left( \begin{pmatrix} \widehat{\overline{\mbox{perf}}}(s; \beta) \\ \widehat{\underline{\mbox{perf}}}(s; \beta)\end{pmatrix} - \begin{pmatrix} \overline{\mbox{perf}}(s; \beta) \\ \underline{\mbox{perf}}(s; \beta) \end{pmatrix} \right) \xrightarrow{d} N\left( 0, \Sigma \right)
\end{equation*}
for $\Sigma = Cov\left( (\overline{\mbox{perf}}_i(\eta), \underline{\mbox{perf}}_i(\eta))^\prime \right)$.
\end{proposition}
\noindent The estimators converge at rate $O_{\mathbb{P}}(1/\sqrt{n})$ plus a term that captures the nuisance estimation error.
Nuisance errors enter multiplicatively --- through the product of propensity score and outcome regression errors --- rather than additively.
The rate condition required for the estimators to be jointly asymptotically normal is satisfied under the standard condition that all nuisance function estimators converge at a rate faster than $O_\mathbb{P}(n^{-1/4})$.
This is the familiar condition from debiased machine learning, satisfied by a wide range of nonparametric regression and modern machine learning methods.
Online Appendix \ref{section: variance estimation for overall performance bounds} constructs a consistent estimator of the asymptotic covariance matrix in Proposition \ref{prop: overall performance estimator, nonparametric outcome bounds}.
This implies that confidence intervals for the identified set $\mathcal{H}(\mbox{perf}(s; \beta))$ can be constructed using standard methods \citep[e.g.,][]{ImbensManski(04), Stoye(09)--PICI}.
\subsection{Bounds on Positive Class Performance}\label{section: estimating pos class bounds}
We next estimate bounds on positive class performance $\mbox{perf}_{+}(s; \beta)$.
Unlike overall performance, this is a ratio of expectations, and the bounds are given by the optimal values of linear fractional programs as shown in Lemma \ref{lem: bounds under MOSM}.
This creates additional challenges in estimation.
In this section, we will develop estimators for these bounds and establish their convergence properties.
\subsubsection{Estimation via Linear Fractional Programming:}
We describe the construction of our estimators using sample-splitting across two folds; extending to cross-fitting with multiple folds is straightforward.
We randomly split the data into two folds (we assume $n$ is divisible by $2$ and each fold contains $n/2$ observations for simplicity).
In the first fold, we estimate the nuisance functions $\hat \eta := (\widehat{\mu}_1(\cdot), \widehat{\pi}_1(\cdot))$ and the bounding functions $(\widehat{\underline{\delta}}(\cdot), \widehat{\overline{\delta}}(\cdot))$ using one of the strategies in Section \ref{subsec: estimating bounds}.
On the second fold, we construct $\phi_{\mu,i}(\hat \eta)$ and solve the sample analogue
\begin{equation}\label{equation: sample pos class LFP}
\widehat{\overline{\mbox{perf}}}_{+}(s; \beta) = \max_{\tilde{\delta} \in \widehat{\Delta}_n} \frac{\mathbb{E}_n\left[ \beta_{0,i} \phi_{\mu, i}(\hat \eta) + \beta_{0,i} \tilde{\delta}_i \right]}{\mathbb{E}_n\left[ \phi_{\mu, i}(\hat \eta) + \tilde{\delta}_i \right]},
\end{equation}
where $\widehat{\Delta}_n := \left\{ \tilde{\delta} \colon (1 - D_i) \widehat{\underline{\delta}}(X_i) \leq \tilde{\delta}_i \leq (1 - D_i) \widehat{\overline{\delta}}(X_i) \right\}$ and $\mathbb{E}_{n}[\cdot]$ denotes the sample average over the second fold.
The optimization problem in Equation \eqref{equation: sample pos class LFP} can be solved efficiently using standard techniques.
It can be written as a linear program by applying the Charnes-Cooper transformation \citep[][]{CharnesCooper(62)}.
See Online Appendix \ref{section: LP reduction for positive class bounds} for details.
The construction of this estimator involves two key ideas.
First, we replace $\mu_1(\cdot)$ in the objective with its estimated influence function $\phi_{\mu,i}(\hat \eta)$. This orthogonalizes against first-stage estimation errors in nuisance functions that only enter the objective.
Second, we plug in estimated bounding functions into the constraint set.
We provide an error decomposition for the positive class performance bounds.
This result characterizes how errors in estimating the nuisance parameters $(\mu_1(\cdot), \pi_1(\cdot))$ and the bounding functions $(\underline{\delta}(\cdot), \overline{\delta}(\cdot))$ propagate to the final estimator.
\begin{proposition}\label{proposition: partial double robustness for pos class estimator}
Define $R_{1,n} = \|\hat{\mu}_{1}(\cdot) - \mu_{1}(\cdot)\| \|\hat{\pi}_{1}(\cdot) - \pi_{1}(\cdot)\|$.
Assume (i) there exists $M < \infty$ such that $\| \beta_0(\cdot) \| \leq M$;
(ii) there exists $\delta > 0$ such that $\mathbb{P}(\pi_1(X_i) \geq \delta) = 1$;
(iii) there exists $\epsilon > 0$ such that $\mathbb{P}(\hat{\pi}_{1}(X_i) \geq \epsilon) = 1$;
(iv) $\|\hat{\mu}_{1}(\cdot) - \mu_{1}(\cdot)\| = o_P(1)$ and $\|\hat{\pi}_{1}(\cdot) - \pi_{1}(\cdot)\| = o_P(1)$;
and (v) $\widehat{\Delta}_n$ is non-empty with probability one.
Then,
\begin{align*}
\left\| \widehat{\overline{\mbox{perf}}}_{+}(s; \beta) - \overline{\mbox{perf}}_{+}(s; \beta) \right\| = O_\mathbb{P}\left( 1/\sqrt{n} + R_{1,n} + \sqrt{\mathbb{E}_{n}[(\widehat{\underline{\delta}}(X_i) - \underline{\delta}(X_i))^2]} + \sqrt{\mathbb{E}_{n}[(\widehat{\overline{\delta}}(X_i) - \overline{\delta}(X_i))^2]} \right).
\end{align*}
\end{proposition}
Proposition \ref{proposition: partial double robustness for pos class estimator} reveals the structure of how estimation errors affect our estimator.
The first term $O_{\mathbb{P}}(1/\sqrt{n})$ reflects sampling error that would arise even if the nuisance functions $(\mu_1(\cdot), \pi_1(\cdot))$ and the bounding functions $(\underline{\delta}(\cdot), \overline{\delta}(\cdot))$ were known.
The second term $O_{\mathbb{P}}(R_{1,n})$ reflects errors from estimating the nuisances $(\mu_1(\cdot), \pi_1(\cdot))$ that appear in the numerator and denominator of the linear fractional program.
These enter through a doubly robust product, and this behavior is inherited from the influence function structure underlying the numerator and denominator.
The final term captures errors from estimating the bounding functions $(\underline{\delta}(\cdot), \overline{\delta}(\cdot))$ that define the constraint set $\Delta$.
This error enters through the root mean square error.
This error can be controlled using standard nonparametric regression techniques, and researchers can trade off bias and variance in familiar ways.
For this reason, Proposition \ref{proposition: partial double robustness for pos class estimator} establishes a ``partial'' double robustness:
given estimated constraints, our estimator is doubly robust to errors in the objective nuisances $(\mu_1(\cdot), \pi_1(\cdot))$.
Constraint estimation errors enter additively through RMSE.
Achieving full double robustness --- where constraint estimation errors would also enter with product structure --- remains an open problem, which we discuss in Remark \ref{remark: path to full double robustness}.
Remark \ref{remark: double robustness for nonsharp bounds} describes a simple alternative approach that obtains fully doubly-robust estimators for non-sharp bounds on positive class performance.
Proposition \ref{proposition: partial double robustness for pos class estimator} and our estimator are practically valuable.
Researchers can use state-of-the-art machine learning methods for estimating the objective nuisances $(\mu_1(\cdot), \pi_1(\cdot))$, while separately controlling the constraint estimation error through their choice of estimator for the bounding functions.
The error decomposition applies to \textit{any} estimator of the bounding functions with favorable $L_2(\mathbb{P})$ properties, and we discuss alternative strategies and their tradeoffs in Section \ref{subsec: estimating bounds}.
This modularity is a key practical advantage of our framework.
\begin{remark}[Connection to existing work]\label{remark: connection to LFP literature}
Our analysis extends existing work on linear fractional programs, including \citet{AronowLee(2013)}, \citet{MiratrixWagerZubizaretta(18)}, \citet{KallusZhou(18)-confoundingrobust}, and \citet{ZhaoSmallBhattacharya(19)}.
These papers analyze linear fractional programs that arise when bounding population means or treatment effects under alternative sensitivity analysis models.
Our analysis advances in two key respects:
(i) we establish partial double robustness for objective nuisances -- a property not achieved in prior work; and
(ii) we provide a unified framework that accommodates multiple bounding strategies (observed outcome bounds, proxy outcome bounds, instrumental variable bounds) through a modular approach to constraint estimation.
\end{remark}
\begin{remark}[Paths toward full double robustness]\label{remark: path to full double robustness}
The proof of Proposition \ref{proposition: partial double robustness for pos class estimator} reveals that the optimizer in the linear fractional program can be characterized as having a threshold structure in $\beta_{0,i}$.
The twist is that the cutoff determining the threshold is itself not fixed but rather a functional of the estimated constraints.
Consequently, estimation error in the constraints affects not only which functions are feasible but also where the optimal threshold lies.
While related, recent work on conditional linear programs \citep[][]{benmichael2025partialidentificationconditionallinear} and aggregated intersection bounds \citep[][]{semenova2023adaptive} are not directly applicable.
While the linear fractional program can be reduced to a linear program through the Charnes-Cooper transformation, the transformation uses a global normalization that couples all observations together, whereas existing work exploits conditional separability across $X$-cells. We leave full resolution of this question for future work and view our partial double robustness result as an important first step.
\end{remark}
\begin{remark}[Estimation of non-sharp bounds]\label{remark: double robustness for nonsharp bounds}
Our estimator targets sharp bounds on positive class performance, but an alternative approach that restores double robustness at the cost of sharpness is possible.
The idea is straightforward: separately estimate the numerator and the denominator using overall performance estimators (Section \ref{section: estimating overall performance bounds}) and then combine via the delta method.
Concretely, define
$$
\widetilde{\overline{\mbox{perf}}}_{+}(s; \beta) = \frac{\sup_{\delta \in \Delta} \mathbb{E}[ \beta_{0,i} \mu_1(X_i) + \beta_{0,i} \pi_0(X_i) \delta(X_i) ] }{\inf_{\delta \in \Delta} \mathbb{E}[\mu_1(X_i) + \pi_0(X_i) \delta(X_i)]}.
$$
Clearly, $\widetilde{\overline{\mbox{perf}}}_{+}(s; \beta) \geq \overline{\mbox{perf}}_{+}(s; \beta)$. The numerator and the denominator can be separately estimated using the methods in Section \ref{section: estimating overall performance bounds}.
This delivers: (i) full double robustness with respect to all nuisances, and (ii) standard inference with straightforward variance estimation.
The cost is that we target non-sharp bounds.
Whether this trade-off is attractive depends on the application.
\end{remark}
\begin{remark}[Interpretation of non-empty constraint sets]\label{remark: interpertation of non-empty constraint sets}
Proposition \ref{proposition: partial double robustness for pos class estimator} assumes the estimated constraint set $\widehat{\Delta}_n$ is non-empty.
If finite-sample variability leads to $\widehat{\underline{\delta}}(x) > \widehat{\overline{\delta}}(x)$ for some $x \in \mathcal{X}$, researchers can either:
(i) clip the bounds to enforce feasibility, which might introduce bias captured by the root mean square terms in Proposition \ref{proposition: partial double robustness for pos class estimator}; or (ii) restrict attention to the feasible covariate region $\tilde{X}_n = \{x : \widehat{\underline{\delta}}(x) \leq \widehat{\overline{\delta}}(x)\}$, thus redefining the estimand.
The latter approach is analogous to trimming observations with extreme estimated propensity scores in causal inference, though this may complicate analysis \citep[e.g.,][]{HotzEtAl(09), Sasaki_Ura_2022}.
\end{remark}
\subsubsection{Estimation of the Bounding Functions}\label{subsec: estimating bounds}
The error decomposition in Proposition \ref{proposition: partial double robustness for pos class estimator} applies to any estimator of the bounding functions with favorable mean square error properties. This modularity gives researchers flexibility.
We discuss two natural strategies and their relative merits, focusing on observed outcome bounds for concreteness.
The simplest approach directly substitutes nuisance estimates into the bounding function formulas.
For observed outcome bounds, this means $\widehat{\underline{\delta}}(x) = (\underline{\Gamma} - 1) \widehat{\mu}_1(x)$ and $\overline{\delta}(x) = (\overline{\Gamma} - 1) \widehat{\mu}_1(x)$.
The mean square error properties of the estimated bounding functions are then inherited from those of $\widehat{\mu}_1(x)$.
This approach requires no additional estimation beyond the first-stage nuisances.
An alternative approach is to construct influence-function-based pseudo-outcomes and regress them on covariates using the full sample (i.e., our pseudo-outcome regression procedure in Section \ref{section: DR Learner, main text}).
We develop this approach below as it may offer bias-variance advantages in practice.
Ultimately, both strategies are valid under our framework.
To understand why pseudo-outcome regression may be advantageous, consider that $\widehat{\mu}_1(x)$ is typically constructed using only observations with $D_i = 1$.
Consequently, standard regression procedures target the mean square error conditional on $D_i = 1$
However, the error bound in Proposition \ref{proposition: partial double robustness for pos class estimator} depends on the \textit{unconditional} mean square error.
These two objectives could lead to different bias-variance tradeoffs when the selected population ($D_i=1$) is small relative to the full population or there is covariate shift between the selected and unselected populations.
In such settings, an estimator that directly targets the unconditional MSE may achieve better performance. Pseudo-outcome regression accomplishes this by constructing unbiased pseudo-outcomes using observations from both the selected and unselected populations, then regressing these pseudo-outcomes on the full sample.
\vspace{-1em}
\paragraph{Estimation Procedure for Observed Outcome Bounds:}
To illustrate how Proposition \ref{proposition: partial double robustness for pos class estimator} may be applied, we estimate bounding functions using our pseudo-outcome regression procedure. This illustrates how our partial double robustness result (Proposition \ref{proposition: partial double robustness for pos class estimator}) can be combined with analyses of nonparametric regression estimators --- in particular, leveraging our oracle inequality for pseudo-outcome regression (Appendix \ref{section: oracle inequality}).
We now randomly split the data into three folds (we assume $n$ is divisible by $3$ and each fold contains $n/3$ observations for simplicity).
We construct nuisance function estimates $\hat \eta$ using the first fold.
On the second fold, we construct the pseudo-outcomes $\phi_{\mu,i}(\hat \eta)$ and regress the estimated pseudo-outcomes on the covariates $X_i$ using a researcher-specified nonparametric regression procedure, yielding $\widehat{\delta}(\cdot)$.
On the third fold, we solve the sample analogue in Equation \eqref{equation: sample pos class LFP}.
Our analysis in Section \ref{section: DR Learner, main text} immediately allows us to analyze the errors from the estimated bounding functions.
\begin{corollary}\label{cor: DR learner rate for pos class}
Under the same conditions as Proposition \ref{proposition: partial double robustness for pos class estimator}, if further $\widehat{\mathbb{E}}_n[\cdot \mid X_i = x]$ satisfies the $L_2(\mathbb{P})$-stability condition (Assumption \ref{asm: L2 stability condition}), then
\begin{align*}
& \sqrt{\mathbb{E}_{n}[(\widehat{\delta}(X_i) - \mu_1(X_i))^2]} = O_{\mathbb{P}}\left( 1/\sqrt{n} + \epsilon^{-1} \| \tilde{R}(\cdot) \| + R_{oracle} \right),
\end{align*}
where $\tilde{R}(x) = \widehat{\mathbb{E}}_n[ (\pi_1(X_i) - \hat{\pi}_1(X_i)) (\mu_{1}(X_i) - \hat{\mu}_1(X_i)) \mid X_i = x ]$, and now $R_{oracle}^2 = \mathbb{E}[ \| \widehat{\delta}_{oracle}(\cdot) - \mu_1(\cdot) \|^2]$ for the oracle pseudo-outcome regression procedure.
\end{corollary}
\section{Empirical Application to Credit Risk Scores}\label{section: empirical application, main text}
We validate the finite-sample performance of our estimators through Monte Carlo simulations in Appendix \ref{section: monte carlo simulations}.
We now apply our framework to credit risk prediction, where financial institutions use scores $s(\cdot)$ to predict default risk.
We observe data on past loan applications but only observe defaults $Y_i = D_i Y_i^*$ for funded applications ($D_i = 1$).
It is implausible to assume unconfounded funding decisions. Applicants who accept loan offers may differ systematically from those who decline in ways that affect their default risk.\footnote{Lending institutions face two distinct inference problems. First, "policy rejects" with $\mathbb{P}(D_i = 1 \mid X_i) = 0$ (i.e., overlap violations) due to underwriting rules require extrapolation-based methods that are beyond our scope. Second, even among fundable applicants with $\mathbb{P}(D_i = 1 \mid X_i) > 0$, those accepting versus declining offers may differ in unobserved ways --- the selective labels problem we address. Extending our framework to handle policy rejects is an important direction for future research.}
We analyze 372,346 personal loan applications submitted to Commonwealth Bank of Australia from July 2017 to July 2019. These loans are repaid in monthly installments and used to purchase vehicles, refinance debt, or cover major expenses \citep[][]{CostonEtAl(21)-Rashomon}. Loan amounts ranged up to AU\$50,000 (median AU\$10,000) with a median offered interest rate of 14.9\% per annum. We observe rich application-level covariates including reported income, occupation, and credit history at CommBank. We only observe whether an applicant defaulted within 5 months $Y_i^* \in \{0, 1\}$ if the application was funded $D_i = 1$.
Approximately one-third of applications were funded, with a 2.0\% default rate among funded loans.
We consider two exercises.
First, we construct credit risk scores robust to varying confounding assumptions. The resulting algorithms produce substantially different risk rankings that could affect credit access.
Second, following recent work in consumer finance \citep[][]{BlattnerNelson(21), FusterEtAl(22)predictably, DiMaggioEtAl(22)-Fintech}, we evaluate how a benchmark score's accuracy varies across income groups.
While an analysis ignoring confounding suggests better performance for higher-income applicants, this disparity disappears under plausible confounding levels.
\subsection{Bounding Individual Default Risk}\label{section: DR-Learner, CBA application}
We construct robust credit risk scores using our estimator from Section \ref{section: DR Learner, main text}, exploring how varying confounding assumptions affect individual default risk predictions.
We estimate the upper bound on the conditional default probability with observed outcome bounds, $\overline{\mu}(\cdot)$.
We split applications into two folds: the first estimates $\pi_1(\cdot)$, $\mu_1(\cdot)$ using random forests; the second regresses estimated pseudo-outcomes on application-level covariates $X_i$ using cross-validated Lasso regression.
We set $\underline{\Gamma} = 1$ and vary $\overline{\Gamma} \in \{1, 2, 3\}$.
We compare our robust estimator (with $\underline{\Gamma} = 1, \overline{\Gamma} = 2$) against a benchmark score trained only on funded applications.
Figure \ref{figure: DR-Learner, CBA application}'s left panel shows the joint distribution of their predictions.
For each benchmark decile, the heatmap displays the percentage of applications falling into each decile of our robust predictions. Applications exhibit substantial reranking: among those in the benchmark's 5th decile, 18.7\% shift to extreme deciles (1-3 or 8-10) under the robust score.
These differences meaningfully affect credit access. For instance, if only the least risky third of applications receive funding, 10.1\% would have their decision reversed when comparing benchmark versus robust scores.
Figure \ref{figure: DR-Learner, CBA application}'s right panel examines how estimated coefficients vary with $\overline{\Gamma}$ (variable definitions in Online Appendix Table \ref{table: DR-Learner, CBA detailed variable description}).
Confounding assumptions affect different covariates heterogeneously.
Some coefficients—total net income and credit bureau score—remain zero across all $\overline{\Gamma}$ values.
Others, like the number of recent credit card applications, receive zero weight in the benchmark but non-zero weight for all $\overline{\Gamma}$>1.
Still others, such as occupation type and maximum recent delinquency, exhibit large magnitude changes as $\overline{\Gamma}$ varies, highlighting how confounding assumptions shape which applicant characteristics the algorithm prioritizes.
\subsection{Bounding the Predictive Performance of a Credit Risk Score}\label{section: robust audit, CBA application}
We now evaluate an existing credit risk score's performance.
We split applications into training and evaluation data.
On the training data, we construct a benchmark score predicting 5-month default risk among funded applications using random forests.
On the evaluation data, we analyze mean square error and ROC curves, following recent mortgage lending studies \citep[][]{BlattnerNelson(21), FusterEtAl(22)predictably}.
Using observed outcome bounds with $\underline{\Gamma} = 1$, we assess how the benchmark's estimated performance varies as we increase $\overline{\Gamma}$, allowing unfunded applicants to have progressively higher default rates than funded applicants with similar characteristics.
Figure \ref{figure: full population evaluation, CBA application}(a) shows how estimated MSE bounds vary with $\overline{\Gamma}$. When $\overline{\Gamma} = 1$ (the unconfounded case), both bounds equal the benchmark's MSE on funded applications (dashed red line; $0.143$). As $\overline{\Gamma}$ increases, allowing unfunded applicants higher default rates, bounds widen but remain informative.
At $\overline{\Gamma} = 2$ --- permitting unfunded applicants to default at twice the rate of funded applicants with similar characteristics --- the estimated bounds are $[0.138, 0.182]$.
Figure \ref{figure: full population evaluation, CBA application}(b) examines how estimated ROC curves vary with $\overline{\Gamma}$.
We summarize ROC curves using Area Under the Curve (AUC), computed via the trapezoidal rule.
Under unconfoundedness ($\overline{\Gamma} = 1$), the benchmark's AUC is $0.637$.
At $\overline{\Gamma} = 1.5$, AUC bounds are $[0.575, 0.705]$.
As with MSE, bounds widen with $\overline{\Gamma}$ but remain informative at empirically plausible confounding levels.
We finally investigate the benchmark credit risk score's predictive disparities across income groups.
We define $G_i \in \{0, 1\}$ to be whether an application is below or above the median personal income of all submitted applications.
We examine how the benchmark credit risk score's mean square error and ROC curve vary across these income groups, thereby studying whether it is a ``noisier'' signal of 5-month default risk for lower income applications, as in \cite{BlattnerNelson(21), FusterEtAl(22)predictably}.
Figure \ref{figure: evaluation by income group, CBA application}(a) compares the mean square error of the benchmark credit risk across income groups.
Interestingly, for $\overline{\Gamma} = 1$, the mean square error of the benchmark credit risk score is significantly larger on applications below the median income (0.169) than those above it (0.136).
This difference persists as we vary assumptions on unobserved confounding.
The upper bound on the mean square error for applications above the median income is only larger than the lower bound on the mean square error for applications below the median income for values $\overline{\Gamma} \geq 1.69$.
Ruling out this mean square error disparity across income groups would require that unfunded applications are at least 1.69 times as likely to default within 5 months as funded applications conditional on covariates.\footnote{Rather than examining disparity bounds across a range of $\overline{\Gamma}$ values, researchers could compute the ``identification breakdown point'' $\overline{\Gamma}^*$: the minimal value of $\overline{\Gamma}$ (holding $\underline{\Gamma}=1$) such that $0 \in H(\text{disp}(s;\beta))$.
A natural estimator is $\widehat{\overline{\Gamma}}^*$, the minimal $\overline{\Gamma}$ such that zero is contained in the estimated $(1-\alpha)$ confidence set $\widehat{H}(\text{disp}(s;\beta))$. Since our confidence sets have valid coverage by the results in Section \ref{section: estimating overall performance bounds} and Appendix \ref{section: identification and estimation of predictive disparities}, $P(\widehat{\overline{\Gamma}}^* \geq \overline{\Gamma}^*) \geq 1-\alpha$.}
Figure \ref{figure: evaluation by income group, CBA application}(b) also compares the ROC curves of the benchmark credit risk score across income groups.
Intriguingly, there exists effectively no disparity across income groups based on the ROC curves.
When $\overline{\Gamma} = 1$, the implied AUC for applications above the median income is $0.603$, whereas the implied AUC is $0.616$ for applications below the median income.
Any disparity in terms of ROC curves is already ruled out at $\overline{\Gamma} = 1.25$ as the bounds on the AUC for applications above the median income are $[0.567, 0.653]$ and $[0.581, 0.654]$ for applications below the median income.
\section{Conclusion}
We developed a framework for robust design and evaluation of predictive algorithms under unobserved confounding.
Our approach nests popular empirical strategies --- observed outcome bounds, proxy outcomes, and instrumental variables --- by bounding outcome differences between selected and unselected units (Assumption \ref{asm: bounding assumption}).
We provide estimators that work well across bounding strategies and performance measures: conditional likelihoods for algorithm design, overall metrics like MSE and calibration, and class-specific measures like TPR and FPR. This breadth makes our framework valuable across domains where predictive algorithms face selective labels—pretrial release, consumer lending, medical diagnosis, hiring, and child welfare.
Several open questions present promising research directions. First, achieving full double robustness for positive class performance bounds remains open, though our technical results suggest possible paths forward. Second, uniform inference for breakdown frontiers would let researchers formally test how strong confounding must be to overturn conclusions \citep[][]{MastenPoirier(20)}.
Third, we evaluate fixed algorithms; extending our results to settings where the same data is used to train and evaluate an algorithm is an important direction.
Finally, extending to multi-valued or continuous outcomes would broaden applicability, though defining class-specific performance in these settings requires care.
As algorithms increasingly drive high-stakes decisions in credit, criminal justice, and healthcare, robust evaluation under realistic confounding becomes essential.
Our framework and these open questions offer pathways to strengthen algorithmic accountability where the stakes are highest.
{\singlespacing
\bibliographystyle{chicago}
\bibliography{ref}}
\newpage \clearpage
\section*{Main Figures}
\begin{figure}[htbp]
\centering
\subfloat[Benchmark risk score vs. robust risk score.]{\includegraphics[width = 3in]{Figures/CBA-Figures/drlearner_conditional_heatmap_Gamma2.jpeg}} \hspace{2.5mm}
\subfloat[Coefficients in risk score as $\overline{\Gamma}$ varies.]{\includegraphics[width=3.25in]{Figures/CBA-Figures/drlearner_overall_coefplot.jpeg}}
\caption{Estimated personal loan credit risk scores as assumptions on unobserved confounding vary.}
\floatfoot{\textit{Notes}:
The left panel summarizes the joint distribution of a benchmark credit risk score's predictions of default risk against our robust predictions of the default risk.
Among applications in each decile of the benchmark credit score's predicted risk distribution, the left panel plots the percentage of applications at each decile of our robust estimator's predicted risk distribution.
Our robust estimator is constructed assuming $\underline{\Gamma} = 1$, $\overline{\Gamma} = 2$, and the benchmark credit score predicts default risk among only funded applications.
The right panel summarizes how the coefficients on a subset of application-level characteristics vary with assumptions on unobserved confounding over $\overline{\Gamma} \in \{1, 2, 3\}$.
The value $\overline{\Gamma} = 1$ corresponds to the benchmark credit score. See Section \ref{section: empirical application, main text} for further discussion.
}
\label{figure: DR-Learner, CBA application}
\end{figure}
\newpage
\begin{figure}[htbp!]
\centering
\subfloat[Mean square error]{\includegraphics[width=3in]{Figures/CBA-Figures/mse_overall_perf.jpg}} \hspace{2.5mm}
\subfloat[ROC curve]{\includegraphics[width=3in]{Figures/CBA-Figures/roc_overall_perf_v2.jpg}}
\caption{Bounds on mean square error and ROC curve of benchmark risk score as we vary the assumption on unobserved confounding.}
\floatfoot{\textit{Notes}: This figure summarizes how the bounds on the benchmark risk score's mean square error and ROC curve vary with our assumptions on unobserved confounding.
In Panel (A), the mean square error among only funded applications is plotted in the red dashed line.
The bounds (black) are constructed using observed outcome bounds, varying $\overline{\Gamma} \in [1, 2]$ and setting $\underline{\Gamma} = 1$.
In Panel (B), the ROC curve among only funded applications is plotted in black.
The bounds on the ROC curve for alternative choices $\overline{\Gamma} \in \{1.25, 1.33, 1.5\}$ are plotted in different colors. See Section \ref{section: empirical application, main text} for further discussion.
}
\label{figure: full population evaluation, CBA application}
\end{figure}
\newpage
\begin{figure}[htbp]
\centering
\subfloat[Mean square error]{\includegraphics[width=3in]{Figures/CBA-Figures/mse_perf_byIncome.jpg}} \hspace{2.5mm}
\subfloat[ROC curve]{\includegraphics[width=3in]{Figures/CBA-Figures/income_roc_overall_perf_v2.jpg}}
\caption{Bounds on mean square error and ROC curve of benchmark credit risk score across income groups as we vary the assumption on unobserved confounding.}
\floatfoot{\textit{Notes}: This figure summarizes how the bounds on the benchmark credit risk score's mean square error and ROC curve vary across income groups and as our assumptions on observed confounding vary.
In Panel (A), the mean square among only funded applications is plotted in the red dashed line. The bounds (black) are constructed using observed outcome bounds, varying $\overline{\Gamma} \in [0, 1]$ and setting $\underline{\Gamma} = 1$.
In Panel (B), the ROC curve among only funded applications is plotted in black.
The bounds on the ROC curve for alternative choices $\overline{\Gamma} \in \{1.25, 1.33, 1.5\}$ are plotted in different colors. See Section \ref{section: empirical application, main text} for further discussion.
}
\label{figure: evaluation by income group, CBA application}
\end{figure}
\clearpage
\newpage
\subfile{Supplement.tex}