1 Introduction↩︎

Machine learning models that achieve strong performance on a held-out test set routinely degrade when deployed in new environments. This degradation is well-documented across high-stakes settings: medical imaging models underdiagnose underrepresented demographic groups [1]–[5], clinical risk scores exhibit performance drops across hospitals and demographic groups [6]–[8], and natural language processing systems underperform on linguistic minorities and dialects [9]–[12]. These failures are often attributed to distribution shift [13] where the training and deployment populations differ systematically. However, in many decision-making settings, distribution shift is only part of the problem: labels are often observed only for individuals selected by a prior human, institutional, or algorithmic decision, which can make the evaluation data an inaccurate reflection of the deployment population.

We study the intersection of these evaluation challenges. The first is covariate shift, where the input distribution of features changes while the relationship between features and outcomes remains stable [14]. Under covariate shift, considering a model’s average performance on a held-out test set can obscure deployment failures; if a model’s error is higher in feature regions that are more common in the deployment setting, then its target risk may be much worse than indicated by the test error [15].

The second is the selective labels problem [16]. In domains such as lending, hiring, and healthcare, outcomes are observed only for individuals selected by a decision policy: approved borrowers, hired applicants, or patients who receive diagnostic medical tests. Such a mechanism induces selection bias because the observed labels are drawn from a decision-dependent subset of the population [17]. As a result, evaluating performance only on labeled samples from the training data can be misleading.

The challenges of covariate shift and selective labels often coexist in precisely the high-stakes settings where reliable pre-deployment evaluation matters most. A salient example is the use of Clinical Decision-making Instruments (CDIs), predictive models that use patient demographics, symptoms, and test results to aid diagnosis, triage, and treatment decisions. Many CDIs are developed using data from patient populations that differ systematically from the populations in which they are ultimately deployed– for example, models trained on predominantly white patient populations are often later used in demographically diverse populations [18]. This practice induces a covariate shift. At the same time, many clinically relevant outcomes are observed only for patients who receive follow-up testing, monitoring, or treatment, and those decisions are shaped by existing clinical policies or the CDI itself. Thus, to accurately determine whether a CDI is suitable for deployment in a new hospital, a practitioner must address not only that the target population may differ systematically from the population in which the model was evaluated but also that outcomes are observed only for patients selected by prior clinical decisions.

This work studies pre-deployment evaluation of an arbitrary black-box prediction model under the joint presence of covariate shift and selective labels. Our first contribution is to formalize the problem of evaluating black-box predictors under simultaneous covariate shift and selective labeling. We frame evaluation around a target risk functional that captures the risk the model would incur in the target population under a fixed general loss function. This estimand answers the practically important question of how a model will perform in its actual deployment population, even when that population looks different from the training data and when outcome labels are systematically missing due to historical decisions.

Second, we establish identification of the target risk using observable quantities in the data under standard assumptions. In particular, we establish identification supposing that the evaluator observes (i) selectively-labeled data from the source distribution (used to train the model) where selection is on observables and (ii) unlabeled samples from the target deployment environment.

Third, we derive the efficient influence function of the target risk estimand and use it to construct a double machine learning (DML) estimator that is asymptotically normal and semiparametrically efficient. Finally, we empirically validate our framework through controlled synthetic, semi-synthetic, and real-world experiments in healthcare settings; we show that our estimator tracks deployment performance more reliably than standard plug-in approaches or methods designed for evaluation under either covariate shift or selective labels in isolation.

1.1 Background and Related Work↩︎

Our work studies model evaluation in a setting that jointly addresses dataset shift and selective outcome observability, two challenges that are typically treated separately. We present a framework for pre-deployment evaluation of a fixed predictor when the deployment population differs from the training population and, in the source data, labels are observed only for a feature-dependent biased subset. Our work builds on the literatures of evaluation under distribution shift, learning and evaluation under selectively observed outcomes, and semiparametric, doubly robust estimation.

1.1.0.1 Covariate shift:

Our work addresses covariate shift [14], in which the marginal feature distribution \(P(X)\) differs between the training and deployment environments while the conditional distribution of outcomes given features \(P(Y | X)\) remains stable (see, e.g., [13], [19], [20] for broader surveys). Classical approaches to addressing covariate shift rely on importance weighting or density ratio estimation [14], [21], [22]. More recent work develops doubly robust procedures for estimation and evaluation under covariate shift [23]–[27].

Beyond methods for correcting covariate shift, a growing body of work addresses the problem of evaluating models under covariate shift [28]–[30] and assessing whether a given shift will degrade performance [31]–[34]. However, even seemingly benign covariate shifts at the population level among labeled samples may mask significant performance degradation for subpopulations whose labels are systematically unobserved in training data [35]. Our work addresses these evaluation and equity concerns that emerge when covariate shift and selective labels occur simultaneously.

1.1.0.2 Multiple forms of distribution shift:

In contrast to covariate shift, conditional shift occurs when the relationship between features and outputs \(P(Y | X)\) changes across environments [36]. We note that the selective labels problem can induce an apparent conditional shift by systematically hiding certain labels, causing the observed distribution to differ from the true conditional law \(P(Y | X)\) [16], [17]. However, our framework imposes additional structure that distinguishes our setting from the joint presence of covariate and conditional shifts [29], [37]: we assume selection on observables, which implies that the unlabeled samples remain informative about both the source covariate distribution and the selection process. This assumption marks a key methodological distinction relative to settings that simultaneously treat covariate and conditional shifts.

1.1.0.3 Selective labels and sample selection:

The selective labels problem [16], [17] arises when outcome labels are observed only for individuals selected by a historical policy or a human decision-maker. This occurs in settings such as credit scoring, criminal justice, and clinical medicine, where outcomes are only revealed after approval, release, or testing decisions [6], [38]. In such settings, the labeled data are not representative of the full population on which the model will ultimately be used, which complicates both training and evaluation [39].

Existing approaches include counterfactual risk estimation using causal inference [40], [41], correcting selection bias using heterogeneity across decision-makers [17], [42], and targeted data acquisition or augmentation for subpopulations underrepresented in the labeled data [35], [43]. Our setting is also connected to questions about feedback loops and performativity, where interventions guided by model predictions change the process by which future labeled data are generated [44], [45].

1.1.0.4 Missing data:

The selective labels problem is a special case of the broader challenge of missing data [46]. Under the assumption of selection on observables, selective labeling can be viewed as a missing at random (MAR) mechanism [47]. Classical approaches to estimation under MAR include inverse probability weighting [48], [49], outcome regression and augmentation [50], [51], and doubly robust estimation [52], [53]. Our approach is similar to the semiparametric missing-data literature in that it combines outcome regression with inverse-probability weighting to construct a doubly robust estimator. However, our setting departs from the standard MAR formulation by aiming to evaluate model performance in a target population with no observed outcomes. This distinction requires a procedure that simultaneously corrects for the selection and covariate transport mechanisms. More broadly, the challenge of missing data also connects to concerns about outcome measurement quality, where the ground-truth outcome is contested or incomplete [51], [54], [55].

1.1.0.5 Variance Reduction:

Our approach relates to methods that use auxiliary predictions or unlabeled samples to improve statistical efficiency. One similar perspective is that of control variates which leverages the known mean of a correlated random variable to reduce the variance of an estimator [56]–[58].

Next, in combining unlabeled samples with regression-based estimates for unlabeled samples, our work similarly relates to prediction-powered inference (PPI) [59]–[62]. PPI shares with our setting the goal of making use of unlabeled target data, though it does not address outcome missingness induced by historical decision policies or covariate shift. Our approach therefore relates to recent work that considers generalized PPI estimators to settings of missingness and distribution shift [63], [64].

1.1.0.6 Double Machine Learning:

Our evaluation framework builds on the literature on doubly robust (DR) and double machine learning (DML) estimators [50], [52], [53]. These methods are well-suited to settings with incomplete or biased data due to the fact that the resulting estimators remain consistent in parametric settings if either the propensity score or outcome model are correctly specified [48], [51], [65]. They enjoy fast rates of convergence in nonparametric settings [66] and have been extended to debiased estimation of more general functionals using deep learning methods [67]. These methods have also been used for estimation in settings closely related to selective labels and dataset shift, e.g., policy learning [68], distribution shift [25], [28]–[30], [69], and data missingness [70], [71].

1.1.0.7 Algorithmic Fairness:

Our setting is also connected to the literature on algorithmic fairness, since both covariate shift and selective outcome observability can systematically misrepresent a model’s performance for subpopulations which are underrepresented in the labeled training data [6], [15], [72]. For instance, models trained on populations from well-resourced institutions and deployed in under-resourced settings may systematically underperform for patient populations with different characteristics [1]–[3], [73].

Moreover, when outcome observability is shaped by historical decisions, labeled data may systematically exclude the subpopulations for whom fairness harms are most significant [16], [17], [39], [40]. This gap presents an important fairness concern: without pre-deployment evaluation methods, practitioners cannot assess a priori whether their models will perpetuate or exacerbate inequities when deployed.

