lava (version 1.9.3)

sim.default: Monte Carlo simulation

Description

Applies a function repeatedly for a specified number of replications or over a list/data.frame with plot and summary methods for summarizing the Monte Carlo experiment. Can be parallelized via the future package (use the future::plan() function).

Usage

# S3 method for default
sim(
  x = NULL,
  R = 100,
  f = NULL,
  colnames = NULL,
  seed = NULL,
  args = list(),
  iter = FALSE,
  mc.cores,
  progressr.message = NULL,
  estimate.index = 1:2,
  ...
)

Arguments

x

function or 'sim' object

R

Number of replications or data.frame with parameters

f

Optional function (i.e., if x is a matrix)

colnames

Optional column names

seed

(optional) Seed (needed with cl=TRUE)

args

(optional) list of named arguments passed to (mc)mapply

iter

If TRUE the iteration number is passed as first argument to (mc)mapply

mc.cores

Optional number of cores. Will use parallel::mcmapply instead of future

progressr.message

Optional message for the progressr progress-bar

estimate.index

If return object inherits from estimate then only these column indices are extracted (estimate, se, lower, upper, p-val)

...

Additional arguments to future.apply::future_mapply()

Details

To parallelize the calculation use the future::plan() function (e.g., future::plan(multisession()) to distribute the calculations over the R replications on all available cores). The output is controlled via the progressr package (e.g., progressr::handlers(global=TRUE) to enable progress information).

See Also

summary.sim() plot.sim() sim.lvm()

Examples

m <- lvm(y~x+e)
distribution(m,~y) <- 0
distribution(m,~x) <- dist_uniform(a=-1.1,b=1.1)
transform(m,e~x) <- function(x) (1*x^4)*rnorm(length(x),sd=1)

onerun <- function(iter=NULL,...,n=2e3,b0=1,idx=2) {
    d <- sim(m,n,p=c("y~x"=b0))
    l <- lm(y~x,d)
    res <- c(coef(summary(l))[idx,1:2],
             confint(l)[idx,],
             coef(estimate(l), mat=TRUE)[idx,2:4])
    names(res) <- c("Estimate","Model.se","Model.lo","Model.hi",
                    "Sandwich.se","Sandwich.lo","Sandwich.hi")
    res
}
val <- sim(onerun,R=10,b0=1)
val

val <- sim(val,R=40,b0=1) ## append results
summary(val,estimate=c(1,1),confint=c(3,4,6,7),true=c(1,1))

summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))
summary(val,estimate=c(1,1),se=c(2,5),true=c(1,1),
        names=c("Model","Sandwich"),confint=TRUE)

if (interactive()) {
    plot(val,estimate=1,c(2,5),true=1,
         names=c("Model","Sandwich"),polygon=FALSE)
    plot(val,estimate=c(1,1),se=c(2,5),main=NULL,
         true=c(1,1),names=c("Model","Sandwich"),
         line.lwd=1,col=c("gray20","gray60"),
         rug=FALSE)
    plot(val,estimate=c(1,1),se=c(2,5),true=c(1,1),
         names=c("Model","Sandwich"))
}

f <- function(a=1, b=1) {
  rep(a*b, 5)
}
R <- Expand(a=1:3, b=1:3)
sim(f, R)
sim(function(a,b) f(a,b), 3, args=c(a=5,b=5))
sim(function(iter=1,a=5,b=5) iter*f(a,b), iter=TRUE, R=5)

## Returning estimate objects with extra per-iteration information
onerun2 <- function(...) {
  x <- rnorm(50)
  y <- 0.5*x + rnorm(50)
  e <- estimate(lm(y ~ x))
  c(e, converged = 1, niter = sample(5:20, 1), drop.ic = TRUE)
  # works identically with summary.estimate objects
  # c(summary(e), converged = 1, niter = sample(5:20, 1), drop.ic = TRUE)
}
val2 <- sim(onerun2, R = 10)
val2
# extra columns are stored but excluded from default summary statistics
idx <- c("(Intercept)", "x", "niter")
summary(val2, estimate=idx, true=c(0,0.5,NA))