simanalyse is an R package to facilitate model comparisons and simulation studies.
To install the latest development version from GitHub
#install.packages("remotes")
remotes::install_github("audrey-b/simanalyse")Simulate 5 datasets using the sims package (here we use only a small number of datasets for the sake of illustration).
library(simanalyse)
#> Loading required package: nlist
#> Registered S3 method overwritten by 'rjags':
#> method from
#> as.mcmc.list.mcarray mcmcr
set.seed(10L)
params <- list(sigma = 2)
constants <- list(mu = 0)
code <- "for(i in 1:10){
y[i] ~ dnorm(mu, 1/sigma^2)}"
sims <- sims::sims_simulate(code,
parameters = params,
constants = constants,
nsims = 5,
silent = TRUE)
print(sims)
#> $y
#> [1] 1.29495655 -0.63919833 0.07602842 -1.55116546 0.89066792
#> [6] -0.82298676 1.03237779 1.69966601 -1.33777215 -0.77232720
#>
#> $mu
#> [1] 0
#>
#> an nlists object of 5 nlist objects each with 2 natomic elementsAnalyse the 5 datasets in “report” mode. This mode runs iterations until convergence, based on r.hat >1.05 and an effective sample size >400. See ?sma_set_mode for other choices of analysis mode.
prior <- "sigma ~ dunif(0, 6)"
results <- sma_analyse_bayesian(sims = sims,
code = code,
code.add = prior,
mode = sma_set_mode("report"))
#> module dic loaded
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Module dic unloadedDerive posterior samples for new parameters.
results.derived <- sma_derive(results, "var=sigma^2", monitor="var")
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
print(results.derived)
#> $mcmcr1
#> $var
#> [1] 2.656139
#>
#> nchains: 3
#> niters: 4000
#>
#>
#> $mcmcr2
#> $var
#> [1] 3.37433
#>
#> nchains: 3
#> niters: 4000
#>
#>
#> $mcmcr3
#> $var
#> [1] 5.681759
#>
#> nchains: 3
#> niters: 4000
#>
#>
#> $mcmcr4
#> $var
#> [1] 2.50968
#>
#> nchains: 3
#> niters: 4000
#>
#>
#> $mcmcr5
#> $var
#> [1] 7.397053
#>
#> nchains: 3
#> niters: 4000The same transformation must be applied to the true parameter values for eventually evaluating the performance (e.g. bias) of the method for those new parameters,
params.derived <- sma_derive(params, "var=sigma^2", monitor="var")
print(params.derived)
#> $var
#> [1] 4
#>
#> nchains: 1
#> niters: 1Evaluate the performance of the model using the 3 analyses
sma_evaluate(results.derived, parameters=params.derived)
#> $bias.var
#> , , 1
#>
#> [,1]
#> [1,] 1.145122
#>
#>
#> $cp.quantile.var
#> , , 1
#>
#> [,1]
#> [1,] 1
#>
#>
#> $mse.var
#> , , 1
#>
#> [,1]
#> [1,] 6.439299
#>
#>
#> an nlist object with 3 natomic elementsSeveral more performance measures are available and can be specified using the measures argument (see ?sma_evaluate for details). You may also create custom performance measures. The example below shows how to reproduce the results above with custom code.
sma_evaluate(results.derived,
measures = "",
parameters = params.derived,
custom_funs = list(estimator = mean,
cp.low = function(x) quantile(x, 0.025),
cp.upp = function(x) quantile(x, 0.975)),
custom_expr_b = "bias = estimator - parameters
mse = (estimator - parameters)^2
cp.quantile = ifelse((parameters >= cp.low) & (parameters <= cp.upp), 1, 0)")
#> $bias.var
#> , , 1
#>
#> [,1]
#> [1,] 1.145122
#>
#>
#> $cp.quantile.var
#> , , 1
#>
#> [,1]
#> [1,] 1
#>
#>
#> $mse.var
#> , , 1
#>
#> [,1]
#> [1,] 6.439299
#>
#>
#> an nlist object with 3 natomic elementsWhen running simulation studies, it is often preferable to save all the results to disk. By default, when the path argument is not specified, results are saved in your working directory.
set.seed(10L)
sims::sims_simulate(code,
parameters = params,
constants = constants,
nsims = 5,
save=TRUE,
exists = NA)
#> [1] TRUE
sma_analyse_bayesian(code = code,
code.add = prior,
mode = sma_set_mode("report"))
#> module dic loaded
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#>
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 10
#> Unobserved stochastic nodes: 1
#> Total graph size: 18
#>
#> Initializing model
#> v data0000001.rds [00:00:00.251]
#> v data0000002.rds [00:00:00.255]
#> v data0000003.rds [00:00:00.245]
#> v data0000004.rds [00:00:00.245]
#> v data0000005.rds [00:00:00.254]
#> Success: 5
#> Failure: 0
#> Remaining: 0
#>
#> Module dic unloaded
sma_derive(code="var=sigma^2", monitor="var")
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> Warning: The following parameters were not in expr and so were dropped from
#> object: 'deviance'.
#> v results0000001.rds [00:00:00.913]
#> v results0000002.rds [00:00:00.968]
#> v results0000003.rds [00:00:01.003]
#> v results0000004.rds [00:00:01.066]
#> v results0000005.rds [00:00:01.071]
#> Success: 5
#> Failure: 0
#> Remaining: 0
#>
sma_evaluate()You may show the files created with
files <- list.files(getwd(), recursive=TRUE, all.files=TRUE)
print(files)
#> [1] ".sims.rds"
#> [2] "analysis0000001/.seeds.rds"
#> [3] "analysis0000001/derived/.parameters.rds"
#> [4] "analysis0000001/derived/deriv0000001.rds"
#> [5] "analysis0000001/derived/deriv0000002.rds"
#> [6] "analysis0000001/derived/deriv0000003.rds"
#> [7] "analysis0000001/derived/deriv0000004.rds"
#> [8] "analysis0000001/derived/deriv0000005.rds"
#> [9] "analysis0000001/results/results0000001.rds"
#> [10] "analysis0000001/results/results0000002.rds"
#> [11] "analysis0000001/results/results0000003.rds"
#> [12] "analysis0000001/results/results0000004.rds"
#> [13] "analysis0000001/results/results0000005.rds"
#> [14] "data0000001.rds"
#> [15] "data0000002.rds"
#> [16] "data0000003.rds"
#> [17] "data0000004.rds"
#> [18] "data0000005.rds"
#> [19] "performance/performance.rds"and read a particular file, e.g.
readRDS(file.path(getwd(), files[19]))
#> $bias.var
#> [1] 1.145122
#>
#> $cp.quantile.var
#> [1] 1
#>
#> $mse.var
#> [1] 6.439299
#>
#> an nlist object with 3 natomic elementsParallelization is achieved using the future package.
To use all available cores on the local machine simply execute the following code before calling any of the package’s functions.
library(future)
plan(multisession)
Please report any issues.
Pull requests are always welcome.
Please note that this project is released with a Contributor Code of Conduct. By contributing, you agree to abide by its terms.