Bayesian ROC Analysis: Binormal, Bibeta, and Nonparametric (Bayesian Bootstrap)
2025-10-03
Introduction
This document brings together three Bayesian approaches to ROC analysis with full posterior uncertainty:
- Parametric binormal model (scores for negatives and positives are Normal with different mean and variance). The ROC and AUC are computed inside the JAGS model using the Normal CDF.
- Parametric bibeta model (scores lie in [0,1] and follow Beta distributions by class). The ROC and AUC are computed in R from posterior draws.
- Nonparametric Bayesian bootstrap (no distributional assumptions). We put Dirichlet(1,…,1) priors on empirical weights and propagate uncertainty to ROC and AUC.
Convention: higher scores imply the positive class. If your direction is reversed, negate the scores or swap definitions accordingly.
Data and Notation
Let \(y_0 = (y_{0j})_{j=1}^{N_0}\) be scores for the negative class and \(y_1 = (y_{1i})_{i=1}^{N_1}\) for the positive class. For threshold \(t\):
- \(\mathrm{TPR}(t) = \Pr(y_1 > t)\)
- \(\mathrm{FPR}(t) = \Pr(y_0 > t)\)
The ROC is the set of pairs \((\mathrm{FPR}(t), \mathrm{TPR}(t))\) over all \(t\). The AUC equals \(\Pr(y_1 > y_0) + \tfrac12 \Pr(y_1 = y_0)\) (the Mann–Whitney probability with half credit for ties).
1) Binormal model in JAGS (ROC & AUC inside the model)
Assume
- \(y_0 \sim N(\mu_0, \sigma_0^2)\)
- \(y_1 \sim N(\mu_1, \sigma_1^2)\)
The binormal ROC has a closed-form AUC under this model.
JAGS model file
Save as bayes_roc_binormal.jags.
binormal_model <- 'model {
# Likelihood
for (i in 1:N0) { y0[i] ~ dnorm(mu0, tau0) }
for (j in 1:N1) { y1[j] ~ dnorm(mu1, tau1) }
# Priors (weakly-informative)
mu0 ~ dnorm(0, 1.0E-6)
mu1 ~ dnorm(0, 1.0E-6)
sigma0 ~ dunif(0, 1000)
sigma1 ~ dunif(0, 1000)
tau0 <- pow(sigma0, -2)
tau1 <- pow(sigma1, -2)
# ROC on a supplied grid of thresholds t[1:K]
for (k in 1:K) {
tpr[k] <- 1 - phi((t[k] - mu1) * sqrt(tau1))
fpr[k] <- 1 - phi((t[k] - mu0) * sqrt(tau0))
}
# Binormal AUC (closed form)
auc <- phi( (mu1 - mu0) / sqrt( (1 / tau0) + (1 / tau1) ) )
# OPTIONAL: Youden\'s J argmax (no which_max in JAGS)
for (k in 1:K) { youden[k] <- tpr[k] - fpr[k] }
best[1] <- youden[1]
idx[1] <- 1
for (k in 2:K) {
better[k] <- step(youden[k] - best[k-1])
best[k] <- better[k]*youden[k] + (1-better[k])*best[k-1]
idx[k] <- better[k]*k + (1-better[k])*idx[k-1]
}
k_star <- idx[K]
t_star <- t[k_star]
}'
writeLines(binormal_model, "bayes_roc_binormal.jags")
cat("Wrote bayes_roc_binormal.jags\n")
#> Wrote bayes_roc_binormal.jagsR driver using jagsUI
library(jagsUI)
set.seed(123)
# Simulate binormal data (edit as needed)
N0 <- 150; N1 <- 120
y0 <- rnorm(N0, mean = 0, sd = 1.0)
y1 <- rnorm(N1, mean = 1.0, sd = 1.2)
# Threshold grid covering the data
all_y <- c(y0, y1)
K <- 101
tgrid <- seq(min(all_y) - 0.5*sd(all_y), max(all_y) + 0.5*sd(all_y), length.out = K)
# JAGS data & settings
jags_data <- list(N0=length(y0), N1=length(y1), y0=y0, y1=y1, K=K, t=tgrid)
params <- c("mu0","mu1","sigma0","sigma1","auc","tpr","fpr","t_star","k_star")
inits <- function() list(
mu0 = mean(y0), mu1 = mean(y1),
sigma0 = sd(y0) + 0.1, sigma1 = sd(y1) + 0.1
)
fit_binormal <- jags(
data = jags_data, inits = inits,
parameters.to.save = params,
model.file = "bayes_roc_binormal.jags",
n.chains = 3, n.adapt = 1000, n.iter = 6000,
n.burnin = 2000, n.thin = 2, parallel = TRUE
)
#>
#> Processing function input.......
#>
#> Done.
#>
#> Beginning parallel processing using 3 cores. Console output will be suppressed.
#>
#> Parallel processing completed.
#>
#> Calculating statistics.......
#>
#> Done.
print(fit_binormal, digits = 3)
#> JAGS output for model 'bayes_roc_binormal.jags', generated by jagsUI.
#> Estimates based on 3 chains of 6000 iterations,
#> adaptation = 1000 iterations (sufficient),
#> burn-in = 2000 iterations and thin rate = 2,
#> yielding 6000 total samples from the joint posterior.
#> MCMC ran in parallel for 0.007 minutes at time 2025-10-03 10:50:41.842751.
#>
#> mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
#> mu0 -0.025 0.078 -0.179 -0.025 0.127 TRUE 0.618 1.000 6000
#> mu1 1.039 0.105 0.836 1.038 1.245 FALSE 1.000 1.001 4033
#> sigma0 0.959 0.055 0.857 0.956 1.073 FALSE 1.000 1.001 6000
#> sigma1 1.147 0.075 1.011 1.143 1.304 FALSE 1.000 1.001 2654
#> auc 0.761 0.029 0.703 0.762 0.815 FALSE 1.000 1.001 2429
#> tpr[1] 1.000 0.000 0.999 1.000 1.000 FALSE 1.000 1.001 2338
#> tpr[2] 0.999 0.000 0.998 1.000 1.000 FALSE 1.000 1.001 2301
#> tpr[3] 0.999 0.001 0.998 1.000 1.000 FALSE 1.000 1.001 2267
#> tpr[4] 0.999 0.001 0.997 0.999 1.000 FALSE 1.000 1.001 2236
#> tpr[5] 0.999 0.001 0.997 0.999 1.000 FALSE 1.000 1.001 2207
#> tpr[6] 0.999 0.001 0.996 0.999 1.000 FALSE 1.000 1.001 2180
#> tpr[7] 0.998 0.001 0.995 0.999 1.000 FALSE 1.000 1.001 2156
#> tpr[8] 0.998 0.001 0.994 0.998 1.000 FALSE 1.000 1.001 2133
#> tpr[9] 0.997 0.002 0.993 0.998 0.999 FALSE 1.000 1.001 2113
#> tpr[10] 0.997 0.002 0.992 0.997 0.999 FALSE 1.000 1.001 2094
#> tpr[11] 0.996 0.002 0.990 0.997 0.999 FALSE 1.000 1.001 2076
#> tpr[12] 0.995 0.003 0.988 0.996 0.999 FALSE 1.000 1.001 2060
#> tpr[13] 0.994 0.003 0.986 0.995 0.998 FALSE 1.000 1.001 2045
#> tpr[14] 0.993 0.004 0.984 0.994 0.998 FALSE 1.000 1.001 2031
#> tpr[15] 0.991 0.004 0.981 0.992 0.997 FALSE 1.000 1.001 2017
#> tpr[16] 0.989 0.005 0.978 0.990 0.996 FALSE 1.000 1.001 2005
#> tpr[17] 0.987 0.006 0.974 0.988 0.995 FALSE 1.000 1.001 1993
#> tpr[18] 0.985 0.007 0.970 0.986 0.994 FALSE 1.000 1.001 1982
#> tpr[19] 0.982 0.007 0.965 0.983 0.993 FALSE 1.000 1.001 1972
#> tpr[20] 0.978 0.008 0.959 0.980 0.991 FALSE 1.000 1.001 1962
#> tpr[21] 0.975 0.009 0.953 0.976 0.989 FALSE 1.000 1.001 1953
#> tpr[22] 0.970 0.010 0.946 0.971 0.986 FALSE 1.000 1.001 1944
#> tpr[23] 0.965 0.012 0.939 0.966 0.983 FALSE 1.000 1.001 1936
#> tpr[24] 0.959 0.013 0.930 0.960 0.980 FALSE 1.000 1.001 1929
#> tpr[25] 0.952 0.014 0.921 0.954 0.975 FALSE 1.000 1.001 1923
#> tpr[26] 0.945 0.015 0.911 0.946 0.970 FALSE 1.000 1.001 1917
#> tpr[27] 0.936 0.017 0.900 0.938 0.964 FALSE 1.000 1.001 1913
#> tpr[28] 0.927 0.018 0.887 0.928 0.958 FALSE 1.000 1.001 1910
#> tpr[29] 0.916 0.020 0.874 0.918 0.950 FALSE 1.000 1.001 1909
#> tpr[30] 0.904 0.021 0.860 0.906 0.941 FALSE 1.000 1.001 1910
#> tpr[31] 0.892 0.023 0.844 0.893 0.931 FALSE 1.000 1.001 1913
#> tpr[32] 0.878 0.024 0.828 0.879 0.920 FALSE 1.000 1.001 1919
#> tpr[33] 0.862 0.025 0.810 0.864 0.907 FALSE 1.000 1.001 1929
#> tpr[34] 0.846 0.027 0.791 0.847 0.894 FALSE 1.000 1.001 1943
#> tpr[35] 0.828 0.028 0.771 0.829 0.879 FALSE 1.000 1.001 1963
#> tpr[36] 0.809 0.029 0.749 0.810 0.863 FALSE 1.000 1.001 1988
#> tpr[37] 0.788 0.030 0.726 0.790 0.845 FALSE 1.000 1.001 2021
#> tpr[38] 0.767 0.031 0.703 0.768 0.825 FALSE 1.000 1.001 2062
#> tpr[39] 0.744 0.032 0.679 0.745 0.805 FALSE 1.000 1.001 2114
#> tpr[40] 0.720 0.033 0.653 0.721 0.783 FALSE 1.000 1.001 2178
#> tpr[41] 0.695 0.034 0.628 0.696 0.759 FALSE 1.000 1.001 2256
#> tpr[42] 0.669 0.035 0.601 0.670 0.734 FALSE 1.000 1.001 2353
#> tpr[43] 0.642 0.035 0.573 0.643 0.709 FALSE 1.000 1.001 2471
#> tpr[44] 0.615 0.035 0.544 0.615 0.682 FALSE 1.000 1.001 2614
#> tpr[45] 0.586 0.036 0.516 0.586 0.656 FALSE 1.000 1.001 2789
#> tpr[46] 0.558 0.036 0.487 0.558 0.628 FALSE 1.000 1.001 3002
#> tpr[47] 0.529 0.036 0.458 0.529 0.599 FALSE 1.000 1.001 3262
#> tpr[48] 0.500 0.036 0.429 0.500 0.571 FALSE 1.000 1.001 3580
#> tpr[49] 0.471 0.036 0.400 0.471 0.542 FALSE 1.000 1.001 3970
#> tpr[50] 0.442 0.036 0.372 0.442 0.513 FALSE 1.000 1.001 4450
#> tpr[51] 0.413 0.036 0.344 0.413 0.484 FALSE 1.000 1.001 5043
#> tpr[52] 0.385 0.036 0.317 0.385 0.457 FALSE 1.000 1.001 5780
#> tpr[53] 0.358 0.035 0.291 0.357 0.429 FALSE 1.000 1.001 6000
#> tpr[54] 0.331 0.035 0.265 0.330 0.401 FALSE 1.000 1.001 6000
#> tpr[55] 0.305 0.034 0.241 0.304 0.374 FALSE 1.000 1.000 6000
#> tpr[56] 0.280 0.033 0.218 0.279 0.348 FALSE 1.000 1.000 6000
#> tpr[57] 0.256 0.032 0.196 0.255 0.322 FALSE 1.000 1.000 6000
#> tpr[58] 0.233 0.031 0.175 0.232 0.298 FALSE 1.000 1.000 6000
#> tpr[59] 0.212 0.030 0.156 0.210 0.274 FALSE 1.000 1.000 6000
#> tpr[60] 0.191 0.029 0.138 0.190 0.252 FALSE 1.000 1.000 6000
#> tpr[61] 0.172 0.028 0.121 0.171 0.231 FALSE 1.000 1.000 6000
#> tpr[62] 0.154 0.027 0.106 0.153 0.211 FALSE 1.000 1.000 6000
#> tpr[63] 0.138 0.025 0.093 0.136 0.192 FALSE 1.000 1.000 6000
#> tpr[64] 0.122 0.024 0.080 0.121 0.173 FALSE 1.000 1.000 6000
#> tpr[65] 0.108 0.023 0.069 0.107 0.157 FALSE 1.000 1.000 6000
#> tpr[66] 0.096 0.021 0.059 0.094 0.141 FALSE 1.000 1.000 6000
#> tpr[67] 0.084 0.020 0.050 0.082 0.127 FALSE 1.000 1.000 6000
#> tpr[68] 0.073 0.018 0.042 0.072 0.113 FALSE 1.000 1.000 6000
#> tpr[69] 0.064 0.017 0.036 0.062 0.101 FALSE 1.000 1.000 6000
#> tpr[70] 0.055 0.016 0.030 0.054 0.090 FALSE 1.000 1.000 6000
#> tpr[71] 0.048 0.014 0.025 0.046 0.080 FALSE 1.000 1.000 6000
#> tpr[72] 0.041 0.013 0.021 0.040 0.070 FALSE 1.000 1.000 6000
#> tpr[73] 0.035 0.012 0.017 0.034 0.062 FALSE 1.000 1.000 6000
#> tpr[74] 0.030 0.011 0.014 0.029 0.054 FALSE 1.000 1.000 6000
#> tpr[75] 0.025 0.009 0.011 0.024 0.048 FALSE 1.000 1.000 6000
#> tpr[76] 0.022 0.008 0.009 0.020 0.041 FALSE 1.000 1.000 6000
#> tpr[77] 0.018 0.007 0.007 0.017 0.036 FALSE 1.000 1.000 6000
#> tpr[78] 0.015 0.007 0.006 0.014 0.031 FALSE 1.000 1.000 6000
#> tpr[79] 0.013 0.006 0.005 0.012 0.027 FALSE 1.000 1.000 6000
#> tpr[80] 0.011 0.005 0.004 0.010 0.023 FALSE 1.000 1.000 6000
#> tpr[81] 0.009 0.004 0.003 0.008 0.020 FALSE 1.000 1.000 6000
#> tpr[82] 0.007 0.004 0.002 0.006 0.017 FALSE 1.000 1.000 5716
#> tpr[83] 0.006 0.003 0.002 0.005 0.014 FALSE 1.000 1.000 5398
#> tpr[84] 0.005 0.003 0.001 0.004 0.012 FALSE 1.000 1.000 5120
#> tpr[85] 0.004 0.002 0.001 0.003 0.010 FALSE 1.000 1.000 4875
#> tpr[86] 0.003 0.002 0.001 0.003 0.009 FALSE 1.000 1.000 4660
#> tpr[87] 0.003 0.002 0.001 0.002 0.007 FALSE 1.000 1.001 4470
#> tpr[88] 0.002 0.001 0.000 0.002 0.006 FALSE 1.000 1.001 4302
#> tpr[89] 0.002 0.001 0.000 0.001 0.005 FALSE 1.000 1.001 4154
#> tpr[90] 0.001 0.001 0.000 0.001 0.004 FALSE 1.000 1.001 4024
#> tpr[91] 0.001 0.001 0.000 0.001 0.003 FALSE 1.000 1.001 3909
#> tpr[92] 0.001 0.001 0.000 0.001 0.003 FALSE 1.000 1.001 3808
#> tpr[93] 0.001 0.001 0.000 0.001 0.002 FALSE 1.000 1.001 3720
#> tpr[94] 0.001 0.000 0.000 0.000 0.002 FALSE 1.000 1.001 3644
#> tpr[95] 0.000 0.000 0.000 0.000 0.002 FALSE 1.000 1.001 3577
#> tpr[96] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.002 3520
#> tpr[97] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.002 3472
#> tpr[98] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.002 3431
#> tpr[99] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.002 3397
#> tpr[100] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.003 3370
#> tpr[101] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.003 3349
#> fpr[1] 0.998 0.001 0.996 0.999 1.000 FALSE 1.000 1.001 5630
#> fpr[2] 0.998 0.001 0.995 0.998 0.999 FALSE 1.000 1.001 6000
#> fpr[3] 0.997 0.002 0.993 0.998 0.999 FALSE 1.000 1.001 6000
#> fpr[4] 0.996 0.002 0.992 0.997 0.999 FALSE 1.000 1.001 6000
#> fpr[5] 0.995 0.002 0.990 0.996 0.999 FALSE 1.000 1.001 6000
#> fpr[6] 0.994 0.003 0.987 0.995 0.998 FALSE 1.000 1.001 6000
#> fpr[7] 0.993 0.003 0.984 0.993 0.997 FALSE 1.000 1.001 6000
#> fpr[8] 0.991 0.004 0.981 0.991 0.997 FALSE 1.000 1.001 6000
#> fpr[9] 0.988 0.005 0.977 0.989 0.995 FALSE 1.000 1.001 6000
#> fpr[10] 0.986 0.006 0.972 0.986 0.994 FALSE 1.000 1.000 6000
#> fpr[11] 0.982 0.006 0.967 0.983 0.992 FALSE 1.000 1.000 6000
#> fpr[12] 0.978 0.007 0.961 0.979 0.990 FALSE 1.000 1.000 6000
#> fpr[13] 0.973 0.009 0.953 0.974 0.987 FALSE 1.000 1.000 6000
#> fpr[14] 0.967 0.010 0.945 0.969 0.983 FALSE 1.000 1.000 6000
#> fpr[15] 0.961 0.011 0.936 0.962 0.979 FALSE 1.000 1.000 6000
#> fpr[16] 0.953 0.012 0.925 0.954 0.974 FALSE 1.000 1.000 6000
#> fpr[17] 0.944 0.014 0.913 0.945 0.967 FALSE 1.000 1.000 6000
#> fpr[18] 0.933 0.015 0.899 0.935 0.960 FALSE 1.000 1.000 6000
#> fpr[19] 0.921 0.017 0.884 0.923 0.951 FALSE 1.000 1.000 6000
#> fpr[20] 0.908 0.018 0.868 0.909 0.941 FALSE 1.000 1.000 6000
#> fpr[21] 0.893 0.020 0.850 0.894 0.928 FALSE 1.000 1.000 6000
#> fpr[22] 0.876 0.021 0.830 0.877 0.915 FALSE 1.000 1.000 6000
#> fpr[23] 0.857 0.023 0.809 0.859 0.899 FALSE 1.000 1.000 6000
#> fpr[24] 0.837 0.024 0.786 0.838 0.882 FALSE 1.000 1.000 6000
#> fpr[25] 0.814 0.026 0.762 0.816 0.862 FALSE 1.000 1.000 6000
#> fpr[26] 0.790 0.027 0.736 0.792 0.841 FALSE 1.000 1.000 6000
#> fpr[27] 0.764 0.028 0.708 0.766 0.817 FALSE 1.000 1.000 6000
#> fpr[28] 0.737 0.029 0.679 0.738 0.792 FALSE 1.000 1.000 6000
#> fpr[29] 0.708 0.030 0.648 0.709 0.765 FALSE 1.000 1.000 6000
#> fpr[30] 0.677 0.031 0.616 0.678 0.736 FALSE 1.000 1.000 6000
#> fpr[31] 0.645 0.031 0.583 0.646 0.705 FALSE 1.000 1.000 6000
#> fpr[32] 0.612 0.032 0.549 0.613 0.673 FALSE 1.000 1.000 6000
#> fpr[33] 0.578 0.032 0.515 0.579 0.639 FALSE 1.000 1.000 6000
#> fpr[34] 0.544 0.032 0.480 0.544 0.605 FALSE 1.000 1.000 6000
#> fpr[35] 0.509 0.032 0.446 0.509 0.572 FALSE 1.000 1.000 6000
#> fpr[36] 0.475 0.032 0.411 0.475 0.537 FALSE 1.000 1.000 6000
#> fpr[37] 0.440 0.032 0.377 0.440 0.503 FALSE 1.000 1.000 6000
#> fpr[38] 0.406 0.032 0.343 0.405 0.469 FALSE 1.000 1.000 6000
#> fpr[39] 0.372 0.031 0.311 0.372 0.434 FALSE 1.000 1.000 6000
#> fpr[40] 0.340 0.031 0.280 0.339 0.401 FALSE 1.000 1.000 6000
#> fpr[41] 0.309 0.030 0.251 0.308 0.368 FALSE 1.000 1.000 6000
#> fpr[42] 0.279 0.029 0.223 0.278 0.337 FALSE 1.000 1.000 6000
#> fpr[43] 0.250 0.028 0.197 0.249 0.307 FALSE 1.000 1.000 6000
#> fpr[44] 0.223 0.027 0.172 0.222 0.279 FALSE 1.000 1.000 6000
#> fpr[45] 0.198 0.026 0.149 0.197 0.251 FALSE 1.000 1.000 6000
#> fpr[46] 0.175 0.025 0.129 0.174 0.226 FALSE 1.000 1.000 6000
#> fpr[47] 0.154 0.024 0.110 0.152 0.202 FALSE 1.000 1.000 6000
#> fpr[48] 0.134 0.022 0.094 0.133 0.180 FALSE 1.000 1.000 6000
#> fpr[49] 0.116 0.021 0.079 0.115 0.159 FALSE 1.000 1.000 6000
#> fpr[50] 0.100 0.019 0.066 0.099 0.140 FALSE 1.000 1.000 6000
#> fpr[51] 0.086 0.018 0.055 0.084 0.123 FALSE 1.000 1.000 6000
#> fpr[52] 0.073 0.016 0.045 0.072 0.107 FALSE 1.000 1.000 6000
#> fpr[53] 0.062 0.015 0.037 0.060 0.093 FALSE 1.000 1.000 6000
#> fpr[54] 0.052 0.013 0.030 0.051 0.080 FALSE 1.000 1.000 6000
#> fpr[55] 0.043 0.012 0.024 0.042 0.069 FALSE 1.000 1.000 6000
#> fpr[56] 0.036 0.010 0.019 0.035 0.059 FALSE 1.000 1.000 6000
#> fpr[57] 0.030 0.009 0.015 0.029 0.050 FALSE 1.000 1.000 6000
#> fpr[58] 0.024 0.008 0.012 0.023 0.042 FALSE 1.000 1.000 6000
#> fpr[59] 0.020 0.007 0.009 0.019 0.036 FALSE 1.000 1.000 6000
#> fpr[60] 0.016 0.006 0.007 0.015 0.030 FALSE 1.000 1.001 6000
#> fpr[61] 0.013 0.005 0.005 0.012 0.025 FALSE 1.000 1.001 6000
#> fpr[62] 0.010 0.004 0.004 0.010 0.021 FALSE 1.000 1.001 6000
#> fpr[63] 0.008 0.004 0.003 0.008 0.017 FALSE 1.000 1.001 6000
#> fpr[64] 0.007 0.003 0.002 0.006 0.014 FALSE 1.000 1.001 5689
#> fpr[65] 0.005 0.003 0.002 0.005 0.011 FALSE 1.000 1.001 5363
#> fpr[66] 0.004 0.002 0.001 0.004 0.009 FALSE 1.000 1.001 5055
#> fpr[67] 0.003 0.002 0.001 0.003 0.007 FALSE 1.000 1.002 4763
#> fpr[68] 0.002 0.001 0.001 0.002 0.006 FALSE 1.000 1.002 4488
#> fpr[69] 0.002 0.001 0.000 0.002 0.005 FALSE 1.000 1.002 4229
#> fpr[70] 0.001 0.001 0.000 0.001 0.004 FALSE 1.000 1.002 3987
#> fpr[71] 0.001 0.001 0.000 0.001 0.003 FALSE 1.000 1.003 3760
#> fpr[72] 0.001 0.001 0.000 0.001 0.002 FALSE 1.000 1.003 3547
#> fpr[73] 0.001 0.000 0.000 0.000 0.002 FALSE 1.000 1.004 3349
#> fpr[74] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.004 3165
#> fpr[75] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.005 2994
#> fpr[76] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.005 2836
#> fpr[77] 0.000 0.000 0.000 0.000 0.001 FALSE 1.000 1.006 2689
#> fpr[78] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.007 2553
#> fpr[79] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.008 1
#> fpr[80] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.009 1
#> fpr[81] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.010 1
#> fpr[82] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.011 1
#> fpr[83] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.013 1
#> fpr[84] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.014 1
#> fpr[85] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.016 1
#> fpr[86] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.018 1
#> fpr[87] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.021 1
#> fpr[88] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.024 1
#> fpr[89] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.027 1
#> fpr[90] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.030 1
#> fpr[91] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.034 1
#> fpr[92] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.038 1
#> fpr[93] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.043 1
#> fpr[94] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.048 1
#> fpr[95] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.054 1
#> fpr[96] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.060 1
#> fpr[97] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.066 1
#> fpr[98] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.073 1
#> fpr[99] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.081 1
#> fpr[100] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.088 1
#> fpr[101] 0.000 0.000 0.000 0.000 0.000 FALSE 1.000 1.096 1
#> t_star 0.638 0.093 0.454 0.621 0.788 FALSE 1.000 1.000 6000
#> k_star 43.204 1.110 41.000 43.000 45.000 FALSE 1.000 1.000 6000
#> deviance 783.318 2.793 779.827 782.686 790.199 FALSE 1.000 1.003 1083
#>
#> Successful convergence based on Rhat values (all < 1.1).
#> Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> For each parameter, n.eff is a crude measure of effective sample size.
#>
#> overlap0 checks if 0 falls in the parameter's 95% credible interval.
#> f is the proportion of the posterior with the same sign as the mean;
#> i.e., our confidence that the parameter is positive or negative.
#>
#> DIC info: (pD = var(deviance)/2)
#> pD = 3.9 and DIC = 787.211
#> DIC is an estimate of expected predictive error (lower is better).
# Summaries for ROC
TPR_sims <- fit_binormal$sims.list$tpr
FPR_sims <- fit_binormal$sims.list$fpr
AUC_sims <- fit_binormal$sims.list$auc
# Direction check: if AUC<0.5, reflect ROC and AUC
if (mean(AUC_sims) < 0.5) {
message("Detected reversed score orientation; reflecting ROC and AUC.")
TPR_old <- TPR_sims; FPR_old <- FPR_sims
TPR_sims <- 1 - FPR_old
FPR_sims <- 1 - TPR_old
AUC_sims <- 1 - AUC_sims
}
TPR_mean <- apply(TPR_sims, 2, mean)
TPR_lo <- apply(TPR_sims, 2, quantile, 0.025)
TPR_hi <- apply(TPR_sims, 2, quantile, 0.975)
FPR_mean <- apply(FPR_sims, 2, mean)
# Plot ROC with credible band
ord <- order(FPR_mean)
plot(FPR_mean[ord], TPR_mean[ord], type="l",
xlab="False Positive Rate", ylab="True Positive Rate")
polygon(c(FPR_mean[ord], rev(FPR_mean[ord])),
c(TPR_lo[ord], rev(TPR_hi[ord])),
border=NA, col=adjustcolor("gray80", 0.8))
lines(FPR_mean[ord], TPR_mean[ord])
abline(0,1,lty=3)
cat(sprintf("Binormal posterior mean AUC = %.3f (95%% CrI: %.3f–%.3f)\n",
mean(AUC_sims), quantile(AUC_sims,0.025), quantile(AUC_sims,0.975)))
#> Binormal posterior mean AUC = 0.761 (95% CrI: 0.703–0.815)2) Bibeta model in JAGS (ROC & AUC computed in R)
Assume scores are in \([0,1]\) and follow Beta distributions by class. Use the mean–precision reparameterization: \(\alpha = m k\), \(\beta = (1-m)k\) with mean \(m \in (0,1)\) and precision \(k > 0\).
JAGS lacks a Beta CDF, so we compute ROC and AUC in R for each
posterior draw using pbeta and dbeta.
JAGS model file
Save as bibeta_roc.jags.
bibeta_model <- 'model {
for (i in 1:N0) { y0[i] ~ dbeta(alpha0, beta0) }
for (j in 1:N1) { y1[j] ~ dbeta(alpha1, beta1) }
m0 ~ dbeta(1,1)
m1 ~ dbeta(1,1)
k0 ~ dgamma(1.0, 0.1)
k1 ~ dgamma(1.0, 0.1)
alpha0 <- m0 * k0
beta0 <- (1 - m0) * k0
alpha1 <- m1 * k1
beta1 <- (1 - m1) * k1
}'
writeLines(bibeta_model, "bibeta_roc.jags")
cat("Wrote bibeta_roc.jags\n")
#> Wrote bibeta_roc.jagsR driver using jagsUI
# Simulate bibeta data (edit as needed)
set.seed(123)
N0 <- 200; N1 <- 180
m0_true <- 0.35; k0_true <- 20
m1_true <- 0.70; k1_true <- 25
a0_true <- m0_true * k0_true; b0_true <- (1 - m0_true) * k0_true
a1_true <- m1_true * k1_true; b1_true <- (1 - m1_true) * k1_true
y0 <- rbeta(N0, a0_true, b0_true)
y1 <- rbeta(N1, a1_true, b1_true)
# If real data contain exact 0/1, jitter slightly into (0,1)
# eps <- 1e-6; y0 <- pmin(pmax(y0, eps), 1-eps); y1 <- pmin(pmax(y1, eps), 1-eps)
jags_data2 <- list(N0=N0, N1=N1, y0=y0, y1=y1)
params2 <- c("m0","m1","k0","k1","alpha0","beta0","alpha1","beta1")
m0_start <- pmin(pmax(mean(y0), 1e-6), 1-1e-6)
m1_start <- pmin(pmax(mean(y1), 1e-6), 1-1e-6)
inits2 <- function() list(
m0 = plogis(rnorm(1, qlogis(m0_start), 0.5)),
m1 = plogis(rnorm(1, qlogis(m1_start), 0.5)),
k0 = rgamma(1, 1, 0.1),
k1 = rgamma(1, 1, 0.1)
)
fit_bibeta <- jags(
data = jags_data2, inits = inits2,
parameters.to.save = params2,
model.file = "bibeta_roc.jags",
n.chains = 3, n.adapt = 1000, n.iter = 6000,
n.burnin = 2000, n.thin = 2, parallel = TRUE
)
#>
#> Processing function input.......
#>
#> Done.
#>
#> Beginning parallel processing using 3 cores. Console output will be suppressed.
#>
#> Parallel processing completed.
#>
#> Calculating statistics.......
#>
#> Done.
print(fit_bibeta, digits = 3)
#> JAGS output for model 'bibeta_roc.jags', generated by jagsUI.
#> Estimates based on 3 chains of 6000 iterations,
#> adaptation = 1000 iterations (sufficient),
#> burn-in = 2000 iterations and thin rate = 2,
#> yielding 6000 total samples from the joint posterior.
#> MCMC ran in parallel for 0.08 minutes at time 2025-10-03 10:50:42.571661.
#>
#> mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
#> m0 0.348 0.007 0.334 0.348 0.362 FALSE 1 1 6000
#> m1 0.698 0.007 0.685 0.698 0.711 FALSE 1 1 6000
#> k0 20.129 1.963 16.526 20.050 24.176 FALSE 1 1 6000
#> k1 24.078 2.498 19.377 23.967 29.088 FALSE 1 1 6000
#> alpha0 7.002 0.680 5.750 6.980 8.380 FALSE 1 1 6000
#> beta0 13.127 1.307 10.698 13.070 15.792 FALSE 1 1 6000
#> alpha1 16.804 1.774 13.453 16.727 20.393 FALSE 1 1 6000
#> beta1 7.275 0.750 5.869 7.241 8.793 FALSE 1 1 6000
#> deviance -698.997 2.850 -702.557 -699.654 -691.665 FALSE 1 1 6000
#>
#> Successful convergence based on Rhat values (all < 1.1).
#> Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> For each parameter, n.eff is a crude measure of effective sample size.
#>
#> overlap0 checks if 0 falls in the parameter's 95% credible interval.
#> f is the proportion of the posterior with the same sign as the mean;
#> i.e., our confidence that the parameter is positive or negative.
#>
#> DIC info: (pD = var(deviance)/2)
#> pD = 4.1 and DIC = -694.934
#> DIC is an estimate of expected predictive error (lower is better).
# Build ROC/AUC from posterior draws
S <- length(fit_bibeta$sims.list$m0)
tgrid <- seq(0, 1, length.out = 201)
TPR <- FPR <- matrix(NA_real_, nrow = S, ncol = length(tgrid))
AUC <- numeric(S)
x_int <- seq(0, 1, length.out = 1001)
dx <- x_int[2] - x_int[1]
for (s in 1:S) {
a0 <- fit_bibeta$sims.list$alpha0[s]; b0 <- fit_bibeta$sims.list$beta0[s]
a1 <- fit_bibeta$sims.list$alpha1[s]; b1 <- fit_bibeta$sims.list$beta1[s]
FPR[s, ] <- 1 - pbeta(tgrid, a0, b0)
TPR[s, ] <- 1 - pbeta(tgrid, a1, b1)
AUC[s] <- sum(dbeta(x_int, a1, b1) * (1 - pbeta(x_int, a0, b0))) * dx
}
# Direction check: if mean AUC < 0.5, reflect ROC and AUC
if (mean(AUC) < 0.5) {
message("Detected reversed score orientation; reflecting ROC and AUC.")
TPR_old <- TPR; FPR_old <- FPR
TPR <- 1 - FPR_old
FPR <- 1 - TPR_old
AUC <- 1 - AUC
}
#> Detected reversed score orientation; reflecting ROC and AUC.
TPR_mean <- apply(TPR, 2, mean)
TPR_lo <- apply(TPR, 2, quantile, 0.025)
TPR_hi <- apply(TPR, 2, quantile, 0.975)
FPR_mean <- apply(FPR, 2, mean)
ord <- order(FPR_mean)
plot(FPR_mean[ord], TPR_mean[ord], type="l",
xlab="False Positive Rate", ylab="True Positive Rate")
polygon(c(FPR_mean[ord], rev(FPR_mean[ord])),
c(TPR_lo[ord], rev(TPR_hi[ord])),
border=NA, col=adjustcolor("gray80", 0.8))
lines(FPR_mean[ord], TPR_mean[ord])
abline(0,1,lty=3)
cat(sprintf("Bibeta posterior mean AUC = %.3f (95%% CrI: %.3f–%.3f)\n",
mean(AUC), quantile(AUC,0.025), quantile(AUC,0.975)))
#> Bibeta posterior mean AUC = 0.991 (95% CrI: 0.985–0.996)3) Nonparametric Bayesian ROC (Bayesian bootstrap)
This approach makes no parametric assumptions on the score distributions. We place independent Dirichlet(1,…,1) priors on class-specific empirical weights. Each posterior draw yields a weighted empirical ROC and AUC.
R implementation
# Utilities
rdirichlet1 <- function(n) { x <- rexp(n); x / sum(x) }
roc_step_weighted <- function(y0, y1, w0, w1) {
stopifnot(length(y0) == length(w0), length(y1) == length(w1))
w0 <- w0 / sum(w0); w1 <- w1 / sum(w1)
y <- c(y0, y1)
lab <- c(rep(0L, length(y0)), rep(1L, length(y1)))
w <- c(w0, w1)
o <- order(y, decreasing = TRUE)
y <- y[o]; lab <- lab[o]; w <- w[o]
runs <- rle(y)
idx_end <- cumsum(runs$lengths)
idx_start <- c(1L, head(idx_end, -1L) + 1L)
tp <- fp <- numeric(length(runs$lengths))
for (g in seq_along(runs$lengths)) {
sl <- idx_start[g]:idx_end[g]
fp[g] <- sum(w[sl][lab[sl] == 0L])
tp[g] <- sum(w[sl][lab[sl] == 1L])
}
FPR <- c(0, cumsum(fp))
TPR <- c(0, cumsum(tp))
FPR[length(FPR)] <- 1; TPR[length(TPR)] <- 1
list(fpr = FPR, tpr = TPR)
}
auc_weighted <- function(y0, y1, w0, w1) {
w0 <- w0 / sum(w0); w1 <- w1 / sum(w1)
s_all <- sort(unique(c(y0, y1)))
w0_at <- tapply(w0, factor(match(y0, s_all), levels = seq_along(s_all)), sum)
w0_at[is.na(w0_at)] <- 0
w1_at <- tapply(w1, factor(match(y1, s_all), levels = seq_along(s_all)), sum)
w1_at[is.na(w1_at)] <- 0
w0_cum_less <- c(0, head(cumsum(w0_at), -1))
sum(w1_at * (w0_cum_less + 0.5 * w0_at))
}
# Example data (replace with your own)
set.seed(1)
N0 <- 200; N1 <- 180
y0 <- rnorm(N0, 0.0, 1.0)
y1 <- rnorm(N1, 0.7, 1.2)
# Posterior sampling
S <- 2000
fpr_grid <- seq(0, 1, length.out = 201)
TPR_mat <- matrix(NA_real_, nrow = S, ncol = length(fpr_grid))
AUC_vec <- numeric(S)
for (s in 1:S) {
w0 <- rdirichlet1(length(y0))
w1 <- rdirichlet1(length(y1))
roc <- roc_step_weighted(y0, y1, w0, w1)
TPR_mat[s, ] <- approx(x = roc$fpr, y = roc$tpr, xout = fpr_grid,
method = "linear", ties = "ordered", rule = 2)$y
AUC_vec[s] <- auc_weighted(y0, y1, w0, w1)
}
# Direction check: if mean AUC < 0.5, reflect ROC and AUC
if (mean(AUC_vec) < 0.5) {
message("Detected reversed score orientation; reflecting ROC and AUC.")
TPR_mat <- 1 - TPR_mat[, ncol(TPR_mat):1, drop=FALSE]
AUC_vec <- 1 - AUC_vec
}
TPR_mean <- apply(TPR_mat, 2, mean)
TPR_lo <- apply(TPR_mat, 2, quantile, 0.025)
TPR_hi <- apply(TPR_mat, 2, quantile, 0.975)
plot(fpr_grid, TPR_mean, type="l",
xlab="False Positive Rate", ylab="True Positive Rate")
polygon(c(fpr_grid, rev(fpr_grid)),
c(TPR_lo, rev(TPR_hi)),
border = NA, col = adjustcolor("gray80", 0.8))
lines(fpr_grid, TPR_mean)
abline(0,1,lty=3)
cat(sprintf("Bayesian-bootstrap posterior mean AUC = %.3f (95%% CrI: %.3f–%.3f)\n",
mean(AUC_vec), quantile(AUC_vec, 0.025), quantile(AUC_vec, 0.975)))
#> Bayesian-bootstrap posterior mean AUC = 0.688 (95% CrI: 0.636–0.739)DAGs for JAGS models
Binormal model DAG
library(DiagrammeR)
grViz('
digraph {
graph [rankdir=LR, fontsize=10]
node [shape=ellipse, style=filled, fontname=Helvetica]
mu0 [fillcolor="#FDE0DC"]
sigma0 [fillcolor="#FDE0DC"]
tau0 [fillcolor="#FFF8CE"]
y0 [label="y0[ ]", fillcolor="#D0E6FF"]
mu1 [fillcolor="#FDE0DC"]
sigma1 [fillcolor="#FDE0DC"]
tau1 [fillcolor="#FFF8CE"]
y1 [label="y1[ ]", fillcolor="#D0E6FF"]
t [label="t[ ]", fillcolor="#D0E6FF"]
tpr [label="tpr[ ]", fillcolor="#FFF8CE"]
fpr [label="fpr[ ]", fillcolor="#FFF8CE"]
auc [label="auc", fillcolor="#FFF8CE"]
mu0 -> y0
sigma0 -> tau0
tau0 -> y0
mu1 -> y1
sigma1 -> tau1
tau1 -> y1
mu0 -> fpr
tau0 -> fpr
mu1 -> tpr
tau1 -> tpr
t -> fpr
t -> tpr
mu0 -> auc
mu1 -> auc
tau0 -> auc
tau1 -> auc
}
')Bibeta model DAG
library(DiagrammeR)
grViz('
digraph {
graph [rankdir=LR, fontsize=10]
node [shape=ellipse, style=filled, fontname=Helvetica]
m0 [fillcolor="#FDE0DC"]
k0 [fillcolor="#FDE0DC"]
m1 [fillcolor="#FDE0DC"]
k1 [fillcolor="#FDE0DC"]
alpha0 [fillcolor="#FFF8CE"]
beta0 [fillcolor="#FFF8CE"]
alpha1 [fillcolor="#FFF8CE"]
beta1 [fillcolor="#FFF8CE"]
y0 [label="y0[ ]", fillcolor="#D0E6FF"]
y1 [label="y1[ ]", fillcolor="#D0E6FF"]
m0 -> alpha0
k0 -> alpha0
m0 -> beta0
k0 -> beta0
alpha0 -> y0
beta0 -> y0
m1 -> alpha1
k1 -> alpha1
m1 -> beta1
k1 -> beta1
alpha1 -> y1
beta1 -> y1
}
')Reproducibility
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sonoma 14.6.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/Los_Angeles
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] DiagrammeR_1.0.11 jagsUI_1.6.2
#>
#> loaded via a namespace (and not attached):
#> [1] cli_3.6.5 knitr_1.50 rlang_1.1.6 xfun_0.53
#> [5] jsonlite_2.0.0 glue_1.8.0 rjags_4-17 htmltools_0.5.8.1
#> [9] sass_0.4.10 rmarkdown_2.30 grid_4.5.1 evaluate_1.0.5
#> [13] jquerylib_0.1.4 visNetwork_2.1.4 fastmap_1.2.0 yaml_2.3.10
#> [17] lifecycle_1.0.4 compiler_4.5.1 coda_0.19-4.1 RColorBrewer_1.1-3
#> [21] htmlwidgets_1.6.4 rstudioapi_0.17.1 lattice_0.22-7 digest_0.6.37
#> [25] R6_2.6.1 parallel_4.5.1 magrittr_2.0.4 bslib_0.9.0
#> [29] tools_4.5.1 cachem_1.1.0