R/plot.metadiag.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 Charachteristics) space with the
#' posterior predictive contours. The predictive curves are approximated using a non-parametric smoother
#' or with a parametric model. For the parametric model the current implementation supports only a
#' logistic link function. The marginal posterior predictive distributions are ploted outside the ROC space.
#'
#'
#'
#' @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 kde2d.n The number of grid points in each direction for the non-parametric density estimation. Can be scalar or a length-2 inter vector.
#'
#' @param color.line Color of the predictive contour line.
#'
#' @param title Optional parameter for setting a title in the plot.
#'
#' @param marginals Plot the posterior marginal predictive histograms.
#'
#' @param bin.hist Number of bins of the marginal histograms.
#'
#' @param color.hist Color of the histograms.
#'
#' @param S Number of predictive rates to be plotted.
#'
#' @param color.pred.points Color of the posterior predictive rates.
#'
#' @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
#'
#'
#'  plot(glas.m1,                    # Fitted model
#'       level = c(0.5, 0.75, 0.95), # Credibility levels
#'       parametric.smooth = TRUE)   # Parametric curve
#'
#'# Plot results: based on a non-parametric smoother of the posterior predictive rates .......
#'
#' plot(glas.m1,                    # Fitted model
#'      level = c(0.5, 0.75, 0.95), # Credibility levels
#'      parametric.smooth = FALSE)  # Non-parametric curve
#'
#'# Using the pipe command in the package dplyr and changing some colors .......
#'
#'library(dplyr)
#'
#'glas.t %>%
#'  metadiag(re = "normal", re.model ="SeSp") %>%
#'    plot(parametric.smooth = FALSE,
#'           S = 100,
#'           color.data.points = "green",
#'           color.pred.points = "blue",
#'           color.line = "black")
#'
#' }

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


plot.metadiag <- function(x,
         parametric.smooth = TRUE,
                     level = c(0.5, 0.75, 0.95),
                     limits.x = c(0, 1),
                     limits.y = c(0, 1),
                   kde2d.n = 25,
                color.line = "red",
                     title = paste("Posterior Predictive Contours (50%, 75% and 95%)"),
                 marginals = TRUE,
                  bin.hist = 30,
                color.hist = "lightblue",
                         S = 500,
         color.pred.points = "lightblue",
         color.data.points = "blue",
         ...)
  {


  #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.pred <- data.frame(fpr.new = 1 - x$BUGSoutput$sims.list$sp.new,
                         tpr.new = x$BUGSoutput$sims.list$se.new)

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

  dat.sample <- dat.pred[sample(1:S), ]

  #dat.sample$Post.pred <- "sample"


  # Base plot ...............................................................................
  baseplot  <- ggplot(dat.sample, aes_string(x = "fpr.new", y = "tpr.new"))+
    geom_point(aes(color = "Predictions")) +
    scale_color_manual(name='', values = c("Predictions" = color.pred.points),
                       guide = "legend") +
    scale_x_continuous(name = "FPR (1 - Specificity)", limits = limits.x) +
    scale_y_continuous(name = "TPR (Sensitivity)", limits = limits.y) +
    geom_point(data = dat.hat, 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)


  # Non-parametric ..........................................................................
  # Estimation of the nonparametric density ...
  x.nopar <- dat.pred[ ,1]
  y.nopar <- dat.pred[ ,2]

  dens <- kde2d(x.nopar, y.nopar, n = kde2d.n)
  dx <- diff(dens$x[1:2])
  dy <- diff(dens$y[1:2])
  sz <- sort(dens$z)
  c1 <- cumsum(sz) * dx * dy
  Levels.nonpar <- approx(c1, sz, xout = 1 - level)$y


  #nparplot  <- baseplot + stat_density2d()

  densdf <- data.frame(expand.grid(fpr.new = dens$x,
                                   tpr.new = dens$y),
                                         dens.z = as.vector(dens$z))

  # Non-parametric ..........................................................................

  # The NULLing out strategy to pach the R CMD  check:
  #                         "no visible binding for global variable 'dens.z'"

  dens.z <- NULL
  nparplot <- baseplot + geom_contour(data = densdf,
                                       aes(z = dens.z),
                                      colour = color.line,
                                      breaks = Levels.nonpar)




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

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


  if(parametric.smooth == TRUE) finalplot <- parplot
  else finalplot <- nparplot

  if(marginals == TRUE) finalplot <- ggMarginal(finalplot, type= "histogram",
                                                fill = color.hist, bins = bin.hist)

  return(finalplot)
}