R/PlotRsmOperatingCharacteristics.R

In RJafroc: Artificial Intelligence Systems and Observer Performance

#' RSM predicted operating characteristics, ROC pdfs and AUCs
#' 
#' @description Visualize RSM predicted ROC, AFROC, wAFROC and FROC curves, and ROC pdfs, 
#'    given \bold{equal-length arrays} of search model parameters: mu, lambda, nu and 
#'    zeta1.
#' 
#' @param mu Array: the RSM mu parameter.
#' 
#' @param lambda Array: the RSM lambda parameter.
#' 
#' @param nu Array: the RSM nu parameter.
#' 
#' @param zeta1 Array, the lowest reporting threshold; if missing the default 
#'    is an array of -3s. 
#'
#' @param lesDistr Array: the probability mass function of the 
#'    lesion distribution for diseased cases. The default is 1. See \link{UtilLesionDistrVector}. 
#' 
#' @param relWeights The relative weights of the lesions; a vector of 
#'    length equal to \code{length(maxLL)}. The default is zero, in which case
#'    equal weights are assumed.
#' 
#' @param  OpChType The type of operating characteristic desired: can be "\code{ROC}", 
#'    "\code{AFROC}", "\code{wAFROC}", "\code{FROC}" or "\code{pdfs}" or "\code{ALL}". 
#'    The default is "\code{ALL}".
#' 
#' @param legendPosition The positioning of the legend: "\code{right}", "\code{left}", 
#'    "\code{top}" or "\code{bottom}". Use "\code{none}" to suppress the legend.
#' 
#' @param legendDirection Allows control on the direction of the legend; 
#'    \code{"horizontal"}, the default, or \code{"vertical"}
#' 
#' @param legendJustification Where to position the legend, default 
#' is bottom right corner c(0,1)
#' 
#' @param nlfRange \bold{This applies to FROC plot only}. The x-axis range, e.g., c(0,2), 
#'    for FROC plot. Default is "\code{NULL}", which means the maximum NLF range, 
#'    as determined by the data.
#' 
#' @param llfRange \bold{This applies to FROC plot only}. The y-axis range, e.g., c(0,1), 
#'    for FROC plot. Default is "\code{NULL}", which means the maximum LLF range,
#'    as determined by the data.
#' 
#' @param nlfAlpha Upper limit of the integrated area under the FROC plot. 
#'    Default is "\code{NULL}", which means the maximum NLF range is used 
#'    (i.e., lambda). Attempt to integrate outside the maximum NLF will 
#'    generate an error.
#' 
#' 
#' @return A list of elements containing five \pkg{ggplot2} objects 
#'    (ROCPlot, AFROCPlot wAFROCPlot, FROCPlot and PDFPlot) and two area measures 
#'    (each of which can have up to two elements), the area under the search 
#'    model predicted ROC curves in up to two treatments, the area under the search 
#'    model predicted AFROC curves in up to two treatments, the area under the 
#'    search model predicted wAFROC curves in up to two treatments, the area under 
#'    the search model predicted FROC curves in up to two treatments.
#' \itemize{
#' \item{\code{ROCPlot}}   {The predicted ROC plots}
#' \item{\code{AFROCPlot}}     {The predicted AFROC plots}
#' \item{\code{wAFROCPlot}}    {The predicted wAFROC plots}
#' \item{\code{FROCPlot}}  {The predicted FROC plots}
#' \item{\code{PDFPlot}}   {The predicted ROC pdf plots, highest rating generated}
#' \item{\code{aucROC}}    {The predicted ROC AUCs, highest rating generated}
#' \item{\code{aucAFROC}}  {The predicted AFROC AUCs}
#' \item{\code{aucwAFROC}} {The predicted wAFROC AUCs}
#' \item{\code{aucFROC}}   {The predicted FROC AUCs}
#' }
#' 
#' @details RSM is the Radiological Search Model described in the book. This 
#'    function is vectorized with respect to the first 4 arguments. For 
#'    \code{lesDistr} the sum must be one. To indicate that all dis. cases
#'    contain 4 lesions, set lesDistr = c(0,0,0,1).  
#' 
#'   
#' 
#' 
#' @import ggplot2
#' 
#' @importFrom stats integrate dbinom dnorm pnorm qnorm
#' 
#' @references
#' Chakraborty DP (2006) A search model and figure of merit for observer data acquired according to the free-response 
#' paradigm, Phys Med Biol 51, 3449-3462.
#' 
#' Chakraborty DP (2006) ROC Curves predicted by a model of visual search, Phys Med Biol 51, 3463--3482.
#'
#' Chakraborty, DP, Yoon, HJ (2008) Operating characteristics predicted by models for diagnostic tasks involving lesion localization, Med Phys, 35:2, 435.
#' 
#' Chakraborty DP (2017) \emph{Observer Performance Methods for Diagnostic Imaging - Foundations, 
#' Modeling, and Applications with R-Based Examples} (CRC Press, Boca Raton, FL). 
#' \url{https://www.routledge.com/Observer-Performance-Methods-for-Diagnostic-Imaging-Foundations-Modeling/Chakraborty/p/book/9781482214840}
#' 
#' @examples
#' ## Following example is for mu = 2, lambda = 1, nu = 0.6, in one treatment and   
#' ## mu = 3, lambda = 1.5, nu = 0.8, in the other treatment. 20% of the diseased 
#' ## cases have a single lesion, 40% have two lesions, 10% have 3 lesions, 
#' ## and 30% have 4 lesions.  
#' 
#' PlotRsmOperatingCharacteristics(mu = c(2, 3), lambda = c(1, 1.5), nu = c(0.6, 0.8),
#'    lesDistr = c(0.2, 0.4, 0.1, 0.3), legendPosition = "bottom")
#' 
#' @export
#' 
PlotRsmOperatingCharacteristics <- function(mu, 
                                            lambda, 
                                            nu,
                                            zeta1,
                                            lesDistr = 1, 
                                            relWeights = 0, 
                                            OpChType = "ALL", 
                                            legendPosition = c(1,0), 
                                            legendDirection = "horizontal", 
                                            legendJustification = c(0,1),
                                            nlfRange = NULL, 
                                            llfRange = NULL, 
                                            nlfAlpha = NULL){
  if (missing(zeta1)) zeta1 <- array(-3, dim = length(mu))
      
  if (!all(c(length(mu) == length(lambda), length(mu) == length(nu), length(mu) == length(zeta1))))
    stop("Parameters mu, lambda, nu and zeta1 have different lengths.")
  
  lesWghtDistr <- UtilLesionWeightsMatrixLesDistr(lesDistr, relWeights)

  plotStep <- 0.01
  # begin bug fix 12/7/21
  # plotStep <- 0.1 # delete after debug
  zeta <- array(list(), dim = length(mu))
  for (i in 1:length(mu)) {
    if (zeta1[i] == -Inf) temp1 <- -3 else temp1 <- zeta1[i]  
    zeta[[i]] <- seq(from = temp1, to = max(mu)+5, by = plotStep)
  }
  # end bug fix 12/7/21
  
  ROCPlot <- NA
  FROCPlot <- NA
  AFROCPlot <- NA
  wAFROCPlot <- NA
  PDFPlot <- NA
  ROCPoints <- data.frame(FPF = NULL, TPF = NULL, Treatment = NULL, stringsAsFactors = FALSE)
  ROCDashes <- data.frame(FPF = NULL, TPF = NULL, Treatment = NULL, stringsAsFactors = FALSE)
  FROCPoints <- data.frame(NLF = NULL, LLF = NULL, Treatment = NULL, stringsAsFactors = FALSE)
  AFROCPoints <- data.frame(FPF = NULL, LLF= NULL, Treatment = NULL, stringsAsFactors = FALSE)
  AFROCDashes <- data.frame(FPF = NULL, LLF= NULL, Treatment = NULL, stringsAsFactors = FALSE)
  wAFROCPoints <- data.frame(FPF = NULL, wLLF= NULL, Treatment = NULL, stringsAsFactors = FALSE)
  wAFROCDashes <- data.frame(FPF = NULL, wLLF= NULL, Treatment = NULL, stringsAsFactors = FALSE)
  abnPDFPoints <- data.frame(pdf = NULL, highestZSample = NULL, Treatment = NULL, stringsAsFactors = FALSE)
  norPDFPoints <- data.frame(pdf = NULL, highestZSample = NULL, Treatment = NULL, stringsAsFactors = FALSE)
  aucROC <- rep(NA, length(mu));aucAFROC <- aucROC;aucwAFROC <- aucROC;aucFROC <- aucROC
  
  for (i in 1:length(mu)){
    if (mu[i] <= 0 ) stop("mu must be greater than zero")
    if (lambda[i] < 0 ) stop("lambda must be greater than zero")
    if (nu[i] < 0 ) stop("nu must be greater than zero")
    
    FPF <- sapply(zeta[[i]], xROC, lambda[i]) # C++ function uses lambda
    TPF <- sapply(zeta[[i]], yROC, mu = mu[i], lambda[i], nu[i], lesDistr)
    NLF <- sapply(zeta[[i]], RSM_NLF, lambda[i])
    LLF <- sapply(zeta[[i]], RSM_LLF, mu = mu[i], nu[i]) # 9/22/22
# begin bug fix
# found bug 11/24/21 two places, here and as indicated by # bug fix 11/24/21
# found ROC AUC did not change with zeta1 as it should   
#    maxFPF <- xROC(-20, lambda[i]) # this is wrong
    maxFPF <- xROC(zeta1[i], lambda[i]) # this is correct
    if( OpChType == "ALL" ||  OpChType == "ROC"){
      ROCPoints <- rbind(ROCPoints, data.frame(FPF = FPF, 
                                               TPF = TPF, 
                                               Treatment = as.character(i), 
                                               stringsAsFactors = FALSE))
      ROCDashes <- rbind(ROCDashes, data.frame(FPF = c(FPF[1], 1), 
                                               TPF = c(TPF[1], 1), 
                                               Treatment = as.character(i), 
                                               stringsAsFactors = FALSE))
#       maxTPF <- yROC(-20, mu[i], lambda, nu, lesDistr)
      maxTPF <- yROC(zeta1[i], mu[i], lambda[i], nu[i], lesDistr)
# end bug fix
      AUC <- integrate(y_ROC_FPF, 0, maxFPF, mu = mu[i], lambda[i], nu[i], lesDistr =lesDistr)$value
      aucROC[i] <- AUC + (1 + maxTPF) * (1 - maxFPF) / 2
    }
    
    if( OpChType == "ALL" ||  OpChType == "FROC"){
      FROCPoints <- rbind(FROCPoints, data.frame(NLF = NLF, 
                                                 LLF = LLF, 
                                                 Treatment = as.character(i), 
                                                 stringsAsFactors = FALSE))
      if (is.null(nlfAlpha)){
        maxNLF <- max(NLF)
        aucFROC[i] <- integrate(integrandFROC, 0, maxNLF, mu= mu[i], lambda[i], nu[i])$value
      }else{
        maxNLF <- max(NLF)
        if (nlfAlpha <= maxNLF){
          aucFROC[i] <- integrate(integrandFROC, 0, nlfAlpha, mu= mu[i], lambda[i], nu[i])$value
        }else{
          stop("nlfAlpha cannot be greater than the maximum of NLF.")
        }
      }
    }
    
    if( OpChType == "ALL" ||  OpChType == "AFROC"){
      AFROCPoints <- rbind(AFROCPoints, data.frame(FPF = FPF, 
                                                   LLF = LLF, 
                                                   Treatment = as.character(i), 
                                                   stringsAsFactors = FALSE))
      AFROCDashes <- rbind(AFROCDashes, data.frame(FPF = c(FPF[1], 1), 
                                                   LLF = c(LLF[1], 1), 
                                                   Treatment = as.character(i), 
                                                   stringsAsFactors = FALSE))
      maxLLF <- RSM_LLF(zeta1[i], mu[i], nu[i]) # bug fix 9/22/22
      AUC <- integrate(y_AFROC_FPF, 0, maxFPF, mu = mu[i], lambda[i], nu[i])$value
      aucAFROC[i] <- AUC + (1 + maxLLF) * (1 - maxFPF) / 2
    }
    
    if( OpChType == "ALL" ||  OpChType == "wAFROC"){
      wLLF <- sapply(zeta[[i]], ywAFROC, mu[i], nu[i], lesDistr, lesWghtDistr) # cpp code requires lesWghtDistr
      wAFROCPoints <- rbind(wAFROCPoints, data.frame(FPF = FPF, 
                                                     wLLF = wLLF, 
                                                     Treatment = as.character(i), 
                                                     stringsAsFactors = FALSE))
      wAFROCDashes <- rbind(wAFROCDashes, data.frame(FPF = c(FPF[1], 1), 
                                                     wLLF = c(wLLF[1], 1), 
                                                     Treatment = as.character(i), 
                                                     stringsAsFactors = FALSE))
      # maxWLLF <- ywAFROC(zeta1[i], mu[i], nu[i], lesDistr, lesWghtDistr) # bug fix 11/24/21
      maxWLLF <- RSM_wLLF(zeta1[i], mu[i], nu[i], lesDistr, relWeights) # bug fix 11/24/21
      AUC <- integrate(y_wAFROC_FPF, 0, maxFPF, mu = mu[i], lambda[i], nu[i], lesDistr, relWeights)$value
      aucwAFROC[i] <- AUC + (1 + maxWLLF) * (1 - maxFPF) / 2
    }
    
    if( OpChType == "ALL" ||  OpChType == "pdfs"){
      deltaFPF <- FPF[1:(length(FPF) - 1)] - FPF[2:length(FPF)]  
      if( OpChType == "ALL" ||  OpChType == "pdfs"){     
        pdfNor <- deltaFPF / plotStep
        norPDFPoints <- rbind(norPDFPoints, 
                              data.frame(pdf = pdfNor[pdfNor > 1e-6], 
                                         highestZSample = zeta[[i]][-1][pdfNor > 1e-6], 
                                         Treatment = as.character(i), 
                                         class = "non-diseased", 
                                         stringsAsFactors = FALSE))
        deltaTPF <- TPF[1:(length(TPF) - 1)] - TPF[2:length(TPF)]
        pdfAbn <- deltaTPF / plotStep
        abnPDFPoints <- rbind(abnPDFPoints, 
                              data.frame(pdf = pdfAbn[pdfAbn > 1e-6], 
                                         highestZSample = zeta[[i]][-1][pdfAbn > 1e-6], 
                                         Treatment = as.character(i), 
                                         class = "diseased", 
                                         stringsAsFactors = FALSE))
      }
    }
  }
  
  if( OpChType == "ALL" ||  OpChType == "ROC") {
    ROCPlot <- with(ROCPoints, {
      ggplot(data = ROCPoints) + 
        geom_line(aes(x = FPF, y = TPF, color = Treatment))  +       
        geom_line(data = ROCDashes, aes(x = FPF, y = TPF, color = Treatment), linetype = 2) +       
        theme(legend.position = legendPosition, legend.direction = legendDirection, legend.justification = c(1, 0)) 
    })
  }
  
  if( OpChType == "ALL" ||  OpChType == "FROC"){
    FROCPlot <- with(FROCPoints, {
      ggplot(data = FROCPoints) + 
        geom_line(aes(x = NLF, y = LLF, color = Treatment))  +       
        scale_x_continuous(expand = c(0, 0), limits = nlfRange) + 
        scale_y_continuous(expand = c(0, 0), limits = llfRange) + 
        theme(legend.position = legendPosition, legend.direction = legendDirection, legend.justification = c(1, 0)) 
    })
  }
  
  if( OpChType == "ALL" ||  OpChType == "AFROC"){
    AFROCPlot <- with(AFROCPoints, {
      ggplot(data = AFROCPoints) + 
        geom_line(aes(x = FPF, y = LLF , color = Treatment)) + 
        geom_line(data = AFROCDashes, aes(x = FPF, y = LLF, color = Treatment), linetype = 2) +       
        theme(legend.position = legendPosition, legend.direction = legendDirection, legend.justification = c(1, 0)) 
    }
    )
  }
  
  if( OpChType == "ALL" ||  OpChType == "wAFROC"){
    wAFROCPlot <- with(wAFROCPoints, {
      ggplot(data = wAFROCPoints) + 
        geom_line(aes(x = FPF, y = wLLF , color = Treatment)) + 
        geom_line(data = wAFROCDashes, aes(x = FPF, y = wLLF, color = Treatment), linetype = 2) +       
         theme(legend.position = legendPosition, legend.direction = legendDirection, legend.justification = c(1, 0)) 
    })
  }
  
  if( OpChType == "ALL" ||  OpChType == "pdfs"){
    if (legendPosition == "top" || legendPosition == "bottom"){
      legendDirection = "horizontal"
    }else{
      legendDirection = "vertical"
    }
    PDFPoints <- rbind(norPDFPoints, abnPDFPoints)
    PDFPlot <- with(PDFPoints, {
      ggplot(data = PDFPoints, 
             aes(x = highestZSample, y = pdf, color = Treatment, linetype = class)) + 
        geom_line()  + theme(legend.position = legendPosition, legend.box = legendDirection) + 
        labs(x = "Highest Z Sample") 
    })
  }
  
  return(list(
    ROCPlot = ROCPlot,
    AFROCPlot = AFROCPlot,
    wAFROCPlot = wAFROCPlot,
    FROCPlot = FROCPlot,
    PDFPlot = PDFPlot,
    aucROC = aucROC,
    aucAFROC = aucAFROC,
    aucwAFROC = aucwAFROC,
    aucFROC = aucFROC
  ))
}

is.wholenumber <- function(x)  round(x) == x