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).
function or 'sim' object
Number of replications or data.frame with parameters
Optional function (i.e., if x is a matrix)
Optional column names
(optional) Seed (needed with cl=TRUE)
(optional) list of named arguments passed to (mc)mapply
If TRUE the iteration number is passed as first argument to (mc)mapply
Optional number of cores. Will use parallel::mcmapply instead of future
Optional message for the progressr progress-bar
Additional arguments to future.apply::future_mapply
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).
summary.sim plot.sim print.sim
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.9950228 0.0134733 0.9685996 1.0214461 0.0192726 0.9572492 1.0327965
#> 2 0.9951989 0.0194513 0.9570519 1.0333459 0.0275779 0.9411471 1.0492506
#> 3 1.0008989 0.0049036 0.9912822 1.0105157 0.0068589 0.9874557 1.0143422
#> 4 1.0112677 0.0176211 0.9767101 1.0458254 0.0253846 0.9615148 1.0610207
#> 5 1.0154954 0.0124853 0.9910097 1.0399810 0.0176186 0.9809635 1.0500272
#> 6 1.0065575 0.0062972 0.9942077 1.0189074 0.0088908 0.9891319 1.0239832
#> 7 1.0095829 0.0056597 0.9984834 1.0206825 0.0080439 0.9938171 1.0253488
#> 8 1.0006804 0.0018267 0.9970979 1.0042630 0.0025955 0.9955934 1.0057675
#> 9 0.9922563 0.0046011 0.9832329 1.0012798 0.0064526 0.9796095 1.0049032
#> 10 0.9992298 0.0008412 0.9975800 1.0008797 0.0011707 0.9969354 1.0015243
#>
#> Estimate Model.se Model.lo Model.hi Sandwich.se Sandwich.lo Sandwich.hi
#> Mean 1.0026191 0.0087161 0.985526 1.019713 0.0123866 0.978342 1.02690
#> SD 0.0077671 0.0065618 0.014008 0.015989 0.0093869 0.018863 0.02102
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.081s
#>
#> Estimate Estimate.1
#> Mean 0.9961042 0.9961042
#> SD 0.0173560 0.0173560
#> Coverage 0.8200000 0.9200000
#>
#> Min 0.9240902 0.9240902
#> 2.5% 0.9576494 0.9576494
#> 50% 0.9993181 0.9993181
#> 97.5% 1.0198415 1.0198415
#> Max 1.0240726 1.0240726
#>
#> Missing 0.0000000 0.0000000
#>
#> True 1.0000000 1.0000000
#> Bias -0.0038958 -0.0038958
#> RMSE 0.0177878 0.0177878
#>
summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))
#> 50 replications Time: 1.081s
#>
#> Model Sandwich
#> Mean 0.9961042 0.996104
#> SD 0.0173560 0.017356
#> SE 0.0099984 0.014101
#> SE/SD 0.5760808 0.812483
#>
#> Min 0.9240902 0.924090
#> 2.5% 0.9576494 0.957649
#> 50% 0.9993181 0.999318
#> 97.5% 1.0198415 1.019841
#> Max 1.0240726 1.024073
#>
#> 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.081s
#>
#> Model Sandwich
#> Mean 0.9961042 0.9961042
#> SD 0.0173560 0.0173560
#> SE 0.0099984 0.0141014
#> SE/SD 0.5760808 0.8124830
#> Coverage 0.8200000 0.9200000
#>
#> Min 0.9240902 0.9240902
#> 2.5% 0.9576494 0.9576494
#> 50% 0.9993181 0.9993181
#> 97.5% 1.0198415 1.0198415
#> Max 1.0240726 1.0240726
#>
#> Missing 0.0000000 0.0000000
#>
#> True 1.0000000 1.0000000
#> Bias -0.0038958 -0.0038958
#> RMSE 0.0177878 0.0177878
#>
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