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
# Default S3 method
sim(
x = NULL,
R = 100,
f = NULL,
colnames = NULL,
seed = NULL,
args = list(),
iter = FALSE,
mc.cores,
progressr.message = NULL,
...
)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
- ...
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).
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.981936 0.010066 0.962195 1.001677 0.014230 0.954045 1.009827
#> 2 0.985010 0.018630 0.948473 1.021546 0.026497 0.933076 1.036943
#> 3 0.993945 0.006745 0.980717 1.007173 0.009338 0.975642 1.012248
#> 4 0.987181 0.009413 0.968721 1.005640 0.013156 0.961396 1.012966
#> 5 0.990049 0.006639 0.977028 1.003069 0.009411 0.971603 1.008495
#> 6 0.980337 0.019269 0.942548 1.018127 0.027398 0.926638 1.034036
#> 7 0.997944 0.001032 0.995920 0.999968 0.001460 0.995082 1.000806
#> 8 1.000065 0.001981 0.996179 1.003951 0.002797 0.994583 1.005546
#> 9 1.006180 0.006264 0.993894 1.018465 0.008815 0.988903 1.023456
#> 10 1.002229 0.006264 0.989943 1.014514 0.008572 0.985427 1.019030
#>
#> Estimate Model.se Model.lo Model.hi Sandwich.se Sandwich.lo Sandwich.hi
#> Mean 0.9924875 0.0086305 0.975562 1.0094131 0.0121675 0.968640 1.016335
#> SD 0.0089484 0.0061147 0.019608 0.0079547 0.0087226 0.024533 0.011951
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: 0.974s
#>
#> Estimate Estimate.1
#> Mean 0.9961265 0.9961265
#> SD 0.0138953 0.0138953
#> Coverage 0.8200000 0.9800000
#>
#> Min 0.9673948 0.9673948
#> 2.5% 0.9722390 0.9722390
#> 50% 0.9976991 0.9976991
#> 97.5% 1.0298590 1.0298590
#> Max 1.0402295 1.0402295
#>
#> Missing 0.0000000 0.0000000
#>
#> True 1.0000000 1.0000000
#> Bias -0.0038735 -0.0038735
#> RMSE 0.0144251 0.0144251
#>
summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))
#> 50 replications Time: 0.974s
#>
#> Model Sandwich
#> Mean 0.9961265 0.996127
#> SD 0.0138953 0.013895
#> SE 0.0091896 0.012960
#> SE/SD 0.6613435 0.932713
#>
#> Min 0.9673948 0.967395
#> 2.5% 0.9722390 0.972239
#> 50% 0.9976991 0.997699
#> 97.5% 1.0298590 1.029859
#> Max 1.0402295 1.040229
#>
#> 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: 0.974s
#>
#> Model Sandwich
#> Mean 0.9961265 0.9961265
#> SD 0.0138953 0.0138953
#> SE 0.0091896 0.0129603
#> SE/SD 0.6613435 0.9327125
#> Coverage 0.8200000 0.9800000
#>
#> Min 0.9673948 0.9673948
#> 2.5% 0.9722390 0.9722390
#> 50% 0.9976991 0.9976991
#> 97.5% 1.0298590 1.0298590
#> Max 1.0402295 1.0402295
#>
#> Missing 0.0000000 0.0000000
#>
#> True 1.0000000 1.0000000
#> Bias -0.0038735 -0.0038735
#> RMSE 0.0144251 0.0144251
#>
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