R/bsroc.R

In bamdit: Bayesian Meta-Analysis of Diagnostic Test Data

#' bsroc
#'
#' This function plots the observed data in the ROC (Receiving Operating Charachteristics) space with the
#' Bayesian SROC (Summary ROC) curve. The predictive curves are approximated using a parametric model.
#'
#' @param m                    The object generated by metadiag.
#' @param level                Credibility levels of the predictive curve.
#' @param title                Optional parameter for setting a title in the plot.
#' @param fpr.x                Grid of values where the conditional distribution is calculated.
#' @param partial.AUC          Automatically calculate the AUC for the observed range of FPRs, default is TRUE.
#'                             When TRUE, the AUC is normalized by the width of the observed FPR interval.
#' @param xlim.bsroc           Graphical limits of the x-axis for the BSROC curve plot.
#' @param ylim.bsroc           Graphical limits of the y-axis for the BSROC curve plot.
#' @param lower.auc            Lower limit of the AUC.
#' @param upper.auc            Upper limit of the AUC.
#' @param col.fill.points      Color used to fill points, default is blue.
#' @param results.bauc         Print results of the Bayesian Area Under the Curve, default value is TRUE.
#' @param results.bsroc        Print results of the Bayesian SROC curve, default value is TRUE.
#' @param plot.outputs         Logical. If TRUE, plots are produced. If FALSE, only numerical results are printed.
#' @param plot.post.bauc       The BSROC and the posterior of the BAUC are ploted in the same page, default is FALSE.
#' @param bins                 Histograms' bins.
#' @param scale.size.area      Scale area for the ploted points, default = 10.
#' @seealso \code{\link{metadiag}}.
#' @keywords file
#'
#' @references Verde P. E. (2010). Meta-analysis of diagnostic test data: A
#' bivariate Bayesian modeling approach. Statistics in Medicine. 29(30):3088-102.
#' doi: 10.1002/sim.4055.
#'
#' @references Verde P. E. (2018). bamdit: An R Package for Bayesian Meta-Analysis
#' of Diagnostic Test Data. Journal of Statistical Software. Volume 86, issue 10, pages 1--32.
#'
#' @examples
#' \dontrun{
#' data(glas)
#' glas.t <- glas[glas$marker == "Telomerase", 1:4]
#' glas.m1 <- metadiag(glas.t, re = "normal", link = "logit")
#' bsroc(glas.m1)
#' bsroc(glas.m1, plot.post.bauc = TRUE)
#' bsroc(glas.m1, plot.outputs = FALSE)
#' bsroc(glas.m1, results.bsroc = FALSE)
#'
#' data(mri)
#' mri.m <- metadiag(mri)
#' bsroc(mri.m)
#' }
#' @import ggplot2
#' @import grid
#' @import gridExtra
#' @import R2jags
#' @import rjags
#' @export
#' @importFrom stats integrate
bsroc <- function(m,
                  level = c(0.05, 0.5, 0.95),
                  title = "Bayesian SROC Curve",
                  fpr.x = seq(0.001, 0.999, by = 0.001),
                  partial.AUC = TRUE,
                  xlim.bsroc = c(0,1),
                  ylim.bsroc = c(0,1),
                  lower.auc = 0,
                  upper.auc = 0.95,
                  col.fill.points = "blue",
                  results.bauc = TRUE,
                  results.bsroc = TRUE,
                  plot.outputs = TRUE,
                  plot.post.bauc = FALSE,
                  bins = 30,
                  scale.size.area = 10)
{
  x <- y.med <- y.up <- y.low <- NULL

  link <- m$link
  link.test <- link == "logit"
  if(!link.test) stop("This link function is not implemented (choose link = logit in metadiag)")

  for(i in seq_along(level)) {
    if(!(0 < level[i] & level[i] < 1)) {
      stop("At least one credibility level is not correct. Use values between 0 and 1!")
    }
  }

  if(!plot.outputs && plot.post.bauc) {
    warning("plot.post.bauc = TRUE ignored because plot.outputs = FALSE")
  }

  # Formatting data
  data <- m$data
  two.by.two <- m$two.by.two

  if(two.by.two == FALSE) {
    tp <- data[,1]
    n1 <- data[,2]
    fp <- data[,3]
    n2 <- data[,4]
  } else {
    tp <- data$TP
    fp <- data$FP
    fn <- data$FN
    tn <- data$TN
    n1 <- tp + fn
    n2 <- fp + tn
  }

  # Data errors
  if(any(tp > n1) || any(fp > n2)) stop("the data is inconsistent")

  if(!missing(data)) {
    tpr <- tp / n1
    fpr <- fp / n2
    n <- n1 + n2
  } else {
    stop("NAs are not alow in this plot function")
  }

  # Stop if BSROC makes no sense
  delta.fpr <- max(fpr) - min(fpr)
  if(delta.fpr < 0.2) {
    warning("The range of the observed FPR is less than 20%. Calculating the BSROC curve makes no sense!")
  }

  # Partial AUC: use observed FPR interval for plotting grid
  if(partial.AUC == TRUE) {
    fpr.x <- seq(min(fpr), max(fpr), by = 0.001)
  }

  # data frame of observed rates
  dat.hat <- data.frame(
    tpr = tpr,
    fpr = fpr,
    n = n
  )

  # Parametric
  param <- m$re.model
  param.test <- param == "DS"

  if(param.test == TRUE) {
    A <- m$BUGSoutput$sims.list$mu.D
    sigma.D <- m$BUGSoutput$sims.list$sigma.D
    sigma.S <- m$BUGSoutput$sims.list$sigma.S
    rho <- m$BUGSoutput$sims.list$rho
    B <- rho * sigma.D / sigma.S
  } else {
    mu.Se <- m$BUGSoutput$sims.list$mu.Se
    mu.Sp <- m$BUGSoutput$sims.list$mu.Sp
    rho <- m$BUGSoutput$sims.list$rho
    sigma.Se <- m$BUGSoutput$sims.list$sigma.Se
    sigma.Sp <- m$BUGSoutput$sims.list$sigma.Sp

    B <- -rho * sigma.Se / sigma.Sp
    A <- mu.Se + B * mu.Sp
  }

  # BSROC
  if(param.test == TRUE) {
    f <- function(FPR, A, B) {
      x <- A / (1 - B) + (B + 1) / (1 - B) * log(FPR / (1 - FPR))
      f.value <- exp(x) / (1 + exp(x))
      return(f.value)
    }
  } else {
    f <- function(FPR, A, B) {
      x <- A + B * log(FPR / (1 - FPR))
      f.value <- exp(x) / (1 + exp(x))
      return(f.value)
    }
  }

  # Simple graph for the BSROC
  BSROC <- Vectorize(f, vectorize.args = c("A", "B"))
  BSROC.out <- BSROC(A, B, FPR = fpr.x)
  bsrocCI <- apply(BSROC.out, 1, quantile, prob = sort(level), na.rm = TRUE)
  bsrocCI <- t(bsrocCI)

  dat.bsroc <- data.frame(
    x = fpr.x,
    y.low = bsrocCI[,1],
    y.med = bsrocCI[,2],
    y.up = bsrocCI[,3]
  )

  bsroc.plot <- ggplot(data = dat.bsroc, aes(x = x, y = y.med)) +
    geom_line(aes(x = x, y = y.up)) +
    geom_line(aes(x = x, y = y.med)) +
    geom_line(aes(x = x, y = y.low)) +
    geom_point(data = dat.hat, aes(size = n, x = fpr, y = tpr),
               shape = 21, alpha = 0.35, fill = col.fill.points) +
    scale_x_continuous(name = "FPR (1 - Specificity)", limits = xlim.bsroc) +
    scale_y_continuous(name = "TPR (Sensitivity)", limits = ylim.bsroc) +
    scale_size_area(max_size = scale.size.area) +
    ggtitle(title)

  # Bayesian Area Under the Curve (BAUC)
  if(partial.AUC == TRUE) {
    lower.auc <- min(fpr)
    upper.auc <- max(fpr)
    auc.width <- upper.auc - lower.auc
    auc.normalized <- TRUE
  } else {
    lower.auc <- min(fpr.x)
    upper.auc <- max(fpr.x)
    auc.width <- upper.auc - lower.auc
    auc.normalized <- FALSE
  }

  # BAUC function
  area <- function(A, B, lower.auc, upper.auc) {
    val <- integrate(f,
                     lower = lower.auc,
                     upper = upper.auc,
                     A = A,
                     B = B)$value
    return(val)
  }

  # Calculations restricted for the slope -1 < B < 1
  index.B.range <- B < 0.99 & B > -0.99
  A.auc <- A[index.B.range]
  B.auc <- B[index.B.range]

  v.area <- Vectorize(area)
  bauc <- v.area(A.auc, B.auc, lower.auc, upper.auc)

  # Normalize only for partial AUC on observed FPR interval
  if(auc.normalized) {
    if(auc.width <= 0) {
      stop("The observed FPR interval has zero width, so the normalized partial AUC cannot be computed.")
    }
    bauc <- bauc / auc.width
  }

  # Posteriors
  A.mean <- round(mean(A.auc), 3)
  A.sd <- round(sd(A.auc), 3)
  A.ci <- round(quantile(A.auc, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), na.rm = TRUE), 3)
  A.post <- cbind(A.mean, A.sd, t(A.ci))

  B.mean <- round(mean(B.auc), 3)
  B.sd <- round(sd(B.auc), 3)
  B.ci <- round(quantile(B.auc, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), na.rm = TRUE), 3)
  B.post <- cbind(B.mean, B.sd, t(B.ci))

  post.AB <- rbind(A.post, B.post)
  rownames(post.AB) <- c("A", "B")
  colnames(post.AB) <- c("Mean", "SD", "2.5%", "25%", "50%", "75%", "97.5%")

  bauc.mean <- round(mean(bauc), 3)
  bauc.sd <- round(sd(bauc), 3)
  bauc.ci <- round(quantile(bauc,
                            probs = c(0.025, 0.25, 0.5, 0.75, 0.975),
                            na.rm = TRUE), 3)
  bauc.post <- cbind(bauc.mean, bauc.sd, t(bauc.ci))
  colnames(bauc.post) <- c("Mean", "SD", "2.5%", "25%", "50%", "75%", "97.5%")
  rownames(bauc.post) <- "BAUC"

  # Print model information only when at least one block of results is requested
  if(results.bsroc || results.bauc) {
    cat("\n")
    cat("------------------------------------------------------------", "\n")
    cat("These results are based on the following random effects model:", "\n")
    cat("------------------------------------------------------------", "\n")
    cat("Link function: ", m$link, "\n")
    cat("Random Effects distribution: ",
        ifelse(m$re == "normal", "Bivariate Normal", "Bivariate Scale Mixture"), "\n")
    cat("Parametrization: ",
        ifelse(m$re.model == "DS", "Differences and Sums", "Sensitivities and Specificities"), "\n", sep = "")
    cat("Splitting study weights: ", ifelse(m$split.w, "Yes", "No"), "\n")
    cat("------------------------------------------------------------", "\n")
    cat("\n")
  }

  # Print BSROC results
  if(results.bsroc == TRUE) {
    cat("------------------------------------------------------------", "\n")
    cat("Posteriors for the parameters of the Bayesian SROC Curve", "\n")
    cat("------------------------------------------------------------", "\n")
    print(post.AB)
    cat("------------------------------------------------------------", "\n")
    cat("\n")
  }

  # Print BAUC results
  if(results.bauc == TRUE) {
    cat("------------------------------------------------------------", "\n")
    cat("Summary results for the Bayesian Area Under the Curve (BAUC)", "\n")

    if(auc.normalized) {
      cat("(Partial AUC restricted to the observed FPR range and normalized)", "\n")
      cat("Observed FPR interval: [", round(lower.auc, 3), ", ", round(upper.auc, 3), "]", sep = "", "\n")
    } else {
      cat("(AUC computed on the full selected FPR range)", "\n")
      cat("Integration interval: [", round(lower.auc, 3), ", ", round(upper.auc, 3), "]", sep = "", "\n")
    }

    cat("------------------------------------------------------------", "\n")
    print(bauc.post)
    cat("------------------------------------------------------------", "\n")
  }

  post.bauc <- ggplot(data.frame(bauc), aes(x = bauc)) +
    geom_histogram(colour = "black", fill = "dodgerblue", bins = bins) +
    xlab(ifelse(auc.normalized,
                "Normalized Bayesian Area Under the Curve",
                "Bayesian Area Under the Curve")) +
    geom_vline(xintercept = median(bauc)) +
    geom_vline(xintercept = quantile(bauc, prob = c(0.025, 0.975)),
               linetype = "longdash") +
    ggtitle("Posterior Distribution")

  if(plot.outputs) {
    if(plot.post.bauc == TRUE) {
      grid.newpage()
      pushViewport(viewport(layout = grid.layout(1, 2)))

      vplayout <- function(x, y) {
        viewport(layout.pos.row = x, layout.pos.col = y)
      }

      print(bsroc.plot, vp = vplayout(1, 1))
      print(post.bauc, vp = vplayout(1, 2))
    } else {
      print(bsroc.plot)
    }
  }
}