R/plotcredibility.R

In bamdit: Bayesian Meta-Analysis of Diagnostic Test Data

#' Generic plot function for metadiag object in bamdit
#'
#' This function plots the observe data in the ROC (Receiving Operating Characteristics) space with the
#' posterior credibility contours.
#'
#'
#'
#' @param x The object generated by the metadiag function.
#'
#' @param parametric.smooth Indicates if the predictive curve is a parametric or non-parametric.
#'
#' @param limits.x Numeric vector of length 2 specifying the x-axis limits. The default value is c(0, 1).
#'
#' @param limits.y Numeric vector of length 2 specifying the x-axis limits. The default value is c(0, 1).
#'
#' @param level Credibility levels of the predictive curve. If parametric.smooth = FALSE, then the probability levels are estimated from the nonparametric surface.
#'
#'
#' @param color.line Color of the predictive contour line.
#'
#' @param title Optional parameter for setting a title in the plot.
#'
#' @param color.data.points Color of the data points.
#'
#' @param ... \dots
#'
#'
#' @seealso \code{\link{metadiag}}.
#'
#'
#' @examples
#'
#'
#'
#' \dontrun{
#' library(bamdit)
#' data("glas")
#' glas.t <- glas[glas$marker == "Telomerase", 1:4]
#' glas.m1 <- metadiag(glas.t,                # Data frame
#'                     re = "normal",         # Random effects distribution
#'                     re.model = "DS",       # Random effects on D and S
#'                     link = "logit",        # Link function
#'                     sd.Fisher.rho   = 1.7, # Prior standard deviation of correlation
#'                     nr.burnin = 1000,      # Iterations for burnin
#'                     nr.iterations = 10000, # Total iterations
#'                     nr.chains = 2,         # Number of chains
#'                     r2jags = TRUE)         # Use r2jags as interface to jags
#'
#'
#'  plotcredibility(glas.m1,         # Fitted model
#'       level = c(0.5, 0.75, 0.95), # Credibility levels
#'       parametric.smooth = TRUE)   # Parametric curve
#' }

#' @importFrom stats pnorm cor qchisq sd approx
#' @importFrom MASS kde2d
#' @import ggplot2
#' @import R2jags
#' @import rjags
#' @import ggExtra
#' @export


plotcredibility <- function(x,
                          parametric.smooth = TRUE,
                          level = c(0.5, 0.75, 0.95),
                          limits.x = c(0, 1),
                          limits.y = c(0, 1),
                          color.line = "red",
                          color.data.points = "blue",
                          title = paste("Posterior Credibility Contours (50%, 75% and 95%)"),
                          ...)
{


  #if(class(m)!="rjags")stop("You have to provide a valid rjags object as fitted model.")

  link <- x$link

  link.test <- link %in% c("logit", "cloglog", "probit")
  if(!link.test)stop("This link function is not implemented.")

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

  # Predictive curve for model m1 ............................................................
  # Patch for the "note" no-visible-binding-for-global-variable

  dat.cred <- data.frame(fpr.new = 1 - x$BUGSoutput$sims.list$sp,
                         tpr.new = x$BUGSoutput$sims.list$se)

  # Formatting data
  data <- x$data
  two.by.two <- x$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")

  dat.hat <- data.frame(tpr = tpr,
                        fpr = fpr,
                        n = n)

  # Base plot ...............................................................................
  baseplot  <- ggplot(dat.hat, aes_string(x = "fpr", y = "tpr"))+
    scale_x_continuous(name = "FPR (1 - Specificity)", limits = limits.x) +
    scale_y_continuous(name = "TPR (Sensitivity)", limits = limits.y) +
    geom_point(aes_string(x = "fpr", y = "tpr", size = "n"),
               shape = 21, alpha = 0.35, fill = color.data.points) +
    scale_size_area(max_size = 10) +
    ggtitle(title)



  # Parametric ..............................................................................
  if(link == "logit"){
    x <- with(dat.cred, log(fpr /(1-fpr )))
    y <- with(dat.cred, log(tpr /(1-tpr )))
  } else if(link == "cloglog")
  {
    x <- with(dat.cred, log(-log(1 - fpr )))
    y <- with(dat.cred, log(-log(1 - tpr )))
  } else if(link == "probit")
  {
    x <- with(dat.cred, qnorm(fpr ))
    y <- with(dat.cred, qnorm(tpr ))
  }

  mean.invlink.x <- mean(x)
  mean.invlink.y <- mean(y)
  sd.invlink.x <- sd(x)
  sd.invlink.y <- sd(y)
  rho.invlink.x.y <- cor(x,y)

  parplot <- baseplot
  for(k in level)
  {
    # Set up the curve 95%...
    cc <- sqrt(qchisq(k, 2))
    tt <- seq(0, 2*pi, len=100)

    mu.x <-  mean(mean.invlink.x) + mean(sd.invlink.x) * cc * cos(tt)
    mu.y <-  mean(mean.invlink.y) + mean(sd.invlink.y) * cc * cos( tt + acos(mean(rho.invlink.x.y)))

    # Pooled summaries
    if(link == "logit"){
      TPR1 <- 1 / ( 1 + exp(-1*mu.y) )  # with logit link
      FPR1 <- 1 / ( 1 + exp(-1*mu.x) )  # with logit link
    } else if(link == "cloglog")
    {
      TPR1 <- 1 - exp(-1*exp(mu.y))      # with cloglog link
      FPR1 <- 1 - exp(-1*exp(mu.x))      # with cloglog link

    } else if(link == "probit")
    {
      TPR1 <- pnorm(mu.y)               # with logit link
      FPR1 <- pnorm(mu.x)               # with logit link
    }

    curve.alpha <- data.frame(FPR1, TPR1)
    parplot <- parplot + geom_path(data = curve.alpha, aes_string(x = "FPR1", y = "TPR1"),
                                   colour = color.line)
  }


 finalplot <- parplot

return(finalplot)
}