1.1.0.8 Notation:

We use \(X \indep Y\) to denote that two random variables \(X\) and \(Y\) are independent. For a measurable function \(h\), we write \(\norm{h}_{L_2(P)}\coloneq\left(\int h(x)^2 dP(x) \right)^{1/2}\) for its \(L_2(P)\) norm. We also write \([n] \coloneq \{1, \ldots, n\}\). For sequences of random variables, \(o_p(\cdot)\) and \(O_p(\cdot)\) denote convergence in probability and stochastic boundedness, respectively. Let \(\sigma(t) = (1 + \exp(-t))^{-1}\) denote the sigmoid function. We denote by \(\mathbb{1}\{A\}\) the indicator of the event \(A\). \(\norm{\cdot}_p\) is the \(\ell_p\) norm. \(I_m\) is the \(m \times m\) identity matrix, and \(\mathbb{1}_m\) is the \(m\)-dimensional all ones vector. Finally, \(\overset{d}{\to}\) denotes convergence in distribution.

2 Problem Setting↩︎

We study the problem of evaluating a fixed prediction model under the joint presence of covariate shift and selective labeling. We observe \(n\) independent and identically distributed (i.i.d.) draws \[\label{eqn:sample} Z \mathrel{\vcenter{:}}= (X, R, R \cdot D, R \cdot D \cdot Y)\tag{1}\] where \(X \in \mathcal{X} \subset \RR^m\) is a covariate vector, \(R \in \{0, 1\}\) is a domain indicator, and \(D \in \{0, 1\}\) indicates whether the outcome of interest is observed. The domain indicator \(R\) satisfies \(R = 1\) for units from the source population and \(R = 0\) for units from the target population. Due to selective labeling, the latent outcome \(Y \in \RR\) is observed only when \(R = 1\) and \(D =1\).

Let \(p_S\) and \(p_T\) denote the densities of \(X|R = 1\) and \(X|R=0\), respectively. For an integrable random variable \(Q\), we write \(\EE_S[Q] \coloneq \EE[Q | R = 1], \EE_T[Q] \coloneq \EE[Q | R = 0]\). When \(Q\) is a function of only \(X\), these reduce to expectations under the covariate distributions \(P_S\) and \(P_T\), respectively.

We aim to evaluate how well a fixed prediction model \(f:\mathcal{X} \to \RR\) predicts \(Y\) in the target environment. For a given loss function \(\ell: \RR \times \RR \to \RR_{\geq 0}\) (e.g., squared loss), our estimand is the target risk \[\label{eqn:target95risk} \psi \mathrel{\vcenter{:}}= \EE_T\left[ \ell \left( f(X), Y \right) \right].\tag{2}\] Note that \(f\) may be an arbitrary fixed predictor (e.g., a deep network, an ensemble model, an LLM scorer), and \(\ell\) may be an arbitrary fixed loss function. Our procedure treats both as black boxes.

3 Identification and Estimation of the Target Risk↩︎

In this section, we identify \(\psi\) under the joint presence of covariate shift and selective labels, and we then derive a novel estimator based on this identification. We use \(L \mathrel{\vcenter{:}}= \ell(f(X), Y)\) as shorthand notation for the loss under \(f\) on observed outcomes. Also define the following nuisance functions: \[\label{eqn:nuisance951} \rho \mathrel{\vcenter{:}}= \PP(R=0), \quad g(X) \mathrel{\vcenter{:}}= \PP(R=0 \mid X),\quad \mu(X) \mathrel{\vcenter{:}}= \EE\!\left[\,L \mid X, R=1, D=1\,\right],\tag{3}\] and the pooled distribution and source distribution propensity scores, respectively, as \[\pi(X) \coloneq \PP(D = 1, R = 1 \mid X), \quad \pi_S(X) \mathrel{\vcenter{:}}= \PP(D=1 \mid X, R = 1).\] We assume throughout that both the source and target domains are represented, i.e., that \[\label{eqn:rho95bound} 0 < \rho < 1.\tag{4}\]

Estimation of \(\psi\) is complicated both by the fact that \(Y\) is unobserved in the target domain (\(R = 0\)) and by the fact that, in the source domain, \(Y\) is selectively observed (when \(D = 1, R = 1\)). To ensure identifiability, we impose the following assumptions, which are standard in the causal inference and transfer learning literatures.

These assumptions identify the target risk from observed data. Intuitively, these conditions require that (i) label observability in the source population is as good as random conditional on \(X\), (ii) the conditional relationship between covariates and outcomes is invariant across the source and target domains, (iii) every covariate profile admits a positive probability of being observed in the source population, and (iv) the source and target distributions’ supports overlap sufficiently. ²

[prop:identification95target95risk] is the central identification result of the paper. To our knowledge, it is the first identification of an evaluation functional that simultaneously transports across populations and corrects for a selective labeling process; existing covariate shift identifications presume fully labeled source data, and existing identifications for selectively labeled data presume a single population. This result lays the foundation for our proposed evaluation procedure by connecting the unobservable target risk to quantities that can be learned from observable, pre-deployment data.

The first equality shows that the target risk can be written as an expectation over the target covariate distribution of the nuisance function \(\mu(X)\), which summarizes the model’s conditional expected loss after accounting for selective labels. We note that \(\mu(X)\) depends only on covariates, i.e., it does not require observing target outcome labels. This expression demonstrates that pre-deployment evaluation is feasible in our setting; if we could estimate \(\mu(\cdot)\) from selectively labeled source data, we could estimate \(\psi\) by averaging \(\mu(X)\) over the target population.

The second expression provides an alternative identification as an expectation with respect to the source population: the first factor \(\tfrac{p_T(X)}{p_S(X)}\) transports the loss estimate from the source to the target population to account for covariate shift, and the second factor \(\tfrac{D}{\pi_S(X)}\) is an inverse propensity weighted (IPW) correction for selective labels. Together, these equalities formalize the insight that, even when deployment labels are unavailable and source labels are observed only through a selective labeling process, we can still evaluate a model’s target risk from pre-deployment data. A proof is provided in Appendix 8.1.

3.1 Efficiency Theory for Target Risk Estimation↩︎

Next, we provide a novel theoretical result that characterizes the efficient influence function (EIF) of \(\psi\). The EIF enables us to construct an estimator that can achieve \(\sqrt{n}\) rates while allowing for flexible nuisance estimation. Let \(\mathcal{P}\) denote the nonparametric model defined by Assumptions 1-4. The following result identifies the efficient influence function of the target risk \(\psi\) and gives the corresponding von Mises expansion. This result demonstrates the bias-correction step that stabilizes the naı̈ve plug-in estimators.

The EIF given in [prop:eif] incorporates two complementary signals from the source and target distributions, respectively, in each term. The first term of 6 is a residual correction computed on observed labeled source units (\(RD = 1\)). Here, \(RD/\pi(X)\) is the Horvitz-Thompson (inverse probability) weight for appearing as a labeled source observation in the pooled sample. This factor corrects the selective-labeling mechanism. The additional factor \(g(X) /\rho\) transports the correction toward the target covariate law. Indeed, if \(p(x) = \rho \cdot p_T(x) + (1 - \rho) \cdot p_S(x)\) denotes the pooled covariate density, then Bayes’ rule gives \(\frac{g(X)}{\rho} = \frac{p_T(X)}{p(X)}\).

The residual \(L - \mu(X)\) centers the labeled source contribution conditional on \(X\), which reduces variance relative to the uncentered representation. The second term is the target domain plug-in component. Since labels (and, therefore, \(L\)) are not observed in the target domain (\(R = 0\)), this term averages the conditional loss regression \(\mu(X)\) over unlabeled target covariates. The factor \((1-R)/\rho\) converts the pooled expectation into an expectation in the target domain, and subtraction of \(\psi(\mathbb{P})\) centers the influence function (since \(\psi = \EE_T[\mu(X)]\)). Thus, the EIF combines a transported labeled-source residual correction with a centered target-domain plug-in term.

The proof of [prop:eif] applies semiparametric calculus to the first identified representation given in  [prop:identification95target95risk]. Differentiating this functional yields the two terms of 6 . The von Mises remainder is obtained by comparing \(\psi(\overline{\PP})-\psi(\PP)\) to the first-order expansion, then collecting and controlling leftover terms, which are the products of nuisance functions. We then invoke [75] to conclude efficiency. A full proof is provided in 8.2.1.

3.2 A novel double machine learning estimator of target risk↩︎

We use the efficient influence function of \(\psi\) (Proposition [prop:eif]) to construct a double machine learning (DML) estimator of the target risk [50], [52], [53]. Our construction uses an EIF-based bias correction that makes the DML estimator robust to errors in the nuisance functions \(\mu, \pi\) and \(g\). The DML approach combines bias correction with sample splitting (or cross-fitting) to remove the bias that would otherwise emerge due to using the same data for nuisance estimation and evaluation.

