run.scenarios: 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, ...)

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).

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 sim.popn (optional)
CH.function: character name of function to simulate capthist
det.args: list of named arguments to sim.capthist (optional)
fit: logical or character; if TRUE a model is fitted with fit.function, otherwise data are generated but no model is fitted
(see also Multi-model fit and Design-only statistics in Details)
fit.function: character name of function to use for model fitting
fit.args: list of named arguments to secr.fit or ipsecr.fit (optional)
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 join is applied to multisession capthist
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 run.scenarios
scen: integer vector of scenario subscripts
repl: integer vector of subscripts in range 1:nrepl

Author

Murray Efford

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 which should be of the same length as `trapset' 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.

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.

Examples

Run this code


## 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)

if (FALSE) {

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'))
}

Run the code above in your browser using DataLab