R/StCadVsRad.R
In dpc10ster/rjafroc: Artificial Intelligence Systems and Observer Performance
Defines functions
LrocOperatingPointsFromRatings
CadVsRadPlots
DualModalityRRRC
SingleModalityRRFC
StCadVsRad
Documented in
StCadVsRad
#' Significance testing: standalone CAD vs. radiologists
#'
#' @description Comparing standalone CAD vs. at least two radiologists interpreting
#' the same cases; \strong{standalone CAD} means that \bold{all} the
#' \bold{designer-level} mark-rating pairs generated by the CAD algorithm
#' are available to the analyst, not just the one or two marks per case
#' displayed to the radiologist (the latter are marks whose ratings exceed a
#' pre-selected threshold). At the very minimum, location-level information,
#' such as in
#' the LROC paradigm, should be used. Ideally, the FROC paradigm should be used.
#' A severe statistical power penalty is paid if one uses the ROC paradigm.
#' See Standalone CAD vs Radiologists chapter, available via \emph{download}
#' link at site
#' \url{https://github.com/dpc10ster/RJafrocBook/blob/gh-pages/RJafrocBook.pdf}
#'
#'
#' @param dataset \strong{The dataset to be analyzed; must be single-modality
#' at least three readers, where the first reader is CAD. }
#' @param FOM The desired FOM; for ROC data it must be \code{"Wilcoxon"}, for FROC data
#' it can be any valid FOM, e.g., \code{"HrAuc"}, \code{"wAFROC"}, etc;
#' for LROC data it must be \code{"Wilcoxon"}, or \code{"PCL"} or \code{"ALROC"}.
#' @param method The desired analysis: "1T-RRFC","1T-RRRC" (the default) or "2T-RRRC",
#' see manuscript for details.
#' @param alpha Significance level of the test, defaults to 0.05.
#' @param FPFValue Only needed for \code{LROC} data \strong{and} FOM = "PCL" or "ALROC";
#' where to evaluate a partial curve based figure of merit. The default is 0.2.
#' @param plots Flag, default is FALSE, i.e., a plot is not displayed.
#' If TRUE, it displays the appropriate operating characteristic for all
#' readers and CAD.
#'
#'
#' @details
#' \itemize{
#' \item{\strong{PCL} is the probability of a correct localization.}
#' \item{The LROC is the plot of PCL (ordinate) vs. FPF.}
#' \item{For LROC data, FOM = "PCL" means the interpolated PCL value
#' at the specified \code{FPFValue}.}
#' \item{For FOM = "ALROC" the trapezoidal area under the LROC
#' from FPF = 0 to FPF = \code{FPFValue} is used.}
#' \item{If \code{method = "1T-RRRC"} the first \strong{reader} is assumed to be CAD.}
#' \item{If \code{method = "2T-RRRC"} the first \strong{modality} is assumed to be CAD.}
#' \item{The NH is that the FOM of CAD equals the average of the readers.}
#' \item{The \code{method = "1T-RRRC"} analysis uses an adaptation of the
#' single-modality multiple-reader Obuchowski Rockette (OR) model described in a
#' paper by Hillis (2007), section 5.3. It is characterized by 3 parameters
#' \code{VarR}, \code{Var} and \code{Cov2}, where the latter two are estimated
#' using the jackknife.}
#' \item{For \code{method = "2T-RRRC"} the analysis replicates the CAD data as many times as
#' necessary so as to form one "modality" of an MRMC pairing, the other
#' "modality" being the radiologists. Then standard ORH analysis is applied. The
#' method is described in Kooi et al. It gives exactly the same final results
#' (F-statistic, ddf and p-value) as \code{"1T-RRRC"} but the intermediate quantities
#' are meaningless.}
#' }
#'
#' @return If \code{method = "1T-RRRC"} the return value is a
#' list with the following elements:
#' @return \item{fomCAD}{The observed FOM for CAD.}
#' @return \item{fomRAD}{The observed FOM array for the readers.}
#' @return \item{avgRadFom}{The average FOM of the readers.}
#' @return \item{avgDiffFom}{The mean of the difference FOM, RAD - CAD.}
#' @return \item{ciAvgDiffFom}{The 95-percent CI of the average difference, RAD - CAD.}
#' @return \item{varR}{The variance of the radiologists.}
#' @return \item{varError}{The variance of the error term in the single-modality
#' multiple-reader OR model.}
#' @return \item{cov2}{The covariance of the error term.}
#' @return \item{tstat}{The observed value of the t-statistic; it's square is
#' equivalent to an F-statistic.}
#' @return \item{df}{The degrees of freedom of the t-statistic.}
#' @return \item{pval}{The p-value for rejecting the NH.}
#' @return \item{Plots}{If argument plots = TRUE, a \pkg{ggplot} object
#' containing empirical operating characteristics
#' corresponding to specified FOM. For example, if \code{FOM} =
#' \code{"Wilcoxon"} an ROC plot object
#' is produced where reader 1 is CAD. If an LROC FOM is selected, an LROC
#' plot is displayed.}
#'
#'
#' @return If \code{method = "2T-RRRC"} the return value is a list
#' with the following elements:
#' @return \item{fomCAD}{The observed FOM for CAD.}
#' @return \item{fomRAD}{The observed FOM array for the readers.}
#' @return \item{avgRadFom}{The average FOM of the readers.}
#' @return \item{avgDiffFom}{The mean of the difference FOM, RAD - CAD.}
#' @return \item{ciDiffFom}{A data frame containing the statistics associated
#' with the average difference, RAD - CAD.}
#' @return \item{ciAvgRdrEachTrt}{A data frame containing the statistics
#' associated with the average FOM in each "modality".}
#' @return \item{varR}{The variance of the pure reader term in the OR model.}
#' @return \item{varTR}{The variance of the modality-reader term error
#' term in the OR model.}
#' @return \item{cov1}{The covariance1 of the error term - same reader,
#' different treatments.}
#' @return \item{cov2}{The covariance2 of the error term -
#' different readers, same modality.}
#' @return \item{cov3}{The covariance3 of the error term - different readers,
#' different treatments.}
#' @return \item{varError}{The variance of the pure error term in the OR model.}
#' @return \item{FStat}{The observed value of the F-statistic.}
#' @return \item{ndf}{The numerator degrees of freedom of the F-statistic.}
#' @return \item{df}{The denominator degrees of freedom of the F-statistic.}
#' @return \item{pval}{The p-value for rejecting the NH.}
#' @return \item{Plots}{see above.}
#'
#'
#'
#' @examples
#' ret1M <- StCadVsRad (dataset09,
#' FOM = "Wilcoxon", method = "1T-RRRC")
#'
#' StCadVsRad(datasetCadLroc,
#' FOM = "Wilcoxon", method = "1T-RRFC")
#'
#' retLroc1M <- StCadVsRad (datasetCadLroc,
#' FOM = "PCL", method = "1T-RRRC", FPFValue = 0.05)
#'
#' ## test with fewer readers
#' dataset09a <- DfExtractDataset(dataset09, rdrs = seq(1:7))
#' ret1M7 <- StCadVsRad (dataset09a,
#' FOM = "Wilcoxon", method = "1T-RRRC")
#'
#' datasetCadLroc7 <- DfExtractDataset(datasetCadLroc, rdrs = seq(1:7))
#' ret1MLroc7 <- StCadVsRad (datasetCadLroc7,
#' FOM = "PCL", method = "1T-RRRC", FPFValue = 0.05)
#'
#' \donttest{
#' ## takes longer than 5 sec on OSX
#' ## retLroc2M <- StCadVsRad (datasetCadLroc,
#' ## FOM = "PCL", method = "2T-RRRC", FPFValue = 0.05)
#'
#' ## ret2MLroc7 <- StCadVsRad (datasetCadLroc7,
#' ## FOM = "PCL", method = "2T-RRRC", FPFValue = 0.05)
#' }
#'
#' @references
#' Hillis SL (2007) A comparison of denominator degrees of freedom methods
#' for multiple observer ROC studies, Statistics in Medicine. 26:596-619.
#'
#' 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}
#'
#' Hupse R, Samulski M, Lobbes M, et al (2013) Standalone computer-aided detection compared to radiologists
#' performance for the detection of mammographic masses, Eur Radiol. 23(1):93-100.
#'
#' Kooi T, Gubern-Merida A, et al. (2016) A comparison between a deep convolutional
#' neural network and radiologists for classifying regions of interest in mammography.
#' Paper presented at: International Workshop on Digital Mammography, Malmo, Sweden.
#'
#'
#' @import ggplot2
#' @importFrom stats var t.test cov
#' @export
StCadVsRad <- function(dataset, FOM, FPFValue = 0.2, method = "1T-RRRC",
alpha = 0.05, plots = FALSE)
{
options(stringsAsFactors = FALSE)
if (length(dataset$ratings$NL[,1,1,1]) != 1) stop("dataset has to be single-modality and at least three readers with CAD as the first reader")
if ((dataset$descriptions$type == "ROC") && (FOM %in% c("PCL", "ALROC"))) stop("Cannot use LROC FOM with ROC data")
if (length(dataset$ratings$NL[1,,1,1]) <= 2) stop("dataset has to have at least two radiologist readers")
if (method == "1T-RRFC") {
ret <- SingleModalityRRFC(dataset, FOM, FPFValue, alpha)
} else if (method == "1T-RRRC") {
ret <- SingleModalityRRRC(dataset, FOM, FPFValue, alpha)
} else if (method == "2T-RRRC") {
ret <- DualModalityRRRC (dataset, FOM, FPFValue, alpha)
} else stop("incorrect method specified")
# if ((dataset$descriptions$type != "LROC") && (method == "2T-RRRC")) {
# stop("2T-RRRC for non LROC data (not implemented) is unnecessary as 1T-method should be used instead.")
# }
if (plots) {
genericPlot <- CadVsRadPlots (dataset, FOM)
retNames <- names(ret)
retNames <- c(retNames, "Plots")
len <- length(ret)
ret1 <- vector("list", len+1)
for (i in 1:len){
ret1[[i]] <- ret[[i]]
}
ret1[[len+1]] <- genericPlot
names(ret1) <- retNames
} else ret1 <- ret
return(ret1)
}
# Handles all dataTypes
SingleModalityRRFC <- function(dataset, FOM, FPFValue, alpha) {
# `as.matrix` is absolutely necessary if following `mean()` function is to work
thetajc <- UtilFigureOfMerit(dataset, FOM, FPFValue)
Psijc <- thetajc[-1] - thetajc[1]
ret <- t.test(Psijc, conf.level = 1-alpha)
Tstat <- as.numeric(ret$statistic)
df <- as.numeric(ret$parameter)
pval <- ret$p.value
CIAvgRadFom <- as.numeric(ret$conf.int)+thetajc[1]
avgDiffFom <- mean(Psijc)
CIAvgDiffFom <- as.numeric(ret$conf.int)
return (list (
fomCAD = thetajc[1],
fomRAD = thetajc[-1],
avgRadFom = mean(thetajc[-1]),
CIAvgRadFom = CIAvgRadFom,
avgDiffFom = avgDiffFom,
CIAvgDiffFom = CIAvgDiffFom,
varR = var(thetajc[-1]),
Tstat = Tstat,
df = df,
pval = pval
))
}
# Formerly
# Anal2007Hillis53
# Anal2007Hillis53
# Anal2007Hillis53
# Handles all datasets
SingleModalityRRRC <- function (dataset, FOM, FPFValue, alpha)
{
ret <- DiffFomVarCov2(dataset, FOM, FPFValue) # VarCov2 (subtract first reader FOMs before getting covariance)
varError <- ret$var; Cov2 <- ret$cov2
J <- length(dataset$ratings$NL[1,,1,1]) - 1 # number of radiologists minus CAD reader
# `as.matrix` is absolutely necessary if following `mean()` function is to work
thetajc <- UtilFigureOfMerit(dataset, FOM, FPFValue)
Psijc <- thetajc[1,2:(J+1)] - thetajc[1,1] # subract CAD from RAD, my Eqn. 13
MSR <- 0 # 1st un-numbered equation on page 607
avgDiffFom <- mean(Psijc)
for (j in 1:J){
MSR <- MSR + (Psijc[j] - avgDiffFom)^2
}
MSR <- MSR / (J - 1)
# Compared to equations in 2013 Hillis paper, in paragraph following Table I
# OK; 10/14/19
MSdenOR_single <- MSR + max(J * Cov2, 0) # # 2nd un-numbered equation on page 607
DdfhSingle <- MSdenOR_single^2 / (MSR^2 / (J - 1)) # 3rd un-numbered equation on page 607
TstatStar <- avgDiffFom / sqrt(MSdenOR_single/J) # in-text equation on line 1 on page 607
# BUT with theta0 = 0
pval <- 1 - pt(abs(TstatStar),DdfhSingle) + pt(-abs(TstatStar),DdfhSingle) # two tailed probability
CIAvgDiffFom <- array(dim=2)
CIAvgDiffFom[1] <- qt(alpha/2,df = DdfhSingle) # Equation 25 on page 607
CIAvgDiffFom[2] <- qt(1-alpha/2,df = DdfhSingle)
CIAvgDiffFom <- CIAvgDiffFom * sqrt(MSdenOR_single/J)
CIAvgDiffFom <- CIAvgDiffFom + avgDiffFom
CIAvgRad <- mean(thetajc[2:(J+1)]) + CIAvgDiffFom - mean(Psijc)
return (list (
fomCAD = thetajc[1],
fomRAD = thetajc[-1],
avgRadFom = mean(thetajc[-1]),
CIAvgRad = CIAvgRad,
avgDiffFom = avgDiffFom,
CIAvgDiffFom = CIAvgDiffFom,
varR = var(as.numeric(thetajc[-1])),
varError = varError,
cov2 = Cov2,
Tstat = TstatStar,
df = DdfhSingle,
pval = pval
))
}
DualModalityRRRC <- function(dataset, FOM, FPFValue, alpha)
{
K <- length(dataset$ratings$NL[1,1,,1])
type <- dataset$descriptions$type
if ((type == "LROC") && (FOM %in% c("PCL", "ALROC")))
{
ret1 <- dataset2ratings(dataset, FOM)
TP <- ret1$zjk2
zjk2Il <- ret1$zjk2Il
K2 <- length(TP[1,])
K1 <- K - K2
FP <- ret1$zjk1[,1:K1]
J <- length(FP[,1]) - 1
combinedNL <- array(-Inf, dim=c(2,J,K,1))
for (j in 1:J){
combinedNL[1,j,1:K1,1] <- FP[1,]
}
combinedNL[2,,1:K1,1] <- FP[2:(J+1),]
combinedLL <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLL[1,j,,1] <- TP[1,]
}
combinedLL[2,,,1] <- TP[2:(J+1),]
combinedLLIl <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLLIl[1,j,,1] <- zjk2Il[1,]
}
combinedLLIl[2,,,1] <- zjk2Il[2:(J+1),]
} else if ((type == "LROC") && (FOM == "Wilcoxon")) {
datasetRoc <- DfLroc2Roc(dataset)
type <- datasetRoc$descriptions$type
NL <- datasetRoc$ratings$NL
LL <- datasetRoc$ratings$LL
K <- length(NL[1,1,,1])
K2 <- length(LL[1,1,,1])
K1 <- K - K2
FP <- NL[1,,1:K1,1]
TP <- LL[1,,,1]
J <- length(FP[,1]) - 1
combinedNL <- array(-Inf, dim=c(2,J,K,1))
for (j in 1:J){
combinedNL[1,j,1:K1,1] <- FP[1,]
}
combinedNL[2,,1:K1,1] <- FP[2:(J+1),]
combinedLL <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLL[1,j,,1] <- TP[1,]
}
combinedLL[2,,,1] <- TP[2:(J+1),]
} else if ((type == "ROC") && (FOM == "Wilcoxon")) {
NL <- dataset$ratings$NL
LL <- dataset$ratings$LL
K <- length(NL[1,1,,1])
K2 <- length(LL[1,1,,1])
K1 <- K - K2
FP <- NL[1,,1:K1,1]
TP <- LL[1,,,1]
J <- length(FP[,1]) - 1
combinedNL <- array(-Inf, dim=c(2,J,K,1))
for (j in 1:J){
combinedNL[1,j,1:K1,1] <- FP[1,]
}
combinedNL[2,,1:K1,1] <- FP[2:(J+1),]
combinedLL <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLL[1,j,,1] <- TP[1,]
}
combinedLL[2,,,1] <- TP[2:(J+1),]
} else stop("Incorrect FOM with LROC data")
IDs <- rep(1, K2)
dim(IDs) <- c(K2,1)
weights <- rep(0, K2)
dim(weights) <- c(K2,1)
if (type == "LROC") {
NL <- combinedNL
LL = combinedLL
LL_IL = combinedLLIl
fileName <- paste0("DualModalityRRRC(", dataset$descriptions$fileName, ")")
name <- NA
design <- "FCTRL"
truthTableStr <- array(dim = c(2,J,K,2))
truthTableStr[,,1:K1,1] <- 1
truthTableStr[,,(K1+1):K,2] <- 1
type <- "LROC"
perCase <- rep(1,K2)
IDs <- dataset$lesions$IDs
weights <- dataset$lesions$weights
modalityID = as.character(c(1,2))
readerID = as.character(1:J)
datasetCombined <- convert2dataset(NL, LL, LL_IL,
perCase, IDs, weights,
fileName, type, name, truthTableStr, design,
modalityID, readerID)
} else {
NL <- combinedNL
LL = combinedLL
LL_IL = NA
fileName <- paste0("DualModalityRRRC(", dataset$descriptions$fileName, ")")
name <- NA
design <- "FCTRL"
truthTableStr <- array(dim = c(2,J,K,2))
truthTableStr[,,1:K1,1] <- 1
truthTableStr[,,(K1+1):K,2] <- 1
type <- "ROC"
perCase <- rep(1,K2)
IDs <- dataset$lesions$IDs
weights <- dataset$lesions$weights
modalityID = as.character(c(1,2))
readerID = as.character(1:J)
datasetCombined <- convert2dataset(NL, LL, LL_IL,
perCase, IDs, weights,
fileName, type, name, truthTableStr, design,
modalityID, readerID)
}
stats1 <- St(datasetCombined, FOM = FOM, method = "OR", alpha = alpha, analysisOption = "RRRC", FPFValue = FPFValue)
thetajc <- stats1$FOMs$foms
thetajc <- as.matrix(thetajc)
fomCAD <- thetajc[1,1]
fomRAD <- thetajc[2,]
avgRadFom <- mean(fomRAD)
varDiffFom <- var(fomRAD)
avgDiffFom <- avgRadFom - fomCAD
FStat <- stats1$RRRC$FTests$FStat[1]
ddf <- stats1$RRRC$FTests$DF[2]
ndf <- stats1$RRRC$FTests$DF[1]
pval <- stats1$RRRC$FTests$p[1]
varR <- stats1$ANOVA$VarCom["VarR", "Estimates"] # corrected 11/26/2020
varTR <- stats1$ANOVA$VarCom["VarTR", "Estimates"] # do:
varError <- stats1$ANOVA$VarCom["Var", "Estimates"] # do:
cov1 <- stats1$ANOVA$VarCom["Cov1", "Estimates"] # do:
cov2 <- stats1$ANOVA$VarCom["Cov2", "Estimates"] # do:
cov3 <- stats1$ANOVA$VarCom["Cov3", "Estimates"] # do:
# corrected 11/26/2020; added minus sign in following ...
ciDiffFom <- stats1$RRRC$ciDiffTrt
# correction needed
# as St returns CAD - RAD (trt1 - trt2)
# not RAD - CAD
t <- ciDiffFom
rowNames <- rownames(t)
x1 <- strsplit(rowNames, split = "-")[[1]]
rowNames <- paste0(x1[2], "-", x1[1])
rownames(t) <- rowNames
t$Estimate <- -ciDiffFom$Estimate
t$t <- -ciDiffFom$t
t$CILower <- -ciDiffFom$CIUpper
t$CIUpper <- -ciDiffFom$CILower
ciDiffFom <- t
# end corrections
ciAvgRdrEachTrt <- stats1$RRRC$ciAvgRdrEachTrt
return (list (
fomCAD = fomCAD,
fomRAD = fomRAD,
avgRadFom = avgRadFom,
avgDiffFom = avgDiffFom,
varDiffFom = varDiffFom,
ciDiffFom = ciDiffFom,
ciAvgRdrEachTrt = ciAvgRdrEachTrt,
varR = varR,
varTR = varTR,
cov1 = cov1,
cov2 = cov2,
cov3 = cov3,
varError = varError,
FStat = FStat,
ndf = ndf,
df = ddf, # for consistency with 1T-RRFC
pval = pval
))
}
CadVsRadPlots <- function(dataset, FOM) {
ret1 <- dataset2ratings(dataset, FOM)
zjk1 <- ret1$zjk1
zjk2 <- ret1$zjk2
type <- dataset$descriptions$type
if (type == "ROC") {
# fixed one of hard coding errors noticed by Alejandro
genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "ROC")$Plot
} else if ((type == "LROC") && ((FOM == "PCL") || (FOM == "ALROC"))) {
if (type == "LROC") genericPlot <- LrocPlots (zjk1, zjk2, seq(1,length(zjk1[,1])-1))$lrocPlot
} else if ((type == "LROC") && (FOM == "Wilcoxon")) {
if (type == "LROC") {
datasetRoc <- DfLroc2Roc(dataset)
genericPlot <- PlotEmpiricalOperatingCharacteristics(datasetRoc, rdrs = 1:length(zjk1[,1]), opChType = "ROC")$Plot
}
} else if ((type == "FROC") && (FOM %in% c("HrAuc", "AFROC", "wAFROC"))) {
if (FOM == "HrAuc") genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "ROC")$Plot
if (FOM == "AFROC") genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "AFROC")$Plot
if (FOM == "wAFROC") genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "wAFROC")$Plot
}else stop("data type has to be ROC, FROC or LROC")
return(genericPlot)
}
##
## this uses the jackknife to estimate var & cov2
## Handles all datasets
DiffFomVarCov2 <- function (dataset, FOM, FPFValue) # for difference FOM, radiologist minus CAD
{
#if (dataset$descriptions$type == "LROC") stop("Dataset must NOT be LROC")
J <- length(dataset$descriptions$readerID)
K <- length(dataset$ratings$NL[1,1,,1])
dsCad <- DfExtractDataset(dataset, trts = 1, rdrs = 1)
dsRad <- DfExtractDataset(dataset, trts = 1, rdrs = c(2:J))
jkFomValuesCad <- UtilPseudoValues(dsCad, FOM, FPFValue)$jkFomValues
jkFomValuesRad <- UtilPseudoValues(dsRad, FOM, FPFValue)$jkFomValues
jkDiffFomValues <- array(dim = c(J-1,K))
# this does not work wasted much time
# this does not work wasted much time
# jkDiffFomValues <- jkFomValuesRad[1,,] - jkFomValuesCad[1,1,] # this does not work
# this does not work wasted much time
# this does not work wasted much time
for (j in 1:(J-1)) jkDiffFomValues[j,] <- jkFomValuesRad[1,j,] - jkFomValuesCad[1,1,]
J <- J - 1 # number of human readers
Covariance <- array(dim = c(J, J))
for (j in 1:J){
for (jp in 1:J){
Covariance[j, jp] <- cov(jkDiffFomValues[j, ], jkDiffFomValues[jp, ])
}
}
varError <- 0;count <- 0
for (j in 1:J){
varError <- varError + Covariance[j, j]
count <- count + 1
}
varError <- varError / count
varError <- varError *(K-1)^2/K
Cov2 <- 0;count <- 0
for (j in 1:J){
for (jp in 1:J){
if (jp != j){
Cov2 <- Cov2 + Covariance[j, jp]
count <- count + 1
}
}
}
Cov2 <- Cov2 / count
Cov2 <- Cov2 *(K-1)^2/K
return (list (
var = varError,
cov2 = Cov2
))
}
LrocOperatingPointsFromRatings <- function( zk1, zk2Cl )
{
if (FALSE) { # this is to check the method, see ChkLrocFoms.xlsx
zk1 <- c(rep(0, 4), 1.424443, rep(-50,5))
zk2Cl <- c(0.5542791, 1.2046176, -50, 3.4596787, 2.732657)
K2 <- length(zk2Cl)
K1 <- length(zk1) - K2
zk1 <- zk1[zk1 != -50]
zk2Cl <- zk2Cl[zk2Cl != -50]
} else {
K2 <- length(zk2Cl)
K1 <- length(zk1) - K2
zk1 <- zk1[is.finite(zk1)]
zk2Cl <- zk2Cl[is.finite(zk2Cl)]
}
FPF <- 1
PCL <- NULL
zk1Temp <- zk1;zk2ClTemp <- zk2Cl
while(1) {
cutoff <- min( c( zk1Temp, zk2ClTemp ) )
zk1Temp <- zk1[ zk1 > cutoff ]
zk2ClTemp <- zk2Cl[ zk2Cl > cutoff ]
FPF1 <- length( zk1Temp ) / K1
PCL1 <- length( zk2ClTemp ) / K2
FPF <- c( FPF, FPF1 )
if (length(PCL) == 0) PCL <- c(PCL1, PCL1) else PCL <- c( PCL, PCL1 )
if( FPF1 == 0 && PCL1 == 0 ) {
break
}
}
FPF <- rev(FPF)
PCL <- rev(PCL)
return( list(
FPF = FPF,
PCL = PCL
) )
}
LrocPlots <- function (zjk1, zjk2, doJ)
{
J <- length(zjk1[,1])
j <- 1;zjk1Temp <- zjk1[j,];zk2Temp <- zjk2[j,]
lroc <- LrocOperatingPointsFromRatings( zjk1Temp, zk2Temp )
FPF <- lroc$FPF;PCL <- lroc$PCL
lrocPlotData <- data.frame(FPF = FPF, PCL = PCL, reader = "R-CAD")
for (j in 2:J) {
if ((j - 1) %in% doJ) {
zjk1Temp <- zjk1[j,]
zk2Temp <- zjk2[j,]
lroc <- LrocOperatingPointsFromRatings( zjk1Temp, zk2Temp )
FPF <- lroc$FPF
PCL <- lroc$PCL
reader = paste0("R-", as.character(j - 1))
lrocPlotData <- rbind(lrocPlotData, data.frame(FPF = FPF, PCL = PCL, reader = reader))
}
}
lrocPlot <- ggplot(data = lrocPlotData, aes(x = FPF, y = PCL, color = reader)) + geom_line()
g <- ggplot_build(lrocPlot)
colors <- as.character(unique(g$data[[1]]$colour))
colors[1] <- "#000000"
sizes <- c(2, rep(1, length(doJ)))
lrocPlot <- ggplot(data = lrocPlotData, aes(x = FPF, y = PCL, color = reader)) + geom_line(aes(linewidth = reader)) +
scale_color_manual(values = colors) + scale_size_manual(values = sizes) +
theme(legend.title = element_blank(), legend.position = legendPosition <- c(1, 0), legend.justification = c(1, 0))
return(list(
lrocPlot = lrocPlot,
afrocPlot = NULL)
)
}
dpc10ster/rjafroc documentation built on Jan. 18, 2024, 4:37 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions LrocOperatingPointsFromRatings CadVsRadPlots DualModalityRRRC SingleModalityRRFC StCadVsRad
Documented in StCadVsRad
#' Significance testing: standalone CAD vs. radiologists
#'
#' @description Comparing standalone CAD vs. at least two radiologists interpreting
#' the same cases; \strong{standalone CAD} means that \bold{all} the
#' \bold{designer-level} mark-rating pairs generated by the CAD algorithm
#' are available to the analyst, not just the one or two marks per case
#' displayed to the radiologist (the latter are marks whose ratings exceed a
#' pre-selected threshold). At the very minimum, location-level information,
#' such as in
#' the LROC paradigm, should be used. Ideally, the FROC paradigm should be used.
#' A severe statistical power penalty is paid if one uses the ROC paradigm.
#' See Standalone CAD vs Radiologists chapter, available via \emph{download}
#' link at site
#' \url{https://github.com/dpc10ster/RJafrocBook/blob/gh-pages/RJafrocBook.pdf}
#'
#'
#' @param dataset \strong{The dataset to be analyzed; must be single-modality
#' at least three readers, where the first reader is CAD. }
#' @param FOM The desired FOM; for ROC data it must be \code{"Wilcoxon"}, for FROC data
#' it can be any valid FOM, e.g., \code{"HrAuc"}, \code{"wAFROC"}, etc;
#' for LROC data it must be \code{"Wilcoxon"}, or \code{"PCL"} or \code{"ALROC"}.
#' @param method The desired analysis: "1T-RRFC","1T-RRRC" (the default) or "2T-RRRC",
#' see manuscript for details.
#' @param alpha Significance level of the test, defaults to 0.05.
#' @param FPFValue Only needed for \code{LROC} data \strong{and} FOM = "PCL" or "ALROC";
#' where to evaluate a partial curve based figure of merit. The default is 0.2.
#' @param plots Flag, default is FALSE, i.e., a plot is not displayed.
#' If TRUE, it displays the appropriate operating characteristic for all
#' readers and CAD.
#'
#'
#' @details
#' \itemize{
#' \item{\strong{PCL} is the probability of a correct localization.}
#' \item{The LROC is the plot of PCL (ordinate) vs. FPF.}
#' \item{For LROC data, FOM = "PCL" means the interpolated PCL value
#' at the specified \code{FPFValue}.}
#' \item{For FOM = "ALROC" the trapezoidal area under the LROC
#' from FPF = 0 to FPF = \code{FPFValue} is used.}
#' \item{If \code{method = "1T-RRRC"} the first \strong{reader} is assumed to be CAD.}
#' \item{If \code{method = "2T-RRRC"} the first \strong{modality} is assumed to be CAD.}
#' \item{The NH is that the FOM of CAD equals the average of the readers.}
#' \item{The \code{method = "1T-RRRC"} analysis uses an adaptation of the
#' single-modality multiple-reader Obuchowski Rockette (OR) model described in a
#' paper by Hillis (2007), section 5.3. It is characterized by 3 parameters
#' \code{VarR}, \code{Var} and \code{Cov2}, where the latter two are estimated
#' using the jackknife.}
#' \item{For \code{method = "2T-RRRC"} the analysis replicates the CAD data as many times as
#' necessary so as to form one "modality" of an MRMC pairing, the other
#' "modality" being the radiologists. Then standard ORH analysis is applied. The
#' method is described in Kooi et al. It gives exactly the same final results
#' (F-statistic, ddf and p-value) as \code{"1T-RRRC"} but the intermediate quantities
#' are meaningless.}
#' }
#'
#' @return If \code{method = "1T-RRRC"} the return value is a
#' list with the following elements:
#' @return \item{fomCAD}{The observed FOM for CAD.}
#' @return \item{fomRAD}{The observed FOM array for the readers.}
#' @return \item{avgRadFom}{The average FOM of the readers.}
#' @return \item{avgDiffFom}{The mean of the difference FOM, RAD - CAD.}
#' @return \item{ciAvgDiffFom}{The 95-percent CI of the average difference, RAD - CAD.}
#' @return \item{varR}{The variance of the radiologists.}
#' @return \item{varError}{The variance of the error term in the single-modality
#' multiple-reader OR model.}
#' @return \item{cov2}{The covariance of the error term.}
#' @return \item{tstat}{The observed value of the t-statistic; it's square is
#' equivalent to an F-statistic.}
#' @return \item{df}{The degrees of freedom of the t-statistic.}
#' @return \item{pval}{The p-value for rejecting the NH.}
#' @return \item{Plots}{If argument plots = TRUE, a \pkg{ggplot} object
#' containing empirical operating characteristics
#' corresponding to specified FOM. For example, if \code{FOM} =
#' \code{"Wilcoxon"} an ROC plot object
#' is produced where reader 1 is CAD. If an LROC FOM is selected, an LROC
#' plot is displayed.}
#'
#'
#' @return If \code{method = "2T-RRRC"} the return value is a list
#' with the following elements:
#' @return \item{fomCAD}{The observed FOM for CAD.}
#' @return \item{fomRAD}{The observed FOM array for the readers.}
#' @return \item{avgRadFom}{The average FOM of the readers.}
#' @return \item{avgDiffFom}{The mean of the difference FOM, RAD - CAD.}
#' @return \item{ciDiffFom}{A data frame containing the statistics associated
#' with the average difference, RAD - CAD.}
#' @return \item{ciAvgRdrEachTrt}{A data frame containing the statistics
#' associated with the average FOM in each "modality".}
#' @return \item{varR}{The variance of the pure reader term in the OR model.}
#' @return \item{varTR}{The variance of the modality-reader term error
#' term in the OR model.}
#' @return \item{cov1}{The covariance1 of the error term - same reader,
#' different treatments.}
#' @return \item{cov2}{The covariance2 of the error term -
#' different readers, same modality.}
#' @return \item{cov3}{The covariance3 of the error term - different readers,
#' different treatments.}
#' @return \item{varError}{The variance of the pure error term in the OR model.}
#' @return \item{FStat}{The observed value of the F-statistic.}
#' @return \item{ndf}{The numerator degrees of freedom of the F-statistic.}
#' @return \item{df}{The denominator degrees of freedom of the F-statistic.}
#' @return \item{pval}{The p-value for rejecting the NH.}
#' @return \item{Plots}{see above.}
#'
#'
#'
#' @examples
#' ret1M <- StCadVsRad (dataset09,
#' FOM = "Wilcoxon", method = "1T-RRRC")
#'
#' StCadVsRad(datasetCadLroc,
#' FOM = "Wilcoxon", method = "1T-RRFC")
#'
#' retLroc1M <- StCadVsRad (datasetCadLroc,
#' FOM = "PCL", method = "1T-RRRC", FPFValue = 0.05)
#'
#' ## test with fewer readers
#' dataset09a <- DfExtractDataset(dataset09, rdrs = seq(1:7))
#' ret1M7 <- StCadVsRad (dataset09a,
#' FOM = "Wilcoxon", method = "1T-RRRC")
#'
#' datasetCadLroc7 <- DfExtractDataset(datasetCadLroc, rdrs = seq(1:7))
#' ret1MLroc7 <- StCadVsRad (datasetCadLroc7,
#' FOM = "PCL", method = "1T-RRRC", FPFValue = 0.05)
#'
#' \donttest{
#' ## takes longer than 5 sec on OSX
#' ## retLroc2M <- StCadVsRad (datasetCadLroc,
#' ## FOM = "PCL", method = "2T-RRRC", FPFValue = 0.05)
#'
#' ## ret2MLroc7 <- StCadVsRad (datasetCadLroc7,
#' ## FOM = "PCL", method = "2T-RRRC", FPFValue = 0.05)
#' }
#'
#' @references
#' Hillis SL (2007) A comparison of denominator degrees of freedom methods
#' for multiple observer ROC studies, Statistics in Medicine. 26:596-619.
#'
#' 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}
#'
#' Hupse R, Samulski M, Lobbes M, et al (2013) Standalone computer-aided detection compared to radiologists
#' performance for the detection of mammographic masses, Eur Radiol. 23(1):93-100.
#'
#' Kooi T, Gubern-Merida A, et al. (2016) A comparison between a deep convolutional
#' neural network and radiologists for classifying regions of interest in mammography.
#' Paper presented at: International Workshop on Digital Mammography, Malmo, Sweden.
#'
#'
#' @import ggplot2
#' @importFrom stats var t.test cov
#' @export
StCadVsRad <- function(dataset, FOM, FPFValue = 0.2, method = "1T-RRRC",
alpha = 0.05, plots = FALSE)
{
options(stringsAsFactors = FALSE)
if (length(dataset$ratings$NL[,1,1,1]) != 1) stop("dataset has to be single-modality and at least three readers with CAD as the first reader")
if ((dataset$descriptions$type == "ROC") && (FOM %in% c("PCL", "ALROC"))) stop("Cannot use LROC FOM with ROC data")
if (length(dataset$ratings$NL[1,,1,1]) <= 2) stop("dataset has to have at least two radiologist readers")
if (method == "1T-RRFC") {
ret <- SingleModalityRRFC(dataset, FOM, FPFValue, alpha)
} else if (method == "1T-RRRC") {
ret <- SingleModalityRRRC(dataset, FOM, FPFValue, alpha)
} else if (method == "2T-RRRC") {
ret <- DualModalityRRRC (dataset, FOM, FPFValue, alpha)
} else stop("incorrect method specified")
# if ((dataset$descriptions$type != "LROC") && (method == "2T-RRRC")) {
# stop("2T-RRRC for non LROC data (not implemented) is unnecessary as 1T-method should be used instead.")
# }
if (plots) {
genericPlot <- CadVsRadPlots (dataset, FOM)
retNames <- names(ret)
retNames <- c(retNames, "Plots")
len <- length(ret)
ret1 <- vector("list", len+1)
for (i in 1:len){
ret1[[i]] <- ret[[i]]
}
ret1[[len+1]] <- genericPlot
names(ret1) <- retNames
} else ret1 <- ret
return(ret1)
}
# Handles all dataTypes
SingleModalityRRFC <- function(dataset, FOM, FPFValue, alpha) {
# `as.matrix` is absolutely necessary if following `mean()` function is to work
thetajc <- UtilFigureOfMerit(dataset, FOM, FPFValue)
Psijc <- thetajc[-1] - thetajc[1]
ret <- t.test(Psijc, conf.level = 1-alpha)
Tstat <- as.numeric(ret$statistic)
df <- as.numeric(ret$parameter)
pval <- ret$p.value
CIAvgRadFom <- as.numeric(ret$conf.int)+thetajc[1]
avgDiffFom <- mean(Psijc)
CIAvgDiffFom <- as.numeric(ret$conf.int)
return (list (
fomCAD = thetajc[1],
fomRAD = thetajc[-1],
avgRadFom = mean(thetajc[-1]),
CIAvgRadFom = CIAvgRadFom,
avgDiffFom = avgDiffFom,
CIAvgDiffFom = CIAvgDiffFom,
varR = var(thetajc[-1]),
Tstat = Tstat,
df = df,
pval = pval
))
}
# Formerly
# Anal2007Hillis53
# Anal2007Hillis53
# Anal2007Hillis53
# Handles all datasets
SingleModalityRRRC <- function (dataset, FOM, FPFValue, alpha)
{
ret <- DiffFomVarCov2(dataset, FOM, FPFValue) # VarCov2 (subtract first reader FOMs before getting covariance)
varError <- ret$var; Cov2 <- ret$cov2
J <- length(dataset$ratings$NL[1,,1,1]) - 1 # number of radiologists minus CAD reader
# `as.matrix` is absolutely necessary if following `mean()` function is to work
thetajc <- UtilFigureOfMerit(dataset, FOM, FPFValue)
Psijc <- thetajc[1,2:(J+1)] - thetajc[1,1] # subract CAD from RAD, my Eqn. 13
MSR <- 0 # 1st un-numbered equation on page 607
avgDiffFom <- mean(Psijc)
for (j in 1:J){
MSR <- MSR + (Psijc[j] - avgDiffFom)^2
}
MSR <- MSR / (J - 1)
# Compared to equations in 2013 Hillis paper, in paragraph following Table I
# OK; 10/14/19
MSdenOR_single <- MSR + max(J * Cov2, 0) # # 2nd un-numbered equation on page 607
DdfhSingle <- MSdenOR_single^2 / (MSR^2 / (J - 1)) # 3rd un-numbered equation on page 607
TstatStar <- avgDiffFom / sqrt(MSdenOR_single/J) # in-text equation on line 1 on page 607
# BUT with theta0 = 0
pval <- 1 - pt(abs(TstatStar),DdfhSingle) + pt(-abs(TstatStar),DdfhSingle) # two tailed probability
CIAvgDiffFom <- array(dim=2)
CIAvgDiffFom[1] <- qt(alpha/2,df = DdfhSingle) # Equation 25 on page 607
CIAvgDiffFom[2] <- qt(1-alpha/2,df = DdfhSingle)
CIAvgDiffFom <- CIAvgDiffFom * sqrt(MSdenOR_single/J)
CIAvgDiffFom <- CIAvgDiffFom + avgDiffFom
CIAvgRad <- mean(thetajc[2:(J+1)]) + CIAvgDiffFom - mean(Psijc)
return (list (
fomCAD = thetajc[1],
fomRAD = thetajc[-1],
avgRadFom = mean(thetajc[-1]),
CIAvgRad = CIAvgRad,
avgDiffFom = avgDiffFom,
CIAvgDiffFom = CIAvgDiffFom,
varR = var(as.numeric(thetajc[-1])),
varError = varError,
cov2 = Cov2,
Tstat = TstatStar,
df = DdfhSingle,
pval = pval
))
}
DualModalityRRRC <- function(dataset, FOM, FPFValue, alpha)
{
K <- length(dataset$ratings$NL[1,1,,1])
type <- dataset$descriptions$type
if ((type == "LROC") && (FOM %in% c("PCL", "ALROC")))
{
ret1 <- dataset2ratings(dataset, FOM)
TP <- ret1$zjk2
zjk2Il <- ret1$zjk2Il
K2 <- length(TP[1,])
K1 <- K - K2
FP <- ret1$zjk1[,1:K1]
J <- length(FP[,1]) - 1
combinedNL <- array(-Inf, dim=c(2,J,K,1))
for (j in 1:J){
combinedNL[1,j,1:K1,1] <- FP[1,]
}
combinedNL[2,,1:K1,1] <- FP[2:(J+1),]
combinedLL <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLL[1,j,,1] <- TP[1,]
}
combinedLL[2,,,1] <- TP[2:(J+1),]
combinedLLIl <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLLIl[1,j,,1] <- zjk2Il[1,]
}
combinedLLIl[2,,,1] <- zjk2Il[2:(J+1),]
} else if ((type == "LROC") && (FOM == "Wilcoxon")) {
datasetRoc <- DfLroc2Roc(dataset)
type <- datasetRoc$descriptions$type
NL <- datasetRoc$ratings$NL
LL <- datasetRoc$ratings$LL
K <- length(NL[1,1,,1])
K2 <- length(LL[1,1,,1])
K1 <- K - K2
FP <- NL[1,,1:K1,1]
TP <- LL[1,,,1]
J <- length(FP[,1]) - 1
combinedNL <- array(-Inf, dim=c(2,J,K,1))
for (j in 1:J){
combinedNL[1,j,1:K1,1] <- FP[1,]
}
combinedNL[2,,1:K1,1] <- FP[2:(J+1),]
combinedLL <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLL[1,j,,1] <- TP[1,]
}
combinedLL[2,,,1] <- TP[2:(J+1),]
} else if ((type == "ROC") && (FOM == "Wilcoxon")) {
NL <- dataset$ratings$NL
LL <- dataset$ratings$LL
K <- length(NL[1,1,,1])
K2 <- length(LL[1,1,,1])
K1 <- K - K2
FP <- NL[1,,1:K1,1]
TP <- LL[1,,,1]
J <- length(FP[,1]) - 1
combinedNL <- array(-Inf, dim=c(2,J,K,1))
for (j in 1:J){
combinedNL[1,j,1:K1,1] <- FP[1,]
}
combinedNL[2,,1:K1,1] <- FP[2:(J+1),]
combinedLL <- array(-Inf, dim=c(2,J,K2,1))
for (j in 1:J){
combinedLL[1,j,,1] <- TP[1,]
}
combinedLL[2,,,1] <- TP[2:(J+1),]
} else stop("Incorrect FOM with LROC data")
IDs <- rep(1, K2)
dim(IDs) <- c(K2,1)
weights <- rep(0, K2)
dim(weights) <- c(K2,1)
if (type == "LROC") {
NL <- combinedNL
LL = combinedLL
LL_IL = combinedLLIl
fileName <- paste0("DualModalityRRRC(", dataset$descriptions$fileName, ")")
name <- NA
design <- "FCTRL"
truthTableStr <- array(dim = c(2,J,K,2))
truthTableStr[,,1:K1,1] <- 1
truthTableStr[,,(K1+1):K,2] <- 1
type <- "LROC"
perCase <- rep(1,K2)
IDs <- dataset$lesions$IDs
weights <- dataset$lesions$weights
modalityID = as.character(c(1,2))
readerID = as.character(1:J)
datasetCombined <- convert2dataset(NL, LL, LL_IL,
perCase, IDs, weights,
fileName, type, name, truthTableStr, design,
modalityID, readerID)
} else {
NL <- combinedNL
LL = combinedLL
LL_IL = NA
fileName <- paste0("DualModalityRRRC(", dataset$descriptions$fileName, ")")
name <- NA
design <- "FCTRL"
truthTableStr <- array(dim = c(2,J,K,2))
truthTableStr[,,1:K1,1] <- 1
truthTableStr[,,(K1+1):K,2] <- 1
type <- "ROC"
perCase <- rep(1,K2)
IDs <- dataset$lesions$IDs
weights <- dataset$lesions$weights
modalityID = as.character(c(1,2))
readerID = as.character(1:J)
datasetCombined <- convert2dataset(NL, LL, LL_IL,
perCase, IDs, weights,
fileName, type, name, truthTableStr, design,
modalityID, readerID)
}
stats1 <- St(datasetCombined, FOM = FOM, method = "OR", alpha = alpha, analysisOption = "RRRC", FPFValue = FPFValue)
thetajc <- stats1$FOMs$foms
thetajc <- as.matrix(thetajc)
fomCAD <- thetajc[1,1]
fomRAD <- thetajc[2,]
avgRadFom <- mean(fomRAD)
varDiffFom <- var(fomRAD)
avgDiffFom <- avgRadFom - fomCAD
FStat <- stats1$RRRC$FTests$FStat[1]
ddf <- stats1$RRRC$FTests$DF[2]
ndf <- stats1$RRRC$FTests$DF[1]
pval <- stats1$RRRC$FTests$p[1]
varR <- stats1$ANOVA$VarCom["VarR", "Estimates"] # corrected 11/26/2020
varTR <- stats1$ANOVA$VarCom["VarTR", "Estimates"] # do:
varError <- stats1$ANOVA$VarCom["Var", "Estimates"] # do:
cov1 <- stats1$ANOVA$VarCom["Cov1", "Estimates"] # do:
cov2 <- stats1$ANOVA$VarCom["Cov2", "Estimates"] # do:
cov3 <- stats1$ANOVA$VarCom["Cov3", "Estimates"] # do:
# corrected 11/26/2020; added minus sign in following ...
ciDiffFom <- stats1$RRRC$ciDiffTrt
# correction needed
# as St returns CAD - RAD (trt1 - trt2)
# not RAD - CAD
t <- ciDiffFom
rowNames <- rownames(t)
x1 <- strsplit(rowNames, split = "-")[[1]]
rowNames <- paste0(x1[2], "-", x1[1])
rownames(t) <- rowNames
t$Estimate <- -ciDiffFom$Estimate
t$t <- -ciDiffFom$t
t$CILower <- -ciDiffFom$CIUpper
t$CIUpper <- -ciDiffFom$CILower
ciDiffFom <- t
# end corrections
ciAvgRdrEachTrt <- stats1$RRRC$ciAvgRdrEachTrt
return (list (
fomCAD = fomCAD,
fomRAD = fomRAD,
avgRadFom = avgRadFom,
avgDiffFom = avgDiffFom,
varDiffFom = varDiffFom,
ciDiffFom = ciDiffFom,
ciAvgRdrEachTrt = ciAvgRdrEachTrt,
varR = varR,
varTR = varTR,
cov1 = cov1,
cov2 = cov2,
cov3 = cov3,
varError = varError,
FStat = FStat,
ndf = ndf,
df = ddf, # for consistency with 1T-RRFC
pval = pval
))
}
CadVsRadPlots <- function(dataset, FOM) {
ret1 <- dataset2ratings(dataset, FOM)
zjk1 <- ret1$zjk1
zjk2 <- ret1$zjk2
type <- dataset$descriptions$type
if (type == "ROC") {
# fixed one of hard coding errors noticed by Alejandro
genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "ROC")$Plot
} else if ((type == "LROC") && ((FOM == "PCL") || (FOM == "ALROC"))) {
if (type == "LROC") genericPlot <- LrocPlots (zjk1, zjk2, seq(1,length(zjk1[,1])-1))$lrocPlot
} else if ((type == "LROC") && (FOM == "Wilcoxon")) {
if (type == "LROC") {
datasetRoc <- DfLroc2Roc(dataset)
genericPlot <- PlotEmpiricalOperatingCharacteristics(datasetRoc, rdrs = 1:length(zjk1[,1]), opChType = "ROC")$Plot
}
} else if ((type == "FROC") && (FOM %in% c("HrAuc", "AFROC", "wAFROC"))) {
if (FOM == "HrAuc") genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "ROC")$Plot
if (FOM == "AFROC") genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "AFROC")$Plot
if (FOM == "wAFROC") genericPlot <- PlotEmpiricalOperatingCharacteristics(dataset, rdrs = 1:length(zjk1[,1]), opChType = "wAFROC")$Plot
}else stop("data type has to be ROC, FROC or LROC")
return(genericPlot)
}
##
## this uses the jackknife to estimate var & cov2
## Handles all datasets
DiffFomVarCov2 <- function (dataset, FOM, FPFValue) # for difference FOM, radiologist minus CAD
{
#if (dataset$descriptions$type == "LROC") stop("Dataset must NOT be LROC")
J <- length(dataset$descriptions$readerID)
K <- length(dataset$ratings$NL[1,1,,1])
dsCad <- DfExtractDataset(dataset, trts = 1, rdrs = 1)
dsRad <- DfExtractDataset(dataset, trts = 1, rdrs = c(2:J))
jkFomValuesCad <- UtilPseudoValues(dsCad, FOM, FPFValue)$jkFomValues
jkFomValuesRad <- UtilPseudoValues(dsRad, FOM, FPFValue)$jkFomValues
jkDiffFomValues <- array(dim = c(J-1,K))
# this does not work wasted much time
# this does not work wasted much time
# jkDiffFomValues <- jkFomValuesRad[1,,] - jkFomValuesCad[1,1,] # this does not work
# this does not work wasted much time
# this does not work wasted much time
for (j in 1:(J-1)) jkDiffFomValues[j,] <- jkFomValuesRad[1,j,] - jkFomValuesCad[1,1,]
J <- J - 1 # number of human readers
Covariance <- array(dim = c(J, J))
for (j in 1:J){
for (jp in 1:J){
Covariance[j, jp] <- cov(jkDiffFomValues[j, ], jkDiffFomValues[jp, ])
}
}
varError <- 0;count <- 0
for (j in 1:J){
varError <- varError + Covariance[j, j]
count <- count + 1
}
varError <- varError / count
varError <- varError *(K-1)^2/K
Cov2 <- 0;count <- 0
for (j in 1:J){
for (jp in 1:J){
if (jp != j){
Cov2 <- Cov2 + Covariance[j, jp]
count <- count + 1
}
}
}
Cov2 <- Cov2 / count
Cov2 <- Cov2 *(K-1)^2/K
return (list (
var = varError,
cov2 = Cov2
))
}
LrocOperatingPointsFromRatings <- function( zk1, zk2Cl )
{
if (FALSE) { # this is to check the method, see ChkLrocFoms.xlsx
zk1 <- c(rep(0, 4), 1.424443, rep(-50,5))
zk2Cl <- c(0.5542791, 1.2046176, -50, 3.4596787, 2.732657)
K2 <- length(zk2Cl)
K1 <- length(zk1) - K2
zk1 <- zk1[zk1 != -50]
zk2Cl <- zk2Cl[zk2Cl != -50]
} else {
K2 <- length(zk2Cl)
K1 <- length(zk1) - K2
zk1 <- zk1[is.finite(zk1)]
zk2Cl <- zk2Cl[is.finite(zk2Cl)]
}
FPF <- 1
PCL <- NULL
zk1Temp <- zk1;zk2ClTemp <- zk2Cl
while(1) {
cutoff <- min( c( zk1Temp, zk2ClTemp ) )
zk1Temp <- zk1[ zk1 > cutoff ]
zk2ClTemp <- zk2Cl[ zk2Cl > cutoff ]
FPF1 <- length( zk1Temp ) / K1
PCL1 <- length( zk2ClTemp ) / K2
FPF <- c( FPF, FPF1 )
if (length(PCL) == 0) PCL <- c(PCL1, PCL1) else PCL <- c( PCL, PCL1 )
if( FPF1 == 0 && PCL1 == 0 ) {
break
}
}
FPF <- rev(FPF)
PCL <- rev(PCL)
return( list(
FPF = FPF,
PCL = PCL
) )
}
LrocPlots <- function (zjk1, zjk2, doJ)
{
J <- length(zjk1[,1])
j <- 1;zjk1Temp <- zjk1[j,];zk2Temp <- zjk2[j,]
lroc <- LrocOperatingPointsFromRatings( zjk1Temp, zk2Temp )
FPF <- lroc$FPF;PCL <- lroc$PCL
lrocPlotData <- data.frame(FPF = FPF, PCL = PCL, reader = "R-CAD")
for (j in 2:J) {
if ((j - 1) %in% doJ) {
zjk1Temp <- zjk1[j,]
zk2Temp <- zjk2[j,]
lroc <- LrocOperatingPointsFromRatings( zjk1Temp, zk2Temp )
FPF <- lroc$FPF
PCL <- lroc$PCL
reader = paste0("R-", as.character(j - 1))
lrocPlotData <- rbind(lrocPlotData, data.frame(FPF = FPF, PCL = PCL, reader = reader))
}
}
lrocPlot <- ggplot(data = lrocPlotData, aes(x = FPF, y = PCL, color = reader)) + geom_line()
g <- ggplot_build(lrocPlot)
colors <- as.character(unique(g$data[[1]]$colour))
colors[1] <- "#000000"
sizes <- c(2, rep(1, length(doJ)))
lrocPlot <- ggplot(data = lrocPlotData, aes(x = FPF, y = PCL, color = reader)) + geom_line(aes(linewidth = reader)) +
scale_color_manual(values = colors) + scale_size_manual(values = sizes) +
theme(legend.title = element_blank(), legend.position = legendPosition <- c(1, 0), legend.justification = c(1, 0))
return(list(
lrocPlot = lrocPlot,
afrocPlot = NULL)
)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.