Let \(\mathcal{S}_n \coloneq \{Z_i\}_{i=1}^n\) denote an evaluation sample of \(n\) i.i.d. draws from \(P\), and let \(\widehat{\pi}\), \(\widehat{g}\), \(\widehat{\mu}\) be nuisance estimators fit on an auxiliary sample \(\mathcal{S}'_n\) consisting of \(n\) samples and independent of \(\mathcal{S}_n\). Motivated by the EIF in 6 , we propose the DML estimator of \(\psi\) as \[\label{eqn:dml95estimator} \widehat{\psi} = \frac{1}{n} \frac{1}{\widehat{\rho}} \sum_{i=1}^n \left[ \frac{R_i D_i }{\widehat{\pi}(X_i)} \widehat{g}(X_i) \left( L_i - \widehat{\mu}(X_i) \right) + (1 - R_i) \widehat{\mu}(X_i)\right]\tag{7}\] where \(\widehat{\rho}=n^{-1}\sum_{i=1}^n\mathbb{1}\{R_i = 0\}\) is the empirical estimator of \(\rho\). The first summand corrects labeled source units (\(R_i = 1, D_i = 1\)) for the selection and covariate shift mechanisms; the second is a regression-based plug-in estimate for target units (\(R_i =0\)).

1 given below establishes that \(\widehat{\psi}\) is \(\sqrt{n}\)-consistent, asymptotically linear, and semiparametrically efficient under standard regularity conditions and mild convergence-rate requirements for nuisance function estimation. For simplicity, we state the results under single splitting and note that analogous guarantees for cross-fitting follow under analogous foldwise conditions.

A more general statement of 1, which formalizes conditions (i)-(iv), is given in 8.4, along with the proof. The main implication of this result is that our proposed estimator supports valid pre-deployment inference, even when the nuisance functions are estimated with flexible machine learning methods. In particular, the conditional loss model \(\mu\), the pooled labeled-source propensity \(\pi\), and the domain propensity \(g\) need not be estimated at parametric rates. It suffices that their errors are small enough for the second-order remainder to be negligible; a sufficient condition is \(o_p(n^{-1/4})\) convergence in \(L_2(P)\) for each nuisance function.

The theorem also highlights why the influence function correction is essential for \(\sqrt{n}\) inference in our setting. Naı̈ve plug-in or IPW estimators generally inherit first-order bias from errors in \(\widehat{\mu}\), \(\widehat{\pi}\), or \(\widehat{g}\), making valid inference difficult. By contrast, the DML correction removes first-order sensitivity to nuisance estimation error: after sample splitting, the remaining bias is governed by products of nuisance errors.

The asymptotic linear expansion in ?? also yields a practical uncertainty quantification procedure. Since the leading term is the average of the influence scores \(\varphi(Z_i; P)\), the asymptotic variance can be estimated by the empirical variance of the estimated influence scores. This enables the construction of standard Wald-type confidence intervals for the target risk. Therefore, the theorem justifies the use of modern prediction methods to estimate the nuisance components while retaining classical \(\sqrt{n}\) inference for the target risk.

4 Experiments↩︎

We present two sets of clinical experiments. First, we perform a semi-synthetic experiment with real-world covariates and simulated outcome and label-observability mechanisms. This setting preserves realistic covariate shifts, which are known to be more complex and difficult to address than shifts simulated in controlled, synthetic experiments [15]. Simulating the label and outcome mechanisms enables comparison of each estimator to the true target risk. Second, we evaluate the estimator on real-world clinical outcomes, where both label observability and outcomes are taken from observed eICU data. Synthetic experiments are deferred to 9. The eICU Collaborative Research Database [76] includes de-identified individual-level electronic health records from over 200,000 admissions to ICUs across multiple hospitals in the United States. We use admission-level patient demographics, vitals, clinical unit type, and hospital ID.

4.1 Real-world Covariate Shifts↩︎

We leverage the fact that the data include multiple hospitals to capture the dynamics of real-world distribution shifts where a model is trained on a population that differs in demographic makeup from the population on which it is deployed. We identify 3 candidate target hospitals that look systematically different from the general population in age, race, and ethnicity and one baseline source hospital with demographics that are representative of the overall eICU sample. We train a model on the source hospital, then evaluate its performance in each candidate, covariate-shifted target hospital.

In particular, our covariate-shifted hospitals include hospital 443 (which tends to have younger patients and more African American patients), hospital 199 (with a typical age profile and more Caucasian patients), hospital 283 (which skews older with a larger share of patients labeled as unknown or “other” ethnicity). We use hospital 208 as a no-covariate-shift benchmark since it is approximately average across these demographics. We illustrate these covariate shifts in Figure 6.

4.2 Semi-Synthetic Design↩︎

We generate outcomes from a sparse logistic model in which older patients have a higher outcome risk. We then train a fixed prediction rule \(f\) on simulated labels from the source training split and hold this rule fixed throughout. We simulate label observability for source patients. The probability that a source label is observed depends on patient covariates, where the selection mechanism is chosen so that selection strength values increasingly distort label observability away from older patients. We vary feature-dependence in the selection mechanism while holding the overall labeling rate approximately constant. As a result, the experiment varies the degree of feature dependence in the selective labeling mechanism while holding the amount of labeled source data roughly constant. The target hospital outcomes are treated as unobserved by the estimators, which matches the pre-deployment evaluation setting in which the evaluator accesses unlabeled target covariates but not target labels. See 11.3 for details.

4.3 Estimators↩︎

For each selection strength value, we compare our DML estimator 7 to three baselines. First, a natural benchmark is the Plug-in Estimator that reweighs observed labeled source losses by the estimated density ratio \(\widehat{w}(x)\) and the inverse propensity weights \(1/\widehat{\pi}(x)\) to account for both covariate shift and selective labels. Second, we compare against a recent approach for evaluating models under covariate shift using double machine learning that we denote Covariate Shift (CS)-only [26]. Finally, we compare against a Selective Labels (SL)-only estimator from [40] that uses double machine learning to address selective labels. The CS-only and SL-only baselines represent leading approaches for their respective settings: each uses double machine learning to leverage flexible nuisance estimation, but each addresses only one of the two challenges considered here. See 11.1 for more details.

4.4 Nuisance Estimation↩︎

We estimate three nuisance components using cross-fitting: (i) a domain classifier for covariate shift reweighting, (ii) a label observability propensity for selective labels, and (iii) an outcome regression model for the conditional expected loss.

First, we estimate the domain propensity \(g(x)\) for covariate shift reweighting via logistic regression on the pooled source and target covariates. Second, we estimate the source label observability propensity \(\pi_S(x)\) via logistic regression. Then, noting that \(\pi(X) = (1 - g(X)) \cdot \pi_S(X)\), we estimate the source outcome-labeling propensity via \(\widehat{\pi}(x) = (1-\widehat{g}(x)) \widehat{\pi}_S(x)\). Finally, we estimate the conditional loss \(\mu(x)\) by noting that for absolute loss with binary outcomes, we can write \(\mu(x) = \eta(x) |1-f(x)| + (1-\eta(x))|f(x)|\) where \(\eta(x) \coloneq \PP(Y = 1|X=x, R = 1, D = 1)\). We use a cross-fitted random forest regression model to estimate \(\eta(x)\) with labeled source units, then use this estimate to construct \(\widehat{\mu}(x)\).

4.5 Evaluation and Uncertainty↩︎

To evaluate each estimator, we compute a Monte Carlo estimate of the ground-truth target risk by simulating an oracle dataset of target samples. For each target covariate vector, we generate an outcome using the same outcome model as the experiment. The ground-truth risk estimate is then computed as the mean absolute error (MAE) between the model prediction and the oracle outcomes. We evaluate each estimator by comparing its estimated target risk to the oracle estimate of ground truth. We report uncertainty via 95% bootstrap confidence intervals. While our proposed DML estimator is an asymptotically linear estimator of the target risk under the conditions of 1, the same need not be true in general for the other estimators. Consequently, we acknowledge that the bootstrap intervals may be inaccurate for such estimators; we nevertheless use the same uncertainty quantification procedure across estimators to provide a consistent comparison.

4.6 Semi-Synthetic Results↩︎

2 reports the bias and RMSE of each estimator relative to the ground-truth target risk as the feature-dependence of the selective labeling mechanism increases on the x-axis. Across the shifted target hospitals, our DML estimator exhibits the least bias. The CS-only estimator is positively biased in all target hospitals, with bias increasing as selection becomes more feature-dependent.

As expected, the SL-only estimator is approximately unbiased in the no-shift case (hospital 208), but exhibits bias in the other shifted settings. The plug-in estimator also exhibits significant bias. Together, these results show that methods addressing only one source of bias, as well as naı̈ve attempts to combine corrections, can systematically misrepresent deployment performance when covariate shift and selective labels occur simultaneously.

The RMSE results show that our DML method consistently has among the lowest RMSE. In some settings (hospital 443), CS-only performs comparably to our method, which has significantly lower error than the plug-in and SL-only estimators. In others (hospital 199), SL-only yields similar error to our method whereas the CS-only and plug-in estimators exhibit significantly higher error. For hospital 283, due to wide uncertainty intervals, we cannot determine which method has the lowest error. The RMSE results suggest that our DML estimator provides a reliable estimate of performance under covariate shift and selective labels.

Overall, the combined bias and RMSE results illustrate that, when covariate shift and selective labeling are both present, correcting for either source of bias in isolation can give a misleading estimate of deployment risk.

4.7 Real-world eICU Experiments↩︎

We next illustrate our estimator in a fully real-world setting using eICU data. As above, we use hospital membership to observe naturally occurring covariate shift. However, now the labeling mechanism and the outcomes are taken directly from the observed clinical data.

We define binary outcomes from laboratory measurements obtained during the first 12–24 hours of ICU admission. Each outcome corresponds to a clinically meaningful threshold value that reflects an underlying complication, such as a heart condition. The prediction model uses patient demographics, admission characteristics, and ICU unit types as covariates. We train a fixed logistic regression prediction rule on labeled observations from hospital 208. We estimate its target risk in three covariate-shifted hospitals, 443, 199, and 283.

The key distinction from the semi-synthetic experiments is that the selective labeling is no longer simulated. Instead, label observability is based on real clinical decision-making: we only observe a given outcome for patients whose clinicians ordered the relevant lab in a specified admission window. Thus, the labeling indicator captures real variation in diagnostic and monitoring practices across hospitals, and, hence, heterogeneity in the selective labeling mechanism. At the same time, hospital membership captures naturally occurring covariate shifts in patient severity, characteristics, and demographics.

We do not observe labels for patients whose labs are not ordered. Thus, we cannot report bias or RMSE against an oracle risk as in the semi-synthetic setting. Instead, we report estimated target risks and interpret the implications of using each estimator in a clinical context.

We evaluate our estimator on four lab-dependent clinical outcomes in the eICU data: indicators of elevated B-type natriuretic peptide (BNP), elevated troponin-I, elevated creatine kinase (CK), and hypoxemia. For each outcome, the label is observed only among ICU patients for whom the relevant laboratory test was ordered within a given time window (24 hours for BNP, troponin-I, CK, and 12 hours for hypoxemia). This observability mechanism exactly reflects the real-world selective labels problem where the outcome is only observed when a clinician orders the corresponding diagnostic test. See 11.4 for further discussion on these outcomes.

4.8 Real-world Results↩︎

3 reports target-risk estimates for the selectively-labeled clinical outcomes across target hospitals and a benchmark no-shift hospital. First, we observe that for some outcomes (BNP, hypoxemia), the estimated target risk changes substantially across hospitals while, for others (troponin-I, CK), the estimates are more stable.

In some instances, the DML, CS-only, and SL-only estimators approximately agree (e.g., BNP 208, troponin-I 199, CK 208 and 283, hypoxemia 283). In others, however, the DML aligns with the CS-only estimator and diverges from the SL-only estimator (e.g., hypoxemia 199), indicating that the covariate shift dominates the risk adjustment. At the same time, the instances of disagreement between the CS-only and DML estimators (e.g., CK 443, hypoxemia 443 and 283) reflect that accounting for selection can meaningfully change the practitioner’s expectation of and confidence about the risk of deploying a clinical model in a new hospital.

We also note the relatively wide intervals for the plug-in estimator across hospitals and outcomes. This underscores the efficiency gains offered by our DML estimator over the naı̈ve direct inverse weighting approach, as reflected in both 1 and the semi-synthetic outcome experiments.

5 Robustness Checks↩︎

The identification result in [prop:identification95target95risk] relies on selection on observables, covariate shift, positivity, and bounded source-target overlap (1- 4). We therefore run an additional set of controlled synthetic experiments to evaluate how each estimator behaves as these assumptions are increasingly violated.

In each experiment, we perturb a single assumption at a time while holding fixed the rest of the data generating process. We consider four violations: (i) the presence of an unobserved variable that jointly affects label observability and outcomes, (ii) conditional outcome distribution shift, (iii) loss of label positivity in a tail region of the source distribution, and (iv) increasing source-target likelihood ratio instability. Full details of the perturbation constructions are given in 11.5.

4 reports bias and RMSE as the severity of each assumption’s violation increases. When selection (and, therefore, label observability) depend on an unobserved variable that also affects outcomes, all estimators become increasingly biased; The DML and plug-in estimators perform best with respect to both bias and RMSE for mild to moderate violations, while larger violations eventually degrade all estimators. Under conditional outcome shift in the target distribution (i.e., where \(P_S(Y | X) \neq P_T(Y | X)\)), all estimators become increasingly biased; the DML and plug-in estimators perform similarly, and outperform all other estimators across the parameter range. Under positivity stress, the DML, plug-in, and CS-only estimators remain unbiased, whereas the SL-only estimator becomes increasingly biased. The DML and CS-only estimators also have the lowest RMSE. Finally, under likelihood ratio stress, the DML estimator remains the least biased and has the lowest RMSE across the parameter range.

Overall, these robustness checks suggest that our proposed estimator is reasonably robust to positivity and likelihood ratio stress. Moreover, while the DML estimator exhibits increasing bias under selection on an unobserved variable and under conditional outcome shift, it is competitive with respect to bias and RMSE in comparison to all other benchmark estimators.

6 Conclusion↩︎

This paper proposes a new framework for pre-deployment evaluation that, to our knowledge, is the first to account jointly for covariate shift and selective labels. Our approach formalizes target risk as the model’s expected deployment performance, establishes conditions for its identification under observed data, and derives an influence-function-based double machine learning estimator. We prove that the proposed estimator is asymptotically linear, which yields asymptotic normality and supports Wald-type inference. In semi-synthetic eICU experiments, our estimator more accurately tracks target risk than standard plug-in procedures and evaluation procedures designed for each problem in isolation.

These results highlight the importance of developing tools that can account for multiple coexisting data challenges. In particular, the combination of covariate shift and selective labels, each of which has been studied extensively in isolation, poses compounded difficulties and is likely to arise in high-stakes domains such as healthcare.

These results come with important qualifications. Identification relies on substantive assumptions, including selection on observables and sufficient source-target overlap. These may fail, for instance, when clinical testing decisions depend on unmeasured severity indicators, or when target patients lie in regions of the covariate space that are impossible to observe in labeled source samples. Our robustness checks in 5 suggest some tolerance to small violations, but future work could develop methods that relax these assumptions or use adaptive data collection to improve overlap.

Accounting for co-occurring data challenges is especially important in fairness and safety assessments. By applying our method to subgroup risks, evaluators can assess fairness or safety properties under deployments characterized by covariate shift and selective labeling.

Finally, our framework has broader applications beyond clinical prediction. Many foundation model and LLM evaluation settings involve selectively observed feedback while the deployment distribution of prompts, users, and downstream tasks often differs substantially from that of the training and testing. Extending our approach to these settings could support more principled evaluation of modern ML systems under the real-world distortions that often arise together.

7 Acknowledgments↩︎

We are grateful to Edward Kennedy for insightful comments. Annie Ulichney’s work is supported in part by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 2146752 and in part by the Cooperative AI PhD Fellowship from the Cooperative AI Foundation. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the National Science Foundation.

8 Proofs↩︎

8.1 Proof of Proposition [prop:identification95target95risk]↩︎

Recall the notation \(L \coloneq \ell(f(X), Y)\). By law of iterated expectations, \[\label{eq:identification95proof951} \psi = \EE_T\left[ L \right] = \EE_T \left[ \EE\left[ L | R = 0, X \right]\right].\tag{8}\] By Assumption 2, \[\label{eq:identification95proof952} \EE_T \left[ \EE\left[ L | R = 0, X \right]\right] = \EE_T \left[ \EE\left[ L | R = 1, X \right]\right].\tag{9}\] By Assumption 1, \[\label{eq:identification95proof953} \EE_T \left[ \EE\left[ L | R = 1, X \right]\right] = \EE_T \left[ \EE\left[ L | R = 1, D = 1, X \right]\right].\tag{10}\] Combining 8 10 and recalling the definition of \(\mu\) yields the first representation \[\label{eq:identification95proof954} \psi = \EE_T\left[ \mu(X) \right].\tag{11}\]

To show the second representation, we start from 11 . Using the change of measure justified by Assumption 4, we obtain: \[\label{eq:identification95proof955} \psi = \int_{x \in \mathcal{X}} \mu(x) p_T(x) dx = \int_{x \in \mathcal{X}} \mu(x) \frac{p_T(x)}{p_S(x)} p_S(x)dx = \EE_S\left[ \frac{p_T(X)}{p_S(X)} \mu(X) \right].\tag{12}\] Under the source law, conditioning on \(R = 1\) is implicit by construction. Hence, \[\label{eq:identification95proof956} \mu(X)=\EE[L \mid X, R=1, D=1]=\EE_S[L \mid X, D=1].\tag{13}\] Since \(D \in \{0, 1\}\), we have \[\label{eq:identification95proof957} \EE_S\left[ D L | X \right] = \PP(D = 1 | X, R = 1) \EE_S[L | X, D = 1] = \pi_S(X) \mu(X).\tag{14}\] By Assumption 3, rearrange 14 to get \[\label{eq:identification95proof958} \mu(X) = \frac{\EE_S\left[ D L | X \right] }{\pi_S(X)}.\tag{15}\] Substituting 15 into 12 yields \[\psi = \EE_S\left[ \frac{p_T(X)}{p_S(X)} \frac{\EE_S\left[ D L | X \right]}{\pi_S(X)} \right].\] Applying the tower property yields the representation \[\psi = \EE_S\left[ \frac{p_T(X)}{p_S(X)} \frac{ D L }{\pi_S(X)} \right]\] which proves the claim. \(\qed\)

8.2 Candidate Influence Function Derivation↩︎

The following lemma recalls well-known results characterizing the influence functions of conditional expectation and density functions. See, e.g., [75].

8.2.1 Proof of Proposition [prop:eif]↩︎

Following the semiparametric calculus of [75], we treat \(\mathcal{X}\) as a discrete set, apply Gateaux differentiation separately to each of \(\mu(\cdot)\) and \(p_T(\cdot)\), and invoke the product rule for influence functions: \[\IF{\psi} = \sum_{x \in \mathcal{X}} \IF{\mu(x)} p_T(x) + \sum_{x \in \mathcal{X}} \mu(x) \IF{p_T(x)}.\]

Applying the building block influence functions ?? and ?? given in Lemma 1 together with Bayes’ Rule, we obtain: \[\IF{\psi} = \frac{R \cdot D }{\pi(X)} \frac{g(X)}{\rho} \left( L - \mu(X) \right) + \frac{(1 - R)}{\rho} \left(\mu(X) - \psi \right)\] as claimed. By 3 detailed below, we see that the von Mises expansion in ?? has a linear term \(\varphi(\cdot; \PP)\) and remainder \(R_2(\PP, \overline{\PP})\) that is the product of nuisance functions. Therefore, by [77], we conclude that \(\varphi\) is the efficient influence function. \(\qed\)

8.3 von Mises Expansion↩︎

8.4 Formal Statement and Proof of Theorem 1↩︎

First, we provide a formal general statement of  1.

8.4.1 Proof of Theorem 2↩︎

Throughout, we suppress the arguments of all nuisance functions for notational convenience. We adopt the notation \[\PP_n f = \frac{1}{n}\sum_{i=1}^n f(Z_i), \quad P f = \int f(z) dP(z) = \EE_P[f(Z)]\] where \(P\) denotes the pooled law of \(Z\). We assume throughout that \(\widehat{\pi}, \widehat{g}, \widehat{\mu}\) are estimated on an auxiliary sample \(\mathcal{S}_n'=\{Z_i'\}_{i=1}^n\), sampled independently from the evaluation sample \(\mathcal{S}_n = \{Z_i\}_{i=1}^n\). Also define \[\widehat{\rho} = \frac{1}{n} \sum_{i=1}^n \mathbb{1}\{R_i=0\}.\]

