Evaluating Treatment Benefit Predictors using Observational Data: Contending with Identification and Confounding Bias

Precision medicine aims to optimize medical care by tailoring treatment decisions to the unique characteristics of each patient. This objective naturally falls in the intersection between predictive analytics and causal inference; the former aims at predicting the outcome of interest given salient patient characteristics, and the latter seeks to answer counterfactual “what if” questions about the outcome. Most progress in clinical prediction models has centered around estimating the risk of adverse outcomes, often guiding treatment decisions by prioritizing high-risk individuals for intervention [1, 2]. However, the success of such guidance relies on an implicit assumption: that those classified as high risk are also those who will benefit the most from the treatment. A more relevant approach for decision-making is to directly predict treatment benefits based on patient characteristics. In this context, treatment benefit refers to the conditional risk reduction for treated individuals given their characteristics, compared to the conditional risk they would face under identical conditions without the treatment. This comparison, which is hypothetical, requires causal reasoning to predict [3]. Such a prediction is often termed “causal prediction” or “counterfactual prediction” [4]. In causal inference, the treatment benefit we are predicting is the conditional average treatment effect (CATE).

We are interested in a function that maps a patient’s characteristics to their putative treatment benefit, which we call a treatment benefit predictor (TBP). Suppose we have a pre-specified TBP developed from a randomized controlled trial (RCT) in one setting, where it demonstrates strong predictive performance. Before being adopted into patient care for the target population with a potentially distinct setting, this TBP needs to be evaluated (validated) [5, 6]. Transportability and generalizability have been studied across various epidemiology-adjacent disciplines, focusing on concepts, assessment methods, and correction techniques, most of which require access to both original and target population data. This includes applications in estimating causal estimands (e.g. average treatment effect (ATE)) [7], identifying causal relations [8], and predicting outcome risks [9] in the target population. However, we do not aim to assess or correct the transportability and generalizability of the pre-specified TBP, as we only have access to data from the target population. Instead, we focus on the predictive performance metrics, though transportability and generalizability may naturally be implied by these results, e.g. if we find the TBP performs well.

Traditionally, performance metrics for risk prediction are categorized into measures of overall fit (e.g. Brier score), discrimination (e.g. c-statistic), and calibration (e.g. calibration plot) [10, 11]. Various performance measures for TBPs have been formulated by extending the concepts from risk prediction to the treatment-benefit paradigm. The overall fit of a TBP is typically evaluated by the discrepancy between predicted treatment benefits and estimated CATE, using metrics like the precision in estimating heterogeneous effects [12, 13]. The discriminative ability reflects the degree to which TBP ranks patients in the correct order of treatment benefit. For instance, the c-for-benefit proposed by van Klaveren et al. is an extension of the c-statistic [14]. However, some possible limitations of the c-for-benefit have been discussed [9, 15, 16]. In particular, it has been shown to lack the properties of a proper scoring rule [17]. Other metrics of discriminatory performance are the concentration of the benefit index ($C_{b}$) [18] and the rank-weighted average treatment effect metrics [19], which are extended from the relative concentration curve (RCC) [20] and Qini curve [21], respectively. Both metrics compare how individuals prioritized by TBP predictions benefit more from treatment than those selected randomly. Note that RCCs and Qini curves are conceptually related to the Lorenz curve and the Gini index [22]. Calibration evaluates the agreement between the predicted and the actual quantities. Van Calster et al. proposed a hierarchical definition in the context of risk prediction [23]. Researchers have extended the ‘weak calibration’ to assess the performance of TBPs via the calibration intercept and slope [9, 16]. They have also extended the ‘moderate calibration’ to examine whether the conditional expected treatment benefit given TBP prediction, equals the TBP prediction [24].

None of these extensions are straightforward, given the unavailability of the counterfactual outcome. Another challenge arises when evaluating TBP performance using observational studies, as these introduce potential confounding bias. This is especially relevant when RCTs are unavailable, underpowered, or not representative of the target population. Additionally, for some interventions where equipoise is not established, conducting a RCT might be unethical.

