I am working on various risk score estimation problems. I assume individual subjects are associated with a true risk $$r_i = f(x_i; \epsilon_i), \quad 0 \leq r_i \leq 1,$$ where $x_i$ is some available information about the subject and $\epsilon_i$ describes aleatoric uncertainty about the risk, which I assume to be heteroskedastic. (I.e., the distribution of $\epsilon_i$ varies as a function of $x_i$.) Furthermore, we are given observations $y_i$ of the (binary) outcome: $$y_i \sim \text{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)$ using a bootstrapping approach. This should include both aleatoric uncertainty (resulting from the action of $\epsilon_i$) and epistemic uncertainty (resulting from uncertainty about the learned model).
Things I have found so far but was unable to transfer to my setting:
- The classical theory for this question seems to fall under the umbrella term "interval estimation for a binomial proportion," see, e.g., Brown et al. (2001). This includes things like Wald CIs etc., and there are also bootstrapping versions of these, see Mantalos et al. (2008) and Wang and Hutson (2013). However, these all consider the setting where a single (population-level) proportion is to be estimated, not where the estimated proportion is a function of some covariates $x$. (I am also unsure whether they consider heteroscedastic aleatoric uncertainty, see my comment below on confidence intervals.)
- There is a branch of research that estimates prediction intervals using bootstrapping methods. This seems exactly like what I want to do (recall that risk estimation is a regression problem and notice that I am not interested in confidence intervals, because these do not cover aleatoric uncertainty). However, the methods I could find so far are 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 (since $r_i$ is unobservable).
What would be / is there a feasible bootstrapping approach to do what I want?
Note: this is a more specific follow-up question to my earlier, broader question.