We first relate the source distribution positivity condition 3 with a positivity condition on the pooled labeled source propensity \(\pi(X)\). Since \[\pi(X) = \PP(R = 1 | X) \pi_S(X),\] by Bayes’ rule and 4, we have that \[\PP(R = 1 | X = x) = \frac{(1-\rho)p_S(x)}{(1-\rho)p_S(x) + \rho p_T(x)} \geq \frac{1 - \rho}{1 - \rho + \rho C}.\] Then, by 4 and 3, we have \[\label{eqn:pooled95positivity} \pi(X) > {\varepsilon }\frac{1 - \rho}{1 - \rho + \rho C}\eqqcolon {\varepsilon}'\tag{20}\] almost surely.

Next, we decompose the left-hand side ?? . Recall the DML estimator 7 : \[\widehat{\psi}=\frac{1}{\widehat{\rho}}\PP_n\left( \frac{RD}{\widehat{\pi}} \widehat{g}(L-\widehat{\mu}) + (1-R) \widehat{\mu} \right).\] Subtracting \(\psi(P)\) from both sides and using \(\widehat{\rho}=\PP_n(1-R)\) gives \[\widehat{\psi}-\psi(P) = \frac{\rho}{\widehat{\rho}} \PP_n \varphi(\cdot; P) + \frac{1}{\widehat{\rho}} \PP_n B.\] where \[B \coloneq RD \left( \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi} \right) (L-\mu) - RD \frac{\widehat{g}}{\widehat{\pi}} (\widehat{\mu} - \mu) + (1-R)(\widehat{\mu}-\mu).\] Rearranging, \[\label{eqn:estimator95theoretical95bounds95rearranged} \widehat{\psi}-\psi(P) = \PP_n \varphi(\cdot; P) + \left( \frac{\rho}{\widehat{\rho}} - 1 \right) \PP_n \varphi(\cdot; P) + \frac{1}{\widehat{\rho}} (\PP_n-P) B + \frac{1}{\widehat{\rho}} PB.\tag{21}\] We proceed by controlling each term of 21 separately.