In this study, we show conceptually how to approach evaluating a TBP from observational data, focusing on a single time-point binary treatment decision. Specifically, for several measures of TBP performance, we re-express the metrics, initially defined in terms of the latent treatment benefit variable, in terms of the distribution of observable quantities. To illustrate, we focus on two existing performance metrics for assessing pre-specified TBPs: $C_{b}$ and the moderate calibration curve. Additionally, we pay close attention to the bias incurred if we are not fully controlling for confounding. In doing so, we extend the understanding of identification and confounding bias from the context of well-understood estimands such as ATE to estimands which describe the predictive performance of a pre-specified TBP [25, 26]. Some applied problems have more complex settings than those of a single time-point binary treatment choice. For instance, a recent study extended the definitions and estimators of calibration and discrimination metrics to time-to-event outcomes and time-varying treatment decisions[27]. However, their focus was not on evaluating predictive performance in scenarios where confounding is not fully controlled or on exploring confounding bias in this context.

Each individual in the target population is described by $(Y^{(0)},Y^{(1)},A,X,Z)$ with joint distribution $\mathbb{P}$. Here $A$ denotes a binary treatment, and $Y^{(a)}$ is the counterfactual outcome that would be observed under treatment $A=a$, where $a\in\{0,1\}$. A larger $Y^{(a)}$ is presumed to be the preferred outcome. The set $X$ consists of $d$ pre-treatment covariates that are observable in routine clinical practice and in observational studies, and are used to predict treatment benefit. In contrast, $Z$ is a distinct set of $p$ covariates that may only be available in the observational study and could be necessary for controlling confounding. Both $X$ and $Z$ may include a mix of continuous and discrete variables. For instance, $X$ could be blood pressure and obesity, available at the point of care, which are used to predict the benefit of statin therapy for cardiovascular diseases. Meanwhile, $Z$ could be socioeconomic status, which is a confounding variable but not often used for predicting benefit from statins.

Individual treatment benefit is quantified as $B=Y^{(1)}-Y^{(0)}$, which is unobservable. For instance, when an individual has received $A=0$, the corresponding outcome $Y^{(1)}$ remains unobserved (and therefore counterfactual). We denote the mean treatment benefit for groups partitioned by both $X$ and $Z$ as $\tau(x,z)=\operatorname{E}[B\mid X=x,Z=z]$. However, in clinical practice, the ideal quantity to inform treatment decisions for an individual with $X=x$ is $\tau_{s}(x)=\operatorname{E}[B\mid X=x]$, which conditions only on $X$. Denote the common causal estimand ATE as $\operatorname{E}[B]$. A TBP predicting the mean treatment benefit based on patient characteristics $X$ is denoted as $h(x)$. This can guide treatment decision-making, for example the care provider offering treatment only to those with $h(x)>0$. We denote the predicted treatment benefit from $h(x)$ as $H:=h(X)$. The best possible TBP is $\tau_{s}(x)$ itself, and the corresponding prediction is $\tau_{s}(X)$.

Observed data from the observational study are realizations of a set of random variables drawn from an underlying probability distribution $\mathbb{P}_{obs}$, denoted as $(Y,A,X,Z)\sim\mathbb{P}_{obs}$, which is a consequence of $\mathbb{P}$. We define $\mu_{a}(x,z)=\operatorname{E}[Y\mid A=a,X=x,Z=z]$ and $\mu_{a}(x)=\operatorname{E}[Y\mid A=a,X=x]$. We denote the propensity score as $e(x,z)=\operatorname{P}(A=1\mid X=x,Z=z)$.

To ensure that $\tau_{s}(X)$ can be identified from the observed data, $\mu_{a}(X,Z)$ must have a causal interpretation. Therefore, the following assumptions are universally required: (1) no interference: between any two individuals, the treatment taken by one does not affect the counterfactual outcomes of the other; (2) consistency: the counterfactual outcome under the observed treatment assignment equals the observed outcome $Y$; and (3) conditional exchangeability: the treatment assignment is independent of the counterfactual outcomes, given the set of variables $X\cup Z$; (4) overlap (positivity): the conditional probability of receiving the active treatment or not is bounded away from $0$ and $1$, i.e., $0<e(x,z)<1$, for all possible $x$ and $z$. The first two assumptions are known as the stable-unit-treatment-value assumption (SUTVA) [28]. The third one assumes no unmeasured confounders given $X$ and $Z$, as illustrated in Figure S1.

