Calculate AUC in a JAGS model

auc_in_jags.R

library(R2jags)

# Create some dummy predictor and response data
set.seed(1)
J <- 2
n <- 5000
X <- cbind(1, matrix(runif(n*J, -1, 1), ncol=2))
b <- 1:3
p <- plogis(X %*% b)
y <- rbinom(n, 1, p)

# Define a JAGS model
M <- function() {
  for(i in 1:n) {
    y[i] ~ dbern(p[i])
    logit(p[i]) <- inprod(X[i, ], b)
  }
  for(j in 1:(J+1)) {
    b[j] ~ dnorm(0, 0.0001)
  }
  for (t in 1:length(thr)) {
    sens[t] <- sum((p > thr[t]) && (y==1))/n1
    spec[t] <- sum((p < thr[t]) && (y==0))/n0
    fpr[t] <- 1 - spec[t]
    fnr[t] <- 1 - sens[t]
  }
  # Calculate area of ROC trapezoid (http://stats.stackexchange.com/a/146174/46761)
  auc <- sum((sens[2:length(sens)]+sens[1:(length(sens)-1)])/2 * 
    -(fpr[2:length(fpr)] - fpr[1:(length(fpr)-1)]))
}

# Fit the JAGS model
j <- jags(list(X=X, y=y, n=n, J=J, n1=sum(y), n0=sum(!y), thr=seq(0, 1, 0.05)), 
          NULL, c('b', 'sens', 'spec', 'fpr', 'fnr', 'auc'), M, 3, 1000)

# Calculate AUC from posterior medians of sensitivity and 1-specificity
library(jagstools) # devtools::install_github('johnbaums/jagstools')
sens <- jagsresults(j, 'sens')[, '50%']
fpr <- jagsresults(j, 'fpr')[, '50%']
plot(fpr, sens, type='o', pch=20, 
     xlab='1 - Specificity', ylab='Sensitivity')
sum((sens[-1]+sens[-length(sens)])/2 * -diff(fpr)) # auc

# Compare to posterior median of AUC
jagsresults(j, 'auc')[, '50%']