Make a simple ROC curve and calculate area under ROC

easyRoc.R

# function
# 'screen' is numeric vector of scores from the screening tool
# 'gs' is a logical vector showing gold standard disease status
# 'vals' is a numeric vector of screening tool values to test (including all values in 'screen')
# returns a ROC plot and a list where the first value is the sensitivity and specificity for all values of 'val', and the second is the area under ROC

easyRoc <- function(screen, gs, vals) { 
  pos <- screen >= rep(vals, each = length(screen))
  pos <- matrix(pos, ncol = length(vals))
  sens <- colSums(gs & pos) / sum(gs)
  spec <- colSums(!gs & !pos) / sum(!gs)
  plot(0, ylim = c(0, 1), xlim = c(0, 1), type = 'n', xlab = 'Specificity', ylab = 'Sensitivity', axes = F)
  lines(1-spec, sens, lwd = 1.5)
  aurx <- c(1, 1, 1-spec)
  aury <- c(0, 1, sens)
  polygon(aurx, aury, col = 'grey93')
  segments(0, 0, 1, 1, lty = 3)
  rect(0, 0, 1, 1)
  axis(1, seq(0, 1, 0.2), seq(1, 0, -0.2), pos = 0)
  axis(2, seq(0, 1, 0.2), pos = 0, las = 2)
  n <- length(aurx)
  a <- sum(aurx[1:(n - 1)] * aury[2:n]) + aurx[n] * aury[1]
  b <- sum(aurx[2:n] * aury[1:(n - 1)]) + aurx[1] * aury[n]
  aur <- 0.5 * (a - b)
  text(0.8, 0.2, paste0('AUR=',round(aur, 2)))
  return(list(res = cbind(cutoff = vals, sens, spec), AUR = aur))
}

# example

dat <- data.frame(screen = c(rbinom(100, 50, 0.2), rbinom(100, 50, 0.5)),
                  disease = c(rbinom(100, 1, 0.1), rbinom(100, 1, 0.6)))
boxplot(screen ~ disease, dat)

par(mar = c(4, 4, 1, 1))
y <- easyRoc(dat$screen, dat$disease, 0:50)