| run.scenarios {secrdesign} | R Documentation |
Simulate Sampling Designs
Description
This function performs simulations to predict the precision of density and other
estimates from simple 1-session SECR designs. Scenarios are specified
via an input dataframe that will usually be constructed with
make.scenarios. Each scenario comprises an index to a detector layout,
the number of sampling occasions, and specified density (D) and detection
parameters (usually g_0 and \sigma).
Detector layouts are provided in a separate list trapset. This
may comprise an actual field design input with read.traps or
‘traps’ objects constructed with make.grid etc., as in the
Examples. Even a single layout must be presented as a component of a
list (e.g., list(make.grid())).
Alternative approaches are offered for predicting precision. Both start
by generating a pseudorandom dataset under the design using the
parameter values for a particular scenario. The first estimates the
parameter values and their standard errors from each dataset by
maximizing the full likelihood, as usual in secr.fit. The second
takes the short cut of computing variances and SE from the Hessian
estimated numerically at the known expected values of the parameters,
without maximizing the likelihood. Set method = "none" in fit.args
for this shortcut.
Usage
run.scenarios(nrepl, scenarios, trapset, maskset, xsigma = 4, nx = 32,
pop.args, CH.function = c("sim.capthist", "simCH"), det.args,
fit = FALSE, fit.function = c("secr.fit", "ipsecr.fit"),
fit.args, chatnsim, extractfn = NULL, multisession = FALSE,
joinsessions = FALSE, ncores = NULL, byscenario = FALSE, seed = 123,
trap.args, prefix = NULL, ...)
fit.models(rawdata, fit = FALSE, fit.function = c("secr.fit", "ipsecr.fit"),
fit.args, chatnsim, extractfn = NULL, ncores = NULL, byscenario = FALSE,
scen, repl, ...)
Arguments
nrepl |
integer number of replicate simulations |
scenarios |
dataframe of simulation scenarios |
trapset |
secr traps object or a list of traps objects or functions |
maskset |
secr mask object or a list of mask objects (optional) |
xsigma |
numeric buffer width as multiple of sigma (alternative to maskset) |
nx |
integer number of cells in mask in x direction (alternative to maskset) |
pop.args |
list of named arguments to
|
CH.function |
character name of function to simulate capthist |
det.args |
list of named arguments to
|
fit |
logical or character; if TRUE a model is fitted with |
fit.function |
character name of function to use for model fitting |
fit.args |
list of named arguments to |
chatnsim |
integer number of simulations for overdispersion of mark-resight models |
extractfn |
function to extract a vector of statistics from secr model |
multisession |
logical; if TRUE groups are treated as additional sessions |
joinsessions |
logical; if TRUE function |
ncores |
integer number of cores for parallel processing or NULL |
byscenario |
logical; if TRUE then each scenario is sent to a different core |
seed |
integer pseudorandom number seed |
trap.args |
list of arguments for trapset components if using function option |
prefix |
character to name files saving output of each scenario |
... |
other arguments passed to extractfn |
rawdata |
‘rawdata’ object from previous call to |
scen |
integer vector of scenario subscripts |
repl |
integer vector of subscripts in range 1:nrepl |
Details
Designs are constructed from the trap layouts in trapset, the
numbers of grids in ngrid, and the numbers of sampling
occasions (secondary sessions) in noccasions. These are
not crossed: the number of designs is the maximum length of any
of these arguments. Any of these arguments whose length is less than
the maximum will be replicated to match.
pop.args is used to customize the simulated population
distribution. It will usually comprise a single list, but may be a
list of lists (one per popindex value in scenarios).
det.args may be used to customize some aspects of the detection
modelling in sim.capthist, but not traps, popn, detectpar,
detectfn, and noccasions, which are controlled directly by the
scenarios. It will usually comprise a single list, but may be a list
of lists (one per detindex value in scenarios).
fit.args is used to customize the fitted model; it will usually
comprise a single list. If you are interested in precision alone, use
fit.args=list(method = 'none') to obtain variance estimates
from the hessian evaluated at the parameter estimates. This is much
faster than a complete model fit, and usually accurate enough.
If no extractfn is supplied then a default is used - see
Examples. Replacement functions should follow this pattern i.e. test
for whether the single argument is an secr object, and if not supply a
named vector of NA values of the correct length.
Using extractfn = summary has the advantage of allowing both model fits and raw statistics to be extracted from one set of simulations. However, this approach requires an additional step to retrieve the desired numeric results from each replicate (see count.summary and predict.summary).
Parallel processing
If byscenario = TRUE then by default each scenario will be run in a separate worker
process using parLapply from parallel (see also Parallel). The number of scenarios should not exceed the available number of cores (set by the 'ncores' argument or a prior call to 'setNumThreads').
If byscenario = FALSE then from secrdesign 2.6.0 onwards the usual multithreading of secr 4.5 is applied. The number of cores should usually be preset with 'setNumThreads'. If ncores is provided then the environment variable RCPP_PARALLEL_NUM_THREADS is reset. The default behaviour of the fitting functions (secr.fit, ipsecr.fit) is to use this value (unless specified in fit.args).
When ‘byscenario = TRUE' the L’Ecuyer pseudorandom generator is used with a separate random
number stream for each core (see clusterSetRNGStream).
For ncores > 1 it pays to keep an eye on the processes from the
Performance page of Windows Task Manager (<ctrl><alt><del>), or ‘top’ in
linux OS. If you interrupt run.scenarios (<Esc> from Windows)
you may occasionally find some processes do not terminate and have to
be manually terminated from the Task Manager - they appear as
Rscript.exe on the Processes page.
Alternate functions for simulation and fitting
The default is to use functions sim.capthist and secr.fit from secr. Either may be substituted by the corresponding function (simCH or ipsecr.fit) from package ipsecr if that has been installed.
Multi-model fit
Multiple models may be fitted to the same simulated data for multi-model inference. This requires both (i) ‘fit = "multifit"’, and (ii) 'fit.args' should be a nested list (fit arguments within models within fit.index) with a separate specification for each model fit. See the vignette for examples.
Design-only statistics
Designs for distance sampling were evaluated by Fewster and Buckland (2004) by computing statistics from simulated detections without fitting a model to estimate the detection parameters. An analogous procedure for SECR is implemented by setting fit = 'design'. A new default extractfn (designextractfn) computes the effective sampling area with the secr function pdot and returns a vector of values -
| n | number of individuals detected | |
| r | number of recaptures | |
| esa | effective sampling area, given the known detection parameters | |
| D | D = n/esa | |
The resulting simulation object is of type 'selectedstatistics' for which the summary method works as usual.
A similar effect may be achieved by providing a custom extractfn and passing arguments to it via the dots argument of run.scenarios.
Miscellaneous
From 2.2.0, two or more rows in scenarios may share the same scenario number. This is used to generate multiple population subclasses (e.g. sexes) differing in density and/or detection parameters. If multisession = TRUE the subclasses become separate sessions in a multi-session capthist object (this may require a custom extractfn). multisession is ignored with a warning if each scenario row has a unique number.
From 2.7.0, each component of ‘trapset’ may be a function that constructs a detector layout. This allows layouts to be constructed dynamically at the time each capthist is generated; arguments of each function are provided in the ‘trap.args’ list. The primary purpose is to allow systematic grids, laceworks etc. to be constructed with a unique random origin for each replicate. The ‘maskset’ argument must be provided - it should cover all potential layouts, regardless of origins.
In fit.models the arguments scen and repl may be used to select a subset of datasets for model fitting.
Mark-resight
chatnsim controls an additional quasi-likelihood model step to adjust for overdispersion of sighting counts. No adjustment happens when chatnsim = 0; otherwise abs(chatnsim) gives the number of simulations to perform to estimate overdispersion. If chatnsim < 0 then the quasilikelihood is used only to re-estimate the variance at the previous MLE (method = "none").
Intermediate output
If 'prefix' is provided than results will be saved for each scenario separately. The filename of scenario 1 is of the form 'prefix1.RDS'. The prefix may include a file path.
Further processing
A summary method is provided (see
summary.secrdesign). It is usually necessary to process
the simulation results further with predict.fittedmodels
and/or select.stats before summarization.
Value
An object of class (x, ‘secrdesign’, ‘list’), where x is one of ‘fittedmodels’, ‘estimatetables’, ‘selectedstatistics’ or ‘rawdata’, with components
call |
function call |
version |
character string including the software version number |
starttime |
character string for date and time of run |
proctime |
processor time for simulations, in seconds |
scenarios |
dataframe as input |
trapset |
list of trap layouts as input |
maskset |
list of habitat masks (input or generated) |
xsigma |
from input |
nx |
from input |
pop.args |
from input |
CH.function |
from input |
det.args |
from input |
fit |
from input |
fit.function |
from input |
fit.args |
from input |
extractfn |
function used to extract statistics from each simulation |
seed |
from input |
nrepl |
from input |
output |
list with one component per scenario |
outputtype |
character code - see vignette |
If fit = FALSE and extractfn = identity the result is of
class (‘rawdata’, ‘secrdesign’, ‘list’). This may be used as input to
fit.models, which interprets each model specification in
fit.args as a new ‘sub-scenario’ of each input scenario
(i.e. all models are fitted to every dataset). The output
possibilities are the same as for run.scenarios.
If subclasses have been defined (i.e. scenarios has multiple rows with the same scenario ID), each simulated capthist object has covariates with a character-valued column named "group" ("1", "2" etc.) (there is also a column "sex" generated automatically by sim.popn).
Note
100 ha = 1 km^2.
fit.function = 'openCR.fit' was deprecated from 2.5.8 and has been removed.
Author(s)
Murray Efford
References
Fewster, R. M. and Buckland, S. T. 2004. Assessment of distance sampling estimators. In: S. T. Buckland, D. R. Anderson, K. P. Burnham, J. L. Laake, D. L. Borchers and L. Thomas (eds) Advanced distance sampling. Oxford University Press, Oxford, U. K. Pp. 281–306.
See Also
Miscellaneous –
secr functions used internally –
To combine output –
Examples
## Simple example: generate and summarise trapping data
## at two densities and for two levels of sampling frequency
scen1 <- make.scenarios(D = c(5,10), sigma = 25, g0 = 0.2, noccasions =
c(5,10))
traps1 <- make.grid() ## default 6 x 6 trap grid
tmp1 <- run.scenarios(nrepl = 20, trapset = traps1, scenarios = scen1,
fit = FALSE)
summary(tmp1)
## Not run:
setNumThreads(7)
##########################################
# new summary method (secrdesign >= 2.8.1)
# assumes fit = TRUE, extractfn = predict
tmp2 <- run.scenarios(nrepl = 10, trapset = traps1, scenarios = scen1,
fit = TRUE, extractfn = predict)
estimateSummary(tmp2, format = "data.frame",
cols = c('scenario', 'noccasions'))
###########################
## 2-phase example
## first make and save rawdata
scen1 <- make.scenarios(D = c(5,10), sigma = 25, g0 = 0.2)
traps1 <- make.grid() ## default 6 x 6 trap grid
tmp1 <- run.scenarios(nrepl = 20, trapset = traps1, scenarios = scen1,
fit = FALSE, extractfn = identity)
## review rawdata
summary(tmp1)
## then fit and summarise models
tmp2 <- fit.models(tmp1, fit.args = list(list(model = g0~1),
list(model = g0~T)), fit = TRUE)
summary(tmp2)
###########################
## Construct a list of detector arrays
## Each is a set of 5 parallel lines with variable between-line spacing;
## the argument that we want to vary (spacey) follows nx, ny and spacex
## in the argument list of make.grid().
spacey <- seq(2000,5000,500)
names(spacey) <- paste('line', spacey, sep = '.')
trapset <- lapply(spacey, make.grid, nx = 101, ny = 5, spacex = 1000,
detector = 'proximity')
## Make corresponding set of masks with constant spacing (1 km)
maskset <- lapply(trapset, make.mask, buffer = 8000, spacing = 1000,
type = 'trapbuffer')
## Generate scenarios
scen <- make.scenarios (trapsindex = 1:length(spacey), nrepeats = 8,
noccasions = 2, D = 0.0002, g0 = c(0.05, 0.1), sigma = 1600, cross = TRUE)
## RSE without fitting model
sim <- run.scenarios (50, scenarios = scen, trapset = trapset, maskset = maskset,
fit = TRUE, fit.args = list(method = 'none'), seed = 123)
## Extract statistics for predicted density
sim <- select.stats(sim, parameter = 'D')
## Plot to compare line spacing
summ <- summary (sim, type='array', fields = c('mean','lcl','ucl'))$OUTPUT
plot(0,0,type='n', xlim=c(1.500,5.500), ylim = c(0,0.36), yaxs = 'i',
xaxs = 'i', xlab = 'Line spacing km', ylab = 'RSE (D)')
xv <- seq(2,5,0.5)
points(xv, summ$mean[,1,'RSE'], type='b', pch=1)
points(xv, summ$mean[,2,'RSE'], type='b', pch=16)
segments(xv, summ$lcl[,1,'RSE'], xv, summ$ucl[,1,'RSE'])
segments(xv, summ$lcl[,2,'RSE'], xv, summ$ucl[,2,'RSE'])
legend(4,0.345, pch=c(1,16), title = 'Baseline detection',
legend = c('g0 = 0.05', 'g0 = 0.1'))
## End(Not run)