# ----------------------
# Author: Andreas Alfons
# KU Leuven
# ----------------------
## for convenience: construct "SimControl" object and re-call function
setMethod("runSimulation",
signature(x = "ANY", setup = "ANY", nrep = "ANY", control = "missing"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
control <- SimControl(contControl=contControl, NAControl=NAControl,
design=design, fun=fun, dots=list(...), SAE=SAE)
if(missing(setup)) {
if(missing(nrep)) runSimulation(x, control=control)
else runSimulation(x, nrep=nrep, control=control)
} else {
if(missing(nrep)) runSimulation(x, setup, control=control)
else runSimulation(x, setup, nrep, control)
}
})
## design-based simulation
setMethod("runSimulation",
signature(x = "data.frame", setup = "VirtualSampleControl",
nrep = "missing", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
setup <- setup(x, setup)
runSimulation(x, setup, control=control)
})
setMethod("runSimulation",
signature(x = "data.frame", setup = "SampleSetup",
nrep = "missing", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
# initializations
nsam <- length(setup)
design <- getDesign(control)
if(nsam == 0) { # nothing to do
return(SimResults(design=design, sampleControl=getControl(setup),
control=control))
}
contControl <- getContControl(control)
epsilon <- if(is.null(contControl)) numeric() else getEpsilon(contControl)
NAControl <- getNAControl(control)
NArate <- if(is.null(NAControl)) numeric() else getNArate(NAControl)
# run the simulations
s <- 1:nsam
tmp <- lapply(s, designSimulation, x, setup, control)
# construct results
getSimResults(tmp, s, epsilon=epsilon, NArate=NArate, design=design,
sampleControl=getControl(setup), control=control)
})
## model-based simulation
setMethod("runSimulation",
signature(x = "VirtualDataControl", setup = "missing",
nrep = "numeric", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
# initializations
if(length(nrep) == 0) stop("'nrep' must be a non-negative integer")
else if(length(nrep) > 1) nrep <- nrep[1]
design <- getDesign(control)
if(nrep == 0) { # nothing to do
return(SimResults(design=design, dataControl=x,
nrep=nrep, control=control))
}
contControl <- getContControl(control)
epsilon <- if(is.null(contControl)) numeric() else getEpsilon(contControl)
NAControl <- getNAControl(control)
NArate <- if(is.null(NAControl)) numeric() else getNArate(NAControl)
# run the simulations
r <- 1:nrep
tmp <- lapply(r, modelSimulation, x, control)
# construct results
getSimResults(tmp, reps=r, epsilon=epsilon, NArate=NArate,
design=design, dataControl=x, control=control)
})
## mixed simulation designs
setMethod("runSimulation",
signature(x = "VirtualDataControl", setup = "VirtualSampleControl",
nrep = "numeric", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
# initializations
if(length(nrep) == 0) stop("'nrep' must be a non-negative integer")
else if(length(nrep) > 1) nrep <- nrep[1]
nsam <- length(setup)
design <- getDesign(control)
if(nrep == 0 || nsam == 0) { # nothing to do
return(SimResults(design=design, dataControl=x,
sampleControl=setup, nrep=nrep, control=control))
}
contControl <- getContControl(control)
epsilon <- if(is.null(contControl)) numeric() else getEpsilon(contControl)
NAControl <- getNAControl(control)
NArate <- if(is.null(NAControl)) numeric() else getNArate(NAControl)
# run the simulations (generate data repeatedly and draw samples)
r <- 1:nrep
tmp <- lapply(r, mixedSimulation, x, setup, control)
# construct results
getSimResults(tmp, 1:nsam, r, epsilon=epsilon, NArate=NArate,
design=design, dataControl=x, sampleControl=setup, control=control)
})
## simulation with repetitions based on (possibly) real data
setMethod("runSimulation",
signature(x = "data.frame", setup = "missing",
nrep = "numeric", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
# initializations
if(length(nrep) == 0) stop("'nrep' must be a non-negative integer")
else if(length(nrep) > 1) nrep <- nrep[1]
design <- getDesign(control)
if(nrep == 0) { # nothing to do
return(SimResults(design=design, nrep=nrep, control=control))
}
contControl <- getContControl(control)
epsilon <- if(is.null(contControl)) numeric() else getEpsilon(contControl)
NAControl <- getNAControl(control)
NArate <- if(is.null(NAControl)) numeric() else getNArate(NAControl)
SAE <- isTRUE(getSAE(control))
# get results (adjustments are needed for parallel computing)
if(length(design)) {
spl <- getStrataSplit(x, design, USE.NAMES=FALSE)
leg <- getStrataLegend(x, design)
if(SAE) {
tmp <- replicate(nrep,
manageSimulationSAE(x, spl, control, leg),
simplify=FALSE)
} else {
xSpl <- lapply(spl, function(s, x) x[s, , drop=FALSE], x)
tmp <- replicate(nrep,
manageSimulationStrata(xSpl, spl, control, leg),
simplify=FALSE)
}
} else {
tmp <- replicate(nrep, manageSimulation(x, control),
simplify=FALSE)
}
# construct results
getSimResults(tmp, reps=1:nrep, epsilon=epsilon,
NArate=NArate, design=design, control=control)
})
## no samples, no repetitions (probably useless, but for completeness)
setMethod("runSimulation",
signature(x = "data.frame", setup = "missing",
nrep = "missing", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ...,
SAE = FALSE) {
runSimulation(x, nrep=1, control=control)
})
setMethod("runSimulation",
signature(x = "VirtualDataControl", setup = "missing",
nrep = "missing", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ..., SAE = FALSE) {
runSimulation(x, nrep=1, control=control)
})
setMethod("runSimulation",
signature(x = "VirtualDataControl", setup = "VirtualSampleControl",
nrep = "missing", control = "SimControl"),
function(x, setup, nrep, control, contControl = NULL,
NAControl = NULL, design = character(), fun, ..., SAE = FALSE) {
runSimulation(x, setup, nrep=1, control=control)
})
## these functions need to be available in the namespace for parallel computing
designSimulation <- function(i, x, setup, control) {
# print(paste(Sys.time(), i, sep=": "))
sam <- drawS3(x, setup, i)
if(nrow(sam) == 0) return(getEmptyResults(control))
design <- getDesign(control)
SAE <- isTRUE(getSAE(control))
if(length(design)) {
spl <- getStrataSplit(sam, design, USE.NAMES=FALSE)
leg <- getStrataLegend(sam, design)
if(SAE) manageSimulationSAE(sam, spl, control, leg)
else {
samSpl <- lapply(spl, function(s, x) x[s, , drop=FALSE], sam)
manageSimulationStrata(samSpl, spl, control, leg)
}
} else manageSimulation(sam, control)
}
# 'i' as first argument is necessary for parallel computing with 'parLapply'
modelSimulation <- function(i, x, control) {
md <- try(generate(x))
if(class(md) == "try-error") return(getEmptyResults(control))
design <- getDesign(control)
SAE <- isTRUE(getSAE(control))
if(length(design)) {
spl <- getStrataSplit(md, design, USE.NAMES=FALSE)
leg <- getStrataLegend(md, design)
if(SAE) manageSimulationSAE(md, spl, control, leg)
else {
mdSpl <- lapply(spl, function(s, x) x[s, , drop=FALSE], md)
manageSimulationStrata(mdSpl, spl, control, leg)
}
} else manageSimulation(md, control)
}
mixedSimulation <- function(i, x, setup, control) {
nsam <- length(setup)
# generate data
md <- try(generate(x))
if(class(md) == "try-error") {
# return empty results for each sample
return(replicate(nsam, getEmptyResults(control), simplify=FALSE))
}
# set up samples (empty vector in case of error in one iteration)
setup <- setup(md, setup)
# run design-based simulations on the generated data with the set up samples
s <- 1:nsam
lapply(s, designSimulation, md, setup, control)
}
manageSimulation <- function(x, control) {
# initializations
contControl <- getContControl(control)
neps <- length(contControl)
NAControl <- getNAControl(control)
nNA <- length(NAControl)
fun <- getFun(control)
useOrig <- "orig" %in% argNames(fun)
dots <- getDots(control)
# get results
if(neps) {
if(nNA) {
# contamination, missings
tmp <- lapply(1:neps,
function(e) {
cx <- try(contaminate(x, contControl, e))
if(class(cx) == "try-error") return(vector("list", nNA))
lapply(1:nNA,
function(n) try({
nx <- setNA(cx, NAControl, n)
ca <- as.call(c(fun, dots))
ca$x <- nx
if(useOrig) ca$orig <- x
getSimResult(eval(ca))
}))
})
do.call("c", tmp)
} else {
# contamination, no missings
lapply(1:neps,
function(e) try({
cx <- contaminate(x, contControl, e)
ca <- as.call(c(fun, dots))
ca$x <- cx
if(useOrig) ca$orig <- x
getSimResult(eval(ca))
}))
}
} else {
if(nNA) {
# no contamination, missings
lapply(1:nNA,
function(n) try({
nx <- setNA(x, NAControl, n)
ca <- as.call(c(fun, dots))
ca$x <- nx
if(useOrig) ca$orig <- x
getSimResult(eval(ca))
}))
} else {
# no contamination, no missings
try({
ca <- as.call(c(fun, dots))
ca$x <- x
if(useOrig) ca$orig <- x
getSimResult(eval(ca))
})
}
}
}
manageSimulationStrata <- function(xs, indices, control, legend) {
# initializations
contControl <- getContControl(control)
neps <- length(contControl)
NAControl <- getNAControl(control)
nNA <- length(NAControl)
fun <- getFun(control)
nam <- argNames(fun)
useOrig <- "orig" %in% nam
useDomain <- "domain" %in% nam
dots <- getDots(control)
# get results
if(neps) {
if(nNA) {
# contamination, missings
tmp <- lapply(1:neps,
function(e) {
cxs <- try(lapply(xs, contaminate, contControl, e))
if(class(cxs) == "try-error") return(vector("list", nNA))
lapply(1:nNA,
function(n) {
try({
tmp <- mapply(
function(cx, x, i) {
nx <- setNA(cx, NAControl, n)
ca <- as.call(c(fun, dots))
ca$x <- nx
if(useOrig) ca$orig <- x
if(useDomain) ca$domain <- i
getSimResult(eval(ca))
}, cxs, xs, indices,
SIMPLIFY=FALSE, USE.NAMES=FALSE)
getSimResultStrata(tmp, legend)
})
})
})
do.call("c", tmp)
} else {
# contamination, no missings
lapply(1:neps,
function(e) {
try({
tmp <- mapply(
function(x, i) {
cx <- contaminate(x, contControl, e)
ca <- as.call(c(fun, dots))
ca$x <- cx
if(useOrig) ca$orig <- x
if(useDomain) ca$domain <- i
getSimResult(eval(ca))
}, xs, indices,
SIMPLIFY=FALSE, USE.NAMES=FALSE)
getSimResultStrata(tmp, legend)
})
})
}
} else {
if(nNA) {
# no contamination, missings
lapply(1:nNA,
function(n) {
try({
tmp <- mapply(
function(x, i) {
nx <- setNA(x, NAControl, n)
ca <- as.call(c(fun, dots))
ca$x <- nx
if(useOrig) ca$orig <- x
if(useDomain) ca$domain <- i
getSimResult(eval(ca))
}, xs, indices,
SIMPLIFY=FALSE, USE.NAMES=FALSE)
getSimResultStrata(tmp, legend)
})
})
} else {
# no contamination, no missings
try({
tmp <- mapply(
function(x, i) {
ca <- as.call(c(fun, dots))
ca$x <- x
if(useOrig) ca$orig <- x
if(useDomain) ca$domain <- i
getSimResult(eval(ca))
}, xs, indices, SIMPLIFY=FALSE, USE.NAMES=FALSE)
getSimResultStrata(tmp, legend)
})
}
}
}
manageSimulationSAE <- function(x, indices, control, legend) {
# initializations
contControl <- getContControl(control)
neps <- length(contControl)
NAControl <- getNAControl(control)
nNA <- length(NAControl)
fun <- getFun(control)
useOrig <- "orig" %in% argNames(fun)
dots <- getDots(control)
# get results
if(neps) {
if(nNA) {
# contamination, missings
tmp <- lapply(1:neps,
function(e) {
cx <- try(contaminate(x, contControl, e))
if(class(cx) == "try-error") return(vector("list", nNA))
lapply(1:nNA,
function(n) {
try({
nx <- setNA(cx, NAControl, n)
ca <- as.call(c(fun, dots))
ca$x <- nx
if(useOrig) ca$orig <- x
tmp <- lapply(indices,
function(i) {
ca$domain <- i
getSimResult(eval(ca))
})
getSimResultStrata(tmp, legend)
})
})
})
do.call("c", tmp)
} else {
# contamination, no missings
lapply(1:neps,
function(e) {
try({
cx <- contaminate(x, contControl, e)
ca <- as.call(c(fun, dots))
ca$x <- cx
if(useOrig) ca$orig <- x
tmp <- lapply(indices,
function(i) {
ca$domain <- i
getSimResult(eval(ca))
})
getSimResultStrata(tmp, legend)
})
})
}
} else {
if(nNA) {
# no contamination, missings
lapply(1:nNA,
function(n) {
try({
nx <- setNA(x, NAControl, n)
ca <- as.call(c(fun, dots))
ca$x <- nx
if(useOrig) ca$orig <- x
tmp <- lapply(indices,
function(i) {
ca$domain <- i
getSimResult(eval(ca))
})
getSimResultStrata(tmp, legend)
})
})
} else {
# no contamination, no missings
try({
ca <- as.call(c(fun, dots))
ca$x <- x
if(useOrig) ca$orig <- x
tmp <- lapply(indices,
function(i) {
ca$domain <- i
getSimResult(eval(ca))
})
getSimResultStrata(tmp, legend)
})
}
}
}
## convenience wrapper
runSim <- function(...) {
res <- runSimulation(...)
call <- match.call()
setCall(res, call)
res
}
## utilities
# check for errors or empty vectors in a list (for simulation results)
checkOK <- function(x) {
# x ... list
if(length(x)) {
sapply(x, function(x) class(x) != "try-error" && length(x$values))
} else logical()
}
## get empty results
getEmptyResults <- function(control) {
neps <- length(getContControl(control))
if(neps == 0) neps <- 1
nNA <- length(getNAControl(control))
if(nNA == 0) nNA <- 1
replicate(neps*nNA, list(values=numeric()))
}
# get result of one simulation run in the correct format
getSimResult <- function(x) {
if(is(x, "numeric")) x <- list(values=x)
else if(is(x, "list")) {
if(!setequal(names(x), c("values", "add"))) {
stop("the components of the list returned by ",
"'fun' must have names 'values' and 'add'")
}
# } else if(is(x, "SimResult")) {
# x <- list(values=getValues(x), add=getAdd(x))
} else {
# stop("'fun' must return a numeric vector ",
# "or an object of class \"SimResult\"")
stop("'fun' must return a numeric vector or a ",
"list with components 'values' and 'add'")
}
x
}
# get result for one stratified simulation run
getSimResultStrata <- function(x, legend) {
values <- lapply(x, function(x) x$values)
values <- do.call("rbind", values)
values <- cbind(legend, values)
add <- lapply(x, function(x) x$add)
list(values=values, add=add)
}
# contruct object to be returned
getSimResults <- function(x, samples = numeric(), reps = numeric(),
epsilon = numeric(), NArate = numeric(), design = character(),
dataControl = NULL, sampleControl = NULL, control) {
nsam <- length(samples)
.nrep <- if(missing(reps)) numeric() else length(reps)
nrep <- length(reps)
neps <- length(epsilon)
origNArate <- NArate
NArate <- convertNArate(NArate)
nNA <- length(NArate)
if(is.na(nNA)) nNA <- 0
# combine results from all runs into one list
if(nsam && nrep) x <- do.call("c", x)
if(neps || nNA) x <- do.call("c", x)
# check for errors or empty results
nruns <- length(x)
ok <- checkOK(x)
x <- x[ok]
if(length(x) == 0) stop("error or empty result in every simulation run")
# get additional information (at least one of 'nrep' or 'nsam' is positive)
ca <- call("expand.grid") # initialize call
if(nNA) ca$NArate <- NArate
if(neps) ca$Epsilon <- epsilon
if(nsam) ca$Sample <- samples
if(nrep) ca$Rep <- reps
info <- eval(ca) # create data.frame with additional information
info <- info[, rev(1:ncol(info)), drop=FALSE] # reverse column order
info <- cbind(Run=1:nruns, info)[ok, , drop=FALSE] # add runs
# elements of 'tmp' are currently lists (faster than S4 objects)
values <- lapply(x, function(x) x$values)
if(length(design)) {
# additional information needs to be adjusted
reps <- sapply(values, getRepetitions)
# info <- info[rep(1:nrow(info), each=reps), , drop=FALSE]
info <- info[rep.int(1:nrow(info), reps), , drop=FALSE]
}
values <- do.call("rbind", values)
# put it all together
values <- cbind(info, values)
ninfo <- ncol(info) + length(design)
nam <- names(values)[-(1:ninfo)]
rownames(values) <- NULL
add <- lapply(x, function(x) x$add)
# return results
SimResults(values=values, add=add, design=design, colnames=nam,
epsilon=epsilon, NArate=origNArate, dataControl=dataControl,
sampleControl=sampleControl, nrep=.nrep, control=control)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.