Nothing
# Initialize the l'Ecuyer random number generator ####
#' @title Initialize the l'Ecuyer random number generator
#'
#' @description See the documentation for the \code{parallel package}.
#' If you require backwards compatibility, please run \code{RNGversion("3.5.0")}.
#'
#' @param seed This determines the position in each stream
#' @param stream This determines the stream
#'
#' @return Nothing
#'
#' @import parallel
#'
#' @export
#'
# @examples
init.lecuyerRNG <- function(seed = 1, stream = 2){
# Set the random number generator
RNGkind("L'Ecuyer-CMRG")
# Set the seed
set.seed(seed)
# Set the streamID
# Skip to the stream corresponding to the taskID
# Each taskID will have its own stream for parallel implementation
# IMPORTANT: .Random.seed must be reassigned in the global environment "<<-" not "<-"
i <- 1
while (i < stream) {
.Random.seed <<- nextRNGStream(.Random.seed)
i <- i + 1
}
}
# Simulate an MRMC data set ####
#' Simulate an MRMC data set
#'
#' @description
#' This program simulates observations from one set of readers scoring one set of cases.
#' It produces one modality and one truth state
#' of ROC data following Roe1997_Acad-Radiol_v4p298 and Roe1997_Acad-Radiol_v4p587.
#' In order to produce an entire ROC data set, please use sim.gRoeMetz.
#'
#' @details
#' The simulation is a linear model with one fixed effect and three
#' normally distributed independent random effects corresponding to readers,
#' cases, and an interaction between the two.
#'
#' L.rc = mu + readerEffect.r + caseEffect.c + readerXcaseEffect.rc
#'
#' @param simMRMC.config [list] of simulation parameters:
#' \itemize{
#' \item modalityID [character] label modalityID
#' \item readerIDs [factor] the ID of each reader
#' \item caseIDs [factor] the ID of each case
#' \item mu [numeric] mean
#' \item var_r [numeric] variance of random reader effect
#' \item var_c [numeric] variance of random case effect
#' \item var_rc [numeric] variance of random reader by case effect
#' }
#'
#' @return L [data.frame] with nC*nR rows of 4 variables
#' \itemize{
#' \item L$modalityID [factor] determined by input modalityID
#' \item L$readerID [factor] determined by input readerIDs
#' \item L$caseID [factor] determined by input caseIDs
#' \item L$score [numeric] R.r + C.c + RC.rc
#' \itemize{
#' \item r = 1,2,...,nR
#' \item c = 1,2,...,nC
#' \item R.r ~ N(0,var_r)
#' \item C.c ~ N(0,var_c)
#' \item RC.rc ~ N(0,var_rc)
#' }
#' }
#'
#' @export
#'
# @examples
# # Create a sample configuration object
# config <- sim.gRoeMetz.config()
# # Simulate an MRMC ROC data set
# dFrame.imrmc <- sim.gRoeMetz(config)
# # Analyze the MRMC ROC data
# result <- doIMRMC(dFrame.imrmc)
simMRMC <- function(simMRMC.config) {
modalityID <- simMRMC.config$modalityID
readerIDs <- simMRMC.config$readerIDs
caseIDs <- simMRMC.config$caseIDs
mu <- simMRMC.config$mu
var_c <- simMRMC.config$var_c
var_r <- simMRMC.config$var_r
var_rc <- simMRMC.config$var_rc
nR <- length(readerIDs)
nC <- length(caseIDs)
# Initialize data array with mu (fixed effect: grand mean)
L_mu <- rep(mu, nC * nR)
# Simulate case effects
L_cases <- stats::rnorm(nC, mean = 0, sd = sqrt(var_c))
# Replicate random case effects for each reader
L_cases <- rep(L_cases, nR)
caseID <- rep(caseIDs, nR)
# Simulate reader effects
L_readers <- stats::rnorm(nR, mean = 0, sd = sqrt(var_r))
# Replicate random reader effects for each case
# rep(nC, nR) indicates that each of the nR elements of readers
# is to be replicated nC times.
L_readers <- rep(L_readers, rep(nC, nR))
readerID <- rep(readerIDs, rep(nC, nR))
# Simulate reader by case interaction
L_reader.case <- stats::rnorm(nC*nR, mean = 0, sd = sqrt(var_rc))
modalityID <- rep(modalityID, nC*nR)
# Put it all together
L <- data.frame(
readerID = readerID,
caseID = caseID,
modalityID = modalityID,
score = L_mu + L_cases + L_readers + L_reader.case,
stringsAsFactors = TRUE
)
}
# Simulate an MRMC data set of an ROC experiment comparing two modalities ####
#' Simulate an MRMC data set of an ROC experiment comparing two modalities
#'
#' @description
#' This procedure simulates an MRMC data set of an ROC experiment comparing two modalities.
#' It is based on Gallas2014_J-Med-Img_v1p031006, which generalizes of the model in
#' Roe1997_Acad-Radiol_v4p298 and Roe1997_Acad-Radiol_v4p587. Specifically, it allows
#' the variance components to depend on the truth and the modality. For the simpler
#' Roe and Metz model, you can enter the smaller set of parameters into
#' sim.gRoeMetz.config and it will return a larger set of parameters that can be used with this function.
#'
#' @details
#' The simulation is a linear model with six fixed effects related to
#' modality and truth and 18 normally distributed independent random effects
#' for readers, cases, and the interaction between the two. Here is the linear model:
#'
#' L.mrct = mu.t + mu.mt \cr
#' + reader.rt + case.ct + readerXcase.rct \cr
#' + modalityXreader.mrt + modalityXcase.mct + modalityXreaderXcase.mrct \cr
#' \itemize{
#' \item m=modality (levels: A and b)
#' \item t=truth (levels: neg and Pos)
#' \item mu.t is the global mean for t=neg and t=pos cases
#' \item mu.mt is the modality specific fixed effects for t=neg and t=pos cases
#' \item the remaining terms are the random effects: all independent normal random variables
#' }
#'
#' @param config [list] of simulation parameters:
#' \itemize{
#' \item Experiment labels and size
#' \itemize{
#' \item modalityID.A: [character] label modality A
#' \item modalityID.B: [character] label modality B
#' \item nR: [numeric] number of readers
#' \item nC.neg: [numeric] number of signal-absent cases
#' \item nC.pos: [numeric] number of signal-present cases
#' }
#' \item There are six fixed effects:
#' \itemize{
#' \item mu.neg: [numeric] signal-absent (neg, global mean)
#' \item mu.pos: [numeric] signal-present (pos, global mean)
#' \item mu.Aneg: [numeric] modality A signal-absent (Aneg, modality effect)
#' \item mu.Bneg: [numeric] modality B signal-absent (Bneg, modality effect)
#' \item mu.Apos: [numeric] modality A signal-present (Apos, modality effect)
#' \item mu.Bpos: [numeric] modality B signal-present (Bpos, modality effect)
#' }
#' \item There are six random effects that are independent of modality
#' \itemize{
#' \item var_r.neg: [numeric] variance of random reader effect
#' \item var_c.neg: [numeric] variance of random case effect
#' \item var_rc.neg: [numeric] variance of random reader by case effect
#' \item var_r.pos: [numeric] variance of random reader effect
#' \item var_c.pos: [numeric] variance of random case effect
#' \item var_rc.pos: [numeric] variance of random reader by case effect
#' }
#' \item There are six random effects that are specific to modality A
#' \itemize{
#' \item var_r.Aneg: [numeric] variance of random reader effect
#' \item var_c.Aneg: [numeric] variance of random case effect
#' \item var_rc.Aneg: [numeric] variance of random reader by case effect
#' \item var_r.Apos: [numeric] variance of random reader effect
#' \item var_c.Apos: [numeric] variance of random case effect
#' \item var_rc.Apos: [numeric] variance of randome reader by case effect
#' }
#' \item There are six random effects that are specific to modality B
#' \itemize{
#' \item var_r.Bneg: [numeric] variance of random reader effect
#' \item var_c.Bneg: [numeric] variance of random case effect
#' \item var_rc.Bneg: [numeric] variance of random reader by case effect
#' \item var_r.Bpos: [numeric] variance of random reader effect
#' \item var_c.Bpos: [numeric] variance of random case effect
#' \item var_rc.Bpos: [numeric] variance of randome reader by case effect
#' }
#' }
#'
#' @return dFrame.imrmc [data.frame] with (nC.neg + nC.pos)*(nR+1) rows including
#' \itemize{
#' \item readerID: [factor] w/ nR levels "reader1", "reader2", ...
#' \item caseID: [factor] w/ nC levels "case1", "case2", ...
#' \item modalityID: [factor] w/ 1 level config$modalityID
#' \item score: [numeric] reader score
#' }
#'
#' Note that the first nC.neg + nC.pos rows specify the truth labels for each case.
#' For these rows, the readerID must be "truth"
#' and the score must be 0 for negative cases and 1 for positive cases.
#'
#' @export
#'
# @examples
# # Create a sample configuration object
# config <- sim.gRoeMetz.config()
# # Simulate an MRMC ROC data set
# dFrame.imrmc <- sim.gRoeMetz(config)
# # Analyze the MRMC ROC data
# result <- doIMRMC(dFrame.imrmc)
#
sim.gRoeMetz <- function(config) {
# Unpack modality labels
modalityID.A <- config$modalityID.A
modalityID.B <- config$modalityID.B
# Unpack experiment size
nR <- config$nR
nC.neg <- config$nC.neg
nC.pos <- config$nC.pos
# Assign readerIDs
readerIDs <- paste("reader", 1:nR, sep = "")
readerIDs <- factor(readerIDs, readerIDs, ordered = TRUE)
# Assign caseIDs
caseIDs.neg <- paste("negCase", 1:nC.neg, sep = "")
caseIDs.neg <- factor(caseIDs.neg, caseIDs.neg, ordered = TRUE)
caseIDs.pos <- paste("posCase", 1:nC.pos, sep = "")
caseIDs.pos <- factor(caseIDs.pos, caseIDs.pos, ordered = TRUE)
caseIDs <- c(as.character(caseIDs.neg), as.character(caseIDs.pos))
caseIDs <- factor(caseIDs, caseIDs, ordered = TRUE)
# Create data frame of truth
dFrame.truth <- data.frame(
readerID = rep("truth", nC.neg + nC.pos),
caseID = c(as.character(caseIDs.neg), as.character(caseIDs.pos)),
modalityID = rep("truth", nC.neg + nC.pos),
score = c(rep(0, nC.neg), rep(1, nC.pos)),
stringsAsFactors = TRUE
)
# Simulate the modality independent random effects, negative cases
neg.config <- list(
nR = nR,
nC = nC.neg,
modalityID = "empty",
readerIDs = readerIDs,
caseIDs = caseIDs.neg,
mu = config$mu.neg,
var_c = config$var_c.neg,
var_r = config$var_r.neg,
var_rc = config$var_rc.neg
)
dFrame.neg <- simMRMC(neg.config)
# Simulate the modality independent random effects, positive cases
pos.config <- list(
nR = nR,
nC = nC.pos,
modalityID = "empty",
readerIDs = readerIDs,
caseIDs = caseIDs.pos,
mu = config$mu.pos,
var_r = config$var_r.pos,
var_c = config$var_c.pos,
var_rc = config$var_rc.pos
)
dFrame.pos <- simMRMC(pos.config)
# Simulate modality A random effects, negative cases
modAneg.config <- list(
nR = nR,
nC = nC.neg,
modalityID = modalityID.A,
readerIDs = readerIDs,
caseIDs = caseIDs.neg,
mu = config$mu.Aneg,
var_r = config$var_r.Aneg,
var_c = config$var_c.Aneg,
var_rc = config$var_rc.Aneg
)
dFrame.modAneg <- simMRMC(modAneg.config)
# Simulate modality A random effects, positive cases
modApos.config <- list(
nR = nR,
nC = nC.pos,
modalityID = modalityID.A,
readerIDs = readerIDs,
caseIDs = caseIDs.pos,
mu = config$mu.Apos,
var_r = config$var_r.Apos,
var_c = config$var_c.Apos,
var_rc = config$var_rc.Apos
)
dFrame.modApos <- simMRMC(modApos.config)
# Simulate modality B random effects, negative cases
modBneg.config <- list(
nR = nR,
nC = nC.neg,
modalityID = modalityID.B,
readerIDs = readerIDs,
caseIDs = caseIDs.neg,
mu = config$mu.Bneg,
var_r = config$var_r.Bneg,
var_c = config$var_c.Bneg,
var_rc = config$var_rc.Bneg
)
dFrame.modBneg <- simMRMC(modBneg.config)
# Simulate modality B random effects, positive cases
modBpos.config <- list(
nR = nR,
nC = nC.pos,
modalityID = modalityID.B,
readerIDs = readerIDs,
caseIDs = caseIDs.pos,
mu = config$mu.Bpos,
var_r = config$var_r.Bpos,
var_c = config$var_c.Bpos,
var_rc = config$var_rc.Bpos
)
dFrame.modBpos <- simMRMC(modBpos.config)
dFrame.modAneg$score <- dFrame.neg$score + dFrame.modAneg$score
dFrame.modApos$score <- dFrame.pos$score + dFrame.modApos$score
dFrame.modBneg$score <- dFrame.neg$score + dFrame.modBneg$score
dFrame.modBpos$score <- dFrame.pos$score + dFrame.modBpos$score
dFrame.imrmc <- rbind(
dFrame.truth,
dFrame.modAneg,
dFrame.modBneg,
dFrame.modApos,
dFrame.modBpos
)
dFrame.imrmc$caseID <- factor(dFrame.imrmc$caseID, as.character(caseIDs), ordered = TRUE)
return(dFrame.imrmc)
}
# Create a configuration object for the sim.gRoeMetz program ####
#' Create a configuration object for the sim.gRoeMetz program
#'
#' @description
#' This function creates a configuration object for the Roe & Metz
#' simulation model to be used as input for the sim.gRoeMetz program.
#' The default model returned when there are no arguments given to the function is
#' the "HH" model from Roe1987_Acad-Radiol_v4p298. Following that paper,
#' The user can specify three parameters related to experiment size (nR, nC.neg, nC.pos)
#' and five parameters parameters specifying a linear model that does not
#' depend on modality or truth (mu.neg, mu.pos, var_r, var_c, var_rc).
#' The mu.pos is set to 1.0, which yields a reader averaged AUC of
#' approximately 0.765.
#'
#' @details If no arguments, this function returns a default simulation
#' configuration for sim.gRoeMetz
#'
#' @param nR Number of readers (default = 5)
#' @param nC.neg Number of signal-absent cases (default = 25)
#' @param nC.pos Number of signal-present cases (default = 25)
#' @param mu.neg Mean fixed effect of signal-absent distribution (default = 0.0) \cr
#' Modality specific parameters are set to zero: mu.Aneg = mu.Bneg = 0
#' @param mu.pos Mean fixed effect of signal-present distribution (default = 1.0) \cr
#' Modality specific parameters are set to zero: mu.Apos = mu.Bpos = 0
#' @param var_r Variance of reader random effect (default = 0.03) \cr
#' var_r.neg = var_r.pos = var_r.Aneg = var_r.Apos = var_r.Bneg = var_r.Bpos = var_r \cr
#' @param var_c Variance of case random effect (default = 0.30) \cr
#' var_c.neg = var_c.pos = var_c.Aneg = var_c.Apos = var_c.Bneg = var_c.Bpos = var_c \cr
#' @param var_rc Variance of reader.by.case random effect (default = 0.20) \cr
#' var_rc.neg = var_rc.pos = var_rc.Aneg = var_rc.Apos = var_rc.Bneg = var_rc.Bpos = var_rc \cr
#'
#' @return config [list] Refer to the sim.gRoeMetz input variable
#'
#' @export
#'
# @examples
# # Create a sample configuration object
# config <- sim.gRoeMetz.config()
# # Simulate an MRMC ROC data set
# dFrame.imrmc <- sim.gRoeMetz(config)
# # Analyze the MRMC ROC data
# result <- doIMRMC(dFrame.imrmc)
sim.gRoeMetz.config <- function(
nR = 5,
nC.neg = 40,
nC.pos = 40,
mu.neg = 0.0,
mu.pos = 1.0,
var_r = 0.03,
var_c = 0.30,
var_rc = 0.20
) {
config <- list(
# Modality labels
modalityID.A = "testA",
modalityID.B = "testB",
# Experiment size
nR = nR,
nC.neg = nC.neg,
nC.pos = nC.pos,
# Model parameters invariant to modality
# Signal-absent global mean
mu.neg = mu.neg,
# Signal-absent variance components
var_r.neg = var_r,
var_c.neg = var_c,
var_rc.neg = var_rc,
# Signal-present global mean
mu.pos = mu.pos,
# Signal-present variance components
var_r.pos = var_r,
var_c.pos = var_c,
var_rc.pos = var_rc,
# Model parameters for modality A
# Signal-present modality effect
mu.Aneg = 0,
# Signal-absent variance components
var_r.Aneg = var_r,
var_c.Aneg = var_c,
var_rc.Aneg = var_rc,
# Signal-absent modality effect
mu.Apos = 0,
# Signal-present variance components
var_r.Apos = var_r,
var_c.Apos = var_c,
var_rc.Apos = var_rc,
# Model parameters for modality B
# Signal-present modality effect
mu.Bneg = 0,
# Signal-absent variance components
var_r.Bneg = var_r,
var_c.Bneg = var_c,
var_rc.Bneg = var_rc,
# Signal-absent modality effect
mu.Bpos = 0,
# Signal-present variance components
var_r.Bpos = var_r,
var_c.Bpos = var_c,
var_rc.Bpos = var_rc
)
return(config)
}
# Simulates a sample MRMC ROC experiment ####
#' Simulates a sample MRMC ROC experiment
#'
#' @return dFrame.imrmc [data.frame] Please refer to the description of the \code{\link{sim.gRoeMetz}} return variable
#' @export
#'
# @examples
# # Simulate a sample MRMC ROC data set
# dFrame.imrmc <- simRoeMetz.example()
# # Analyze the MRMC ROC data
# result <- doIMRMC(dFrame.imrmc)
#
simRoeMetz.example <- function() {
config <- sim.gRoeMetz.config()
# Simulate data
dFrame.imrmc <- sim.gRoeMetz(config)
tempA <- dFrame.imrmc[dFrame.imrmc$modalityID == "testA", ]
tempApos <- tempA[grep("pos", tempA$caseID), ]
tempAneg <- tempA[grep("neg", tempA$caseID), ]
tempB <- dFrame.imrmc[dFrame.imrmc$modalityID == "testB", ]
tempBpos <- tempB[grep("pos", tempA$caseID), ]
tempBneg <- tempB[grep("neg", tempA$caseID), ]
cat("test A pos: mean, var = ", mean(tempApos$score), ",", stats::var(tempApos$score), "\n")
cat("test B pos: mean, var = ", mean(tempBpos$score), ",", stats::var(tempBpos$score), "\n")
cat("test A neg: mean, var = ", mean(tempAneg$score), ",", stats::var(tempAneg$score), "\n")
cat("test B neg: mean, var = ", mean(tempBneg$score), ",", stats::var(tempBneg$score), "\n")
return(dFrame.imrmc)
}
## roeMetzConfigs ####
#' @title roeMetzConfigs
#'
#' @name roeMetzConfigs
#'
#' @description This is a data frame containing the configuration parameters
#' used in Roe1997_Acad-Radiol_v4p298. Each row corresponds to one of the twelve configurations
#' appearing in Table 1 of that paper in a format that can be the input to \code{sim.gRoeMetz}.
#'
#' @details The columns of this data frame are as follows
#' \itemize{
#' \item Experiment labels and size
#' \itemize{
#' \item modalityID.A: [character] label modality A
#' \item modalityID.B: [character] label modality B
#' \item nR: [numeric] number of readers
#' \item nC.neg: [numeric] number of signal-absent cases
#' \item nC.pos: [numeric] number of signal-present cases
#' }
#' \item There are six fixed effects:
#' \itemize{
#' \item mu.neg: [numeric] signal-absent (neg, global mean)
#' \item mu.pos: [numeric] signal-present (pos, global mean)
#' \item mu.Aneg: [numeric] modality A signal-absent (Aneg, modality effect)
#' \item mu.Bneg: [numeric] modality B signal-absent (Bneg, modality effect)
#' \item mu.Apos: [numeric] modality A signal-present (Apos, modality effect)
#' \item mu.Bpos: [numeric] modality B signal-present (Bpos, modality effect)
#' }
#' \item There are six random effects that are independent of modality
#' \itemize{
#' \item var_r.neg: [numeric] variance of random reader effect
#' \item var_c.neg: [numeric] variance of random case effect
#' \item var_rc.neg: [numeric] variance of random reader by case effect
#' \item var_r.pos: [numeric] variance of random reader effect
#' \item var_c.pos: [numeric] variance of random case effect
#' \item var_rc.pos: [numeric] variance of random reader by case effect
#' }
#' \item There are six random effects that are specific to modality A
#' \itemize{
#' \item var_r.Aneg: [numeric] variance of random reader effect
#' \item var_c.Aneg: [numeric] variance of random case effect
#' \item var_rc.Aneg: [numeric] variance of random reader by case effect
#' \item var_r.Apos: [numeric] variance of random reader effect
#' \item var_c.Apos: [numeric] variance of random case effect
#' \item var_rc.Apos: [numeric] variance of randome reader by case effect
#' }
#' \item There are six random effects that are specific to modality B
#' \itemize{
#' \item var_r.Bneg: [numeric] variance of random reader effect
#' \item var_c.Bneg: [numeric] variance of random case effect
#' \item var_rc.Bneg: [numeric] variance of random reader by case effect
#' \item var_r.Bpos: [numeric] variance of random reader effect
#' \item var_c.Bpos: [numeric] variance of random case effect
#' \item var_rc.Bpos: [numeric] variance of randome reader by case effect
#' }
#' }
#'
# @examples
# # Extract one configuration
# config <- roeMetzConfigs[1, ]
# # Create an MRMC ROC data set
# df.iMRMC <- sim.gRoeMetz(config)
# # Analyze the MRMC ROC data set
# result <- doIMRMC(df.iMRMC)
NULL
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.