#' Performs DBM or OR significance testing for factorial dataset Old Code
#'
#' @description (DBM) or Obuchowski-Rockette (OR) significance testing, for
#' specified dataset.
#'
#' @param dataset The dataset to be analyzed.
#'
#' @param FOM The figure of merit
#'
#' @param FPFValue Only needed for \code{LROC} data.
#'
#' @param alpha The significance level of the test; the default is 0.05
#'
#' @param method The testing method to be used: \code{"DBM"} (the default),
#' representing the Dorfman-Berbaum-Metz method or \code{"OR"}, representing
#' the Obuchowski-Rockette method.
#'
#' @param covEstMethod The covariance matrix estimation method
#'
#' @param nBoots The number of bootstraps (defaults to 200)
#'
#' @param analysisOption Determines which factors are regarded as random vs.
#' fixed:
#'
#' @keywords internal
#'
#' @export
#
# Code from version 0.0.1 of RJafroc (see RJafroc_0.0.1.tar)
# with corrections that fixes the Lucy McGowan found bug in version 0.0.1.
# This is intended to check the current code for errors that might creep in as
# I attempt to improve the organization of the code and the output.
#
StOldCode <- function(dataset, FOM, FPFValue = 0.2, alpha = 0.05, method = "DBM",
covEstMethod = "jackknife", nBoots = 200, analysisOption = "ALL")
{
if (dataset$descriptions$design != "FCTRL") stop("old code requires FCTRL study design")
if (method == "DBM"){
return(DBMHAnalysis(dataset, FOM, alpha, analysisOption)) # original code: StOldCode.R
} else if (method == "OR"){
return(ORHAnalysis(dataset, FOM, alpha, covEstMethod, nBoots, analysisOption)) # original code: StOldCode.R
} else {
errMsg <- sprintf("analysis method must be `DBM` or `OR`", method)
stop(errMsg)
}
}
# Old Code from 0.0.1 (the first version on CRAN)
# One error in p-value (noted by Lucy D'Agostino McGowan) was corrected prior to 7/27/2023
# Lots of other similar errors were corrected on 7/27/2023
# Example of error and correction:
################################################################################
# tPr[i] <- 2 * pt(tStat[i], df) 7/27/2023
# Corrected on 7/27/2023 as follows:
# tPr[i] <- 2 * pt(abs(tStat[i]), df, lower.tail = FALSE)
################################################################################
# Otherwise the following codes (DBMHAnalysis and ORHAnalysis) is identical to XZ code
DBMHAnalysis <- function(dataset, FOM, alpha, analysisOption)
{
NL <- dataset$ratings$NL
LL <- dataset$ratings$LL
lesionVector <- dataset$lesions$perCase
lesionID <- dataset$lesions$IDs
lesionWeight <- dataset$lesions$weights
maxNL <- dim(NL)[4]
dataType <- dataset$descriptions$type
modalityID <- dataset$descriptions$modalityID
readerID <- dataset$descriptions$readerID
I <- length(modalityID)
J <- length(readerID)
K <- dim(NL)[3]
K2 <- dim(LL)[3]
K1 <- K - K2
if (!analysisOption %in% c("RRRC", "FRRC", "RRFC", "ALL")){
errMsg <- sprintf("%s is not an available analysisOption.", analysisOption)
stop(errMsg)
}
if (I < 2) {
stop("The analysis requires at least 2 modalities.")
}
fomArray <- UtilFigureOfMerit(dataset, FOM)
trMeans <- rowMeans(fomArray)
if (FOM %in% c("MaxNLF", "HrSp")) {
jkFOMArray <- array(dim = c(I, J, K1))
pseudoValues <- array(dim = c(I, J, K1))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K1) {
nl <- NL[i, j, -k, ]
ll <- LL[i, j, , ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector,
lesionID, lesionWeight, maxNL, FOM)
pseudoValues[i, j, k] <- fomArray[i, j] * K1 - jkFOMArray[i, j, k] * (K1 - 1)
}
pseudoValues[i, j, ] <- pseudoValues[i, j, ] + (fomArray[i, j] - mean(pseudoValues[i, j, ]))
}
}
} else if (FOM %in% c("MaxLLF", "HrSe")) {
jkFOMArray <- array(dim = c(I, J, K2))
pseudoValues <- array(dim = c(I, J, K2))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K2) {
nl <- NL[i, j, -(k + K1), ]
ll <- LL[i, j, -k, ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2 - 1, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector[-k],
lesionID[-k, ], lesionWeight[-k, ], maxNL, FOM)
pseudoValues[i, j, k] <- fomArray[i, j] * K2 - jkFOMArray[i, j, k] * (K2 - 1)
}
pseudoValues[i, j, ] <- pseudoValues[i, j, ] + (fomArray[i, j] - mean(pseudoValues[i, j, ]))
}
}
} else {
jkFOMArray <- array(dim = c(I, J, K))
pseudoValues <- array(dim = c(I, J, K))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K) {
if (k <= K1) {
nl <- NL[i, j, -k, ]
ll <- LL[i, j, , ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector,
lesionID, lesionWeight, maxNL, FOM)
} else {
nl <- NL[i, j, -k, ]
ll <- LL[i, j, -(k - K1), ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2 - 1, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector[-(k - K1)],
lesionID[-(k - K1), ], lesionWeight[-(k - K1), ], maxNL, FOM)
}
pseudoValues[i, j, k] <- fomArray[i, j] * K - jkFOMArray[i, j, k] * (K - 1)
}
pseudoValues[i, j, ] <- pseudoValues[i, j, ] + (fomArray[i, j] - mean(pseudoValues[i, j, ]))
}
}
}
K <- length(pseudoValues[1, 1, ])
msT <- 0
for (i in 1:I) {
msT <- msT + (mean(pseudoValues[i, , ]) - mean(pseudoValues))^2
}
msT <- msT * K * J/(I - 1)
msR <- 0
for (j in 1:J) {
msR <- msR + (mean(pseudoValues[, j, ]) - mean(pseudoValues))^2
}
msR <- msR * K * I/(J - 1)
msC <- 0
for (k in 1:K) {
msC <- msC + (mean(pseudoValues[, , k]) - mean(pseudoValues))^2
}
msC <- msC * I * J/(K - 1)
msTR <- 0
for (i in 1:I) {
for (j in 1:J) {
msTR <- msTR + (mean(pseudoValues[i, j, ]) - mean(pseudoValues[i, , ]) -
mean(pseudoValues[, j, ]) + mean(pseudoValues))^2
}
}
msTR <- msTR * K/((I - 1) * (J - 1))
msTC <- 0
for (i in 1:I) {
for (k in 1:K) {
msTC <- msTC + (mean(pseudoValues[i, , k]) - mean(pseudoValues[i, , ]) -
mean(pseudoValues[, , k]) + mean(pseudoValues))^2
}
}
msTC <- msTC * J/((I - 1) * (K - 1))
msRC <- 0
for (j in 1:J) {
for (k in 1:K) {
msRC <- msRC + (mean(pseudoValues[, j, k]) - mean(pseudoValues[, j, ]) -
mean(pseudoValues[, , k]) + mean(pseudoValues))^2
}
}
msRC <- msRC * I/((J - 1) * (K - 1))
msTRC <- 0
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K) {
msTRC <- msTRC + (pseudoValues[i, j, k] - mean(pseudoValues[i, j, ]) -
mean(pseudoValues[i, , k]) - mean(pseudoValues[, j, k]) +
mean(pseudoValues[i, , ]) + mean(pseudoValues[, j, ]) +
mean(pseudoValues[, , k]) - mean(pseudoValues))^2
}
}
}
msTRC <- msTRC/((I - 1) * (J - 1) * (K - 1))
varR <- (msR - msTR - msRC + msTRC)/(I * K)
varC <- (msC - msTC - msRC + msTRC)/(I * J)
varTR <- (msTR - msTRC)/K
varTC <- (msTC - msTRC)/J
varRC <- (msRC - msTRC)/I
varErr <- msTRC
varComp <- c(varR, varC, varTR, varTC, varRC, varErr)
varCompName <- c("Var(R)", "Var(C)", "Var(T*R)", "Var(T*C)", "Var(R*C)", "Var(Error)")
varComp <- data.frame(varComp, row.names = varCompName)
msArray <- c(msT, msR, msC, msTR, msTC, msRC, msTRC)
dfArray <- c(I - 1, J - 1, K - 1, (I - 1) * (J - 1), (I - 1) * (K - 1), (J - 1) * (K - 1), (I - 1) * (J - 1) * (K - 1))
ssArray <- msArray * dfArray
msArray <- c(msArray, NA)
dfArray <- c(dfArray, sum(dfArray))
ssArray <- c(ssArray, sum(ssArray))
sourceArray <- c("T", "R", "C", "TR", "TC", "RC", "TRC", "Total")
anovaY <- data.frame(Source = sourceArray, SS = ssArray, DF = dfArray, MS = msArray)
msRSingle <- array(0, dim = c(I))
msCSingle <- array(0, dim = c(I))
msRCSingle <- array(0, dim = c(I))
for (i in 1:I) {
for (j in 1:J) {
msRSingle[i] <- msRSingle[i] + (mean(pseudoValues[i, j, ]) - mean(pseudoValues[i, , ]))^2
}
msRSingle[i] <- msRSingle[i] * K/(J - 1)
for (k in 1:K) {
msCSingle[i] <- msCSingle[i] + (mean(pseudoValues[i, , k]) - mean(pseudoValues[i, , ]))^2
}
msCSingle[i] <- msCSingle[i] * J/(K - 1)
for (j in 1:J) {
for (k in 1:K) {
msRCSingle[i] <- msRCSingle[i] +
(mean(pseudoValues[i, j, k]) -
mean(pseudoValues[i, j, ]) - mean(pseudoValues[i, , k]) + mean(pseudoValues[i, , ]))^2
}
}
msRCSingle[i] <- msRCSingle[i]/((J - 1) * (K - 1))
}
sourceArraySingle <- c("R", "C", "RC")
dfArraySingle <- c(J - 1, K - 1, (J - 1) * (K - 1))
msArraySingle <- t(cbind(msRSingle, msCSingle, msRCSingle))
msSingleTable <- data.frame(sourceArraySingle, dfArraySingle, msArraySingle, row.names = NULL)
colnames(msSingleTable) <- c("Source", "DF", modalityID)
diffTRMeans <- array(dim = choose(I, 2))
diffTRName <- array(dim = choose(I, 2))
ii <- 1
for (i in 1:I) {
if (i == I)
break
for (ip in (i + 1):I) {
diffTRMeans[ii] <- trMeans[i] - trMeans[ip]
diffTRName[ii] <- paste(modalityID[i], modalityID[ip], sep = "-")
ii <- ii + 1
}
}
msNum <- msT
if (analysisOption %in% c("RRRC", "ALL")) {
# ************ RRRC ****************
if (J > 1) {
msDenRRRC <- msTR + max(msTC - msTRC, 0)
fRRRC <- msNum/msDenRRRC
ddfRRRC <- msDenRRRC^2/(msTR^2/((I - 1) * (J - 1)))
pRRRC <- 1 - pf(fRRRC, I - 1, ddfRRRC)
stdErrRRRC <- sqrt(2 * msDenRRRC/J/K)
tStat <- vector()
tPr <- vector()
CIRRRC <- array(dim = c(length(diffTRMeans), 2))
for (i in 1:length(diffTRMeans)) {
tStat[i] <- diffTRMeans[i]/stdErrRRRC
################################################################################
# tPr[i] <- 2 * pt(tStat[i], ddfRRRC) 7/27/2023
# corrected 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), ddfRRRC, lower.tail = FALSE)
################################################################################
CIRRRC[i, ] <- sort(c(diffTRMeans[i] - qt(alpha/2, ddfRRRC) * stdErrRRRC,
diffTRMeans[i] + qt(alpha/2, ddfRRRC) * stdErrRRRC))
}
ciDiffTrtRRRC <- data.frame(Treatment = diffTRName,
Estimate = diffTRMeans,
StdErr = rep(stdErrRRRC, choose(I, 2)),
DF = rep(ddfRRRC, choose(I, 2)),
t = tStat,
p = tPr,
CI = CIRRRC)
colnames(ciDiffTrtRRRC) <- c("Treatment",
"Estimate",
"StdErr",
"DF",
"t",
"Pr > t",
"CI Lower",
"CI Upper")
dfSingleRRRC <- array(dim = I)
msDenSingleRRRC <- array(dim = I)
stdErrSingleRRRC <- array(dim = I)
CISingleRRRC <- array(dim = c(I, 2))
for (i in 1:I) {
msDenSingleRRRC[i] <- msRSingle[i] + max(msCSingle[i] - msRCSingle[i], 0)
dfSingleRRRC[i] <- msDenSingleRRRC[i]^2/msRSingle[i]^2 * (J - 1)
stdErrSingleRRRC[i] <- sqrt(msDenSingleRRRC[i]/J/K)
CISingleRRRC[i, ] <- sort(c(trMeans[i] - qt(alpha/2,
dfSingleRRRC[i]) * stdErrSingleRRRC[i],
trMeans[i] + qt(alpha/2,
dfSingleRRRC[i]) * stdErrSingleRRRC[i]))
}
ciAvgRdrEachTrtRRRC <- data.frame(Treatment = modalityID,
Area = trMeans,
StdErr = stdErrSingleRRRC,
DF = dfSingleRRRC,
CI = CISingleRRRC,
row.names = NULL)
colnames(ciAvgRdrEachTrtRRRC) <- c("Treatment",
"Area",
"StdErr",
"DF",
"CI Lower",
"CI Upper")
} else {
fRRRC <- NA
ddfRRRC <- NA
pRRRC <- NA
ciDiffTrtRRRC <- NA
ciAvgRdrEachTrtRRRC <- NA
}
if (analysisOption == "RRRC")
return(list(fomArray = as.data.frame(fomArray),
anovaY = anovaY,
anovaYi = msSingleTable,
varComp = varComp,
fRRRC = fRRRC,
ddfRRRC = ddfRRRC,
pRRRC = pRRRC,
ciDiffTrtRRRC = ciDiffTrtRRRC,
ciAvgRdrEachTrtRRRC = ciAvgRdrEachTrtRRRC))
}
if (analysisOption %in% c("FRRC", "ALL")) {
# ************ FRRC ****************
msDenFRRC <- msTC
fFRRC <- msNum/msDenFRRC
ddfFRRC <- (I - 1) * (K - 1)
pFRRC <- 1 - pf(fFRRC, I - 1, ddfFRRC)
stdErrFRRC <- sqrt(2 * msDenFRRC/J/K)
tStat <- vector()
tPr <- vector()
CIFRRC <- array(dim = c(length(diffTRMeans), 2))
for (i in 1:length(diffTRMeans)) {
tStat[i] <- diffTRMeans[i]/stdErrFRRC
################################################################################
# tPr[i] <- 2 * pt(tStat[i], ddfFRRC) # 7/27/2023
# corrected 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), ddfFRRC, lower.tail = FALSE)
################################################################################
CIFRRC[i, ] <- sort(c(diffTRMeans[i] - qt(alpha/2, ddfFRRC) * stdErrFRRC,
diffTRMeans[i] + qt(alpha/2, ddfFRRC) * stdErrFRRC))
}
ciDiffTrtFRRC <- data.frame(Treatment = diffTRName,
Estimate = diffTRMeans,
StdErr = rep(stdErrFRRC, choose(I, 2)),
DF = rep(ddfFRRC, choose(I, 2)),
t = tStat,
p = tPr,
CI = CIFRRC)
colnames(ciDiffTrtFRRC) <- c("Treatment",
"Estimate",
"StdErr",
"DF",
"t",
"Pr > t",
"CI Lower",
"CI Upper")
dfSingleFRRC <- array(dim = I)
msDenSingleFRRC <- array(dim = I)
stdErrSingleFRRC <- array(dim = I)
CISingleFRRC <- array(dim = c(I, 2))
for (i in 1:I) {
msDenSingleFRRC[i] <- msCSingle[i]
dfSingleFRRC[i] <- (K - 1)
stdErrSingleFRRC[i] <- sqrt(msDenSingleFRRC[i]/J/K)
CISingleFRRC[i, ] <- sort(c(trMeans[i] - qt(alpha/2,
dfSingleFRRC[i]) * stdErrSingleFRRC[i],
trMeans[i] + qt(alpha/2, dfSingleFRRC[i]) * stdErrSingleFRRC[i]))
}
ciAvgRdrEachTrtFRRC <- data.frame(Treatment = modalityID, Area = trMeans, StdErr = stdErrSingleFRRC, DF = dfSingleFRRC, CI = CISingleFRRC, row.names = NULL)
colnames(ciAvgRdrEachTrtFRRC) <- c("Treatment",
"Area",
"StdErr",
"DF",
"CI Lower",
"CI Upper")
ssTFRRC <- array(0, dim = c(J))
ssCFRRC <- array(0, dim = c(J))
ssTCFRRC <- array(0, dim = c(J))
for (j in 1:J) {
for (i in 1:I) {
ssTFRRC[j] <- ssTFRRC[j] + (mean(pseudoValues[i, j, ]) - mean(pseudoValues[, j, ]))^2
}
ssTFRRC[j] <- ssTFRRC[j] * K
for (k in 1:K) {
ssCFRRC[j] <- ssCFRRC[j] + (mean(pseudoValues[, j, k]) - mean(pseudoValues[, j, ]))^2
}
ssCFRRC[j] <- ssCFRRC[j] * I
for (i in 1:I) {
for (k in 1:K) {
ssTCFRRC[j] <- ssTCFRRC[j] + (mean(pseudoValues[i, j, k]) -
mean(pseudoValues[i, j, ]) -
mean(pseudoValues[, j, k]) +
mean(pseudoValues[, j, ]))^2
}
}
}
sourceArrayFRRC <- c("T", "C", "TC")
dfArrayFRRC <- c(I - 1, K - 1, (I - 1) * (K - 1))
ssArrayFRRC <- t(cbind(ssTFRRC, ssCFRRC, ssTCFRRC))
ssTableFRRC <- data.frame(sourceArrayFRRC,
dfArrayFRRC,
ssArrayFRRC,
row.names = NULL)
colnames(ssTableFRRC) <- c("Source", "DF", readerID)
msArrayFRRC <- ssArrayFRRC
for (i in 1:3) msArrayFRRC[i, ] <- ssArrayFRRC[i, ]/dfArrayFRRC[i]
msTableFRRC <- data.frame(sourceArrayFRRC,
dfArrayFRRC,
msArrayFRRC,
row.names = NULL)
colnames(msTableFRRC) <- c("Source", "DF", readerID)
diffTRMeansFRRC <- array(dim = c(J, choose(I, 2)))
for (j in 1:J) {
ii <- 1
for (i in 1:I) {
if (i == I)
break
for (ip in (i + 1):I) {
diffTRMeansFRRC[j, ii] <- fomArray[i, j] - fomArray[ip, j]
ii <- ii + 1
}
}
}
diffTRMeansFRRC <- as.vector(t(diffTRMeansFRRC))
stdErrFRRC <- sqrt(2 * msArrayFRRC[3, ]/K)
stdErrFRRC <- rep(stdErrFRRC, choose(I, 2))
dim(stdErrFRRC) <- c(J, choose(I, 2))
stdErrFRRC <- as.vector(t(stdErrFRRC))
readerNames <- rep(readerID, choose(I, 2))
dim(readerNames) <- c(J, choose(I, 2))
readerNames <- as.vector(t(readerNames))
trNames <- rep(diffTRName, J)
dfReaderFRRC <- rep(K - 1, length(stdErrFRRC))
CIReaderFRRC <- array(dim = c(length(stdErrFRRC), 2))
tStat <- vector()
tPr <- vector()
for (i in 1:length(stdErrFRRC)) {
tStat[i] <- diffTRMeansFRRC[i]/stdErrFRRC[i]
################################################################################
# tPr[i] <- 2 * pt(tStat[i], dfReaderFRRC[i]) # 7/27/2023
# corrected 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), dfReaderFRRC[i], lower.tail = FALSE)
################################################################################
CIReaderFRRC[i, ] <- sort(c(diffTRMeansFRRC[i] - qt(alpha/2, dfReaderFRRC[i]) * stdErrFRRC[i],
diffTRMeansFRRC[i] + qt(alpha/2, dfReaderFRRC[i]) * stdErrFRRC[i]))
}
ciDiffTrtEachRdr <- data.frame(Reader = readerNames,
Treatment = trNames,
Estimate = diffTRMeansFRRC,
StdErr = stdErrFRRC,
DF = dfReaderFRRC,
t = tStat,
p = tPr,
CI = CIReaderFRRC)
colnames(ciDiffTrtEachRdr) <- c("Reader",
"Treatment",
"Estimate",
"StdErr",
"DF",
"t",
"Pr > t",
"CI Lower",
"CI Upper")
if (analysisOption == "FRRC")
return(list(fomArray = as.data.frame(fomArray),
anovaY = anovaY,
anovaYi = msSingleTable,
varComp = varComp,
fFRRC = fFRRC,
ddfFRRC = ddfFRRC,
pFRRC = pFRRC,
ciDiffTrtFRRC = ciDiffTrtFRRC,
ciAvgRdrEachTrtFRRC = ciAvgRdrEachTrtFRRC,
ssAnovaEachRdr = ssTableFRRC,
msAnovaEachRdr = msTableFRRC,
ciDiffTrtEachRdr = ciDiffTrtEachRdr))
}
if (analysisOption %in% c("RRFC", "ALL")) {
# ************ RRFC ****************
if (J > 1) {
msDenRRFC <- msTR
fRRFC <- msNum/msDenRRFC
ddfRRFC <- ((I - 1) * (J - 1))
pRRFC <- 1 - pf(fRRFC, I - 1, ddfRRFC)
stdErrRRFC <- sqrt(2 * msDenRRFC/J/K)
tStat <- vector()
tPr <- vector()
CIRRFC <- array(dim = c(length(diffTRMeans), 2))
for (i in 1:length(diffTRMeans)) {
tStat[i] <- diffTRMeans[i]/stdErrRRFC
################################################################################
# tPr[i] <- 2 * pt(tStat[i], ddfRRFC) # 7/27/2023
# corrected 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), ddfRRFC, lower.tail = FALSE)
################################################################################
CIRRFC[i, ] <- sort(c(diffTRMeans[i] - qt(alpha/2, ddfRRFC) * stdErrRRFC,
diffTRMeans[i] + qt(alpha/2, ddfRRFC) * stdErrRRFC))
}
ciDiffTrtRRFC <- data.frame(Treatment = diffTRName,
Estimate = diffTRMeans,
StdErr = rep(stdErrRRFC, choose(I, 2)),
DF = rep(ddfRRFC, choose(I, 2)),
t = tStat,
p = tPr,
CI = CIRRFC)
colnames(ciDiffTrtRRFC) <- c("Treatment",
"Estimate",
"StdErr",
"DF",
"t",
"Pr > t",
"CI Lower",
"CI Upper")
dfSingleRRFC <- array(dim = I)
msDenSingleRRFC <- array(dim = I)
stdErrSingleRRFC <- array(dim = I)
CISingleRRFC <- array(dim = c(I, 2))
for (i in 1:I) {
msDenSingleRRFC[i] <- msRSingle[i]
dfSingleRRFC[i] <- (J - 1)
stdErrSingleRRFC[i] <- sqrt(msDenSingleRRFC[i]/J/K)
CISingleRRFC[i, ] <- sort(c(trMeans[i] - qt(alpha/2, dfSingleRRFC[i]) * stdErrSingleRRFC[i],
trMeans[i] + qt(alpha/2, dfSingleRRFC[i]) * stdErrSingleRRFC[i]))
}
ciAvgRdrEachTrtRRFC <- data.frame(Treatment = modalityID,
Area = trMeans,
StdErr = stdErrSingleRRFC,
DF = dfSingleRRFC,
CI = CISingleRRFC,
row.names = NULL)
colnames(ciAvgRdrEachTrtRRFC) <- c("Treatment",
"Area",
"StdErr",
"DF",
"CI Lower",
"CI Upper")
} else {
fRRFC <- NA
ddfRRFC <- NA
pRRFC <- NA
ciDiffTrtRRFC <- NA
ciAvgRdrEachTrtRRFC <- NA
}
if (analysisOption == "RRFC")
return(list(fomArray = as.data.frame(fomArray),
anovaY = anovaY,
anovaYi = msSingleTable,
varComp = varComp,
fRRFC = fRRFC,
ddfRRFC = ddfRRFC,
pRRFC = pRRFC,
ciDiffTrtRRFC = ciDiffTrtRRFC,
ciAvgRdrEachTrtRRFC = ciAvgRdrEachTrtRRFC))
}
# analysisOption == "ALL"
return(list(fomArray = as.data.frame(fomArray),
anovaY = anovaY,
anovaYi = msSingleTable,
varComp = varComp,
fRRRC = fRRRC,
ddfRRRC = ddfRRRC,
pRRRC = pRRRC,
ciDiffTrtRRRC = ciDiffTrtRRRC,
ciAvgRdrEachTrtRRRC = ciAvgRdrEachTrtRRRC,
fFRRC = fFRRC,
ddfFRRC = ddfFRRC,
pFRRC = pFRRC,
ciDiffTrtFRRC = ciDiffTrtFRRC,
ciAvgRdrEachTrtFRRC = ciAvgRdrEachTrtFRRC,
ssAnovaEachRdr = ssTableFRRC,
msAnovaEachRdr = msTableFRRC,
ciDiffTrtEachRdr = ciDiffTrtEachRdr,
fRRFC = fRRFC,
ddfRRFC = ddfRRFC,
pRRFC = pRRFC,
ciDiffTrtRRFC = ciDiffTrtRRFC,
ciAvgRdrEachTrtRRFC = ciAvgRdrEachTrtRRFC))
}
# Old Code from 0.0.1 (the first version on CRAN)
# The error in p-value (noted by Lucy A) has been corrected below
# Otherwise the code is identical to XZ code
ORHAnalysis <- function(dataset, FOM, alpha, covEstMethod, nBoots, analysisOption)
{
NL <- dataset$ratings$NL
LL <- dataset$ratings$LL
lesionVector <- dataset$lesions$perCase
lesionID <- dataset$lesions$IDs
lesionWeight <- dataset$lesions$weights
maxNL <- dim(NL)[4]
modalityID <- dataset$descriptions$modalityID
readerID <- dataset$descriptions$readerID
I <- length(modalityID)
J <- length(readerID)
K <- dim(NL)[3]
K2 <- dim(LL)[3]
dim(NL) <- c(I, J, K, maxNL)
dim(LL) <- c(I, J, K2, max(lesionVector))
if (!analysisOption %in% c("RRRC", "FRRC", "RRFC", "ALL")){
errMsg <- sprintf("%s is not an available analysisOption.", analysisOption)
stop(errMsg)
}
if (I < 2) {
stop("The analysis requires at least 2 modalities")
}
if (!covEstMethod %in% c("jackknife", "bootstrap", "DeLong")) {
errMsg <- paste0(covEstMethod, " is not an allowed covariance estimation method.")
stop(errMsg)
}
fomArray <- UtilFigureOfMerit(dataset, FOM)
trMeans <- rowMeans(fomArray)
fomMean <- mean(fomArray)
ret <- EstimateVarCov(fomArray, NL, LL, lesionVector, lesionID, lesionWeight, maxNL, FOM, covEstMethod, nBoots)
var <- ret$var
cov1 <- ret$cov1
cov2 <- ret$cov2
cov3 <- ret$cov3
msT <- 0
for (i in 1:I) {
msT <- msT + (mean(fomArray[i, ]) - fomMean)^2
}
msT <- J * msT/(I - 1)
msR <- 0
for (j in 1:J) {
msR <- msR + (mean(fomArray[, j]) - fomMean)^2
}
msR <- I * msR/(J - 1)
msTR <- 0
for (i in 1:I) {
for (j in 1:J) {
msTR <- msTR + (fomArray[i, j] - mean(fomArray[i, ]) - mean(fomArray[, j]) + fomMean)^2
}
}
msTR <- msTR/((J - 1) * (I - 1))
varTR <- msTR - var + cov1 + max(cov2 - cov3, 0)
varR <- (msR - var - (I - 1) * cov1 + cov2 + (I - 1) * cov3 - varTR)/I
varCovArray <- c(varR, varTR, cov1, cov2, cov3, var)
nameArray <- c("Var(R)", "Var(T*R)", "COV1", "COV2", "COV3", "Var(Error)")
varComp <- data.frame(varCov = varCovArray, row.names = nameArray)
varSingle <- vector(length = I)
cov2Single <- vector(length = I)
for (i in 1:I) {
fomSingle <- fomArray[i, ]
nl <- NL[i, , , ]
ll <- LL[i, , , ]
dim(fomSingle) <- c(1, J)
dim(nl) <- c(1, J, K, maxNL)
dim(ll) <- c(1, J, K2, max(lesionVector))
ret <- EstimateVarCov(fomSingle, nl, ll, lesionVector, lesionID, lesionWeight, maxNL, FOM, covEstMethod, nBoots)
varSingle[i] <- ret$var
if (J > 1) {
cov2Single[i] <- ret$cov2
} else {
cov2Single[i] <- 0
}
}
varEchRder <- vector(length = J)
cov1EchRder <- vector(length = J)
for (j in 1:J) {
fomSingle <- fomArray[, j]
nl <- NL[, j, , ]
ll <- LL[, j, , ]
dim(fomSingle) <- c(I, 1)
dim(nl) <- c(I, 1, K, maxNL)
dim(ll) <- c(I, 1, K2, max(lesionVector))
ret <- EstimateVarCov(fomSingle, nl, ll, lesionVector, lesionID, lesionWeight, maxNL, FOM, covEstMethod, nBoots)
varEchRder[j] <- ret$var
cov1EchRder[j] <- ret$cov1
}
msRSingle <- array(0, dim = c(I))
for (i in 1:I) {
msRSingle[i] <- sum((fomArray[i, ] - trMeans[i])^2)/(J - 1)
}
diffTRMeans <- array(dim = choose(I, 2))
diffTRName <- array(dim = choose(I, 2))
ii <- 1
for (i in 1:I) {
if (i == I)
break
for (ip in (i + 1):I) {
diffTRMeans[ii] <- trMeans[i] - trMeans[ip]
diffTRName[ii] <- paste(modalityID[i], modalityID[ip], sep = "-")
ii <- ii + 1
}
}
msNum <- msT
# ************ RRRC ****************
if (analysisOption %in% c("RRRC", "ALL")) {
if (J > 1) {
msDenRRRC <- msTR + max(J * (cov2 - cov3), 0)
fRRRC <- msNum/msDenRRRC
ddfRRRC <- msDenRRRC^2/(msTR^2/((I - 1) * (J - 1)))
pRRRC <- 1 - pf(fRRRC, I - 1, ddfRRRC)
stdErrRRRC <- sqrt(2 * msDenRRRC/J)
tStat <- vector()
tPr <- vector()
CIRRRC <- array(dim = c(length(diffTRMeans), 2))
for (i in 1:length(diffTRMeans)) {
tStat[i] <- diffTRMeans[i]/stdErrRRRC
################################################################################
# tPr[i] <- 2 * pt(tStat[i], ddfRRRC) # this was only ONE of many such mistakes in the original version of CRAN code
# corrected prior to 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), ddfRRRC, lower.tail = FALSE)
################################################################################
CIRRRC[i, ] <- sort(c(diffTRMeans[i] - qt(alpha/2, ddfRRRC) * stdErrRRRC,
diffTRMeans[i] + qt(alpha/2, ddfRRRC) * stdErrRRRC))
}
ciDiffTrtRRRC <- data.frame(Treatment = diffTRName,
Estimate = diffTRMeans,
StdErr = rep(stdErrRRRC, choose(I, 2)),
DF = rep(ddfRRRC, choose(I, 2)),
t = tStat,
p = tPr,
CI = CIRRRC)
colnames(ciDiffTrtRRRC) <- c("Treatment", "Estimate", "StdErr", "DF", "t", "Pr > t", "CI Lower", "CI Upper")
dfSingleRRRC <- array(dim = I)
msDenSingleRRRC <- array(dim = I)
stdErrSingleRRRC <- array(dim = I)
CISingleRRRC <- array(dim = c(I, 2))
for (i in 1:I) {
msDenSingleRRRC[i] <- msRSingle[i] + max(J * cov2Single[i], 0)
dfSingleRRRC[i] <- msDenSingleRRRC[i]^2/msRSingle[i]^2 * (J - 1)
stdErrSingleRRRC[i] <- sqrt(msDenSingleRRRC[i]/J)
CISingleRRRC[i, ] <- sort(c(trMeans[i] - qt(alpha/2, dfSingleRRRC[i]) * stdErrSingleRRRC[i],
trMeans[i] + qt(alpha/2, dfSingleRRRC[i]) * stdErrSingleRRRC[i]))
}
ciAvgRdrEachTrtRRRC <- data.frame(Treatment = modalityID,
Area = trMeans,
StdErr = stdErrSingleRRRC,
DF = dfSingleRRRC,
CI = CISingleRRRC,
row.names = NULL)
colnames(ciAvgRdrEachTrtRRRC) <- c("Treatment", "Area", "StdErr", "DF", "CI Lower", "CI Upper")
} else {
fRRRC <- NA
ddfRRRC <- NA
pRRRC <- NA
ciDiffTrtRRRC <- NA
ciAvgRdrEachTrtRRRC <- NA
}
if (analysisOption == "RRRC"){
return(list(fomArray = as.data.frame(fomArray),
msT = msT,
msR = msR,
msRSingle = msRSingle,
msTR = msTR,
varComp = varComp,
fRRRC = fRRRC,
ddfRRRC = ddfRRRC,
pRRRC = pRRRC,
ciDiffTrtRRRC = ciDiffTrtRRRC,
ciAvgRdrEachTrtRRRC = ciAvgRdrEachTrtRRRC))
}
}
# ************ FRRC ****************
if (analysisOption %in% c("FRRC", "ALL")) {
if (J > 1) {
msDenFRRC <- var - cov1 + (J - 1) * (cov2 - cov3)
} else {
msDenFRRC <- var - cov1
}
fFRRC <- msNum/msDenFRRC
ddfFRRC <- Inf
pFRRC <- 1 - pf(fFRRC, I - 1, ddfFRRC)
stdErrFRRC <- sqrt(2 * msDenFRRC/J)
tStat <- vector()
tPr <- vector()
CIFRRC <- array(dim = c(length(diffTRMeans), 2))
for (i in 1:length(diffTRMeans)) {
tStat[i] <- diffTRMeans[i]/stdErrFRRC
################################################################################
# tPr[i] <- 2 * pt(tStat[i], ddfFRRC) 7/27/2023
# Corrected on 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), ddfFRRC, lower.tail = FALSE)
################################################################################
CIFRRC[i, ] <- sort(c(diffTRMeans[i] - qt(alpha/2, ddfFRRC) * stdErrFRRC,
diffTRMeans[i] + qt(alpha/2, ddfFRRC) * stdErrFRRC))
}
ciDiffTrtFRRC <- data.frame(Treatment = diffTRName, Estimate = diffTRMeans, StdErr = rep(stdErrFRRC, choose(I, 2)), DF = rep(ddfFRRC, choose(I, 2)), t = tStat, p = tPr, CI = CIFRRC)
colnames(ciDiffTrtFRRC) <- c("Treatment", "Estimate", "StdErr", "DF", "t", "Pr > t", "CI Lower", "CI Upper")
dfSingleFRRC <- array(dim = I)
msDenSingleFRRC <- array(dim = I)
stdErrSingleFRRC <- array(dim = I)
CISingleFRRC <- array(dim = c(I, 2))
for (i in 1:I) {
msDenSingleFRRC[i] <- varSingle[i] + (J - 1) * cov2Single[i]
dfSingleFRRC[i] <- Inf
stdErrSingleFRRC[i] <- sqrt(msDenSingleFRRC[i]/J)
CISingleFRRC[i, ] <- sort(c(trMeans[i] - qt(alpha/2, dfSingleFRRC[i]) * stdErrSingleFRRC[i], trMeans[i] + qt(alpha/2, dfSingleFRRC[i]) * stdErrSingleFRRC[i]))
}
ciAvgRdrEachTrtFRRC <- data.frame(Treatment = modalityID, Area = trMeans, StdErr = stdErrSingleFRRC, DF = dfSingleFRRC, CI = CISingleFRRC, row.names = NULL)
colnames(ciAvgRdrEachTrtFRRC) <- c("Treatment", "Area", "StdErr", "DF", "CI Lower", "CI Upper")
diffTRMeansFRRC <- array(dim = c(J, choose(I, 2)))
for (j in 1:J) {
ii <- 1
for (i in 1:I) {
if (i == I)
break
for (ip in (i + 1):I) {
diffTRMeansFRRC[j, ii] <- fomArray[i, j] - fomArray[ip, j]
ii <- ii + 1
}
}
}
diffTRMeansFRRC <- as.vector(t(diffTRMeansFRRC))
stdErrFRRC <- sqrt(2 * (varEchRder - cov1EchRder))
stdErrFRRC <- rep(stdErrFRRC, choose(I, 2))
dim(stdErrFRRC) <- c(J, choose(I, 2))
stdErrFRRC <- as.vector(t(stdErrFRRC))
readerNames <- rep(readerID, choose(I, 2))
dim(readerNames) <- c(J, choose(I, 2))
readerNames <- as.vector(t(readerNames))
trNames <- rep(diffTRName, J)
dfReaderFRRC <- rep(Inf, length(stdErrFRRC))
CIReaderFRRC <- array(dim = c(length(stdErrFRRC), 2))
tStat <- vector()
tPr <- vector()
for (i in 1:length(stdErrFRRC)) {
tStat[i] <- diffTRMeansFRRC[i]/stdErrFRRC[i]
################################################################################
# tPr[i] <- 2 * pt(tStat[i], dfReaderFRRC[i]) 7/27/2023
# Corrected on 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), dfReaderFRRC[i], lower.tail = FALSE)
################################################################################
CIReaderFRRC[i, ] <- sort(c(diffTRMeansFRRC[i] - qt(alpha/2, dfReaderFRRC[i]) * stdErrFRRC[i],
diffTRMeansFRRC[i] + qt(alpha/2, dfReaderFRRC[i]) * stdErrFRRC[i]))
}
ciDiffTrtEachRdr <- data.frame(Reader = readerNames,
Treatment = trNames,
Estimate = diffTRMeansFRRC,
StdErr = stdErrFRRC,
DF = dfReaderFRRC,
t = tStat,
p = tPr,
CI = CIReaderFRRC)
colnames(ciDiffTrtEachRdr) <- c("Reader", "Treatment", "Estimate", "StdErr", "DF", "t", "Pr > t", "CI Lower", "CI Upper")
varCovEachRdr <- data.frame(readerID, varEchRder, cov1EchRder)
colnames(varCovEachRdr) <- c("Reader", "Var", "Cov1")
if (analysisOption == "FRRC"){
return(list(fomArray = as.data.frame(fomArray),
msT = msT,
msTR = msTR,
varComp = varComp,
fFRRC = fFRRC,
ddfFRRC = ddfFRRC,
pFRRC = pFRRC,
ciDiffTrtFRRC = ciDiffTrtFRRC,
ciAvgRdrEachTrtFRRC = ciAvgRdrEachTrtFRRC,
ciDiffTrtEachRdr = ciDiffTrtEachRdr,
varCovEachRdr = varCovEachRdr
))
}
}
# ************ RRFC ****************
if (analysisOption %in% c("RRFC", "ALL")) {
if (J > 1) {
msDenRRFC <- msTR
fRRFC <- msNum/msDenRRFC
ddfRRFC <- ((I - 1) * (J - 1))
pRRFC <- 1 - pf(fRRFC, I - 1, ddfRRFC)
stdErrRRFC <- sqrt(2 * msDenRRFC/J)
tStat <- vector()
tPr <- vector()
CIRRFC <- array(dim = c(length(diffTRMeans), 2))
for (i in 1:length(diffTRMeans)) {
tStat[i] <- diffTRMeans[i]/stdErrRRFC
################################################################################
# tPr[i] <- 2 * pt(tStat[i], ddfRRFC) 7/27/2023
# Corrected on 7/27/2023 as follows:
tPr[i] <- 2 * pt(abs(tStat[i]), ddfRRFC, lower.tail = FALSE)
################################################################################
CIRRFC[i, ] <- sort(c(diffTRMeans[i] - qt(alpha/2, ddfRRFC) * stdErrRRFC,
diffTRMeans[i] + qt(alpha/2, ddfRRFC) * stdErrRRFC))
}
ciDiffTrtRRFC <- data.frame(Treatment = diffTRName, Estimate = diffTRMeans, StdErr = rep(stdErrRRFC, choose(I, 2)), DF = rep(ddfRRFC, choose(I, 2)), t = tStat, p = tPr, CI = CIRRFC)
colnames(ciDiffTrtRRFC) <- c("Treatment", "Estimate", "StdErr", "DF", "t", "Pr > t", "CI Lower", "CI Upper")
dfSingleRRFC <- array(dim = I)
msDenSingleRRFC <- array(dim = I)
stdErrSingleRRFC <- array(dim = I)
CISingleRRFC <- array(dim = c(I, 2))
for (i in 1:I) {
msDenSingleRRFC[i] <- msRSingle[i]
dfSingleRRFC[i] <- (J - 1)
stdErrSingleRRFC[i] <- sqrt(msDenSingleRRFC[i]/J)
CISingleRRFC[i, ] <- sort(c(trMeans[i] - qt(alpha/2, dfSingleRRFC[i]) * stdErrSingleRRFC[i], trMeans[i] + qt(alpha/2, dfSingleRRFC[i]) * stdErrSingleRRFC[i]))
}
ciAvgRdrEachTrtRRFC <- data.frame(Treatment = modalityID, Area = trMeans, StdErr = stdErrSingleRRFC, DF = dfSingleRRFC, CI = CISingleRRFC, row.names = NULL)
colnames(ciAvgRdrEachTrtRRFC) <- c("Treatment", "Area", "StdErr", "DF", "CI Lower", "CI Upper")
} else {
fRRFC <- NA
ddfRRFC <- NA
pRRFC <- NA
ciDiffTrtRRFC <- NA
ciAvgRdrEachTrtRRFC <- NA
}
if (analysisOption == "RRFC"){
return(list(fomArray = as.data.frame(fomArray),
msT = msT,
msTR = msTR,
varComp = varComp,
fRRFC = fRRFC,
ddfRRFC = ddfRRFC,
pRRFC = pRRFC,
ciDiffTrtRRFC = ciDiffTrtRRFC,
ciAvgRdrEachTrtRRFC = ciAvgRdrEachTrtRRFC))
}
}
return(list(fomArray = as.data.frame(fomArray),
msT = msT,
msTR = msTR,
varComp = varComp,
fRRRC = fRRRC,
ddfRRRC = ddfRRRC,
pRRRC = pRRRC,
ciDiffTrtRRRC = ciDiffTrtRRRC,
ciAvgRdrEachTrtRRRC = ciAvgRdrEachTrtRRRC,
fFRRC = fFRRC,
ddfFRRC = ddfFRRC,
pFRRC = pFRRC,
ciDiffTrtFRRC = ciDiffTrtFRRC,
ciAvgRdrEachTrtFRRC = ciAvgRdrEachTrtFRRC,
ciDiffTrtEachRdr = ciDiffTrtEachRdr,
varCovEachRdr = varCovEachRdr,
fRRFC = fRRFC,
ddfRRFC = ddfRRFC,
pRRFC = pRRFC,
ciDiffTrtRRFC = ciDiffTrtRRFC,
ciAvgRdrEachTrtRRFC = ciAvgRdrEachTrtRRFC))
}
# this function is needed ONLY with OldCode
EstimateVarCov <- function(fomArray, NL, LL, lesionVector, lesionID, lesionWeight, maxNL, fom, covEstMethod, nBoots) {
UNINITIALIZED <- RJafrocEnv$UNINITIALIZED
I <- dim(NL)[1]
J <- dim(NL)[2]
K <- dim(NL)[3]
K2 <- dim(LL)[3]
K1 <- K - K2
if (covEstMethod == "jackknife") {
if (fom %in% c("MaxNLF", "HrSp")) {
jkFOMArray <- array(dim = c(I, J, K1))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K1) {
nl <- NL[i, j, -k, ]
ll <- LL[i, j, , ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector, lesionID, lesionWeight, maxNL, fom)
}
}
}
} else if (fom %in% c("MaxLLF", "HrSe")) {
jkFOMArray <- array(dim = c(I, J, K2))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K2) {
nl <- NL[i, j, -(k + K1), ]
ll <- LL[i, j, -k, ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2 - 1, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector[-k], lesionID[-k, ], lesionWeight[-k, ], maxNL, fom)
}
}
}
} else {
jkFOMArray <- array(dim = c(I, J, K))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:K) {
if (k <= K1) {
nl <- NL[i, j, -k, ]
ll <- LL[i, j, , ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector, lesionID, lesionWeight, maxNL, fom)
} else {
nl <- NL[i, j, -k, ]
ll <- LL[i, j, -(k - K1), ]
dim(nl) <- c(K - 1, maxNL)
dim(ll) <- c(K2 - 1, max(lesionVector))
jkFOMArray[i, j, k] <- MyFOMOldCode(nl, ll, lesionVector[-(k - K1)], lesionID[-(k - K1), ], lesionWeight[-(k - K1), ], maxNL, fom)
}
}
}
}
}
Cov <- FOM2ORVarComp(jkFOMArray, varInflFactor = TRUE)
var <- Cov$Var
cov1 <- Cov$Cov1
cov2 <- Cov$Cov2
cov3 <- Cov$Cov3
}
if (covEstMethod == "bootstrap") {
if (fom %in% c("MaxNLF", "HrSp")) {
fomBsArray <- array(dim = c(I, J, nBoots))
for (b in 1:nBoots) {
kBs <- ceiling(runif(K1) * K1)
for (i in 1:I) {
for (j in 1:J) {
nlBs <- NL[i, j, kBs, ]
llBs <- LL[i, j, , ]
dim(nlBs) <- c(K1, maxNL)
dim(llBs) <- c(K2, max(lesionVector))
fomBsArray[i, j, b] <- MyFOMOldCode(nlBs, llBs, lesionVector, lesionID, lesionWeight, maxNL, fom)
}
}
}
} else if (fom %in% c("MaxLLF", "HrSe")) {
fomBsArray <- array(dim = c(I, J, nBoots))
for (b in 1:nBoots) {
kBs <- ceiling(runif(K2) * K2)
for (i in 1:I) {
for (j in 1:J) {
nlBs <- NL[i, j, c(1:K1, (kBs + K1)), ]
llBs <- LL[i, j, kBs, ]
dim(nlBs) <- c(K, maxNL)
dim(llBs) <- c(K2, max(lesionVector))
fomBsArray[i, j, b] <- MyFOMOldCode(nlBs, llBs, lesionVector[kBs], lesionID[kBs, ], lesionWeight[kBs, ], maxNL, fom)
}
}
}
} else {
fomBsArray <- array(dim = c(I, J, nBoots))
for (b in 1:nBoots) {
kBs1 <- ceiling(runif(K1) * K1)
kBs2 <- ceiling(runif(K2) * K2)
for (i in 1:I) {
for (j in 1:J) {
nlBs <- NL[i, j, c(kBs1, kBs2 + K1), ]
llBs <- LL[i, j, kBs2, ]
dim(nlBs) <- c(K, maxNL)
dim(llBs) <- c(K2, max(lesionVector))
fomBsArray[i, j, b] <- MyFOMOldCode(nlBs, llBs, lesionVector[kBs2], lesionID[kBs2, ], lesionWeight[kBs2, ], maxNL, fom)
}
}
}
}
Cov <- FOM2ORVarComp(fomBsArray, varInflFactor = FALSE)
var <- Cov$var
cov1 <- Cov$cov1
cov2 <- Cov$cov2
cov3 <- Cov$cov3
}
if (covEstMethod == "DeLong") {
if (!fom %in% c("Wilcoxon", "HrAuc", "ROI"))
stop("DeLong\"s method can only be used for trapezoidal figures of merit.")
if (fom == "ROI") {
kI01 <- which(apply((NL[1, 1, , ] != UNINITIALIZED), 1, any))
numKI01 <- rowSums((NL[1, 1, , ] != UNINITIALIZED))
I01 <- length(kI01)
I10 <- K2
N <- sum((NL[1, 1, , ] != UNINITIALIZED))
M <- sum(lesionVector)
V01 <- array(dim = c(I, J, I01, maxNL))
V10 <- array(dim = c(I, J, I10, max(lesionVector)))
for (i in 1:I) {
for (j in 1:J) {
for (k in 1:I10) {
for (el in 1:lesionVector[k]) {
V10[i, j, k, el] <- (sum(as.vector(NL[i, j, , ][NL[i, j, , ] != UNINITIALIZED]) < LL[i, j, k, el]) + 0.5 * sum(as.vector(NL[i, j, , ][NL[i, j, , ] != UNINITIALIZED]) == LL[i, j, k, el]))/N
}
}
for (k in kI01) {
for (el in 1:maxNL) {
if (NL[i, j, k, el] == UNINITIALIZED)
next
V01[i, j, k, el] <- (sum(NL[i, j, k, el] < as.vector(LL[i, j, , ][LL[i, j, , ] != UNINITIALIZED])) + 0.5 * sum(NL[i, j, k, el] == as.vector(LL[i, j, , ][LL[i, j, , ] != UNINITIALIZED])))/M
}
}
}
}
s10 <- array(0, dim = c(I, I, J, J))
s01 <- array(0, dim = c(I, I, J, J))
s11 <- array(0, dim = c(I, I, J, J))
for (i in 1:I) {
for (ip in 1:I) {
for (j in 1:J) {
for (jp in 1:J) {
for (k in 1:I10) {
s10[i, ip, j, jp] <- (s10[i, ip, j, jp] + (sum(V10[i, j, k, !is.na(V10[i, j, k, ])]) - lesionVector[k] * fomArray[i, j]) * (sum(V10[ip, jp, k, !is.na(V10[ip, jp, k, ])]) - lesionVector[k] * fomArray[ip,
jp]))
}
for (k in kI01) {
s01[i, ip, j, jp] <- (s01[i, ip, j, jp] + (sum(V01[i, j, k, !is.na(V01[i, j, k, ])]) - numKI01[k] * fomArray[i, j]) * (sum(V01[ip, jp, k, !is.na(V01[ip, jp, k, ])]) - numKI01[k] * fomArray[ip,
jp]))
}
for (k in 1:K2) {
if (all(NL[ip, jp, k + K1, ] == UNINITIALIZED))
next
s11[i, ip, j, jp] <- (s11[i, ip, j, jp] + (sum(V10[i, j, k, !is.na(V10[i, j, k, ])]) - lesionVector[k] * fomArray[i, j]) * (sum(V01[ip, jp, k + K1, !is.na(V01[ip, jp, k + K1, ])]) - numKI01[K1 +
k] * fomArray[ip, jp]))
}
}
}
}
}
s10 <- s10 * I10/(I10 - 1)/M
s01 <- s01 * I01/(I01 - 1)/N
s11 <- s11 * K/(K - 1)
S <- array(0, dim = c(I, I, J, J))
for (i in 1:I) {
for (ip in 1:I) {
for (j in 1:J) {
for (jp in 1:J) {
S[i, ip, j, jp] <- s10[i, ip, j, jp]/M + s01[i, ip, j, jp]/N + s11[i, ip, j, jp]/(M * N) + s11[ip, i, jp, j]/(M * N)
}
}
}
}
} else {
# ROI
V10 <- array(dim = c(I, J, K2))
V01 <- array(dim = c(I, J, K1))
for (i in 1:I) {
for (j in 1:J) {
nl <- NL[i, j, 1:K1, ]
ll <- cbind(NL[i, j, (K1 + 1):K, ], LL[i, j, , ])
dim(nl) <- c(K1, maxNL)
dim(ll) <- c(K2, maxNL + max(lesionVector))
fp <- apply(nl, 1, max)
tp <- apply(ll, 1, max)
for (k in 1:K2) {
V10[i, j, k] <- (sum(fp < tp[k]) + 0.5 * sum(fp == tp[k]))/K1
}
for (k in 1:K1) {
V01[i, j, k] <- (sum(fp[k] < tp) + 0.5 * sum(fp[k] == tp))/K2
}
}
}
s10 <- array(dim = c(I, I, J, J))
s01 <- array(dim = c(I, I, J, J))
for (i in 1:I) {
for (ip in 1:I) {
for (j in 1:J) {
for (jp in 1:J) {
s10[i, ip, j, jp] <- sum((V10[i, j, ] - fomArray[i, j]) * (V10[ip, jp, ] - fomArray[ip, jp]))/(K2 - 1)
s01[i, ip, j, jp] <- sum((V01[i, j, ] - fomArray[i, j]) * (V01[ip, jp, ] - fomArray[ip, jp]))/(K1 - 1)
}
}
}
}
S <- s10/K2 + s01/K1
}
Cov <- FOM2ORVarComp(S, varInflFactor = FALSE)
var <- Cov$var
cov1 <- Cov$cov1
cov2 <- Cov$cov2
cov3 <- Cov$cov3
}
return(list(
var = var,
cov1 = cov1,
cov2 = cov2,
cov3 = cov3
))
}
# this function is needed ONLY with OldCode
MyFOMOldCode <- function(nl, ll, lesionNum, lesionID, lesionWeight, maxNL, fom) {
UNINITIALIZED <- RJafrocEnv$UNINITIALIZED
K <- length(nl[, 1])
K2 <- length(ll[, 1])
K1 <- K - K2
FOM <- NA
if (fom == "Wilcoxon") {
truth <- c(rep(0, K1), rep(1, K2))
ratings <- c(nl[1:K1], ll)
FOM <- TrapezoidalArea1(ratings, truth)
} else if (fom == "HrAuc") {
truth <- c(rep(0, K1), rep(1, K2))
ratings <- array(dim = c(K1 + K2))
for (k in 1:K1) {
ratings[k] <- max(nl[k, ])
}
for (k in (K1 + 1):(K1 + K2)) {
ratings[k] <- max(c(nl[k, ], ll[k - K1, ]))
}
FOM <- TrapezoidalArea1(ratings, truth)
} else if (fom == "HrSe") {
tpCount <- 0
for (k in 1:K2) {
tempMax <- UNINITIALIZED
for (el in 1:lesionNum[k]) tempMax <- max(tempMax, ll[k, el])
for (el in 1:maxNL) tempMax <- max(tempMax, nl[k + K1, el])
if (tempMax > UNINITIALIZED)
tpCount <- tpCount + 1
}
FOM <- tpCount/K2
} else if (fom == "HrSp") {
fpCount <- 0
for (k in 1:K1) {
tempMax <- UNINITIALIZED
for (el in 1:maxNL) tempMax <- max(tempMax, nl[k, el])
if (tempMax > UNINITIALIZED)
fpCount <- fpCount + 1
}
FOM <- 1 - fpCount/K1
} else if (fom == "SongA1") {
x <- array(UNINITIALIZED, dim = c(K))
y1 <- array(UNINITIALIZED, dim = c(K2))
y2 <- array(UNINITIALIZED, dim = c(K2))
kInd <- which(apply(nl != UNINITIALIZED, 1, any))
for (k in kInd){
x[ k ] <- mean(nl[ k, ( nl[ k, ] != UNINITIALIZED )])
}
y1 <- x[(K1 + 1):K]
x <- x[1:K1]
px0 <- sum(x == UNINITIALIZED)
kInd <- which(apply(ll != UNINITIALIZED, 1, any))
for (k in kInd){
y2[ k ] <- mean(ll[ k, ( ll[ k, ] != UNINITIALIZED )])
}
y <- apply(cbind(y1, y2), 1, max)
py0 <- sum(y == UNINITIALIZED)
FOM <- 0
x <- x[x != UNINITIALIZED]
y <- y[y != UNINITIALIZED]
for ( k1 in 1 : length(x)){
FOM <- FOM + sum(x[k1] < y) + 0.5 * sum(x[k1] == y)
}
FOM <- FOM / (K1 * K2) + (px0 / K1) * ( 1 - 0.5 * py0 / K2)
} else if (fom == "SongA2") {
n0 <- apply(nl[1:K1, ] != UNINITIALIZED, 1, sum)
n1 <- apply(cbind(nl[(K1 + 1):K, ], ll) != UNINITIALIZED, 1, sum)
k1Ind <- which(n0 > 0)
k2Ind <- which(n1 > 0)
px0 <- K1 - length(k1Ind)
py0 <- K2 - length(k2Ind)
v1 <- 0.5
tw <- 0
for (k1 in k1Ind) {
for (k2 in k2Ind) {
v2 <- 0
x <- nl[k1, ][nl[k1, ] != UNINITIALIZED]
y1 <- nl[K1 + k2, ][nl[K1 + k2, ] != UNINITIALIZED]
y2 <- ll[k2, ][ll[k2, ] != UNINITIALIZED]
for (l in 1:length(x)){
v2 <- v2 + sum(x[l] < y1) + 0.5 * sum(x[l] == y1) + sum(x[l] < y2) + 0.5 * sum(x[l] == y2)
}
v2 <- v2/(n0[k1] * n1[k2])
tw <- tw + (v1 < v2) + 0.5 * (v1 == v2)
}
}
FOM <- tw/(K1 * K2) + px0/K1 * (1 - 0.5 * py0/K2)
} else if (fom == "AFROC1") {
numLesTotal <- sum(lesionNum)
les <- as.vector(ll[lesionID != UNINITIALIZED])
fp <- array(dim = c(K))
for (k in 1:K) {
fp[k] <- max(nl[k, ])
}
ratings <- c(fp, les)
truth <- c(rep(0, K), rep(1, numLesTotal))
FOM <- TrapezoidalArea1(ratings, truth)
} else if (fom == "AFROC") {
numLesTotal <- sum(lesionNum)
les <- as.vector(ll[lesionID != UNINITIALIZED])
fp <- array(dim = c(K1))
for (k in 1:K1) {
fp[k] <- max(nl[k, ])
}
ratings <- c(fp, les)
truth <- c(rep(0, K1), rep(1, length(les)))
FOM <- TrapezoidalArea1(ratings, truth)
} else if (fom == "wAFROC1") {
numLesTotal <- sum(lesionNum)
les <- as.vector(ll[lesionID != UNINITIALIZED])
fp <- array(dim = c(K))
for (k in 1:K) {
fp[k] <- max(nl[k, ])
}
ratings <- c(fp, les)
truth <- c(rep(0, K), rep(1, numLesTotal))
weights <- lesionWeight[lesionWeight != UNINITIALIZED]
FOM <- TrapezoidalAreaWeighted(ratings, truth, weights, K2)
} else if (fom == "wAFROC") {
numLesTotal <- sum(lesionNum)
les <- as.vector(ll[lesionID != UNINITIALIZED])
fp <- array(dim = c(K1))
for (k in 1:K1) {
fp[k] <- max(nl[k, ])
}
ratings <- c(fp, les)
truth <- c(rep(0, K1), rep(1, numLesTotal))
weights <- lesionWeight[lesionWeight != UNINITIALIZED]
FOM <- TrapezoidalAreaWeighted(ratings, truth, weights, K2)
} else if (fom == "MaxLLF") {
numMarksTotal <- sum(ll != UNINITIALIZED)
numLesTotal <- sum(lesionNum)
FOM <- numMarksTotal/numLesTotal
} else if (fom == "MaxNLF") {
numMarksTotal <- sum(nl[1:K1, ] != UNINITIALIZED)
FOM <- numMarksTotal/K1
} else if (fom == "MaxNLFAllCases") {
numMarksTotal <- sum(nl != UNINITIALIZED)
FOM <- numMarksTotal/K
} else if (fom == "ROI") {
nn <- 0
ns <- sum(lesionNum)
tw <- 0
for (k in 1:K) {
for (el in 1:maxNL) {
if (nl[k, el] == UNINITIALIZED)
next
nn <- nn + 1
tw <- tw + sum(nl[k, el] < as.vector(ll)) + 0.5 * sum(nl[k, el] == as.vector(ll))
}
}
FOM <- tw/(nn * ns)
} else {
errMsg <- paste0(fom, " is not an available figure of merit.")
stop(errMsg)
}
return(FOM)
}
# this function is needed ONLY with OldCode
TrapezoidalArea1 <- function(rocRatings, truth) {
K2 <- sum(truth)
K1 <- length(truth) - K2
K1Indx <- (1:length(truth))[!as.logical(truth)] #index of incorrect cases; or normal cases
K2Indx <- (1:length(truth))[as.logical(truth)] #index of correct cases; or abnormal cases
S <- 0
for (k in K1Indx) {
S <- S + sum(rocRatings[k] < rocRatings[K2Indx]) + 0.5 * sum(rocRatings[k] == rocRatings[K2Indx])
}
S <- S/(K1 * K2)
return(S)
}
# this function is needed ONLY with OldCode
TrapezoidalAreaWeighted <- function(rocRatings, truth, weights, numAbn) {
K2 <- sum(truth)
K1 <- length(truth) - K2
K1Indx <- (1:length(truth))[!as.logical(truth)] #index of incorrect cases; or normal cases
K2Indx <- (1:length(truth))[as.logical(truth)] #index of correct cases; or abnormal cases
S <- 0
for (k in K1Indx) {
S <- S + sum((rocRatings[k] < rocRatings[K2Indx]) * weights) + 0.5 * sum((rocRatings[k] == rocRatings[K2Indx]) * weights)
}
S <- S/(K1 * numAbn)
return(S)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.