/
bsims03-workflow.Rmd
205 lines (170 loc) · 5.28 KB
/
bsims03-workflow.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
---
title: "Simulation workflow in bSims"
author: "Peter Solymos"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Simulation workflow in bSims}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup,include=FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
par(mar = c(1, 1, 1, 1))
set.seed(429)
suppressPackageStartupMessages(library(bSims))
```
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.
```{r sim1}
library(bSims)
tint <- c(2, 4, 6, 8, 10)
rint <- c(0.5, 1, 1.5, 2, Inf) # unlimited
## no road
b1 <- bsims_all(
road = 0,
density = c(1, 1, 0),
tint = tint,
rint = rint)
## road
b2 <- bsims_all(
road = 0.5,
density = c(1, 1, 0),
tint = tint,
rint = rint)
b1
b2
```
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:
```{r sim2}
b1$new()
b2$new()
```
Run multiple realizations is done as:
```{r sim3,eval=FALSE}
B <- 25 # number of runs
bb1 <- b1$replicate(B)
bb2 <- b2$replicate(B)
```
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.
```{r grid1}
s <- expand_list(
road = c(0, 0.5, 1),
density = list(c(1, 1, 0)),
tint = list(tint),
rint = list(rint))
str(s)
```
We now can use this list of settings to run simulations for each.
The following illustrates the use of multiple cores:
```{r grid2,eval=FALSE}
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):
```{r grid3}
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)
```
## 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:
```{r equal}
bsims_all(
road = 0.5,
density = 1)
bsims_all(
list(
road = 0.5,
density = 1))
bsims_all(
data.frame(
road = 0.5,
density = 1))
```
```{r variables}
# 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)
# cue rate
phi <- drop(exp(model.matrix(~x1+I(x1^2), x) %*% c(-1,-0.25,-1)))
summary(phi)
# 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
```