#' Fit the binormal model to selected modality and reader in an ROC dataset
#'
#' @description Fit the binormal model-predicted ROC curve for a dataset.
#' This is the R equivalent of ROCFIT or RSCORE
#'
#' @param dataset The ROC dataset
#' @param trt The desired modality, default is 1
#' @param rdr The desired reader, default is 1
#'
#' @return The returned value is a list with the following elements:
#' @return \item{a}{The mean of the diseased distribution;
#' the non-diseased distribution is assumed to have zero mean}
#' @return \item{b}{The standard deviation of the non-diseased
#' distribution. The diseased distribution is assumed to have
#' unit standard deviation}
#' @return \item{zetas}{The binormal model cutoffs, zetas or thresholds}
#' @return \item{AUC}{The binormal model fitted ROC-AUC}
#' @return \item{StdAUC}{The standard deviation of AUC}
#' @return \item{NLLIni}{The initial value of negative LL}
#' @return \item{NLLFin}{The final value of negative LL}
#' @return \item{ChisqrFitStats}{The chisquare goodness of fit results}
#' @return \item{covMat}{The covariance matrix of the parameters}
#' @return \item{fittedPlot}{A \pkg{ggplot2} object containing the
#' fitted operating characteristic along with the empirical operating
#' points. Use \code{print()} to display the object}
#'
#' @details
#' In the binormal model ratings (more accurately the latent decision variables)
#' from diseased cases are sampled from \eqn{N(a,1)} while ratings for
#' non-diseased cases are sampled from \eqn{N(0,b^2)}. To avoid clutter error
#' bars are only shown for the lowest and uppermost operating points. An FROC
#' dataset is internally converted to a highest rating inferred ROC dataset. To
#' many bins containing zero counts will cause the algorithm to fail; so be sure
#' to bin the data appropriately to fewer bins, where each bin has at least one
#' count.
#'
#'
#' @examples
#' \donttest{
#' ## Test with an included ROC dataset
#' retFit <- FitBinormalRoc(dataset02);## print(retFit$fittedPlot)
#'
#'
#' ## Test with an included FROC dataset; it needs to be binned
#' ## as there are more than 5 discrete ratings levels
#' binned <- DfBinDataset(dataset05, desiredNumBins = 5, opChType = "ROC")
#' retFit <- FitBinormalRoc(binned);## print(retFit$fittedPlot)
#'
#'
#' ## Test with single interior point data
#' fp <- c(rep(1,7), rep(2, 3))
#' tp <- c(rep(1,5), rep(2, 5))
#' dataset <- Df2RJafrocDataset(fp, tp)
#' retFit <- FitBinormalRoc(dataset);## print(retFit$fittedPlot)
#'
#' ## Test with two interior data points
#' fp <- c(rep(1,7), rep(2, 5), rep(3, 3))
#' tp <- c(rep(1,3), rep(2, 5), rep(3, 7))
#' dataset <- Df2RJafrocDataset(fp, tp)
#' retFit <- FitBinormalRoc(dataset);## print(retFit$fittedPlot)
#'
#' ## Test with TONY data for which chisqr can be calculated
#' ds <- DfFroc2Roc(dataset01)
#' retFit <- FitBinormalRoc(ds, 2, 3);## print(retFit$fittedPlot)
#' ## retFit$ChisqrFitStats
#'
#' ## Test with included degenerate ROC data
#' retFit <- FitBinormalRoc(datasetDegenerate);## print(retFit$fittedPlot)
#' }
#'
#'
#'
#' @references
#' Dorfman DD, Alf E (1969) Maximum-Likelihood Estimation of Parameters of Signal-Detection Theory and Determination of Confidence Intervals -
#' Rating-Method Data, Journal of Mathematical Psychology 6, 487-496.
#'
#' Grey D, Morgan B (1972) Some aspects of ROC curve-fitting: normal and logistic models. Journal of Mathematical Psychology 9, 128-139.
#'
#' @importFrom bbmle mle2
#' @importFrom stats lm
#' @importFrom binom binom.confint
#' @export
#'
FitBinormalRoc <- function(dataset, trt = 1, rdr = 1){
if (dataset$descriptions$type == "FROC") dataset <- DfFroc2Roc(dataset)
minA <- RJafrocEnv$minA
maxA <- RJafrocEnv$maxA
minB <- RJafrocEnv$minB
maxB <- RJafrocEnv$maxB
minZeta <- RJafrocEnv$minZeta
maxZeta <- RJafrocEnv$maxZeta
errorMsg <- "" # keep track of any generated warnings
fp <- dataset$ratings$NL[trt,rdr,,1];fp <- fp[fp != -Inf]
tp <- dataset$ratings$LL[trt,rdr,,1]
lenZetas <- length(unique(c(fp, tp))) - 1
zetas <- array(dim = lenZetas)
plotZeta <- seq(-3, 8, by = 0.1)
ret1 <- RawOpPtsROC2ROC(fp, tp)
fpCounts <- ret1$fpCounts;tpCounts <- ret1$tpCounts;fpf <- ret1$fpf;tpf <- ret1$tpf
temp <- isDataDegenerate (fpf, tpf)
msg <- temp$msg
if (msg != "") errorMsg <- paste0(errorMsg, msg)
if (temp$ret) {
warning("Data is degenerate; binormal model cannot fit it unambiguously: use CBM or RSM method")
return(list(
mu = NA,
alpha = NA,
zetas = NA,
AUC = NA,
StdAUC = NA,
NLLIni = NA,
NLLFin = NA,
covMat = NA,
fittedPlot = NA,
errorMsg = errorMsg
))
}
if (length(fpf) != 1) {
phiInvFpf <- qnorm(fpf)
phiInvTpf <- qnorm(tpf)
phiInvFpf <- qnorm(fpf); phiInvFpf[phiInvFpf == -Inf] <- minZeta; phiInvFpf[phiInvFpf == Inf] <- maxZeta
phiInvTpf <- qnorm(tpf); phiInvTpf[phiInvTpf == -Inf] <- minZeta; phiInvTpf[phiInvTpf == Inf] <- maxZeta
fit <- lm(phiInvTpf ~ phiInvFpf) # straight line fit method of estimating a and b
aIni <- fit$coefficients[[1]] # these is the initial estimate of a
bIni <- fit$coefficients[[2]] # these is the initial estimate of b
# thresholds can be estimated by by applying inverse function to Eqn. xx and solving for zetas
zetaNorIni <- rev(-phiInvFpf) # see Eqn. xx
zetaAbnIni <- rev((aIni - phiInvTpf)/bIni) # see Eqn. xx
zetasIni <- (zetaNorIni + zetaAbnIni) / 2 # average the two estimates
} else {
aIni <- 1
bIni <- 1
zetasIni <- 0
}
a <- aIni; b <- bIni;zetas <- zetasIni
# while (1){
#param <- c(a, b, zetas)
namesVector <- c(c("aFwd", "bFwd"), paste0("zetaFwd", 1:length(zetas)))
aFwd <- ForwardValue(a, minA, maxA)
bFwd <- ForwardValue(b, minB, maxB)
zetaFwd <- ForwardZetas(zetas)
parameters <- c(list(aFwd, bFwd), as.list(zetaFwd))
names(parameters) <- namesVector
BMNLLNew <- AddArguments(BMNLL, length(zetas))
ret <- mle2(BMNLLNew, start = parameters, method = "BFGS", eval.only = TRUE,
data = list(fi = fpCounts, ti = tpCounts))
NLLIni <- ret@min
ret <- mle2(BMNLLNew, start = parameters, method = "BFGS",
data = list(fi = fpCounts, ti = tpCounts))
if (ret@details$convergence != 0) stop("mle2 error in binormal fit")
NLLFin <- ret@min
vcov <- ret@vcov
a <- InverseValue(ret@coef[1], minA, maxA)
b <- InverseValue(ret@coef[2], minB, maxB)
zetas <- InverseZetas(ret@coef[3:length(ret@coef)])
AUC <- pnorm(a / sqrt(1 + b^2))
covMat <- vcov[1:2,1:2]
StdAUC <- StdDevBinormalAuc1(ret@coef[1], ret@coef[2], covMat) ## !!!dpc!!! looks right; can it be proved?
ChisqrFitStats <- ChisqrGoodnessOfFit(fpCounts, tpCounts,
parameters = c(a,b,zetas), model = "BINORMAL")
fpfPred <- 1 - pnorm(plotZeta)
tpfPred <- pnorm(a - b * plotZeta)
fittedPlot <- genericPlotROC (fp, tp, fpfPred, tpfPred)
# calculate covariance matrix using un-transformed variables
namesVector <- c(c("a", "b"), paste0("zeta", 1:length(zetas)))
parameters <- c(list(a, b), as.list(zetas))
names(parameters) <- namesVector
BMNLLNewNoTransf <- AddArguments2(BMNLLNoTransf, length(zetas))
ret <- mle2(BMNLLNewNoTransf, start = parameters, method = "BFGS",
data = list(fi = fpCounts, ti = tpCounts))
# NLLChk <- ret@min
# if (NLLChk < NLLFin) {
# if (abs(NLLChk - NLLFin)/NLLFin > 1e-6) next else break
# } else break
#
# }
covMat <- ret@vcov
return(list(
a = a,
b = b,
zetas = zetas,
AUC = AUC,
StdAUC = StdAUC,
NLLIni = NLLIni,
NLLFin = NLLFin,
ChisqrFitStats = ChisqrFitStats,
covMat = covMat,
fittedPlot = fittedPlot,
errorMsg = errorMsg
))
}
BMNLL <- function(aFwd, bFwd, fi, ti){
minA <- RJafrocEnv$minA
maxA <- RJafrocEnv$maxA
minB <- RJafrocEnv$minB
maxB <- RJafrocEnv$maxB
a <- InverseValue(aFwd, minA, maxA)
b <- InverseValue(bFwd, minB, maxB)
allParameters <- names(formals())
zetaPos <- regexpr("zeta", allParameters)
zetaFwd <- unlist(mget(allParameters[which(zetaPos == 1)]))
zeta <- InverseZetas(zetaFwd)
zetas <- c(-Inf, zeta, Inf)
L <- BMNLLInner(a, b, zetas, fi, ti)
return(-L)
}
# ROC paradigm negative of log likelihood function, without transformations
BMNLLNoTransf <- function (a, b, fi, ti){
allParameters <- names(formals())
zetaPos <- regexpr("zeta", allParameters)
zetas <- unlist(mget(allParameters[which(zetaPos == 1)]))
zetas <- c(-Inf, zetas, Inf)
L <- BMNLLInner(a, b, zetas, fi, ti)
return(-L)
}
# analytical derivatives using Cov derived using no transformations
StdDevBinormalAuc <- function (a, b,Cov)
{
derivWrtA <- dnorm (a/sqrt(1+b^2))/sqrt(1+b^2)
derivWrtB <- dnorm (a/sqrt(1+b^2))*(-a*b*(1+b^2)^(-1.5))
VarAUC <- (derivWrtA)^2*Cov[1,1]+(derivWrtB)^2*Cov[2,2] +
2 * derivWrtA*derivWrtB*Cov[1,2]
return (sqrt(VarAUC))
}
# inputs are transformed parameters; covMat is also wrt trans. parameters
StdDevBinormalAuc1 <- function (aFwd, bFwd,covMatFwd)
{
minA <- RJafrocEnv$minA
maxA <- RJafrocEnv$maxA
minB <- RJafrocEnv$minB
maxB <- RJafrocEnv$maxB
derivsFwd <- jacobian(func = tempAucBinormal, c(aFwd,bFwd))
VarAz <- derivsFwd %*% covMatFwd %*% t(derivsFwd)
if (VarAz < 0 || is.na(VarAz)) return(NA) else return(sqrt(VarAz))
}
tempAucBinormal <- function (fwdParms){
minA <- RJafrocEnv$minA
maxA <- RJafrocEnv$maxA
minB <- RJafrocEnv$minB
maxB <- RJafrocEnv$maxB
a <- InverseValue(fwdParms[1], minA, maxA)
b <- InverseValue(fwdParms[2], minB, maxB)
AUC <- pnorm(a / sqrt(1 + b^2))
return (AUC)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.