secrdesign (version 2.5.5)

run.scenarios: Simulate Sampling Designs

Description

This function performs simulations to predict the precision of abundance 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())).

If ncores > 1 then by default each scenario will be run in a separate worker process using parLapply from parallel (see also Parallel).

If byscenario = FALSE then replicates are split among cores (the default is to split scenarios among cores), which is useful if you have more cores than scenarios. Dividing replicates among cores (byscenario = FALSE) also largely avoids the inefficiency that results when some workers finish much sooner than others (load balancing is not an option in run.scenarios). Setting ncores greater than the number of scenarios causes an error when byscenario = TRUE.

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 of openCR.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" for this shortcut.

Usage

run.scenarios(nrepl, scenarios, trapset, maskset, xsigma = 4, nx = 32,
    pop.args, det.args, fit = FALSE, fit.function = c("secr.fit", "openCR.fit"), 
    fit.args, chatnsim, extractfn = NULL, multisession = FALSE,
    ncores = 1, byscenario = TRUE, seed = 123, ...)
fit.models(rawdata, fit = FALSE, fit.function = c("secr.fit", "openCR.fit"), 
    fit.args, chatnsim, extractfn = NULL, ncores = 1, byscenario = TRUE,
    scen, repl, ...)

Arguments

nrepl

integer number of replicate simulations

scenarios

dataframe of simulation scenarios

trapset

secr traps object or a list of traps objects

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)

det.args

list of named arguments to sim.capthist (optional)

fit

logical; if TRUE a model is fitted with secr.fit, otherwise data are generated but no model is fitted

fit.function

character name of function to use for model fitting

fit.args

list of named arguments to secr.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

byscenario

logical; if TRUE and ncores>1 then scenarios are sent to different cores

seed

integer pseudorandom number seed

…

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

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

det.args

from input

fit

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

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.

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.

The L'Ecuyer pseudorandom generator is used with a separate random number stream for each core (see clusterSetRNGStream).

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.

In fit.models the arguments scen and repl may be used to select a subset of datasets for model fitting.

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

See Also

predict.fittedmodels, scenarioSummary, select.stats, summary.secrdesign, summary.selectedstatistics, sim.popn, sim.capthist, secr.fit

Examples

# NOT RUN {
## 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 {
###########################
## 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, ncores = 4)
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,
    ncores = 8, 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'))$summary
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'))
# }
# NOT RUN {
# }