Given a pre-specified TBP $h(x)$, we intend to evaluate its predictive performance on observational data from the target population. Among all possible predictive performance metrics for evaluating TBP, we exemplify with two metrics by specifying the estimand of each, which helps conceptually in understanding how to measure the predictive performance of $h(\cdot)$ using observational data.

If we treat the observational data as if they arose from a RCT, or if we do not sufficiently control for confounding, confounding bias may emerge. It is essential to investigate this bias and grasp how the lack of full control for confounding might affect the accuracy of our TBP performance evaluations. We focus on the confounding bias that occurs when $X$ alone is insufficient to control for confounding and denote the bias as a function of $X$, which is

where $\tau_{s}(X)$ is obtained by taking the expectation of $\tau(X,Z)$ over $Z$. To illustrate the propagation of $\text{bias}(X)$ to performance measures, we denote the inaccurate $C_{b}$ and $\operatorname{E}[B\mid H]$, calculated without controlling for $Z$, as $\tilde{C}_{b}$ and $\tilde{E}[B\mid H]$, respectively.

The bias of $C_{b}$ and the bias function of the moderate calibration curve, $\operatorname{E}[B\mid H=h]$, will be influenced by $\text{bias}(X)$. For $C_{b}$, this bias affects the identification of both $\operatorname{E}[B]$ and $\operatorname{E}[B\eta(X)]$. The discrepancy from actual $\operatorname{E}[B]$ is $\operatorname{E}[\text{bias}(X)]$, while the deviation from $\operatorname{E}[B\eta(H)]$ is $\operatorname{E}[\text{bias}(X)\eta(H)]$. However, expressing the bias of $C_{b}$ is complex as it involves the difference between two ratios, which is in the form of:

This bias value depends on $\text{bias}(X)$, $B$ and $\eta(H)$.

For $\operatorname{E}[B\mid H=h]$, the deviation from the accurate assessment can be expressed as

which is a function of $h$. It depends on $\text{bias}(X)$ and the association between $H$ and $X$.

According to Equation (2) and Equation (3), the deviations in $C_{b}$ and moderate calibration may yield zero value(s) even with non-zero $\text{bias}(X)$. In particular, zero deviation in $C_{b}$ would occur when the numerator in Equation (2) equals zero, and zero deviations in moderate calibration curve would occur when $\operatorname{E}[\text{bias}(X)\mid H=h]=0$ for some $h$. In the next section, we further investigate these biases in several illustrative examples to demonstrate how $\text{bias}(X)$ influences and propagates to the evaluation measures.

In clinical settings, TBPs may offer valuable guidance for physicians and patients in making informed treatment decisions. Before being transported to a new population, these TBPs need to be evaluated in that population. Observational data, where treatment assignment is not random, might be the only opportunity for such evaluation. Consequently, addressing confounding bias is crucial when assessing TBPs based on observational data from the target population.

This study demonstrated conceptually how to evaluate pre-specified TBPs using observational data. We delved into two specific measures, $C_{b}$ and moderate calibration curve, which are estimands defined in terms of the latent benefit variable $B$. We showed the identification of $C_{b}$ and $\operatorname{E}[B\mid H]$ by re-expressing the respective estimands in terms of the distribution of $(Y,A,X,Z)$. The absence of full confounding control leads to bias in identifying $C_{b}$ and $\operatorname{E}[B\mid H]$. The behavior of these biases are more complex than when targeting commonly encountered estimands such as ATE. We showed how the bias associated with $\tau_{s}(X)$, denoted as $\text{bias}(X)$, propagates to that associated with $C_{b}$ and $\tau_{s}(X)$. We demonstrated the behaviors of these biases in two synthetic populations with $\text{bias}(X)>0$. Positive $\text{bias}(X)$ function will always result in $\tilde{\operatorname{E}}[B\mid H]>\operatorname{E}[B\mid H]$; nevertheless, it might result $\tilde{C}_{b}<C_{b}$ for some choices of $h(x)$. In the two synthetic populations, we showed that positive $\text{bias}(X)$ functions lead to overestimation of $\operatorname{E}[B]$ and $\operatorname{E}[B\eta(H)]$ but underestimation of $C_{b}$ for our choices of $h(x)$. However, it is also possible that the overestimation of $\operatorname{E}[B\eta(H)]$ is greater than that of $\operatorname{E}[B]$, leading to an overestimation of $C_{b}$. When $\text{bias}(X)$ function takes both positive and negative values, behaviors of bias in identifying $C_{b}$ and $\operatorname{E}[B\mid H]$ become more unpredictable, making it difficult to intuit the specific direction of bias.

