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 conditionl distribution is calculated.
#' @param partial.AUC          Automatically calculate the AUC for the observed range of FPRs, default is TRUE.
#' @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 FALSE.
#' @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 Statisticsl Software. Volume 86, issue 10, pages 1--32.
#'
#'
#' @examples
#'
#' ## execute analysis
#' \dontrun{
#' # Example: data from Glas et al. (2003).....................................
#'
#' 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)
#'
#' # Example: data from Scheidler et al. (1997)
#' # In this example the range of the observed FPR is less than 20%.
#' # Calculating the BSROC curve makes no sense! You will get a warning message!
#'
#' 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.01, 0.95, 0.01),
                  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 = FALSE,
                  plot.post.bauc = FALSE,
                  bins = 30,
                  scale.size.area = 10)
{

  link <- m$link
  link.test <- link == "logit"      #%in%  c("logit", "cloglog", "probit")
  if(!link.test)stop("This link function is not implemented (choose link = logit in metadiag)")

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



  # 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
  }


  if(tp>n1 || 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..............................................................................
  if(partial.AUC==TRUE) fpr.x <-  seq(min(fpr), max(fpr), 0.01)


  # 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){
  # Patch for the "note" no-visible-binding-for-global-variable
        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
        {A <- m$BUGSoutput$sims.list$mu.Se
  rho.SeSp <- m$BUGSoutput$sims.list$rho
  sigma.Se <- m$BUGSoutput$sims.list$sigma.Se
  sigma.Sp <- m$BUGSoutput$sims.list$sigma.Sp
         B <- sigma.Se/sigma.Sp*rho.SeSp
         }


  ###############################################################################
  # BSROC
  # Problem here ... it must return a bound value!!
  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 =T)
    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_string(x = "x", y = "y.med")) +
                   geom_line(aes_string(x = "x", y = "y.up")) +
                   geom_line(aes_string(x = "x", y = "y.med")) +
                   geom_line(aes_string(x = "x", y = "y.low")) +
                   geom_point(data = dat.hat, aes_string(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)
  }

  # 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 (Verde 2008 pag 11.).....

  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)

  # 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 results .........................................................

  if(results.bauc==TRUE){
  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")
  cat("Splitting study weights: ", ifelse(m$split.w, "Yes", "No"), "\n")
  cat("------------------------------------------------------------","\n")
  cat("\n")
  cat("------------------------------------------------------------","\n")
  cat("Posteriors for the parameters of the Bayesian SROC Curve","\n")
  cat("------------------------------------------------------------","\n")
   print(post.AB)
  cat("\n")
  cat("------------------------------------------------------------","\n")
  cat("Summary results for the Bayesian Area Under the Curve (BAUC)","\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("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.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)
  }
}