simulation: Simulate Capture Histories

Description

Generate non-spatial or spatial open-population data and fit models.

Usage

sim.nonspatial (N, turnover = list(), p, nsessions, noccasions = 1, intervals = NULL, 
    recapfactor = 1, seed = NULL, savepopn = FALSE, ...)
    
runsim.nonspatial (nrepl = 100, seed = NULL, ncores = NULL, fitargs = list(), 
    extractfn = predict, ...)
runsim.spatial (nrepl = 100, seed = NULL, ncores = NULL, popargs = list(), 
    detargs = list(), fitargs = list(), extractfn = predict, intervals = NULL)
sumsims (sims, parm = 'phi', session = 1, dropifnoSE = TRUE, svtol = NULL, 
    maxcode = 3, true = NULL)
runsim.RMark (nrepl = 100, model = "CJS", model.parameters = NULL, extractfn,
    seed = NULL, ...)

Value

sim.nonspatial --

A capthist object representing an open-population sample

runsim.nonspatial and runsim.spatial --

List with one component (output from extractfn) for each replicate. Each component also has attributes 'eigH' and 'fit' taken from the output of openCR.fit. See Examples to extract convergence codes from 'fit' attribute.

sumsims --

Data.frame with rows `estimate', `SE.estimate', `lcl', `ucl', `RSE', `CI.length' and columns for median, mean, SD and n. If `true' is specified there are additional rows are `Bias' and `RB', and columns for `rRMSE' and `COV'.

Arguments

N: integer population size
turnover: list as described for turnover
p: numeric detection probability
nsessions: number of primary sessions
noccasions: number of secondary sessions per primary session
intervals: intervals between primary sessions (see Details)
recapfactor: numeric multiplier for capture probability after first capture
seed: random number seed see random numbers
savepopn: logical; if TRUE the generated population is saved as an attribute of the capthist object
...: other arguments passed to sim.popn (sim.nonspatial) or sim.nonspatial (runsims)
nrepl: number of replicates
ncores: integer number of cores to be used for parallel processing (see Details)
popargs: list of arguments for sim.popn
detargs: list of arguments for sim.capthist
fitargs: list of arguments for openCR.fit
extractfn: function applied to each fitted openCR model
sims: list output from runsim.nonspatial or runsim.spatial
parm: character name of parameter to summarise
session: integer vector of session numbers to summarise
dropifnoSE: logical; if TRUE then replicates are omitted when SE missing for parm
svtol: numeric; minimum singular value (eigenvalue) considered non-zero
maxcode: integer; maximum accepted value of convergence code
true: true value of requested parm in given session
model: character; RMark model type
model.parameters: list with RMark model specification (see ?mark)

Details

For sim.nonspatial -- If intervals is specified then the number of primary and secondary sessions is inferred from intervals and nsessions and noccasions are ignored. If N and p are vectors of length 2 then subpopulations of the given initial size are sampled with the differing capture probabilities and the resulting capture histories are combined.

runsim.spatial is a relatively simple wrapper for sim.popn, sim.capthist, and openCR.fit. Some arguments are set automatically: the sim.capthist argument 'renumber' is always FALSE; argument 'seed' is ignored within 'popargs' and 'detargs'; if no 'traps' argument is provided in 'detargs' then 'core' from 'popargs' will be used; detargs$popn and fitargs$capthist are derived from the preceding step. The 'type' specified in fitargs may refer to a non-spatial or spatial open-population model ('CJS', 'JSSAsecrfCL' etc.). If the intervals argument is specified it is used to set the intervals attribute of the simulated capthist object; turnover parameters in sim.popn are not scaled by intervals.

Control of parallel processing changed in openCR 1.5.0 to conform to secr. In runsim.nonspatial and runsim.spatial, if ncores is NULL (the default) then the number of cores used for multithreading by openCR.fit is controlled by the environment variable RCPP_PARALLEL_NUM_THREADS. Use the secr function setNumThreads to set RCPP_PARALLEL_NUM_THREADS to a value greater than the default (2, from openCR 1.5 onwards).

Otherwise, (ncores specified in runsim.nonspatial or runsim.spatial) 'ncores' is set to 1 for each replicate and the replicates are split across the specified number of cores.

