R/FitCorCbm.R
In dpc10ster/rjafroc: Artificial Intelligence Systems and Observer Performance
Defines functions
PlotCorCbmFit
LLCorCBM
nLLCorCBMNoTrns
nLLCorCBM
StatsCorCbm
FitCorCbm
Documented in
FitCorCbm
#' Fit CORCBM to a paired ROC dataset
#'
#' @description Fit the Correlated Contaminated Binormal Model (CORCBM)
#' to a paired ROC dataset.
#' The \strong{ROC} dataset has to be formatted as a
#' \strong{single modality}, \strong{two-reader} dataset, even though the actual
#' pairing may be different, see details.
#'
#'
#' @param dataset A \strong{paired ROC} dataset
#'
#'
#' @return A list containing three objects:
#' @return \item{fitCorCbmRet}{list(\code{FPCounts},\code{TPCounts},
#' \code{muX},\code{muY},\code{alphaX},\code{alphaY},\code{rhoNor},
#' \code{rhoAbn2},\code{zetaX},\code{zetaY},\code{covMat},\code{fixParam})}
#' @return \item{stats}{list(\code{aucX},\code{aucX},\code{stdAucX},
#' \code{stdAucY},\code{stdErr},\code{areaStat},\code{areaPval})}
#' @return \item{fittedPlot}{The fitted plot with operating points, error bars,
#' for both conditions}
#'
#' @details The conditions (X, Y) can be two readers interpreting images in the same
#' modality, the same reader interpreting images in different treatments, or
#' different readers interpreting images in 2 different treatments. Function
#' \code{\link{DfExtractCorCbmDataset}} can be used to construct a dataset suitable for
#' \code{FitCorCbm}. With reference to the returned values, and assuming R bins
#' in condition X and L bins in conditon Y,
#' \code{FPCounts} is the R x L matrix containing the counts for non-diseased cases,
#' \code{TPCounts} is the R x L matrix containing the counts for diseased cases;
#' \code{muX},\code{muY},\code{alphaX},\code{alphaY},\code{rhoNor},\code{rhoAbn2} are
#' the CORCBM parameters; \code{aucX},\code{aucX} are the AUCs in the two conditions;
#' \code{stdAucX},\code{stdAucY} are the corresponding standard errors;\code{stdErr}
#' contains the standard errors of the parameters of the model; \code{areaStat},
#' \code{areaPval},\code{covMat} are the area-statistic, the p-value and the covariance
#' matrix of the parameters. If a parameter approaches a limit, e.g., \code{rhoNor}
#' = 0.9999, it is held constant at near the limiting value and the covariance matrix
#' has one less dimension (along each edge) for each parameter that is held constant.
#' The indices of the parameters held fixed are in \code{fitCorCbmRet$fixParam}.
#'
## following examples generate excessive CPU time NOTES on devtools::check_win_xx()
## but it is instructive to execute them independently
## see also the publication referenced below
##
## dataset <- DfExtractCorCbmDataset(dataset05, trts = 1, rdrs = c(4,7))
## ret <- FitCorCbm(dataset)
## print(ret$fitCorCbmRet)
## print(ret$stats)
## print(ret$fittedPlot)
##
##
## ret <- FitCorCbm(datasetBinned123)
## print(ret$fitCorCbmRet)
## print(ret$stats)
## print(ret$fittedPlot)
## Also try two other datasets ending with 124 and 125
#'
#'
#' @references
#' Zhai X, Chakraborty DP (2017) A bivariate contaminated binormal model for robust fitting of proper ROC
#' curves to a pair of correlated, possibly degenerate, ROC datasets. Medical Physics. 44(6):2207--2222.
#'
#'
#' @importFrom bbmle mle2
#' @importFrom Rcpp cppFunction
#' @importFrom mvtnorm pmvnorm
#' @importFrom stats cor
#'
#' @export
#'
#'
FitCorCbm <- function(dataset){
if (dataset$descriptions$type != "ROC") {
stop("This program requires an ROC dataset")
}
minMu <- RJafrocEnv$minMu
maxMu <- RJafrocEnv$maxMu
minAlpha <- RJafrocEnv$minAlpha
maxAlpha <- RJafrocEnv$maxAlpha
minRho <- RJafrocEnv$minRho
maxRho <- RJafrocEnv$maxRho
# `as.matrix` is NOT absolutely necessary as `mean()` function is not used here
aucArray <- UtilFigureOfMerit(dataset, FOM = "Wilcoxon")
maxAUC <- max(aucArray)
while (pnorm(maxMu / sqrt(2)) <= maxAUC){
maxMu <- qnorm(maxAUC) * sqrt(2) + 0.5
}
I <- length(dataset$ratings$NL[,1,1,1])
J <- length(dataset$ratings$NL[1,,1,1])
K <- length(dataset$ratings$NL[1,1,,1])
K2 <- length(dataset$ratings$LL[1,1,,1])
K1 <- K - K2
FP <- dataset$ratings$NL[1,,1:K1,1]
TP <- dataset$ratings$LL[1,,,1]
z1X <- FP[1, ];z1Y <- FP[2, ]
z2X <- TP[1, ];z2Y <- TP[2, ]
rhoNor <- cor(z1X, z1Y)
rhoAbn <- cor(z2X, z2Y)
iniParamX <- FitCbmRoc(dataset, 1, 1) # arm X of pairing
muX <- iniParamX$mu
zetaX <- iniParamX$zetas
alphaX <- iniParamX$alpha
iniParamY <- FitCbmRoc(dataset, 1, 2) # arm Y of pairing
muY <- iniParamY$mu
zetaY <- iniParamY$zetas
alphaY <- iniParamY$alpha
muXFwd <- ForwardValue(muX, minMu, maxMu)
muYFwd <- ForwardValue(muY, minMu, maxMu)
alphaXFwd <- ForwardValue(alphaX, minAlpha, maxAlpha)
alphaYFwd <- ForwardValue(alphaY, minAlpha, maxAlpha)
rhoNorFwd <- ForwardValue(rhoNor, minRho, maxRho)
rhoAbn2Fwd <- ForwardValue(rhoAbn, minRho, maxRho)
zetaXFwd <- ForwardZetas(zetaX)
zetaYFwd <- ForwardZetas(zetaY)
parameters <- c(list(muXFwd, muYFwd, alphaXFwd, alphaYFwd, rhoNorFwd, rhoAbn2Fwd),
as.list(zetaXFwd), as.list(zetaYFwd))
namesVector <- c("muXFwd", "muYFwd", "alphaXFwd", "alphaYFwd","rhoNorFwd", "rhoAbn2Fwd")
namesVector <- c(namesVector, paste0("zetaXFwd", 1:length(zetaXFwd)))
namesVector <- c(namesVector, paste0("zetaYFwd", 1:length(zetaYFwd)))
names(parameters) <- namesVector
nBinsX <- length(unique(c(FP[1,], TP[1,])))
nBinsY <- length(unique(c(FP[2,], TP[2,])))
FPCounts <- array(dim = c(nBinsX, nBinsY))
TPCounts <- array(dim = c(nBinsX, nBinsY))
for (bX in 1:nBinsX){
for (bY in 1:nBinsY){
FPCounts[bX, bY] <- sum((z1X == bX & z1Y == bY))
TPCounts[bX, bY] <- sum((z2X == bX & z2Y == bY))
}
}
nLLCorCBMAdd <- AddZetaXFwd(nLLCorCBM, length(zetaX))
nLLCorCBMAdd <- AddZetaYFwd(nLLCorCBMAdd, length(zetaY))
ret <- mle2(nLLCorCBMAdd, start = parameters, eval.only = TRUE,
method = "BFGS", data = list(FPCounts = FPCounts, TPCounts = TPCounts, maxMu = maxMu))
NLLIni <- as.numeric(ret@min)
ret <- mle2(nLLCorCBMAdd, start = parameters,
method = "BFGS", data = list(FPCounts = FPCounts, TPCounts = TPCounts, maxMu = maxMu))
NLLFin <- as.numeric(ret@min)
allCoef1 <- ret@coef
muX <- InverseValue(allCoef1[1], minMu, maxMu)
muY <- InverseValue(allCoef1[2], minMu, maxMu)
alphaX <- InverseValue(allCoef1[3], minAlpha, maxAlpha)
alphaY <- InverseValue(allCoef1[4], minAlpha, maxAlpha)
rhoNor <- InverseValue(allCoef1[5], minRho, maxRho)
rhoAbn2 <- InverseValue(allCoef1[6], minRho, maxRho)
zetaX <- InverseZetas(allCoef1[7:(6 + length(zetaX))])
zetaY <- InverseZetas(allCoef1[(7 + length(zetaX)):length(allCoef1)])
CorCBMLLNoTrnsAdd <- AddZetaX(nLLCorCBMNoTrns, length(zetaX))
CorCBMLLNoTrnsAdd <- AddZetaY(CorCBMLLNoTrnsAdd, length(zetaY))
param <- c(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2, zetaX, zetaY)
fixParam <- c(which(param[1:2] > (maxMu - 0.01)),
2+which(param[3:4] > (maxAlpha - 0.01)),
4+which(param[5:6] > (maxRho - 0.01)))
if (length(fixParam) > 0){
# some parameters are fixed
fixList <- list()
afterIndex <- c()
for (p in fixParam){
if (p == 1){
muX <- maxMu - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(muXFwd = ForwardValue(muX, minMu, maxMu)))
} else if (p == 2){
muY <- maxMu - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(muYFwd = ForwardValue(muY, minMu, maxMu)))
}else if (p == 3){
alphaX <- maxAlpha - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(alphaXFwd = ForwardValue(alphaX, minAlpha, maxAlpha)))
}else if (p == 4){
alphaY <- maxAlpha - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(alphaYFwd = ForwardValue(alphaY, minAlpha, maxAlpha)))
}else if (p == 5){
rhoNor <- ifelse(test = rhoNor > 0, yes = 0.99, no = -0.99)
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(rhoNorFwd = ForwardValue(rhoNor, minRho, maxRho)))
}else if (p == 6){
rhoAbn2 <- ifelse(test = rhoAbn2 > 0, yes = 0.99, no = -0.99)
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(rhoAbn2Fwd = ForwardValue(rhoAbn2, minRho, maxRho)))
}
}
muXFwd <- ForwardValue(muX, minMu, maxMu)
muYFwd <- ForwardValue(muY, minMu, maxMu)
alphaXFwd <- ForwardValue(alphaX, minAlpha, maxAlpha)
alphaYFwd <- ForwardValue(alphaY, minAlpha, maxAlpha)
rhoNorFwd <- ForwardValue(rhoNor, minRho, maxRho)
rhoAbn2Fwd <- ForwardValue(rhoAbn2, minRho, maxRho)
zetaXFwd <- ForwardZetas(zetaX)
zetaYFwd <- ForwardZetas(zetaY)
parameters <- c(list(muXFwd, muYFwd, alphaXFwd, alphaYFwd, rhoNorFwd, rhoAbn2Fwd),
as.list(zetaXFwd), as.list(zetaYFwd))
namesVector <- c("muXFwd", "muYFwd", "alphaXFwd", "alphaYFwd","rhoNorFwd", "rhoAbn2Fwd")
namesVector <- c(namesVector, paste0("zetaXFwd", 1:length(zetaXFwd)))
namesVector <- c(namesVector, paste0("zetaYFwd", 1:length(zetaYFwd)))
names(parameters) <- namesVector
ret2 <- mle2(nLLCorCBMAdd, start = parameters, method = "BFGS", fixed = fixList,
data = list(FPCounts = FPCounts, TPCounts = TPCounts, maxMu = maxMu), trace = TRUE)
NLLFin <- ret2@min
allCoef2 <- ret2@coef
for (i in 1:length(afterIndex)) {
allCoef2 <- append(allCoef2, fixList[[i]], after = afterIndex[i])
}
names(allCoef2) <- namesVector
muX <- InverseValue(allCoef2[1], minMu, maxMu)
muY <- InverseValue(allCoef2[2], minMu, maxMu)
alphaX <- InverseValue(allCoef2[3], minAlpha, maxAlpha)
alphaY <- InverseValue(allCoef2[4], minAlpha, maxAlpha)
rhoNor <- InverseValue(allCoef2[5], minRho, maxRho)
rhoAbn2 <- InverseValue(allCoef2[6], minRho, maxRho)
zetaX <- InverseZetas(allCoef2[7:(6 + length(zetaX))])
zetaY <- InverseZetas(allCoef2[(7 + length(zetaX)):length(allCoef2)])
parameters <- c(list(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2),
as.list(zetaX), as.list(zetaY))
namesVector <- c("muX", "muY", "alphaX", "alphaY","rhoNor", "rhoAbn2")
namesVector <- c(namesVector, paste0("zetaX", 1:length(zetaX)))
namesVector <- c(namesVector, paste0("zetaY", 1:length(zetaY)))
names(parameters) <- namesVector
ret3 <- mle2(CorCBMLLNoTrnsAdd, start = parameters, method = "BFGS", fixed = parameters[fixParam],
data = list(FPCounts = FPCounts, TPCounts = TPCounts))
} else {
# all parameters are used
allCoef2 <- allCoef1
names(allCoef2) <- namesVector
muX <- InverseValue(allCoef2[1], minMu, maxMu)
muY <- InverseValue(allCoef2[2], minMu, maxMu)
alphaX <- InverseValue(allCoef2[3], minAlpha, maxAlpha)
alphaY <- InverseValue(allCoef2[4], minAlpha, maxAlpha)
rhoNor <- InverseValue(allCoef2[5], minRho, maxRho)
rhoAbn2 <- InverseValue(allCoef2[6], minRho, maxRho)
zetaX <- InverseZetas(allCoef2[7:(6 + length(zetaX))])
zetaY <- InverseZetas(allCoef2[(7 + length(zetaX)):length(allCoef2)])
parameters <- c(list(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2),
as.list(zetaX), as.list(zetaY))
namesVector <- c("muX", "muY", "alphaX", "alphaY","rhoNor", "rhoAbn2")
namesVector <- c(namesVector, paste0("zetaX", 1:length(zetaX)))
namesVector <- c(namesVector, paste0("zetaY", 1:length(zetaY)))
names(parameters) <- namesVector
ret3 <- mle2(CorCBMLLNoTrnsAdd, start = parameters, method = "BFGS",
data = list(FPCounts = FPCounts, TPCounts = TPCounts))
}
covMat <- ret3@vcov
fittedPlot <- PlotCorCbmFit(
list(
FPCounts = FPCounts,
TPCounts = TPCounts,
muX = muX,
alphaX = alphaX,
muY = muY,
alphaY = alphaY,
rhoNor = rhoNor,
rhoAbn2 = rhoAbn2,
zetaX = zetaX,
zetaY = zetaY
)
)
fitCorCbmRet <- list(
FPCounts = FPCounts,
TPCounts = TPCounts,
muX = muX,
alphaX = alphaX,
muY = muY,
alphaY = alphaY,
rhoNor = rhoNor,
rhoAbn2 = rhoAbn2,
zetaX = zetaX,
zetaY = zetaY,
covMat = covMat,
fixParam = fixParam
)
stats <- StatsCorCbm(fitCorCbmRet)
return(list(
fitCorCbmRet = fitCorCbmRet,
stats = stats,
fittedPlot = fittedPlot
))
}
StatsCorCbm <- function(fitCorCbmRet) {
muX <- fitCorCbmRet$muX
muY <- fitCorCbmRet$muY
alphaX <- fitCorCbmRet$alphaX
alphaY <- fitCorCbmRet$alphaY
rhoNor <- fitCorCbmRet$rhoNor
rhoAbn2 <- fitCorCbmRet$rhoAbn2
zetaX <- fitCorCbmRet$zetaX
zetaY <- fitCorCbmRet$zetaY
FPCounts <- fitCorCbmRet$FPCounts
TPCounts <- fitCorCbmRet$TPCounts
covMat <- fitCorCbmRet$covMat
fixParam <- fitCorCbmRet$fixParam
vars <- diag(covMat)
stdErr <- sqrt(vars)
dMuX <- alphaX * dnorm(muX / sqrt(2)) / sqrt(2)
dAlphaX <- -1 / 2 + pnorm(muX / sqrt(2))
stdAucX <-
sqrt(vars[1] * dMuX ^ 2 + vars[3] * dAlphaX ^ 2 + dMuX * dAlphaX * covMat[1, 3])
dMuY <- alphaY * dnorm(muY / sqrt(2)) / sqrt(2)
dAlphaY <- -1 / 2 + pnorm(muY / sqrt(2))
stdAucY <-
sqrt(vars[2] * dMuY ^ 2 + vars[4] * dAlphaY ^ 2 + dMuY * dAlphaY * covMat[2, 4])
stdErr <- c(stdErr, stdAucX, stdAucY)
if (length(fixParam) > 0) {
for (p in fixParam) {
stdErr <- append(stdErr, NA, after = p - 1)
if (p == 3) {
alphaX <- 1
} else if (p == 4) {
alphaY <- 1
} else if (p == 5) {
rhoNor <- 1
} else if (p == 6) {
rhoAbn2 <- 1
}
}
stdErr <- as.numeric(stdErr)
}
stdErr[2:3] <- stdErr[3:2] # switch the position of alphaX and muY
aucX <- (1 - alphaX) * 0.5 + alphaX * pnorm(muX / sqrt(2))
aucY <- (1 - alphaY) * 0.5 + alphaY * pnorm(muY / sqrt(2))
aucDiff <- aucX - aucY
varDiff <- 0
derivs <- c(dMuX,-dMuY, dAlphaX,-dAlphaY)
for (k in 1:4) {
for (l in 1:4) {
varDiff <- derivs[l] * derivs[k] * covMat[k, l] + varDiff
}
}
stdDiff <- sqrt(varDiff)
areaStat <- abs(aucDiff) / stdDiff
areaPval <- 1 - pnorm(areaStat)
return(
list(
aucX = aucX,
aucY = aucY,
stdAucX = as.numeric(stdAucX),
stdAucY = as.numeric(stdAucY),
stdErr = stdErr,
areaStat = areaStat,
areaPval = areaPval
)
)
}
###############################################################################
nLLCorCBM <- function(muXFwd, muYFwd, alphaXFwd, alphaYFwd, rhoNorFwd, rhoAbn2Fwd, FPCounts, TPCounts, maxMu){
minMu <- RJafrocEnv$minMu
#maxMu <- RJafrocEnv$maxMu
minAlpha <- RJafrocEnv$minAlpha
maxAlpha <- RJafrocEnv$maxAlpha
minRho <- RJafrocEnv$minRho
maxRho <- RJafrocEnv$maxRho
allParameters <- names(formals())
zetaXPos <- regexpr("zetaX", allParameters)
zetaXFwd <- unlist(mget(allParameters[which(zetaXPos == 1)]))
zetaYPos <- regexpr("zetaY", allParameters)
zetaYFwd <- unlist(mget(allParameters[which(zetaYPos == 1)]))
muX <- InverseValue(muXFwd, minMu, maxMu)
muY <- InverseValue(muYFwd, minMu, maxMu)
alphaX <- InverseValue(alphaXFwd, minAlpha, maxAlpha)
alphaY <- InverseValue(alphaYFwd, minAlpha, maxAlpha)
rhoNor <- InverseValue(rhoNorFwd, minRho, maxRho)
rhoAbn <- InverseValue(rhoAbn2Fwd, minRho, maxRho)
zetaX <- InverseZetas(zetaXFwd)
zetaY <- InverseZetas(zetaYFwd)
LL <- LLCorCBM (muX, muY, alphaX, alphaY, rhoNor, rhoAbn, zetaX, zetaY, FPCounts, TPCounts)
return(-LL)
}
###############################################################################
nLLCorCBMNoTrns <- function(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2, FPCounts, TPCounts){
allParameters <- names(formals())
zetaXPos <- regexpr("zetaX", allParameters)
zetaX <- unlist(mget(allParameters[which(zetaXPos == 1)]))
zetaYPos <- regexpr("zetaY", allParameters)
zetaY <- unlist(mget(allParameters[which(zetaYPos == 1)]))
LL <- LLCorCBM (muX, muY, alphaX, alphaY, rhoNor, rhoAbn2, zetaX, zetaY, FPCounts, TPCounts)
return(-LL)
}
###############################################################################
LLCorCBM <- function(muX, muY, alphaX, alphaY, rhoNor, rhoAbn, zetaX, zetaY, FPCounts, TPCounts) {
if ((abs(zetaX[1]) > 2) ||
(abs(zetaY[1]) > 2) ||
(zetaX[length(zetaX)] > 10) ||
(zetaY[length(zetaY)] > 10)) {
return(-1e10)
}
zetaX <- c(-Inf, zetaX, Inf)
zetaY <- c(-Inf, zetaY, Inf)
if (is.unsorted(zetaX)) return(-1e10)
if (is.unsorted(zetaY)) return(-1e10)
rhoAbn1 <- rhoNor
rhoAbn2 <- rhoAbn
rhoAbn12 <- mean(c(rhoAbn1, rhoAbn2))
rhoAbn4 <- rhoAbn12
sigmaNor <- rbind(c(1, rhoNor), c(rhoNor, 1))
sigmaAbn1 <- rbind(c(1, rhoAbn1), c(rhoAbn1, 1))
sigmaAbn2 <- rbind(c(1, rhoAbn2), c(rhoAbn2, 1))
sigmaAbn3 <- rbind(c(1, rhoAbn12), c(rhoAbn12, 1))
sigmaAbn4 <- rbind(c(1, rhoAbn4), c(rhoAbn4, 1))
LLNor <- 0
LLAbn <- 0
nBinsX <- dim(FPCounts)[1]
nBinsY <- dim(FPCounts)[2]
for (bX in 1:nBinsX){
for (bY in 1:nBinsY){
if (FPCounts[bX, bY] > 0) {
tempNor <- pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]), sigma = sigmaNor)
if (tempNor >= 0) {
LLNor <- LLNor +
FPCounts[bX, bY] * log(tempNor)
} else {
return(-1e10)
}
}
if (TPCounts[bX, bY] > 0) {
tempAbn <- (1 - alphaX) * (1 - alphaY) *
pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]), mean = c(0, 0), sigma = sigmaAbn1) +
(alphaX) * (1 - alphaY) * pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]),
mean = c(muX, 0), sigma = sigmaAbn3) +
(1 - alphaX) * (alphaY) * pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]),
mean = c(0, muY), sigma = sigmaAbn4) +
alphaX * alphaY * pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]),
mean = c(muX, muY), sigma = sigmaAbn2)
if (tempAbn >= 0) {
LLAbn <- LLAbn +
TPCounts[bX, bY] * log(tempAbn)
} else {
return(-1e10)
}
}
}
}
LL <- LLNor + LLAbn
return(LL)
}
#' @import ggplot2
PlotCorCbmFit <- function(retFitCorCBM){
muX <- retFitCorCBM$muX
muY <- retFitCorCBM$muY
alphaX <- retFitCorCBM$alphaX
alphaY <- retFitCorCBM$alphaY
FPCounts <- retFitCorCBM$FPCounts
TPCounts <- retFitCorCBM$TPCounts
K1 <- sum(FPCounts);K2 <- sum(TPCounts)
plotZeta <- seq(-3, max(muX,muY)+2, by = 0.1)
plotCBM <- NULL
plotOpPnts <- NULL
FPFX <- 1 - pnorm(plotZeta)
TPFX <- (1 - alphaX) * (1 - pnorm(plotZeta)) + alphaX * (1 - pnorm(plotZeta, mean = muX))
plotCBM <- rbind(plotCBM, data.frame(FPF = FPFX, TPF = TPFX, Condition = "X"))
FPFX <- cumsum(rev(rowSums(FPCounts)))/K1
TPFX <- cumsum(rev(rowSums(TPCounts)))/K2
FPFX <- FPFX[-length(FPFX)]
TPFX <- TPFX[-length(TPFX)]
plotOpPnts <- rbind(plotOpPnts, data.frame(FPF = FPFX, TPF = TPFX, Condition = "X"))
FPFY <- 1 - pnorm(plotZeta)
TPFY <- (1 - alphaY) * (1 - pnorm(plotZeta)) + alphaY * (1 - pnorm(plotZeta, mean = muY))
plotCBM <- rbind(plotCBM, data.frame(FPF = FPFY, TPF = TPFY, Condition = "Y"))
FPFY <- cumsum(rev(colSums(FPCounts)))/K1
TPFY <- cumsum(rev(colSums(TPCounts)))/K2
FPFY <- FPFY[-length(FPFY)]
TPFY <- TPFY[-length(TPFY)]
plotOpPnts <- rbind(plotOpPnts, data.frame(FPF = FPFY, TPF = TPFY, Condition = "Y"))
Condition <- NULL;FPF <- NULL;TPF <- NULL
fittedPlot <- ggplot(data = plotCBM, mapping = aes(x = FPF, y = TPF, color = Condition)) +
geom_line(data = plotCBM, linewidth = 1) +
geom_point(data = plotOpPnts, size = 4) +
theme(legend.position = c(1, 0), legend.direction = "horizontal")
fittedPlot <- fittedPlot +
geom_line(data = plotCBM, mapping = aes(linetype = Condition), linewidth = 1) +
geom_point(data = plotOpPnts, mapping = aes(shape = Condition), size = 3) +
theme(legend.title=element_blank(),
legend.position = c(1, 0),
legend.direction = "horizontal",
legend.justification = c(1, 0),
legend.key.size = unit(1, "cm")) #+
# scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))
ciIndxX <- c(1, length(FPFX))
FPF <- FPFX[ciIndxX]
TPF <- TPFX[ciIndxX]
ciX <- binom.confint(x = FPF * K1, n = K1, methods = "exact")
ciY <- binom.confint(x = TPF * K2, n = K2, methods = "exact")
ciXUpper <- ciX$upper
ciXLower <- ciX$lower
ciYUpper <- ciY$upper
ciYLower <- ciY$lower
for (p in 1:length(FPF)){
ciX <- data.frame(FPF = c(ciXUpper[p], ciXLower[p]), TPF = c(TPF[p], TPF[p]))
ciY <- data.frame(FPF = c(FPF[p], FPF[p]), TPF = c(ciYUpper[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = ciY, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = ciX, aes(x = FPF, y = TPF), color = "black")
barRgt <- data.frame(FPF = c(ciXUpper[p], ciXUpper[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barLft <- data.frame(FPF = c(ciXLower[p], ciXLower[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barUp <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYUpper[p], ciYUpper[p]))
barBtm <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYLower[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = barRgt, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barLft, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barUp, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barBtm, aes(x = FPF, y = TPF), color = "black")
}
ciIndxY <- c(1, length(FPFY))
FPF <- FPFY[ciIndxY]
TPF <- TPFY[ciIndxY]
ciX <- binom.confint(x = FPF * K1, n = K1, methods = "exact")
ciY <- binom.confint(x = TPF * K2, n = K2, methods = "exact")
ciXUpper <- ciX$upper
ciXLower <- ciX$lower
ciYUpper <- ciY$upper
ciYLower <- ciY$lower
for (p in 1:length(FPF)){
ciX <- data.frame(FPF = c(ciXUpper[p], ciXLower[p]), TPF = c(TPF[p], TPF[p]))
ciY <- data.frame(FPF = c(FPF[p], FPF[p]), TPF = c(ciYUpper[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = ciY, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = ciX, aes(x = FPF, y = TPF), color = "black")
barRgt <- data.frame(FPF = c(ciXUpper[p], ciXUpper[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barLft <- data.frame(FPF = c(ciXLower[p], ciXLower[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barUp <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYUpper[p], ciYUpper[p]))
barBtm <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYLower[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = barRgt, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barLft, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barUp, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barBtm, aes(x = FPF, y = TPF), color = "black")
}
return(fittedPlot)
}
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 PlotCorCbmFit LLCorCBM nLLCorCBMNoTrns nLLCorCBM StatsCorCbm FitCorCbm
Documented in FitCorCbm
#' Fit CORCBM to a paired ROC dataset
#'
#' @description Fit the Correlated Contaminated Binormal Model (CORCBM)
#' to a paired ROC dataset.
#' The \strong{ROC} dataset has to be formatted as a
#' \strong{single modality}, \strong{two-reader} dataset, even though the actual
#' pairing may be different, see details.
#'
#'
#' @param dataset A \strong{paired ROC} dataset
#'
#'
#' @return A list containing three objects:
#' @return \item{fitCorCbmRet}{list(\code{FPCounts},\code{TPCounts},
#' \code{muX},\code{muY},\code{alphaX},\code{alphaY},\code{rhoNor},
#' \code{rhoAbn2},\code{zetaX},\code{zetaY},\code{covMat},\code{fixParam})}
#' @return \item{stats}{list(\code{aucX},\code{aucX},\code{stdAucX},
#' \code{stdAucY},\code{stdErr},\code{areaStat},\code{areaPval})}
#' @return \item{fittedPlot}{The fitted plot with operating points, error bars,
#' for both conditions}
#'
#' @details The conditions (X, Y) can be two readers interpreting images in the same
#' modality, the same reader interpreting images in different treatments, or
#' different readers interpreting images in 2 different treatments. Function
#' \code{\link{DfExtractCorCbmDataset}} can be used to construct a dataset suitable for
#' \code{FitCorCbm}. With reference to the returned values, and assuming R bins
#' in condition X and L bins in conditon Y,
#' \code{FPCounts} is the R x L matrix containing the counts for non-diseased cases,
#' \code{TPCounts} is the R x L matrix containing the counts for diseased cases;
#' \code{muX},\code{muY},\code{alphaX},\code{alphaY},\code{rhoNor},\code{rhoAbn2} are
#' the CORCBM parameters; \code{aucX},\code{aucX} are the AUCs in the two conditions;
#' \code{stdAucX},\code{stdAucY} are the corresponding standard errors;\code{stdErr}
#' contains the standard errors of the parameters of the model; \code{areaStat},
#' \code{areaPval},\code{covMat} are the area-statistic, the p-value and the covariance
#' matrix of the parameters. If a parameter approaches a limit, e.g., \code{rhoNor}
#' = 0.9999, it is held constant at near the limiting value and the covariance matrix
#' has one less dimension (along each edge) for each parameter that is held constant.
#' The indices of the parameters held fixed are in \code{fitCorCbmRet$fixParam}.
#'
## following examples generate excessive CPU time NOTES on devtools::check_win_xx()
## but it is instructive to execute them independently
## see also the publication referenced below
##
## dataset <- DfExtractCorCbmDataset(dataset05, trts = 1, rdrs = c(4,7))
## ret <- FitCorCbm(dataset)
## print(ret$fitCorCbmRet)
## print(ret$stats)
## print(ret$fittedPlot)
##
##
## ret <- FitCorCbm(datasetBinned123)
## print(ret$fitCorCbmRet)
## print(ret$stats)
## print(ret$fittedPlot)
## Also try two other datasets ending with 124 and 125
#'
#'
#' @references
#' Zhai X, Chakraborty DP (2017) A bivariate contaminated binormal model for robust fitting of proper ROC
#' curves to a pair of correlated, possibly degenerate, ROC datasets. Medical Physics. 44(6):2207--2222.
#'
#'
#' @importFrom bbmle mle2
#' @importFrom Rcpp cppFunction
#' @importFrom mvtnorm pmvnorm
#' @importFrom stats cor
#'
#' @export
#'
#'
FitCorCbm <- function(dataset){
if (dataset$descriptions$type != "ROC") {
stop("This program requires an ROC dataset")
}
minMu <- RJafrocEnv$minMu
maxMu <- RJafrocEnv$maxMu
minAlpha <- RJafrocEnv$minAlpha
maxAlpha <- RJafrocEnv$maxAlpha
minRho <- RJafrocEnv$minRho
maxRho <- RJafrocEnv$maxRho
# `as.matrix` is NOT absolutely necessary as `mean()` function is not used here
aucArray <- UtilFigureOfMerit(dataset, FOM = "Wilcoxon")
maxAUC <- max(aucArray)
while (pnorm(maxMu / sqrt(2)) <= maxAUC){
maxMu <- qnorm(maxAUC) * sqrt(2) + 0.5
}
I <- length(dataset$ratings$NL[,1,1,1])
J <- length(dataset$ratings$NL[1,,1,1])
K <- length(dataset$ratings$NL[1,1,,1])
K2 <- length(dataset$ratings$LL[1,1,,1])
K1 <- K - K2
FP <- dataset$ratings$NL[1,,1:K1,1]
TP <- dataset$ratings$LL[1,,,1]
z1X <- FP[1, ];z1Y <- FP[2, ]
z2X <- TP[1, ];z2Y <- TP[2, ]
rhoNor <- cor(z1X, z1Y)
rhoAbn <- cor(z2X, z2Y)
iniParamX <- FitCbmRoc(dataset, 1, 1) # arm X of pairing
muX <- iniParamX$mu
zetaX <- iniParamX$zetas
alphaX <- iniParamX$alpha
iniParamY <- FitCbmRoc(dataset, 1, 2) # arm Y of pairing
muY <- iniParamY$mu
zetaY <- iniParamY$zetas
alphaY <- iniParamY$alpha
muXFwd <- ForwardValue(muX, minMu, maxMu)
muYFwd <- ForwardValue(muY, minMu, maxMu)
alphaXFwd <- ForwardValue(alphaX, minAlpha, maxAlpha)
alphaYFwd <- ForwardValue(alphaY, minAlpha, maxAlpha)
rhoNorFwd <- ForwardValue(rhoNor, minRho, maxRho)
rhoAbn2Fwd <- ForwardValue(rhoAbn, minRho, maxRho)
zetaXFwd <- ForwardZetas(zetaX)
zetaYFwd <- ForwardZetas(zetaY)
parameters <- c(list(muXFwd, muYFwd, alphaXFwd, alphaYFwd, rhoNorFwd, rhoAbn2Fwd),
as.list(zetaXFwd), as.list(zetaYFwd))
namesVector <- c("muXFwd", "muYFwd", "alphaXFwd", "alphaYFwd","rhoNorFwd", "rhoAbn2Fwd")
namesVector <- c(namesVector, paste0("zetaXFwd", 1:length(zetaXFwd)))
namesVector <- c(namesVector, paste0("zetaYFwd", 1:length(zetaYFwd)))
names(parameters) <- namesVector
nBinsX <- length(unique(c(FP[1,], TP[1,])))
nBinsY <- length(unique(c(FP[2,], TP[2,])))
FPCounts <- array(dim = c(nBinsX, nBinsY))
TPCounts <- array(dim = c(nBinsX, nBinsY))
for (bX in 1:nBinsX){
for (bY in 1:nBinsY){
FPCounts[bX, bY] <- sum((z1X == bX & z1Y == bY))
TPCounts[bX, bY] <- sum((z2X == bX & z2Y == bY))
}
}
nLLCorCBMAdd <- AddZetaXFwd(nLLCorCBM, length(zetaX))
nLLCorCBMAdd <- AddZetaYFwd(nLLCorCBMAdd, length(zetaY))
ret <- mle2(nLLCorCBMAdd, start = parameters, eval.only = TRUE,
method = "BFGS", data = list(FPCounts = FPCounts, TPCounts = TPCounts, maxMu = maxMu))
NLLIni <- as.numeric(ret@min)
ret <- mle2(nLLCorCBMAdd, start = parameters,
method = "BFGS", data = list(FPCounts = FPCounts, TPCounts = TPCounts, maxMu = maxMu))
NLLFin <- as.numeric(ret@min)
allCoef1 <- ret@coef
muX <- InverseValue(allCoef1[1], minMu, maxMu)
muY <- InverseValue(allCoef1[2], minMu, maxMu)
alphaX <- InverseValue(allCoef1[3], minAlpha, maxAlpha)
alphaY <- InverseValue(allCoef1[4], minAlpha, maxAlpha)
rhoNor <- InverseValue(allCoef1[5], minRho, maxRho)
rhoAbn2 <- InverseValue(allCoef1[6], minRho, maxRho)
zetaX <- InverseZetas(allCoef1[7:(6 + length(zetaX))])
zetaY <- InverseZetas(allCoef1[(7 + length(zetaX)):length(allCoef1)])
CorCBMLLNoTrnsAdd <- AddZetaX(nLLCorCBMNoTrns, length(zetaX))
CorCBMLLNoTrnsAdd <- AddZetaY(CorCBMLLNoTrnsAdd, length(zetaY))
param <- c(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2, zetaX, zetaY)
fixParam <- c(which(param[1:2] > (maxMu - 0.01)),
2+which(param[3:4] > (maxAlpha - 0.01)),
4+which(param[5:6] > (maxRho - 0.01)))
if (length(fixParam) > 0){
# some parameters are fixed
fixList <- list()
afterIndex <- c()
for (p in fixParam){
if (p == 1){
muX <- maxMu - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(muXFwd = ForwardValue(muX, minMu, maxMu)))
} else if (p == 2){
muY <- maxMu - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(muYFwd = ForwardValue(muY, minMu, maxMu)))
}else if (p == 3){
alphaX <- maxAlpha - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(alphaXFwd = ForwardValue(alphaX, minAlpha, maxAlpha)))
}else if (p == 4){
alphaY <- maxAlpha - 0.01
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(alphaYFwd = ForwardValue(alphaY, minAlpha, maxAlpha)))
}else if (p == 5){
rhoNor <- ifelse(test = rhoNor > 0, yes = 0.99, no = -0.99)
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(rhoNorFwd = ForwardValue(rhoNor, minRho, maxRho)))
}else if (p == 6){
rhoAbn2 <- ifelse(test = rhoAbn2 > 0, yes = 0.99, no = -0.99)
afterIndex <- c(afterIndex, p - 1)
fixList <- c(fixList, list(rhoAbn2Fwd = ForwardValue(rhoAbn2, minRho, maxRho)))
}
}
muXFwd <- ForwardValue(muX, minMu, maxMu)
muYFwd <- ForwardValue(muY, minMu, maxMu)
alphaXFwd <- ForwardValue(alphaX, minAlpha, maxAlpha)
alphaYFwd <- ForwardValue(alphaY, minAlpha, maxAlpha)
rhoNorFwd <- ForwardValue(rhoNor, minRho, maxRho)
rhoAbn2Fwd <- ForwardValue(rhoAbn2, minRho, maxRho)
zetaXFwd <- ForwardZetas(zetaX)
zetaYFwd <- ForwardZetas(zetaY)
parameters <- c(list(muXFwd, muYFwd, alphaXFwd, alphaYFwd, rhoNorFwd, rhoAbn2Fwd),
as.list(zetaXFwd), as.list(zetaYFwd))
namesVector <- c("muXFwd", "muYFwd", "alphaXFwd", "alphaYFwd","rhoNorFwd", "rhoAbn2Fwd")
namesVector <- c(namesVector, paste0("zetaXFwd", 1:length(zetaXFwd)))
namesVector <- c(namesVector, paste0("zetaYFwd", 1:length(zetaYFwd)))
names(parameters) <- namesVector
ret2 <- mle2(nLLCorCBMAdd, start = parameters, method = "BFGS", fixed = fixList,
data = list(FPCounts = FPCounts, TPCounts = TPCounts, maxMu = maxMu), trace = TRUE)
NLLFin <- ret2@min
allCoef2 <- ret2@coef
for (i in 1:length(afterIndex)) {
allCoef2 <- append(allCoef2, fixList[[i]], after = afterIndex[i])
}
names(allCoef2) <- namesVector
muX <- InverseValue(allCoef2[1], minMu, maxMu)
muY <- InverseValue(allCoef2[2], minMu, maxMu)
alphaX <- InverseValue(allCoef2[3], minAlpha, maxAlpha)
alphaY <- InverseValue(allCoef2[4], minAlpha, maxAlpha)
rhoNor <- InverseValue(allCoef2[5], minRho, maxRho)
rhoAbn2 <- InverseValue(allCoef2[6], minRho, maxRho)
zetaX <- InverseZetas(allCoef2[7:(6 + length(zetaX))])
zetaY <- InverseZetas(allCoef2[(7 + length(zetaX)):length(allCoef2)])
parameters <- c(list(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2),
as.list(zetaX), as.list(zetaY))
namesVector <- c("muX", "muY", "alphaX", "alphaY","rhoNor", "rhoAbn2")
namesVector <- c(namesVector, paste0("zetaX", 1:length(zetaX)))
namesVector <- c(namesVector, paste0("zetaY", 1:length(zetaY)))
names(parameters) <- namesVector
ret3 <- mle2(CorCBMLLNoTrnsAdd, start = parameters, method = "BFGS", fixed = parameters[fixParam],
data = list(FPCounts = FPCounts, TPCounts = TPCounts))
} else {
# all parameters are used
allCoef2 <- allCoef1
names(allCoef2) <- namesVector
muX <- InverseValue(allCoef2[1], minMu, maxMu)
muY <- InverseValue(allCoef2[2], minMu, maxMu)
alphaX <- InverseValue(allCoef2[3], minAlpha, maxAlpha)
alphaY <- InverseValue(allCoef2[4], minAlpha, maxAlpha)
rhoNor <- InverseValue(allCoef2[5], minRho, maxRho)
rhoAbn2 <- InverseValue(allCoef2[6], minRho, maxRho)
zetaX <- InverseZetas(allCoef2[7:(6 + length(zetaX))])
zetaY <- InverseZetas(allCoef2[(7 + length(zetaX)):length(allCoef2)])
parameters <- c(list(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2),
as.list(zetaX), as.list(zetaY))
namesVector <- c("muX", "muY", "alphaX", "alphaY","rhoNor", "rhoAbn2")
namesVector <- c(namesVector, paste0("zetaX", 1:length(zetaX)))
namesVector <- c(namesVector, paste0("zetaY", 1:length(zetaY)))
names(parameters) <- namesVector
ret3 <- mle2(CorCBMLLNoTrnsAdd, start = parameters, method = "BFGS",
data = list(FPCounts = FPCounts, TPCounts = TPCounts))
}
covMat <- ret3@vcov
fittedPlot <- PlotCorCbmFit(
list(
FPCounts = FPCounts,
TPCounts = TPCounts,
muX = muX,
alphaX = alphaX,
muY = muY,
alphaY = alphaY,
rhoNor = rhoNor,
rhoAbn2 = rhoAbn2,
zetaX = zetaX,
zetaY = zetaY
)
)
fitCorCbmRet <- list(
FPCounts = FPCounts,
TPCounts = TPCounts,
muX = muX,
alphaX = alphaX,
muY = muY,
alphaY = alphaY,
rhoNor = rhoNor,
rhoAbn2 = rhoAbn2,
zetaX = zetaX,
zetaY = zetaY,
covMat = covMat,
fixParam = fixParam
)
stats <- StatsCorCbm(fitCorCbmRet)
return(list(
fitCorCbmRet = fitCorCbmRet,
stats = stats,
fittedPlot = fittedPlot
))
}
StatsCorCbm <- function(fitCorCbmRet) {
muX <- fitCorCbmRet$muX
muY <- fitCorCbmRet$muY
alphaX <- fitCorCbmRet$alphaX
alphaY <- fitCorCbmRet$alphaY
rhoNor <- fitCorCbmRet$rhoNor
rhoAbn2 <- fitCorCbmRet$rhoAbn2
zetaX <- fitCorCbmRet$zetaX
zetaY <- fitCorCbmRet$zetaY
FPCounts <- fitCorCbmRet$FPCounts
TPCounts <- fitCorCbmRet$TPCounts
covMat <- fitCorCbmRet$covMat
fixParam <- fitCorCbmRet$fixParam
vars <- diag(covMat)
stdErr <- sqrt(vars)
dMuX <- alphaX * dnorm(muX / sqrt(2)) / sqrt(2)
dAlphaX <- -1 / 2 + pnorm(muX / sqrt(2))
stdAucX <-
sqrt(vars[1] * dMuX ^ 2 + vars[3] * dAlphaX ^ 2 + dMuX * dAlphaX * covMat[1, 3])
dMuY <- alphaY * dnorm(muY / sqrt(2)) / sqrt(2)
dAlphaY <- -1 / 2 + pnorm(muY / sqrt(2))
stdAucY <-
sqrt(vars[2] * dMuY ^ 2 + vars[4] * dAlphaY ^ 2 + dMuY * dAlphaY * covMat[2, 4])
stdErr <- c(stdErr, stdAucX, stdAucY)
if (length(fixParam) > 0) {
for (p in fixParam) {
stdErr <- append(stdErr, NA, after = p - 1)
if (p == 3) {
alphaX <- 1
} else if (p == 4) {
alphaY <- 1
} else if (p == 5) {
rhoNor <- 1
} else if (p == 6) {
rhoAbn2 <- 1
}
}
stdErr <- as.numeric(stdErr)
}
stdErr[2:3] <- stdErr[3:2] # switch the position of alphaX and muY
aucX <- (1 - alphaX) * 0.5 + alphaX * pnorm(muX / sqrt(2))
aucY <- (1 - alphaY) * 0.5 + alphaY * pnorm(muY / sqrt(2))
aucDiff <- aucX - aucY
varDiff <- 0
derivs <- c(dMuX,-dMuY, dAlphaX,-dAlphaY)
for (k in 1:4) {
for (l in 1:4) {
varDiff <- derivs[l] * derivs[k] * covMat[k, l] + varDiff
}
}
stdDiff <- sqrt(varDiff)
areaStat <- abs(aucDiff) / stdDiff
areaPval <- 1 - pnorm(areaStat)
return(
list(
aucX = aucX,
aucY = aucY,
stdAucX = as.numeric(stdAucX),
stdAucY = as.numeric(stdAucY),
stdErr = stdErr,
areaStat = areaStat,
areaPval = areaPval
)
)
}
###############################################################################
nLLCorCBM <- function(muXFwd, muYFwd, alphaXFwd, alphaYFwd, rhoNorFwd, rhoAbn2Fwd, FPCounts, TPCounts, maxMu){
minMu <- RJafrocEnv$minMu
#maxMu <- RJafrocEnv$maxMu
minAlpha <- RJafrocEnv$minAlpha
maxAlpha <- RJafrocEnv$maxAlpha
minRho <- RJafrocEnv$minRho
maxRho <- RJafrocEnv$maxRho
allParameters <- names(formals())
zetaXPos <- regexpr("zetaX", allParameters)
zetaXFwd <- unlist(mget(allParameters[which(zetaXPos == 1)]))
zetaYPos <- regexpr("zetaY", allParameters)
zetaYFwd <- unlist(mget(allParameters[which(zetaYPos == 1)]))
muX <- InverseValue(muXFwd, minMu, maxMu)
muY <- InverseValue(muYFwd, minMu, maxMu)
alphaX <- InverseValue(alphaXFwd, minAlpha, maxAlpha)
alphaY <- InverseValue(alphaYFwd, minAlpha, maxAlpha)
rhoNor <- InverseValue(rhoNorFwd, minRho, maxRho)
rhoAbn <- InverseValue(rhoAbn2Fwd, minRho, maxRho)
zetaX <- InverseZetas(zetaXFwd)
zetaY <- InverseZetas(zetaYFwd)
LL <- LLCorCBM (muX, muY, alphaX, alphaY, rhoNor, rhoAbn, zetaX, zetaY, FPCounts, TPCounts)
return(-LL)
}
###############################################################################
nLLCorCBMNoTrns <- function(muX, muY, alphaX, alphaY, rhoNor, rhoAbn2, FPCounts, TPCounts){
allParameters <- names(formals())
zetaXPos <- regexpr("zetaX", allParameters)
zetaX <- unlist(mget(allParameters[which(zetaXPos == 1)]))
zetaYPos <- regexpr("zetaY", allParameters)
zetaY <- unlist(mget(allParameters[which(zetaYPos == 1)]))
LL <- LLCorCBM (muX, muY, alphaX, alphaY, rhoNor, rhoAbn2, zetaX, zetaY, FPCounts, TPCounts)
return(-LL)
}
###############################################################################
LLCorCBM <- function(muX, muY, alphaX, alphaY, rhoNor, rhoAbn, zetaX, zetaY, FPCounts, TPCounts) {
if ((abs(zetaX[1]) > 2) ||
(abs(zetaY[1]) > 2) ||
(zetaX[length(zetaX)] > 10) ||
(zetaY[length(zetaY)] > 10)) {
return(-1e10)
}
zetaX <- c(-Inf, zetaX, Inf)
zetaY <- c(-Inf, zetaY, Inf)
if (is.unsorted(zetaX)) return(-1e10)
if (is.unsorted(zetaY)) return(-1e10)
rhoAbn1 <- rhoNor
rhoAbn2 <- rhoAbn
rhoAbn12 <- mean(c(rhoAbn1, rhoAbn2))
rhoAbn4 <- rhoAbn12
sigmaNor <- rbind(c(1, rhoNor), c(rhoNor, 1))
sigmaAbn1 <- rbind(c(1, rhoAbn1), c(rhoAbn1, 1))
sigmaAbn2 <- rbind(c(1, rhoAbn2), c(rhoAbn2, 1))
sigmaAbn3 <- rbind(c(1, rhoAbn12), c(rhoAbn12, 1))
sigmaAbn4 <- rbind(c(1, rhoAbn4), c(rhoAbn4, 1))
LLNor <- 0
LLAbn <- 0
nBinsX <- dim(FPCounts)[1]
nBinsY <- dim(FPCounts)[2]
for (bX in 1:nBinsX){
for (bY in 1:nBinsY){
if (FPCounts[bX, bY] > 0) {
tempNor <- pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]), sigma = sigmaNor)
if (tempNor >= 0) {
LLNor <- LLNor +
FPCounts[bX, bY] * log(tempNor)
} else {
return(-1e10)
}
}
if (TPCounts[bX, bY] > 0) {
tempAbn <- (1 - alphaX) * (1 - alphaY) *
pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]), mean = c(0, 0), sigma = sigmaAbn1) +
(alphaX) * (1 - alphaY) * pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]),
mean = c(muX, 0), sigma = sigmaAbn3) +
(1 - alphaX) * (alphaY) * pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]),
mean = c(0, muY), sigma = sigmaAbn4) +
alphaX * alphaY * pmvnorm(c(zetaX[bX], zetaY[bY]), c(zetaX[bX + 1], zetaY[bY + 1]),
mean = c(muX, muY), sigma = sigmaAbn2)
if (tempAbn >= 0) {
LLAbn <- LLAbn +
TPCounts[bX, bY] * log(tempAbn)
} else {
return(-1e10)
}
}
}
}
LL <- LLNor + LLAbn
return(LL)
}
#' @import ggplot2
PlotCorCbmFit <- function(retFitCorCBM){
muX <- retFitCorCBM$muX
muY <- retFitCorCBM$muY
alphaX <- retFitCorCBM$alphaX
alphaY <- retFitCorCBM$alphaY
FPCounts <- retFitCorCBM$FPCounts
TPCounts <- retFitCorCBM$TPCounts
K1 <- sum(FPCounts);K2 <- sum(TPCounts)
plotZeta <- seq(-3, max(muX,muY)+2, by = 0.1)
plotCBM <- NULL
plotOpPnts <- NULL
FPFX <- 1 - pnorm(plotZeta)
TPFX <- (1 - alphaX) * (1 - pnorm(plotZeta)) + alphaX * (1 - pnorm(plotZeta, mean = muX))
plotCBM <- rbind(plotCBM, data.frame(FPF = FPFX, TPF = TPFX, Condition = "X"))
FPFX <- cumsum(rev(rowSums(FPCounts)))/K1
TPFX <- cumsum(rev(rowSums(TPCounts)))/K2
FPFX <- FPFX[-length(FPFX)]
TPFX <- TPFX[-length(TPFX)]
plotOpPnts <- rbind(plotOpPnts, data.frame(FPF = FPFX, TPF = TPFX, Condition = "X"))
FPFY <- 1 - pnorm(plotZeta)
TPFY <- (1 - alphaY) * (1 - pnorm(plotZeta)) + alphaY * (1 - pnorm(plotZeta, mean = muY))
plotCBM <- rbind(plotCBM, data.frame(FPF = FPFY, TPF = TPFY, Condition = "Y"))
FPFY <- cumsum(rev(colSums(FPCounts)))/K1
TPFY <- cumsum(rev(colSums(TPCounts)))/K2
FPFY <- FPFY[-length(FPFY)]
TPFY <- TPFY[-length(TPFY)]
plotOpPnts <- rbind(plotOpPnts, data.frame(FPF = FPFY, TPF = TPFY, Condition = "Y"))
Condition <- NULL;FPF <- NULL;TPF <- NULL
fittedPlot <- ggplot(data = plotCBM, mapping = aes(x = FPF, y = TPF, color = Condition)) +
geom_line(data = plotCBM, linewidth = 1) +
geom_point(data = plotOpPnts, size = 4) +
theme(legend.position = c(1, 0), legend.direction = "horizontal")
fittedPlot <- fittedPlot +
geom_line(data = plotCBM, mapping = aes(linetype = Condition), linewidth = 1) +
geom_point(data = plotOpPnts, mapping = aes(shape = Condition), size = 3) +
theme(legend.title=element_blank(),
legend.position = c(1, 0),
legend.direction = "horizontal",
legend.justification = c(1, 0),
legend.key.size = unit(1, "cm")) #+
# scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))
ciIndxX <- c(1, length(FPFX))
FPF <- FPFX[ciIndxX]
TPF <- TPFX[ciIndxX]
ciX <- binom.confint(x = FPF * K1, n = K1, methods = "exact")
ciY <- binom.confint(x = TPF * K2, n = K2, methods = "exact")
ciXUpper <- ciX$upper
ciXLower <- ciX$lower
ciYUpper <- ciY$upper
ciYLower <- ciY$lower
for (p in 1:length(FPF)){
ciX <- data.frame(FPF = c(ciXUpper[p], ciXLower[p]), TPF = c(TPF[p], TPF[p]))
ciY <- data.frame(FPF = c(FPF[p], FPF[p]), TPF = c(ciYUpper[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = ciY, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = ciX, aes(x = FPF, y = TPF), color = "black")
barRgt <- data.frame(FPF = c(ciXUpper[p], ciXUpper[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barLft <- data.frame(FPF = c(ciXLower[p], ciXLower[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barUp <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYUpper[p], ciYUpper[p]))
barBtm <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYLower[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = barRgt, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barLft, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barUp, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barBtm, aes(x = FPF, y = TPF), color = "black")
}
ciIndxY <- c(1, length(FPFY))
FPF <- FPFY[ciIndxY]
TPF <- TPFY[ciIndxY]
ciX <- binom.confint(x = FPF * K1, n = K1, methods = "exact")
ciY <- binom.confint(x = TPF * K2, n = K2, methods = "exact")
ciXUpper <- ciX$upper
ciXLower <- ciX$lower
ciYUpper <- ciY$upper
ciYLower <- ciY$lower
for (p in 1:length(FPF)){
ciX <- data.frame(FPF = c(ciXUpper[p], ciXLower[p]), TPF = c(TPF[p], TPF[p]))
ciY <- data.frame(FPF = c(FPF[p], FPF[p]), TPF = c(ciYUpper[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = ciY, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = ciX, aes(x = FPF, y = TPF), color = "black")
barRgt <- data.frame(FPF = c(ciXUpper[p], ciXUpper[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barLft <- data.frame(FPF = c(ciXLower[p], ciXLower[p]), TPF = c(TPF[p] - 0.01, TPF[p] + 0.01))
barUp <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYUpper[p], ciYUpper[p]))
barBtm <- data.frame(FPF = c(FPF[p] - 0.01, FPF[p] + 0.01), TPF = c(ciYLower[p], ciYLower[p]))
fittedPlot <- fittedPlot + geom_line(data = barRgt, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barLft, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barUp, aes(x = FPF, y = TPF), color = "black") +
geom_line(data = barBtm, aes(x = FPF, y = TPF), color = "black")
}
return(fittedPlot)
}
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.