Risk score uncertainty quantification

Question

I am working on various risk score estimation problems. I assume individual subjects are associated with a true risk $$ r_i = f(x_i, \varepsilon)$$ where $x_i$ is some available information about the subject and $\varepsilon_i$ describes aleatoric uncertainty about the risk, which I assume to be heteroskedastic. (I.e., the distribution of $\varepsilon_i$ varies as a function of $x_i$.) Furthermore, we are given observations $y_i$ of the (binary) outcome: $$ y_i \sim Ber(r_i).$$ The canonical way of estimating such risk scores is, of course, to use a probability model (say, a logistic regression or xgboost model), perform (log loss) regression, and interpret the model output (which lies in $[0, 1]$) as the risk estimate $\hat{r}_i(x_i)$ of that subject (possibly applying calibration to ensure that the risk scores coincide with the observed incidence, i.e., $\hat{r}_i(x_i) \approx E[y_i \mid X=x_i]$).

What I am now interested in is quantifying the uncertainty of individual risk estimates $\hat{r}_i(x_i)$. This should account for both aleatoric uncertainty (resulting from the action of $\epsilon_i$) and epistemic uncertainty (resulting from uncertainty about the learned model). However, I am having trouble finding canonical (ideally model-independent) methods to do so.

What I have learned so far:

• (Regression) prediction intervals seem closest to what I am interested in (since risk estimation is a regression problem). In my case, I would be interested in prediction intervals concerning the risk estimates $\hat{r}_i$, not the actual outcome $y_i$. Various comments and answers here on stats.SE suggest that this "does not make sense", but I do not see why that would be the case? Why exactly would it not make sense to apply any canonical PI method to the risk model and interpret the resulting interval as one within which a given individual's true risk is contained with likelihood X? (Notice that I am not interested in confidence intervals, because these only cover epistemic uncertainty.) I tried adapting standard bootstrapping-based approaches such as the one described here, but the problem is that that approach is based on resampling the residuals - and that does not seem feasible in my setting, because the "risk residuals" $r_i - \hat{r}_i$ are unobservable.
• There is a branch of work on quantile regression that can also be used to estimate prediction intervals. Again, I am having trouble adapting this line of work to the classification / binary outcome setting.
• The model I wrote down above differs from the standard logistic regression model due to the addition of the aleatoric risk uncertainty $\varepsilon_i$. My model has two independent sources of aleatoric uncertainty: 1) uncertainty in the risk score, and 2) uncertainty in the outcome $y_i$. Are there canonical methods that treat this setting? Or does my framing not make sense for any reason? For instance, I could imagine writing down a generative probabilistic model in the above form and then using MCMC or variational inference. This is something I am currently trying to implement, but I would prefer a model-agnostic approach.

To summarize: what would be a canonical (model-independent) method for quantifying (heteroscedastic) aleatoric and epistemic uncertainty of binary risk scores? (Bonus points given for ways to actually do this in python.)

Answer 1

It seems that what I am asking for is fundamentally impossible: in particular, it is impossible to quantify the aleatoric uncertainty associated with risk scores estimated from binary observations.

The problem can be seen quite simply. Consider two populations, A and B. In population A, every subject has an inherent risk of $r=0.5$ of a binary outcome. In population B, half of the subjects have a risk of $r=0.1$ and the other half has a risk of $r=0.9$. In both cases, the distribution of the binary outcome $y$ is given by $y \sim \mathrm{Ber}(p=0.5)$ – in other words, the two groups (with different levels of aleatoric uncertainty) are completely indistinguishable based on the distribution of the binary outcome. More generally, $$a \cdot \mathrm{Ber}(p_a) + b \cdot \mathrm{Ber}(p_b) = \mathrm{Ber}(a \cdot p_a + b \cdot p_b), \quad a+b=1.$$ Thus, it is impossible to quantify the level of aleatoric uncertainty given just the binary observations. (Without prior information about the underlying distribution of risk scores, that is.)

Notice that (as opposed to what I initially thought and wrote in my question), this is crucially different from the usual regression setting. Consider the following (regression) variation of the model I described in my question: $$ \mu_i = f(x_i, \varepsilon_i) \sim \mathcal{N}(\theta(x_i), \sigma(x_i)^2), \quad \quad y_i = \mathcal{N}(\mu_i, s^2).$$ In this case, the distribution of $y_i$ actually differs between groups with different levels of aleatoric noise (i.e., different values of $\sigma(x_i)^2$): we have $$y_i \sim \mathcal{N}(\theta(x_i), s^2+\sigma(x_i)^2),$$ and thus the aleatoric noise distribution is, at least to some degree, observable from the observations $y_i$.

There are a number of recent theoretical results in the literature pointing in similar directions, e.g., that distribution-free confidence intervals for risk predictions in binary classification are necessarily "wide," i.e, cover the whole interval [0, 1] in most cases (Gupta et al. 2022, Barber 2020).