# Risk score uncertainty quantification

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.)

• "Canonical" might have some (heuristic, informal) meaning in mathematics, but is not relevant in statistics, where the model must be chosen based on circumstances and assumptions. Some things can be said about your approach, though, of which one of the most salient is that an additive error model for a probability like $r_i$ is unlikely to be suitable in most applications. Rather than asking such a generic question, then, please consider describing your actual problem.
– whuber
May 24, 2022 at 16:13
• A prediction interval gives an interval within which an observation is expected to be contained. A confidence interval gives an interval in which the mean of a distribution is expected to be contained. It doesn't make sense for a (binary) prediction interval to be something like $(0.3, 0.7)$, while that is a totally reasonable (even if wide) confidence interval. // What's wrong with having a measure of model performance? If the model tends to do well, then the predicted risk is believable. If the model tends to do poorly, then there is more uncertainty in the risk estimate.
– Dave
May 24, 2022 at 16:15
• @whuber Thanks - you're right of course, the additive model did not make sense and is also not what I am actually assuming. I changed the formulation. However, I am really interested in general approaches, not in ways to do this for one specific model - hence my broad question. Nevertheless, I just asked a follow-up question on bootstrapping approaches specifically, and will probably draw out more sub-questions. I somehow expected the answer to this to be "Sure, you use the standard XYZ something interval, d'uh," hence my broad question. :-) May 28, 2022 at 18:22
• @Dave Well, what would you call an interval that contains the true risk with 95% probability then? I have been grappling with this terminology for a while now; right now it seems to me that it is really a "risk prediction interval" wrt the risk regression problem. This is different from a "prediction set" wrt a classifier built on a risk regression model , which would be either {0}, {0, 1} or {1}. But happy to learn if this is the wrong terminology? May 28, 2022 at 18:25
• @Dave Concerning your second question, I really need uncertainty quantification for individual risk predictions, not an overall model assessment. I edited the question to make this clearer. May 28, 2022 at 18:26

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$$.