run.scenarios: Simulate Sampling Designs

In secrdesign: Sampling Design for Spatially Explicit Capture-Recapture

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,
    ncores = NULL, byscenario = FALSE, seed = 123, trap.args, ...)

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 sim.popn (optional)
CH.function	character name of function to simulate capthist
det.args	list of named arguments to sim.capthist (optional)
fit	logical; if TRUE a model is fitted with fit.function, otherwise data are generated but no model is fitted (see also 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
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
...	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

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.

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

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

select.stats,

summary.secrdesign,

summary.estimatetables,

summary.selectedstatistics,

estimateSummary

Miscellaneous –

predict.fittedmodels,

scenarioSummary,

count.summary,

predict.summary

secr functions used internally –

sim.popn,

sim.capthist,

secr.fit

To combine output –

rbind.estimatetables,

rbind.selectedstatistics,

c.estimatetables,

c.selectedstatistics

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)

secrdesign documentation built on March 31, 2023, 10:25 p.m.