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).

# S3 method for default
sim(
  x = NULL,
  R = 100,
  f = NULL,
  colnames = NULL,
  seed = NULL,
  args = list(),
  iter = FALSE,
  ...
)

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

...

Additional arguments to (mc)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 print.sim

Examples

m <- lvm(y~x+e)
distribution(m,~y) <- 0
distribution(m,~x) <- uniform.lvm(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,],
             estimate(l,only.coef=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
#>    Estimate Model.se Model.lo Model.hi Sandwich.se Sandwich.lo Sandwich.hi
#> 1  0.928718 0.033366 0.863282 0.994155 0.047050    0.836502    1.020935   
#> 2  0.995991 0.001633 0.992788 0.999193 0.002255    0.991571    1.000411   
#> 3  0.998523 0.003280 0.992090 1.004955 0.004726    0.989261    1.007785   
#> 4  1.001350 0.004145 0.993220 1.009480 0.005812    0.989960    1.012741   
#> 5  0.999646 0.002153 0.995424 1.003869 0.003048    0.993673    1.005620   
#> 6  1.004597 0.006829 0.991204 1.017990 0.009495    0.985987    1.023207   
#> 7  0.998076 0.009524 0.979397 1.016754 0.013488    0.971640    1.024511   
#> 8  0.997917 0.003406 0.991237 1.004596 0.004853    0.988405    1.007428   
#> 9  0.999195 0.001333 0.996582 1.001809 0.001906    0.995459    1.002932   
#> 10 1.000199 0.000730 0.998768 1.001631 0.001023    0.998195    1.002204   
#> 
#>      Estimate  Model.se Model.lo  Model.hi Sandwich.se Sandwich.lo Sandwich.hi
#> Mean 0.992421 0.0066400 0.979399 1.0054433   0.0093655    0.974065   1.0107773
#> SD   0.022501 0.0097691 0.041125 0.0074485   0.0137735    0.048867   0.0090643

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))
#> 50 replications					Time: 1.056s
#> 
#>            Estimate Estimate.1
#> Mean      0.9941568  0.9941568
#> SD        0.0152271  0.0152271
#> Coverage  0.8800000  0.9800000
#>                               
#> Min       0.9287185  0.9287185
#> 2.5%      0.9546727  0.9546727
#> 50%       0.9984014  0.9984014
#> 97.5%     1.0140017  1.0140017
#> Max       1.0169750  1.0169750
#>                               
#> Missing   0.0000000  0.0000000
#>                               
#> True      1.0000000  1.0000000
#> Bias     -0.0058432 -0.0058432
#> RMSE      0.0163097  0.0163097
#> 

summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))
#> 50 replications					Time: 1.056s
#> 
#>             Model Sandwich
#> Mean    0.9941568 0.994157
#> SD      0.0152271 0.015227
#> SE      0.0084199 0.011813
#> SE/SD   0.5529578 0.775782
#>                           
#> Min     0.9287185 0.928718
#> 2.5%    0.9546727 0.954673
#> 50%     0.9984014 0.998401
#> 97.5%   1.0140017 1.014002
#> Max     1.0169750 1.016975
#>                           
#> Missing 0.0000000 0.000000
#> 
summary(val,estimate=c(1,1),se=c(2,5),true=c(1,1),names=c("Model","Sandwich"),confint=TRUE)
#> 50 replications					Time: 1.056s
#> 
#>               Model   Sandwich
#> Mean      0.9941568  0.9941568
#> SD        0.0152271  0.0152271
#> SE        0.0084199  0.0118129
#> SE/SD     0.5529578  0.7757816
#> Coverage  0.8800000  0.9800000
#>                               
#> Min       0.9287185  0.9287185
#> 2.5%      0.9546727  0.9546727
#> 50%       0.9984014  0.9984014
#> 97.5%     1.0140017  1.0140017
#> Max       1.0169750  1.0169750
#>                               
#> Missing   0.0000000  0.0000000
#>                               
#> True      1.0000000  1.0000000
#> Bias     -0.0058432 -0.0058432
#> RMSE      0.0163097  0.0163097
#> 

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)
#>   [,1] [,2] [,3] [,4] [,5]
#> 1 1    1    1    1    1   
#> 2 2    2    2    2    2   
#> 3 3    3    3    3    3   
#> 4 2    2    2    2    2   
#> 5 4    4    4    4    4   
#> 6 6    6    6    6    6   
#> 7 3    3    3    3    3   
#> 8 6    6    6    6    6   
#> 9 9    9    9    9    9   
#> 
#>        [,1]   [,2]   [,3]   [,4]   [,5]
#> Mean 4.0000 4.0000 4.0000 4.0000 4.0000
#> SD   2.5495 2.5495 2.5495 2.5495 2.5495
sim(function(a,b) f(a,b), 3, args=c(a=5,b=5))
#>   [,1] [,2] [,3] [,4] [,5]
#> 1 25   25   25   25   25  
#> 2 25   25   25   25   25  
#> 3 25   25   25   25   25  
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> Mean   25   25   25   25   25
#> SD      0    0    0    0    0
sim(function(iter=1,a=5,b=5) iter*f(a,b), iter=TRUE, R=5)
#>   [,1] [,2] [,3] [,4] [,5]
#> 1  25   25   25   25   25 
#> 2  50   50   50   50   50 
#> 3  75   75   75   75   75 
#> 4 100  100  100  100  100 
#> 5 125  125  125  125  125 
#> 
#>        [,1]   [,2]   [,3]   [,4]   [,5]
#> Mean 75.000 75.000 75.000 75.000 75.000
#> SD   39.528 39.528 39.528 39.528 39.528