#' 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 -Inf.
#'
#' @param lesDistr Array: the probability mass function of the lesion
#' distribution for diseased cases. The default is 1. See \link{UtilLesDistr}.
#'
#' @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 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 modality and
#' ## mu = 3, lambda = 1.5, nu = 0.8, in the other modality. 20% of the diseased
#' ## cases have a single lesion, 40% have two lesions, 10% have 3 lesions,
#' ## and 30% have 4 lesions.
#'
#' res <- 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 = "bottom",
legendDirection = "horizontal",
legendJustification = c(0,1),
nlfRange = NULL,
llfRange = NULL,
nlfAlpha = NULL){
if (missing(zeta1)) zeta1 <- array(-Inf, 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 <- UtilLesWghtsLD(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_cpp, lambda[i]) # C++ function uses lambda
TPF <- sapply(zeta[[i]], yROC_cpp, 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_cpp(-20, lambda[i]) # this is wrong
maxFPF <- xROC_cpp(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_cpp(-20, mu[i], lambda, nu, lesDistr)
maxTPF <- yROC_cpp(zeta1[i], mu[i], lambda[i], nu[i], lesDistr)
# end bug fix
AUC <- integrate(y_ROC_FPF_cpp, 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 max(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]], RSM_wLLF_cpp, 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 <- RSM_wLLF_R(zeta1[i], mu[i], nu[i], lesDistr, relWeights)
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 # 5/3/23
}
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.