Confidence intervals for accuracy metrics: pros and cons of the bootstrap method compared to variance estimation

Question

I want to report the CI for metrics like F1 and AUC. I'm a bit confused on when it's better to bootstrap it, or when to use a formula.

For F1 there are several works that estimate the variance and derive a (symmetric) confidence interval, like: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8936911.

Similarly for AUC there is the Delong method, and others.

However for all these metrics, we could bootstrap them by sampling with replacement, computing the metric for every sample, and then taking the percentiles of the distribution, like the 2.5% 97.5% and percentiles.

Which of these approaches is expected to be better?

I guess an issue with the bootstrap method is that it's slow to compute it thousands of times. Is that the only reason there are works to create a closed formula for the CI ?

Answer 1

There's situations where you can bootstrap. If your sample size is large, then bootstrapping when possible is really convenient for a number of reasons:

1. If it works, then it works for almost any metric you can define, while frequentist analytical solutions tend to be derived for one metric at a time (bad luck if your metric of interest is not covered)
2. It's easy to implement (assuming a straightforward bootstrapping scheme is appropriate for the situation and the metric is easy to calculate for a sample).

Where bootstrapping has problems is

1. When sample size is small, then the discreteness of bootstrapping can make it inefficient.
2. When another simple solution is easily available/already implemented, why do it, when it often takes longer to run/is a little bit more complex to implement than running a well-validated package that someone else already wrote?
3. Sometimes you just cannot bootstrap in a valid way. E.g. multiple-rate-multiple-case (MRMC) studies tend to make it very hard to bootstrap (some solutions still have a bootstrap component, but it's nightmarishly complicated*), because everything is in a way correlated with everything (you have both the same raters looking at cases, as well as the same cases being looked).

As the other answer mentioned, one approach is a form of mixed effects model that describes the whole process. We took that approach for a MRMC study, where there were 3 possible diagnoses, recently. The Methods Section of the article outlines the approach:

We analyzed the primary outcome using a Bayesian model that jointly modeled each patient’s true disease status (ie, expert panel diagnosis) as a categorical random variable and the diagnoses given for each patient by each physician or algorithm (determined using multinomial logistic regression). The multinomial logistic regression model included a separate intercept term for each combination of disease (asthma, COPD, and ACO) and group (PCPs, pulmonologists, and the AC/DC tool), as well as a random case and random rater effect.

We took this approach, because we could not find any published methods that applied to the metrics our clinical colleagues wanted to look at (most MRMC literature is about AuROC and as mentioned, simple bootstrapping does not work). We ended up fitting this in a Bayesian way using Stan, because that deals with the uncertainties sensibly. It might be possible to fit such a model with something like PROC NLMIXED in SAS and to get confidence intervals that way (they have a decent automatic procedure for SEs there that uses the delta method) in a large study, but the Bayesian approach deals more naturally with (nearly) empty cells.

* Regarding "nightmarishly complicated", look at the complexity in these papers for getting AuROC estimates with CIs from MRMC studies:

• Brandon D. Gallas, Gene A. Pennello, and Kyle J. Myers. Multireader multicase variance analysis for binary data. J. Opt. Soc. Am. A, 24(12): B70–B80, Dec 2007. doi: 10.1364/JOSAA.24.000B70.
• Brandon D. Gallas, Andriy Bandos, Frank W. Samuelson, and Robert F. Wagner. A framework for random-effects roc analysis: Biases with the boot- strap and other variance estimators. Communications in Statistics - Theory and Methods, 38(15):2586–2603, 2009. doi: 10.1080/03610920802610084.

Answer 2

I tend to have data with inherent structure, such as multiple samples per patient, multiple measurements per sample and the like. In that situation, the statistically independent unit is typically the patient rather than the row of the data matrix. (In some cases, independence is even more complicated, with multiple top-level sources of variance.)

While I can easily bootstrap patients (or more general: set up a resampling procedure so that statistical independence is obeyed*) the analytical c.i. calculations I've met so far assume statistical independence at the row level. I.e., they do not accomodate more complex real-life data structure.

It's not that I would speak against using analytical formulations of variance (though I'd recommend to closely check the "small print" of the respective method) - I just don't meet nice data that would allow me to use them...

However, if your data meets the requirements of the analytical/parametric calculation, I'd expect it to have more power than bootstrapping.

What I do use is e.g. binomial c.i.s for proportions like sensitivity, specificity or the predictive values with patient number as sample size. Often, I can show that model instability is negligible compared to this, and I thus get a sensible approximation. (And the result is often anyways so wide that it does a good job serving as warning not to overinterprete the results...)

Sometimes I use this also as back-of-the-envelope calculation before experiments start: the c.i. won't get any better than that - if that is not satisfactory, more patients are needed right from the beginning.

---

Mixed models that describe the underlying structure for the prediction error may be a way here, but I haven't seen any papers applying this so far. It is, however, an approach that I'm looking into.

I have a sneaking suspicion that situations where the variance on the error estimate can sensibly be approximated as stemming from a single source are actually very rare (and my educated guess is that this is where patient-to-patient variance dominates everything and if that's the case with the final error estimate, it is rarely of practical use).
Of course, one may still sum the contributions of various sources of error into one total variance - but also there, analytical calculation would need to account for the data structure and the corresponding sources of variance.

---

*Since I typically also work with high-dimensional data, I anyways have to set up such splitting procedures to generate suitable test data to produce the generalization error estimates: the fanciest variance estimator is of no use if my validation procedure has a (large) optimistic bias. Also keep in mind that optimistically biased test results on classification tasks often lead to optimistically biased (too narrow) variance estimates as well.

---

Update: back of the envelope calculations

Since I don't use F1, I'll outline mostly what I do for sensitivity (or analogously specificity). Sensitivity is the (observed) fraction of correctly predicted tested cases of the positive class, and we can thus use binomial confidence intervals. The width of these depends on the denominator, i.e. the number of tested cases of the positive class, and on the true value of the sensitivity (with maximum width for 50 % sens and minimum width for 0 and 100 % sens).

For the minimum confidence interval width, we can use the rule of 3: assuming we's observe no errors, the 95 % c.i. for sensitivity would be roughy ranging from 1 to $1 - \frac{3}{n}$

However, we can also calculate minimum and maximum c.i. width depending on the number of tested positive cases before any experimental results are available:

... and compare that to the the requirements of the application. (One can also set up calculation to estimate reasonable sample sizes based on expected performance and application requirements. One scenario I've been using is specifying for an application situation "super nice", "OK" and "not really acceptable" performance. If for the discussed sample size "super nice" performance couldn't be reliably distinguished from "not really acceptable", sample size is insufficient)

To get guesstimates for other measures that are not simply an observed fraction of a pre-specified number of test cases (such as F1), you could either start with the variance of the binomially distributed simple fractions and do error propagation, or set up simulations (which would be quite similar to bootstrapping).