8.4.1.1 Preliminary bounds:

Define the event \[\mathcal{E} \coloneq \left\{\inf_{x \in \mathcal{X}} \widehat{\pi}(x) \geq {\varepsilon}'/2, \widehat{\rho} \geq \rho/2, \widehat{g}(x) \in (0, 1) \forall x \in \mathcal{X} \right\}.\] By Assumption ?? , \(\PP(\mathcal{E}) \to 1\). On \(\mathcal{E}\), we have \[\label{eqn:event95vareps95bounds} \frac{1}{\widehat{\rho}} \leq \frac{2}{\rho}, \quad \left| \frac{\rho}{\widehat{\rho}} - 1 \right| \leq \frac{2}{\rho} |\widehat{\rho}-\rho|, \quad \norm{ \frac{1}{\widehat{\pi}}}_{\infty} \leq \frac{2}{{\varepsilon}'}, \quad \norm{ \frac{\widehat{g}}{\widehat{\pi}}}_{\infty} \leq \frac{2}{{\varepsilon}'}.\tag{22}\] Moreover, since \(\widehat{\rho} = \frac{1}{n} \sum_{i = 1}^n \mathbb{1}\{R_i = 0\}\) is the sample mean of i.i.d. \(\text{Bernoulli}(\rho)\) random variables, we have that \[\label{eqn:estimator-guarantees-rho-bound} |\widehat{\rho}- \rho| = O_p(n^{-1/2}).\tag{23}\]

Also, on \(\mathcal{E}\), we have the bound \[\label{eqn:195minus95pi95hp} \left| 1 - \frac{\pi}{\widehat{\pi}}\right| = \frac{|\pi-\widehat{\pi}|}{\widehat{\pi}} \leq \frac{2}{{\varepsilon}'}|\pi-\widehat{\pi}|.\tag{24}\]

8.4.1.2 First term:

First, we note that we can write \[P \varphi = \frac{1}{\rho} \EE\left[ g(X) \EE\left[ \frac{RD}{\pi} (L-\mu) | X\right] \right] + \frac{1}{\rho} \EE\left[(1-R)(\mu - \psi) \right].\] For the first term, by Law of Iterated Expectation and the definition of \(\mu\), \[\label{eqn:estimator95theoretical95bounds95L95mu95mean95zero} \EE\left[ \frac{RD}{\pi} (L-\mu) | X\right] = \EE[L - \mu | X, R = 1, D = 1] = 0.\tag{25}\] For the second term, by Law of Iterated Expectation and Law of Total Probability \[\EE\left[(1-R)(\mu - \psi) \right] = \EE\left[(\mu(X) - \psi)\EE\left[1-R|X \right] \right] = \EE\left[ g(X) (\mu - \psi) \right] = \rho \EE_T\left[ \mu - \psi \right] = 0.\] Hence, \(P \varphi = 0\) and it follows that \[\PP_n \varphi(\cdot; P) = (\PP_n - P) \varphi(\cdot; P).\] By ?? , \(\EE[\varphi(Z; P)^2] < \infty\), so, by the central limit theorem, \[\label{eq:estimator95guarantees95first95term95final95bound} \PP_n \varphi(\cdot; P) = O_p(n^{-1/2}).\tag{26}\]

8.4.1.3 Second term:

By 26 and 22 , we get \[\left( \frac{\rho}{\widehat{\rho}} - 1 \right) \PP_n \varphi(\cdot; P) = O_p(n^{-1/2}) O_p(n^{-1/2}) = o_p(n^{-1/2}).\]

8.4.1.4 Third term:

First note that, since the nuisance functions are constructed on the independent sample \(\mathcal{S}_n'\), \[\EE\left[ (\PP_n-P) B \mid \mathcal{S}_n'\right] = 0.\] Next, on \(\mathcal{E}\), by 22 and triangle inequality, \[|B| \leq RD \left| \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}\right| |L - \mu| + RD\left| \frac{\widehat{g}}{\widehat{\pi}}\right| |\widehat{\mu}-\mu | + (1-R) |\widehat{\mu}-\mu | \leq RD \left| \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}\right| |L - \mu| + \left( 1 + \frac{2}{{\varepsilon}'}\right) |\widehat{\mu}-\mu |.\] Then, by Cauchy-Schwarz, \[\label{eq:estimator95property95variance95term} P(B^2) \leq 2 P\left( \left( \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}\right)^2 RD (L-\mu)^2 \right) + 2 \left( 1 + \frac{2}{{\varepsilon}'}\right)^2 P\left( (\widehat{\mu}-\mu )^2 \right).\tag{27}\] For the first term of 27 , by iterated expectation and ?? , \[\begin{align} P\left( \left( \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}\right)^2 RD (L-\mu)^2 \right) &= \EE\left[ \left( \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}\right)^2 \EE\left[RD (L-\mu)^2 | X \right]\right]\\ &= \EE\left[ \left( \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}\right)^2 \pi(X)\EE\left[ (L-\mu)^2 | X, R = 1, D = 1 \right]\right] \\ &\leq M \norm{\frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}}_{L_2(P)}^2 \end{align}\] Moreover, by 22 and 20 , on \(\mathcal{E}\), \[\left| \frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi} \right|\leq \frac{2}{{\varepsilon}'}|\widehat{g}-g| + |g|\left|\frac{1}{\widehat{\pi}} - \frac{1}{\pi}\right|\leq \frac{2}{{\varepsilon}'}|\widehat{g}-g| + \frac{2}{{{\varepsilon}'}^2}|\widehat{\pi}-\pi|.\] Therefore, \[\norm{\frac{\widehat{g}}{\widehat{\pi}} - \frac{g}{\pi}}_{L_2(P)} \leq \frac{2}{{\varepsilon}'} \norm{\widehat{g}-g}_{L_2(P)} + \frac{2}{{\varepsilon}'^2} \norm{\widehat{\pi}-\pi}_{L_2(P)} = o_p(1)\] where the final equality applies ?? . We use ?? once again to evaluate the second term of 27 : \[P[(\widehat{\mu}-\mu)^2] = o_p(1).\] All together, we have \[P(B^2) = o_p(1).\] Hence, the conditional variance of the empirical process term on \(\mathcal{E}\) is bounded as \[\operatorname{Var}\left((\PP_n-P) B \mid \mathcal{S}_n' \right) = \frac{1}{n} \operatorname{Var}(B(Z)|\mathcal{S}_n') \leq \frac{1}{n} \EE[B^2|\mathcal{S}_n'] = \frac{1}{n} P[B^2] = o_p(n^{-1}).\] It follows that, by Chebyshev’s inequality, \[(\PP_n-P) B = o_p(n^{-1/2}).\] By 22 we conclude \[\frac{1}{\widehat{\rho}} (\PP_n-P) B = o_p(n^{-1/2})\]

