Wrapper function for mclapply

# S3 method for default
sim(
  x = NULL,
  R = 100,
  f = NULL,
  colnames = NULL,
  messages = lava.options()$messages,
  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

messages

Messages

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

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,messages=0)
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.043s
#> 
#>            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.043s
#> 
#>             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.043s
#> 
#>               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