While we focused on population-level performance of TBPs and used synthetic populations, several findings from this study are relevant for practice. The primary objective of evaluating a pre-determined TBP is to estimate $\tau(x,z)$, for all possible $x$ and $z$, from sample data. Estimating $\tau(x,z)$ is a classical topic in the casual inference literature. For instance, it can be estimated through outcome regression, inverse probability weighting [31], or a combination of both, such as with the doubly robust method [32]. Additionally, Bayesian additive regression trees (BART) have been shown to offer several advantages in estimating CATE [12]. Building on the work presented here, we have initiated a new project that explores BART in TBP assessment using a real-world observational study.

Even with a scheme to estimate $\tau(\cdot)$, further effect will be needed to estimate specific performance metrics, such as $C_{b}$ and the calibration curve. For instance, one possible estimator of $C_{b}$ uses sample averages to estimate expectations in the estimand and estimate the CDF of $H$, and possibly its $f_{H}(h)$, through the empirical distribution of $H$. For the calibration curve, the sample average within groups sharing the same $H$ can be used to estimate the conditional expectation of the estimated $\tau(X,Z)$ given $H=h$ when $H$ is discrete. For continuous $H$, one approach is to discretize $H$ into equally sized bins and use the sample average within each bin to estimate the conditional expectation[24]. Moderate calibration can also be assessed non-parametrically, without the need for grouping or smoothing factors [33]. To our knowledge, this approach, based on the cumulative sum of prediction errors has yet to be applied to TBP assessment. Overall, exploring and tackling the challenges associated with estimating different performance measures is worthwhile.

In summary, we presented a conceptual framework for evaluating the predictive performance of TBPs using observational data. We demonstrated that in the absence of complete confounding control, the unpredictable nature of bias can significantly impact the reliability of TBP assessments. Moreover, this framework can be readily applied to practical settings, enabling the estimation of performance metrics such as $C_{b}$ and calibration curves.

For variables $B$ and $H$, the RCC is $R(p)=\frac{\operatorname{E}[BI(H\leq h)]}{\operatorname{E}[B]}$, where $p$ represents the $p$-th quantile concerning the value of $H$. We assume that $E[B]>0$.

In the first population example, we defined a distribution $\mathbb{P}$ and a linear function $\tau(X,Z)$ with a total of 18 parameters. To begin with, we design moderately and strongly calibrated TBPs. To design a moderately calibrated $h_{2}(x_{1},x_{2})$, we specific a set of coefficients $c=(b_{0},b_{1},b_{1},b_{2})$, where

To design a strongly calibrated $h_{3}(x_{1},x_{2})$, we set $c=(c_{0},c_{1},c_{2},c_{3})$, where

Finally, we look at the discrimination and calibration assessments of $h_{1}$, $h_{2}$, $h_{3}$ and $h_{4}$ when there is no confounding bias. When $\beta_{1}=0$, $\text{bias}(X)=0$, resulting in zero bias for identifying both the $C_{b}$ and $\operatorname{E}[B\mid H]$. Results are shown in Figure S5. Particularly, the red and blue calibration curves for $h_{2}(X_{1},X_{2})$, $h_{3}(X_{1},X_{2})$ and $h_{4}(X_{1},X_{2})$ lie on the 45-degree diagonal line. Note that all values in this section are rounded to three decimal places.

We first compute $\operatorname{E}[B\mid H]$. Given $\operatorname{E}[B\mid H]=\operatorname{E}[\tau_{s}(X)\mid H]$, we need to find the joint PDF of $(\tau_{s}(X),H)$, where $\tau_{s}(X)=\max(X_{1},X_{2})$ and $H=X_{1}+X_{2}$. Since $\tau_{s}$ is not a one-to-one function, we look at the pre-image of the point $(\tau_{s}(x),h)$. When $h\leq 2\tau_{s}(x)$ and $\tau_{s}(x)\leq h$, the pre-image of any point $(\tau_{s}(x),h)$ consists of two points, which are expressed as