8.4.1.5 Fourth term:

As in 25 , \[\EE\left[ \frac{RD}{\pi} (L-\mu) | X\right] = 0.\] Thus, \[PB = - \EE \left[ RD \frac{\widehat{g}}{\widehat{\pi}} (\widehat{\mu} - \mu)\right] + \EE \left[(1-R)(\widehat{\mu}-\mu)\right].\] Using law of iterated expectation combined with \(\EE[RD|X] = \pi(X)\) and \(\EE[1-R|X]=g(X)\), we obtain \[\begin{align} PB &= - \EE \left[ \frac{\widehat{g}}{\widehat{\pi}} \pi (\widehat{\mu} - \mu)\right] + \EE \left[g (\widehat{\mu}-\mu)\right] = \EE \left[ \left( g - \widehat{g} \frac{\pi}{\widehat{\pi}}\right) (\widehat{\mu} - \mu)\right]. \end{align}\]

Therefore, on the event \(\mathcal{E}\), \[\begin{align} |PB| &\leq \norm{\widehat{\mu}-\mu}_{L_2(P)}\norm{g - \widehat{g} \frac{\pi}{\widehat{\pi}}}_{L_2(P)}\\ &\leq \norm{\widehat{\mu}-\mu}_{L_2(P)}\left( \norm{g - \widehat{g}}_{L_2(P)} + \norm{\widehat{g}\left(1 - \frac{\pi}{\widehat{\pi}}\right)}_{L_2(P)} \right)\\ &\leq \norm{\widehat{\mu}-\mu}_{L_2(P)}\left( \norm{g - \widehat{g}}_{L_2(P)} + \frac{2}{\varepsilon'} \norm{\widehat{\pi} - \pi }_{L_2(P)} \right) \end{align}\] where the first inequality uses Cauchy-Schwarz, the second line follows by triangle inequality, and the last line uses ?? , 22 , and 24 . By ?? and 22 , we reach \[\frac{1}{\widehat{\rho}} |PB| = o_p(n^{-1/2}).\]

8.4.1.6 Collecting bounds:

Returning to 21 , we have established so far that \[\left( \frac{\rho}{\widehat{\rho}} - 1 \right) \PP_n \varphi(\cdot; P) = o_p(n^{-1/2}), \quad \frac{1}{\widehat{\rho}} (\PP_n-P) B = o_p(n^{-1/2}), \quad \frac{1}{\widehat{\rho}} |PB| = o_p(n^{-1/2}).\] Hence, \[\widehat{\psi}-\psi(P) = \frac{1}{n} \sum_{i = 1}^n \varphi(Z_i; P) + o_p(n^{-1/2})\] which establishes asymptotic linearity ?? . Finally, since \(\EE[\varphi(Z; P)^2] < \infty\), applying the central limit theorem with Slutsky’s theorem gives ?? as claimed. \(\qed\)

9 Synthetic Experiments↩︎

We demonstrate the performance of our proposed DML estimator and compare it to the existing doubly robust pre-deployment evaluation methods designed to evaluate model performance under either covariate shift or selective labels.

9.0.0.1 Synthetic Data Generation:

In each Monte Carlo iteration, we draw \(n_S\) source covariates and \(n_T\) target Gaussian covariates, each \(m\)-dimensional: \(X_i^{(S)} \sim \NN(\mu_S, \Sigma_S), \quad i \in \{1, \ldots n_S\}\), \(X_j^{(T)} \sim \NN(\mu_T, \Sigma_T), \quad j \in \{1, \ldots n_T\}\). We simulate covariate shift by varying \(\mu_S \neq \mu_T\) and/or \(\Sigma_S \neq \Sigma_T\). We pool the covariates into a single dataset \(\{(X_k, R_k)\}_{k \in [n]}\) where \(n \coloneq n_S + n_T\), \(R_k = 1\) denotes a source unit, and \(R_k = 0\) denotes a target unit.

9.0.0.2 Labeling and Outcome Models:

For source units \((R = 1)\), label observability follows a covariate-dependent Bernoulli mechanism with an overlap floor: \[\label{eqn:labeling95sythetic} D | (X, R = 1) \sim \operatorname{Bern}(\pi_S(X)), \quad \pi_S(X) = \pi_{\min} + (1-2\pi_{\min}) \sigma(-\alpha^\top X).\tag{28}\] The vector \(\alpha\) controls the direction and strength of the selective labels mechanism, where the sign convention in 28 makes labels less likely to be observed in the direction \(\alpha\). In the fixed-selection experiments, we set \(\alpha = c \cdot \mathbb{1}_m\) for a fixed, experiment-dependent scalar \(c\). In experiments that vary selection strength, we sweep over values of \(c\), with larger values corresponding to stronger dependence of label observability on the covariates. The constant \(\pi_{\min} \in (0, 1/2)\) ensures that Assumption 3 is satisfied with \(\varepsilon = \pi_{\min}\).

We generate continuous latent outcomes according to \[\label{eqn:outcome95sythetic} Y = \xi(X) + \varepsilon, \quad \varepsilon \sim \NN(0, \sigma_\varepsilon^2),\tag{29}\] where the conditional mean is defined relative to the fixed prediction model \(f\). In particular, we construct the conditional mean as \(\xi(x) = f(x) + \nu(x),\) such that, under squared loss, \(L = (f(X)-Y)^2\), the conditional expected loss satisfies \(\mu(x) = \EE[L|X = x] = \nu(x)^2 + \sigma_\varepsilon^2\). In other words, this construction allows us to specify where the prediction model \(f\) performs poorly through \(\nu(x)\). To understand the robustness of our framework to misspecification, we take \(\nu(x)\) to be nonlinear throughout our experiments and refer the reader to Appendix 11.2 for details.

9.1 Synthetic Experiment Results↩︎

The bias and Root Mean Squared Error (RMSE) of each estimator across four experiments, each varying either the magnitude of covariate shift, selection strength, or the ratio of samples from the source and target distributions are reported in Figure 5.

In the mean-shift experiment (Figure [fig:synthetic-mean-bias], Figure [fig:synthetic-mean-rmse]), the target distribution mean is shifted increasingly from the source mean. We construct this mean shift so that it is in the direction of increasing prediction error, i.e., in the direction of \(\nu(x)\). This means that, as the mean shift magnitude grows, the target risk increases. We observe that SL-only and CS-only estimators grow increasingly biased for the target risk. The plug-in and DML estimators remain approximately unbiased. However, the DML estimator has lower RMSE than the other estimators, especially for larger shifts.

The covariance-shift experiments yield the same qualitative conclusion, except that the SL-only estimator decreases in bias as the shift increases. However, we note that, as reported in Figure 7, this is merely a product of the fact that the SL-only estimate of risk is a constant, initially upward biased estimate of the increasing true target risk.

As selection increasingly depends on covariates, the CS-only estimator deteriorates because it treats the selectively labeled source sample as representative after covariate transport. The SL-only estimator remains biased because it corrects selective labels only in the source distribution. Both DML and the plug-in remain approximately unbiased, but DML once again achieves uniformly smaller RMSE.

Finally, the target/source ratio experiment varies the amount of unlabeled target data while holding both covariate shift and selective labeling fixed. The CS-only and SL-only estimators exhibit non-vanishing bias as the sample size increases. The DML and plug-in estimators remain approximately unbiased over the interval, while the DML maintains lower RMSE than the plug-in.

10 Additional Figures↩︎

11 Experimentation Details↩︎

11.1 Estimator Details↩︎

Here we provide precise formulae for the estimators with which we compare the DML estimator.

11.1.0.1 DML Estimator:

We also compute the DML estimator for the target risk by implementing 7 . \[\widehat{\psi}_{\rm{DML}} = \frac{1}{n} \sum_{k = 1}^n \left\{ \frac{R_k \cdot D_k }{\widehat{\pi}(X_k)} \cdot \frac{\widehat{g}(X_k)}{\widehat{\rho}} \cdot (L_k - \widehat{\mu}(X_k) ) + \frac{1 - R_k}{\widehat{\rho}} \cdot \widehat{\mu}(X_k) \right\}.\]

11.1.0.2 Naı̈ve (Plug-in) Estimator:

A natural benchmark is the plug-in estimator that directly combines observed labeled source losses. The observed losses are reweighted by the estimated density ratio \(\widehat{w}(x)\) and the inverse propensity weights \(1/\widehat{\pi}(x)\) to account for both the covariate shift and selective labels. We compute: \[\widehat{\psi}_{\text{plug-in}} = \frac{1}{n} \sum_{k = 1}^{n} \frac{R_k \cdot D_k}{\widehat{\pi}(X_k)} \cdot \frac{\widehat{g}(X_k)}{\widehat{\rho}} \cdot L_k.\]

11.1.0.3 Covariate Shift Benchmark:

To isolate the importance of correcting for selective labels, we implement the doubly robust estimator for evaluation under covariate shift from [26]. Their estimator addresses the setting where the source and target covariate distributions differ according to Assumption 2; however, it assumes that outcomes are observed for all source units. We adapt their estimator to our setting where labels are only observed for labeled source units (i.e., where \(R_i = 1, D_i = 1\)): \[\widehat{\psi}_{\rm{CS}} = \frac{1}{n_T} \sum_{i = 1}^{n} \left\{ \frac{\widehat{g}(X_i)}{1 - \widehat{g}(X_i)} \cdot R_i D_i \cdot (L_i - \widehat{\mu}(X_i)) + (1 - R_i) \widehat{\mu}(X_i) \right\}.\] Their estimator corrects for covariate shift via the density ratio \(\tfrac{\widehat{g}(X_i)}{1 - \widehat{g}(X_i)}\), but it does not account for selective labels.

11.1.0.4 Selective Labels Benchmark:

To isolate the importance of correcting for covariate shift, we implement the doubly robust counterfactual risk estimator from [40]. Their estimator addresses the setting where outcome labels are only observed for labeled units and uses inverse propensity weighting to account for the selective labeling mechanism. However, their setting assumes that the training and deployment populations share the same covariate distribution. We use their estimator to evaluate the risk on the source distribution using selectively labeled source samples: \[\widehat{\psi}_{\rm{SL}} = \frac{1}{n_S}\sum_{i = 1}^{n_S} \left\{ \frac{D_i}{\widehat{\pi}_S(X_i)} (L_i - \widehat{\mu}(X_i)) + \widehat{\mu}(X_i) \right\}.\]

11.2 Synthetic Experiments↩︎

This appendix provides the full specification of the synthetic experiments. In all synthetic experiments, we fix the dimension as \(m = 5\). In the mean-shift, covariance shift, and selection-strength shift experiments reported in Figures [fig:synthetic-mean-bias], [fig:synthetic-mean-rmse], [fig:synthetic-cov-bias], [fig:synthetic-cov-rmse], [fig:synthetic-selection-bias], [fig:synthetic-selection-rmse], we fix \(n_S = 1000, n_T = 4000\). We use \(n_{\rm{oracle}} = 50,000\) Monte Carlo iterations to estimate the target risk and implement \(N_{\rm{it}} = 200\) iterations for each setting in the grid. Across all experiments, we evaluate the fixed, non-linear prediction model \(f(x) = 0.30 x^\top u + 0.20 x_1\) where \(x_1\) denotes the first coordinate of \(x\) and \(u \coloneq \left(\frac{\mathbb{1}_m}{\norm{\mathbb{1}_m}_2}\right)\). First, we take \(\nu(x) = 1.80 \sigma(1.50 x^\top u) + 0.50 \sin(x_1) - 0.60\). We generate outcomes according to \(Y = \xi(x) + \varepsilon\) where \(\varepsilon \sim \NN(0, \sigma_\varepsilon^2)\) and \(\xi(x)=f(x) + \nu(x)\). Thus, the gap between the prediction rule and the conditional mean outcome given \(X = x\) is specified by the nonlinear choice \(\nu(x)\). Also note that we take \(\sigma_\varepsilon = 0.15\) and \(\pi_{\min} = 0.005\) across all synthetic experiments. We use \(3\)-fold cross-fitting to fit nuisance parameters across all synthetic experiments.

11.2.0.1 Increasing mean-shift:

Fix \(\mu_S = (0, \ldots, 0) \in \RR^m\), \(\Sigma_T = \Sigma_S = I_m\), \(\alpha = 1 \cdot u\) (i.e., \(c = 1\)), \(\beta = u\). For \(s \in [0.0, 0.222, 0.444, 0.667, 0.889, 1.111, 1.333, 1.556, 1.778, 2.0]\), we take \(\mu_T = \mu_S + s \cdot u\).

11.2.0.2 Increasing covariance-shift:

Fix \(\mu_S = \mu_T = (0, \ldots, 0) \in \RR^m\), \(\Sigma_S = 3I_m\), \(\alpha = 1 \cdot u\). For \(c \in [0.5, 0.667, 0.833, 1.0, 1.167, 1.333, 1.5, 1.667, 1.833, 2.0]\), we take \(\Sigma_T(x) = I_m + (c-1)u u^\top\).

11.2.0.3 Increasing selection strength:

Fix \(\mu_S=(0,\ldots,0)\in\RR^m\), \(\Sigma_S=I_m\), \(\Sigma_T=I_m\), and \(\mu_T = \mu_S + 0.5u\). For \(a \in [0.0, 0.222, 0.444, 0.667, 0.889, 1.111, 1.333, 1.556, 1.778, 2.0]\), take \(\alpha = a\cdot u\).

11.2.0.4 Increasing Target/source sample size ratio:

Fix \(n_S=1000\) and vary \[\frac{n_T}{n_S} \in [0.05, 0.1, 0.25, 0.5, 1, 2, 5, 10, 20, 50, 100].\] For each ratio \(r\), we set \(n_T=\mathrm{round}(1000r)\) and fix \(\mu_S=(0,\ldots,0)\in\RR^m\), \(\mu_T=\mu_S+0.75u\), \(\Sigma_S=I_m\), \(\Sigma_T=I_m+(1.25-1)uu^\top\), and \(\alpha = a\cdot u\).

11.3 Semi-synthetic Experiments↩︎

11.3.0.1 Data:

We construct a covariate vector with \(m = 26\) patient-level demographics and admission variables. Specifically, we include indicators for gender, ethnicity, age bin, and ICU unit type, together with admission height, weight, and body mass index (BMI). Age is binned into categories \(0\)–\(19\), \(20\)–\(39\), \(40\)–\(59\), \(60\)–\(79\), \(80\)–\(89\), and \(90+\). The continuous variables are median imputed to remove non-finite values and then standardized.

The fixed prediction rule \(f\) is trained once on a training split from hospital 208 using 1095 observations and evaluated on 730 test observations from a test split. The oracle risk is computed over a held-out target evaluation sample of sizes 1829, 2120, 1049, 1301 for hospitals 443, 199, 283, and 208, respectively.

11.3.0.2 Target and source hospitals:

Though eICU hospitals are de-identified, the data report hospital-level metadata which we give here for additional context on the real-world shifts of our experiments. All hospitals are in the United States. Hospital 208 is a large non-teaching hospital in the South, hospital 443 is a large teaching hospital in the South, hospital 199 is a large teaching hospital in the Northeast, and hospital 283 is a mid-sized non-teaching hospital in the Midwest.

11.3.0.3 Setup:

Across semi-synthetic experiments, we generate binary outcomes according to the logistic model \(Y | X = x \sim \operatorname{Bernoulli}(\sigma(x^\top \beta))\). The coefficient vector \(\beta \in \RR^{26}\) is sparse. In particular, for the coordinates corresponding to the covariates age \(60\)–\(79\), age \(80\)–\(89\), age \(90+\), BMI, and weight, we set the coordinates of \(\beta\) to be \(2.5, 4.0, 5.0, 0.5, 0.3\), respectively, with all other entries \(0\). This makes it so that older patients are higher risk with some additional variation by weight and BMI.

To train a fixed \(f\), we draw binary labels on the training split held out for the prediction model and fit a random forest classifier with maximum depth \(6\) and minimum leaf size \(10\). This classifier is held fixed across all Monte Carlo iterations and all values of the selection strength sweep. In the shifted hospital experiments, the fixed model is the one trained on Hospital 208, which we take to be the source hospital.

For the labels, we construct \(\alpha \in \RR^{26}\) as follows. First, set the coordinates corresponding to age \(80\)–\(90\) and \(90+\) of \(\alpha'\) to be \(-2, -2.5\), respectively, with all other entries 0, then set \(\alpha\) to be the unit-norm vector in the same direction as \(\alpha'\). Then, set \(s(x) = x^\top \alpha + (x^\top \alpha)^2 - \frac{1}{n_S} \sum_{i: R_i = 1} (X_i^\top \alpha)^2\). For selection strength \(a\), the labeling probability for \(X = x\) is \(e_a(x) = \pi_{\min} + (1-2 \pi_{\min}) \sigma(C_a + a s(x))\) where \(\pi_{\min}=0.05\) and \(C_a\) is an intercept calibrated for each \(a\) such that the labeling probability is approximately held constant at \(0.3\). Then, we draw \(D_i |(X_i, R_i = 1) \sim \operatorname{Bernoulli}(e_a(x))\).

We sweep \(a \in \left[ 0, 0.214,\,0.429,\,0.643,\,0.857, 1.071, 1.286, 1.500 \right]\). We run \(20\) Monte Carlo iterations per selection strength value. The target estimand is mean absolute prediction error. Since the outcome model is known in the semi-synthetic setting, we evaluate the oracle target risk via \(\lambda_{\textrm{true}} = \frac{1}{N_\textrm{oracle}} \sum_{i \in \mathcal{S}_\textrm{oracle}} \left(\sigma(x^\top \beta)|f(X_i) - 1| + (1-\sigma(x^\top \beta))|f(X_i)|\right)\) where \(\mathcal{S}_\textrm{oracle}\) denotes the oracle hold-out set. All nuisance functions are fit with 5-fold cross-fitting. The domain and source classifiers are trained via logistic regression.

For the loss regression, we estimate the conditional outcome probability \(\widehat{p}(x)\) via a random forest regression with 300 trees, max depth 4, minimum leaf size 50. This is trained with labeled source observations and we then plug the estimate into the conditional absolute loss formula \(\widehat{p}(x)|f(x)-1|+(1-\widehat{p}(x))|f(x)|\).

11.3.0.4 Bias and RMSE:

For a fixed setting across all synthetic and semi-synthetic experiments, let \(\widehat{\psi}_{e, i}\) denote the target risk estimate from estimator \(e\) in Monte Carlo iteration \(i \in \{1, \ldots,N_{\textrm{it}}\}\). Let \(\psi^{\textrm{oracle}}\) denote the corresponding oracle estimator of the target risk. We define Bias and RMSE for estimator \(e\) as: \[\textrm{Bias}_e = \frac{1}{N_{\textrm{it}}} \sum_{i = 1}^{N_{\textrm{it}}} (\widehat{\psi}_{e, i} - \psi^{\textrm{oracle}}), \quad \textrm{RMSE}_e = \sqrt{\frac{1}{N_{\textrm{it}}} \sum_{i = 1}^{N_{\textrm{it}}} (\widehat{\psi}_{e, i} - \psi^{\textrm{oracle}})^2}.\]

11.4 Real-world Outcome Experiments↩︎

11.4.0.1 Outcome clinical significance:

Here we provide clinical context for outcomes of interest in our real-world outcome experiments. We select these outcomes because they are clinically meaningful, yet not universally observed for patients in the ICU: the outcomes are unobserved unless, for instance, a clinician suspects a condition and subsequently orders further tests. Elevated BNP is an indication of cardiac stress. Troponin-I is a marker of damage to the heart muscle, and is an indicator of myocardial injury, sepsis, or shock, among other conditions. CK is an enzyme released in skeletal muscle injury, heart damage, infections, rhabdomyolysis, and other serious conditions. PaO2 measures the pressure of oxygen dissolved in the plasma. Low PaO2 characterizes hypoxemia which can reflect underlying conditions such as lung conditions and heart issues.

11.4.0.2 Outcome details:

We consider four binary outcomes corresponding to lab results. Each is defined using measurements taken shortly after ICU admission. In particular, we define the elevated BNP, elevated troponin-I, elevated CK, and hypoxemia outcomes, respectively as: \[\begin{align} Y_{\textrm{BNP}} &= \mathbb{1}\{\max(\mathrm{BNP}) > 500 \text{ pg/mL within 24 hours of ICU admission}\},\\ Y_{\textrm{trop}} &= \mathbb{1}\{\max(\mathrm{troponin-I }) > 1.0 \text{ ng/mL within 24 hours of ICU admission}\},\\ Y_{\textrm{CK}} &= \mathbb{1}\{\max(\mathrm{CK}) > 500 \text{ U/L within 24 hours of ICU admission}\},\\ Y_{\textrm{hyp}} &= \mathbb{1}\{\min(\mathrm{PaO}_2) < 60 \text{ mmHg within 12 hours of ICU admission}\}. \end{align}\]

11.5 Robustness Check Details↩︎

The robustness check experiments use the same notation as the synthetic experiments. The fixed prediction rule to be evaluated is \(f(x) = 0.30 x^\top u + 0.20 x_1\) where \(u = \mathbb{1}_m \in \RR^m\). We take \(\nu(x) = 1.90 \sigma(1.50 x^\top u) + 0.50 \sin(x_1) - 0.60\).

11.5.0.1 Selection on observables violation:

To evaluate the robustness to 1, we introduce a latent variable \(U \sim \NN(0, 1)\) that affects both source label observability \(D\) and the outcome \(Y\). In particular, we draw \(D | (X, U, R = 1) \sim \operatorname{Bernoulli}\left(\pi_{\min} + (1-2\pi_{\min}) \sigma(-\alpha^\top X - cU) \right)\) and \(Y = f(x) + \nu(x) +cU + \varepsilon\) where \(\varepsilon \sim \NN(0, \sigma_\varepsilon^2)\). The scalar \(c\) is the parameter we sweep, with \(c = 0\) corresponding to selection on observables. As \(c\) increases, label observability and outcomes jointly increasingly depend on \(U\).

11.5.0.2 Bounded likelihood-ratio stress:

To stress 4, we modify the target distribution to include a mixture component from a distant region in the covariate space. In particular, we take \(X^{(T)} \sim (1-\lambda) \NN(\mu_T, \Sigma_T) + \lambda \NN(\mu_T + au, \Sigma_T)\), where \(u\) is the unit vector in the direction of increasing loss, \(a\) is a fixed offset parameter, and \(\lambda\) is the robustness parameter that we sweep. Note that \(\lambda =0\) corresponds to the synthetic experiment setting, and, as \(\lambda\) increases, more target mass is placed in a region with low source density. This increases the likelihood ratio.

11.5.0.3 Covariate shift violation:

To evaluate robustness to 2, we incorporate a target-only conditional mean shift. In particular, we suppose that target outcomes experience a mean shift such that the target conditional expected loss is \(\mu(x | R = 0) = \left( \nu(x) + \left( \gamma x^\top u\right) \right)^2 + \sigma_\varepsilon^2\). The scalar \(\gamma\) is the parameter we sweep, where \(\gamma = 0\) corresponds to the covariate shift assumption where \(P_S(Y | X) = P_T(Y | X)\).

11.5.0.4 Positivity stress:

To stress 3, we enforce that a tail region of the source distribution has no observed labels. For a fraction \(q\) of tail, we set \(D = 0\) in an extreme lower-tail region of the source law. We sweep \(q\), where higher values cause more instability for the inverse-propensity weights.

11.5.0.5 Additional risk estimate plot:

Note that 9 reports the corresponding estimated target risk curves for the same four experiments whose results are reported in 4. The red curve denotes the oracle target risk computed from MC estimates as in the synthetic experiments.

References↩︎