sumsims assumes output from runsim.nonspatial and runsim.spatial with `extractfn = predict' or `extractfn = summary'. Missing SE usually reflects non-identifiability of a parameter or failure of maximisation, so these replicates are dropped by default. If svtol is specified then the rank of the Hessian is determined by counting eigenvalues that exceed svtol, and replicates are dropped if the rank is less than the number of beta parameters. A value of 1e-5 is suggested for svtol in AIC.openCR, but smaller values may be appropriate for larger models (MARK has its own algorithm for this threshold).

Replicates are also dropped by sumsims if the convergence code exceeds 'maxcode'. The maximisation functions nlm (used for method = 'Newton-Raphson', the default), and optim (all other methods) return different convergence codes; their help pages should be consulted. The default is to accept code = 3 from nlm, although the status of such maximisations is ambiguous.

Examples

Run this code


if (FALSE) {

cores <- 2   # for CRAN check; increase as available

ch <- sim.nonspatial(100, list(phi = 0.7, lambda = 1.1), p = 0.3, nsessions = 8, noccasions=2)
openCR.fit(ch, type = 'CJS')

turnover <- list(phi = 0.85, lambda = 1.0, recrmodel = 'constantN')
set.seed(123)

## using type = 'JSSAlCL' and extractfn = predict

fitarg <- list(type = 'JSSAlCL', model = list(p~t, phi~t, lambda~t))
out <- runsim.nonspatial(nrepl = 100, N = 100, ncores = cores, turnover = turnover, 
   p = 0.2, recapfactor = 4, nsessions = 10, noccasions = 1, fitargs = fitarg)
sumsims(out, 'lambda', 1:10)

## using type = 'Pradelg' and extractfn = derived
## homogeneous p
fitarg <- list(type = 'Pradelg', model = list(p~t, phi~t, gamma~t))
outg <- runsim.nonspatial(nrepl = 100, N = 100, ncores = cores, turnover = turnover, 
    p = 0.2, recapfactor = 4, nsessions = 10, noccasions = 1, 
    fitargs = fitarg, extractfn = derived)
apply(sapply(outg, function(x) x$estimates$lambda),1,mean)

turnover <- list(phi = 0.85, lambda = 1.0, recrmodel = 'discrete')

## 2-class mixture for p
outg2 <- runsim.nonspatial(nrepl = 100, N = c(50,50), ncores = cores, turnover = turnover, 
    p = c(0.3,0.9), recapfactor = 1, nsessions = 10, noccasions = 1, 
    fitargs = fitarg, extractfn = derived)
outg3 <- runsim.nonspatial(nrepl = 100, N = c(50,50), ncores = cores, turnover = turnover, 
    p = c(0.3,0.3), recapfactor = 1, nsessions = 10, noccasions = 1, 
    fitargs = fitarg, extractfn = derived)
apply(sapply(outg2, function(x) x$estimates$lambda),1,mean)

plot(2:10, apply(sapply(outg2, function(x) x$estimates$lambda),1,mean)[-1],
    type='o', xlim = c(1,10), ylim = c(0.9,1.1))

## RMark

extfn <- function(x) x$results$real$estimate[3:11]
MarkPath <- 'c:/mark'  # customise as needed
turnover <- list(phi = 0.85, lambda = 1.0, recrmodel = 'discrete')
outrm <- runsim.RMark (nrepl = 100, model = 'Pradlambda', extractfn = extfn, 
                       model.parameters = list(Lambda=list(formula=~time)),
                       N = c(200,200), turnover = turnover, p = c(0.3,0.9),
                       recapfactor = 1, nsessions = 10, noccasions = 1)
extout <- apply(do.call(rbind, outrm),1,mean)

## Spatial

grid <- make.grid()
msk <- make.mask(grid, type = 'trapbuffer', nx = 32)
turnover <- list(phi = 0.8, lambda = 1)
poparg <- list(D = 10, core = grid, buffer = 100, Ndist = 'fixed', nsessions = 6, 
    details = turnover)
detarg <- list(noccasions = 5, detectfn = 'HHN', detectpar = list(lambda0 = 0.5, sigma = 20))
fitarg <- list(type = 'JSSAsecrfCL', mask = msk, model = list(phi~1, f~1))
sims <- runsim.spatial (nrepl = 7, ncores = cores, pop = poparg, det = detarg, fit = fitarg)
sumsims(sims)

## extract the convergence code from nlm for each replicate in preceding simulation
sapply(lapply(sims, attr, 'fit'), '[[', 'code')
## if method != 'Newton-Raphson then optim is used and the code is named 'convergence'
# sapply(lapply(sims, attr, 'fit'), '[[', 'convergence')

}