Simulation workflow in bSims

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Peter Solymos

We recommend exploring the simulation settings interactively in the shiny apps using run_app("bsimsH") app for the homogeneous habitat case and the run_app("bsimsHER") app for the stratified habitat case. The apps represent the simulation layers as tabs, the last tab presenting the settings that can be copied onto the clipboard and pasted into the R session or code. In simple situations, comparing results from a few different settings might be enough.

Let us consider the following simple comparison: we want to see how much of an effect does roads have when the only effect is that the road stratum is unsuitable. Otherwise there are no behavioral or detectability effects of the road.

The bsims_all function accepts all the arguments we discussed before for the simulation layers. Unspecified arguments will be taken to be the default value. However, bsims_all does not evaluate these arguments, but it creates a closure with the settings. Realizations can be drawn as:

The replicate function takes an argument for the number of replicates (B) and returns a list of transcript objects with B elements. The cl argument can be used to parallelize the work, it can be a numeric value on Unix/Linux/OSX, or a cluster object on any OS. The recover = TRUE argument allows to run simulations with error catching.

Simulated objects returned by bsims_all will contain different realizations and all the conditionally independent layers. Use a customized layered approach if former layers are meant to be kept identical across runs.

In more complex situations the shiny apps will help identifying corner cases that are used to define a gradient of settings for single or multiple simulation options. You can copy the bsims_all settings from the app to be used in simulations.

Sensitivity analysis

In a sensitivity analysis we evaluate how varying one or more settings affect the estimates. This requires setting up a series of values for a setting (argument) while keeping others constant.

Let us consider the following scenario: we would like to evaluate how the estimates are changing with increasing road width. We will use the expand_list function which creates a list from all combinations of the supplied inputs. Note that we need to wrap vectors inside list() to avoid interpreting those as values to iterate over.

s <- expand_list(
  road = c(0, 0.5, 1),
  density = list(c(1, 1, 0)),
  tint = list(tint),
  rint = list(rint))
str(s)
#> List of 3
#>  $ :List of 4
#>   ..$ road   : num 0
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf
#>  $ :List of 4
#>   ..$ road   : num 0.5
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf
#>  $ :List of 4
#>   ..$ road   : num 1
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf

We now can use this list of settings to run simulations for each. The following illustrates the use of multiple cores:

b <- lapply(s, bsims_all)
nc <- 4 # number of cores to use
library(parallel)
cl <- makeCluster(nc)
bb <- lapply(b, function(z) z$replicate(B, cl=cl))
stopCluster(cl)

In some cases, we want to evaluate crossed effects of multiple settings. This will give us information about how these settings interact. For example, road width and spatial pattern (random vs. clustered):

s <- expand_list(
  road = c(0, 0.5),
  xy_fun = list(
    NULL,
    function(d) exp(-d^2/1^2) + 0.5*(1-exp(-d^2/4^2))),
  density = list(c(1, 1, 0)),
  tint = list(tint),
  rint = list(rint))
str(s)
#> List of 4
#>  $ :List of 5
#>   ..$ road   : num 0
#>   ..$ xy_fun : NULL
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf
#>  $ :List of 5
#>   ..$ road   : num 0.5
#>   ..$ xy_fun : NULL
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf
#>  $ :List of 5
#>   ..$ road   : num 0
#>   ..$ xy_fun :function (d)  
#>   .. ..- attr(*, "srcref")= 'srcref' int [1:8] 5 5 5 53 5 53 5 5
#>   .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x1210c6a68> 
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf
#>  $ :List of 5
#>   ..$ road   : num 0.5
#>   ..$ xy_fun :function (d)  
#>   .. ..- attr(*, "srcref")= 'srcref' int [1:8] 5 5 5 53 5 53 5 5
#>   .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x1210c6a68> 
#>   ..$ density: num [1:3] 1 1 0
#>   ..$ tint   : num [1:5] 2 4 6 8 10
#>   ..$ rint   : num [1:5] 0.5 1 1.5 2 Inf

Varying landscapes

Studying covariate effects on density, cue rates, and detection distances sometimes require that we simulate a series of landscapes that differ.

We exploit the fact that arguments to bsims_all can be supplied as a list, which is the same as a single-row data frame. This should work for all arguments that accept atomic vectors as arguments:

bsims_all(
  road = 0.5,
  density = 1)
#> bSims wrapper object with settings:
#>   road   : 0.5
#>   density: 1


bsims_all(
  list(
    road = 0.5,
    density = 1))
#> bSims wrapper object with settings:
#>   road   : 0.5
#>   density: 1


bsims_all(
  data.frame(
    road = 0.5,
    density = 1))
#> bSims wrapper object with settings:
#>   road   : 0.5
#>   density: 1

# number of stations to visit
n <- 5

# random predictors: continuous and discrete
x <- data.frame(x1=runif(n,-1,2), x2=rnorm(n))

# density
D <- drop(exp(model.matrix(~x2, x) %*% c(0,-0.5)))
summary(D)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.9735  1.1403  1.2043  2.2396  2.3035  5.5766


# cue rate
phi <- drop(exp(model.matrix(~x1+I(x1^2), x) %*% c(-1,-0.25,-1)))
summary(phi)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> 0.004842 0.005648 0.200322 0.172966 0.301399 0.352621


# this data frame collects the columns to be used as arguments
s <- data.frame(
    D=D,
    vocal_rate = phi, 
    duration = 10,
    condition = "det1",
    tau = 1)

# each row from s becomes a simulation settings object
bb <- lapply(1:n, function(i) bsims_all(s[i,]))

# define how you want the data extracted
get_counts <- function(b) {
    o <- b$new() # simulate
    get_table(o)[1,1]
}

x$y <- sapply(bb, get_counts)
x
#>           x1         x2 y
#> 1  0.1158007 -3.4371555 4
#> 2  1.9597121 -0.2626537 0
#> 3  1.9224377 -0.3717355 0
#> 4 -0.5886233  0.0536659 6
#> 5  0.6645918 -1.6688566 4