\renewcommand{\vec}[1]{\mathbf{#1}} \newcommand{\bm}[1]{\mbox{\boldmath $#1$}}
mycache < FALSE Sys.setenv(RCPP_PARALLEL_BACKEND = "tinythread") ## to dodge CRAN ASAN issue
\vspace{12pt}
The R package openCR fits both nonspatial and spatial capturerecapture models to data from open animal populations, where there is turnover during sampling. The interface generally resembles that of secr (Efford 2022a) upon which openCR depends for some functions. This document explains the purpose and general features of openCR. Help pages should be consulted for more detail on particular functions. The vignette openCRkernel.pdf explains the intricacies of movement kernels (see also Efford and Schofield 2022). Worked examples using published datasets are given in another vignette openCRexamples.pdf[^footnote1]. The spatial model was described by Efford and Schofield (2020).
This is still something of a work in progress, so be careful to check results 'make sense' and be aware of limitations.
[^footnote1]: These supplementary vignettes are not included with the package. It is intended to distribute them on the website https://www.otago.ac.nz/density/. Otherwise contact the author.
openCR fits nonspatial openpopulation models of the CormackJollySeber (CJS) and JollySeberSchwarzArnason (JSSA[^footnote2] or 'POPAN') types. JSSA models are offered in both full and conditional likelihood forms, each with several parameterizations of recruitment, and incorporating Pollock's robust design. Conditionallikelihood JSSA models are also called PradelLinkBarker (PLB) models. Pradel analyses are also provided.
[^footnote2]: As far as I know, this abbreviation was first used by Pledger et al. (2010). Recognising the contributions of Crosbie and Manly, Schofield and Barker (2009) and Cowen et al. (2010) referred to it the CrosbieManlyArnasonSchwarz (CMAS) model. Link and Barker (2010) used 'CrosbieManlySchwarzArnason' (CMSA) for the same model. CMSA has since been used by various authors, including Schofield and Barker (2016). JSSA is used in openCR because this highlights its evolution from the widely known JollySeber model. POPAN refers to the software of Schwarz and Arnason (1996), recycled as the name of a data type in MARK.
Spatial versions of the CJS and JSSA model types are also provided[^footnote2a]. The spatial models allow for 'multi', 'proximity' or 'count' detectors as defined in secr. Several functions are implemented for the decline in hazard of detection with distance. Movement between primary sessions may be modelled (cf Ergon and Gardner 2014; Glennie et al. 2019), but particular care is needed, especially with respect to kernel truncation.
[^footnote2a]: A direct spatial implementation of CJS fails because the distribution of detected animals is not uniform at first detection, but rather biased towards the vicinity of the detectors. Set details = list(CJSp1 = TRUE)
to model first detections and get sensible estimates (model type 'CJSsecr').
Data are assumed to be from a robust design. Secondary sampling sessions are nested within primary sessions and all turnover (births, deaths, immigration or emigration) is between primary sessions (Pollock 1982). There may be a single secondary session per primary session (this limits identifiability of some parameters).
Models are specified using formula notation as in secr. Possible predictors include both predefined variables for learned responses, trend over time, etc., and userprovided covariates. Models are fitted by numerically maximizing the log likelihood. The likelihood is formed as a product over capture histories (Pledger at al. 2010) rather than from summary statistics. The fitted model is an object of class 'openCR' for which generic methods are implemented (print
, predict
, AIC
, plot
etc.).
Variation in a parameter between primary sessions is modelled as e.g., model = phi ~ session
[^footnote3]. Withinsession variation in detection parameters may also be modelled (see field vole example in openCRexamples.pdf).
[^footnote3]: This is equivalent of ~t in Lebreton et al. (1992) or ~time in RMark, and openCR recognises ~ t
as a synonym of ~ session
.
A selection of parameterizations is offered for recruitment in JSSA models. Models can also be parameterized in terms of the timespecific population size (nonspatial models) or density (spatial models), avoiding the superpopulation parameter.
Superpopulation size (or density in the case of secr models) may be computed as a derived parameter from 'CL' models with the function derived(), which also computes timespecific population sizes and densities.
openCR has definite limitations that may or may not be addressed in future versions. Important differences between secr and openCR are noted here. Online help is not guaranteed: users should attempt to solve their own problems, or seek help from other users via phidot or secrgroup.
We start with a simple nonspatial example. Lebreton et al. (1992) demonstrated CormackJollySeber methods with a dataset on European Dipper (Cinclus cinclus) collected by Marzolin (1988). The object dipperCH
distributed with openCR provides these data in the secr 'capthist' format. See the Examples section of its help page ?dipperCH
for code to input the data from other sources.
library(openCR) # also loads secr options(digits = 4, width = 90) # for more readable output
Dippers were captured annually over 19811987.
m.array(dipperCH, never.recap = T) # compare Lebreton et al. 1992 Table 10
We can fit a CormackJollySeber model directly with openCR.fit
and display the estimates:
dipper.phi.t < openCR.fit(dipperCH, type = 'CJS', model = phi~t) predict(dipper.phi.t) plot(dipper.phi.t, par = 'phi', ylim = c(0,1), pch = 16, col = 'red')
From this example you can see some of the virtues of openCR
See openCRexamples.pdf for more extensive analyses of this dataset.
There is a large literature on openpopulation capturerecapture modelling. Almost all modern models derive from the CormackJollySeber (CJS) or JollySeber (JS) models (Seber 1982), with refinements by Crosbie and Manly (1985), Schwarz and Arnason (1996), Pradel (1996) and others. The MARK software (White and Burnham 1999) implemented many of these developments and remains the standard. This section describes differences among models as they relate to openCR.
The split between the CJS and JS model lineages is fundamental. CJS models do not model the first capture of each animal; they condition on that capture and model subsequent recapture probabilities $p$ and apparent survival $\phi$. CJS estimates of apparent survival are robust and useful (Lebreton et al. 1992), but CJS models stop short of estimating abundance, recruitment or population trend.
JS models model the first capture of each animal, and lead either directly or indirectly to estimates of abundance and recruitment. The modern development of JS methods rests heavily on Schwarz and Arnason (1996), so openCR follows Pledger et al. (2010) in using the label 'JSSA'. JSSA models were the basis of the POPAN software, which led to the POPAN data type in MARK. JSSA models are the main focus of openCR.
The JSSA model appears in several different forms whose unity is obscured by differing parameterizations of recruitment. The classic POPAN formulation uses entry probabilities: the members of a notional superpopulation enter the population with timespecific probability $\beta_j$ (PENT in MARK), an idea from Crosbie and Manly (1985). Other parameterizations are
Estimates of recruitment or implied recruitment from any one of these six parameterizations can be used to infer the others[^footnote4]. The choice of parameterization rests on which is more natural for the problem in hand (and allows the desired constraints to be applied) and on practicalities (some are more likely to give numerical problems than others).
[^footnote4]: except for some mostly trivial differences relating to removals
Schwarz (2001) is illuminating (see also chapter on JollySeber models by Schwarz and Arnason in the MARK book, Cooch and White 2019). Pradel (1996), Williams, Nichols and Conroy (2002: p.518 et seq.), Pledger et al. (2003, 2010) and Link and Barker (2005) also comment on and compare JS parameterizations. See also the MARK help page on 'Recruitment Parameters in JollySeber models' ('Recruitment Parameters' in the help index).
For each JSSA recruitment parameterization there is a choice between models that include the total number of detected individuals ($u_\cdot$ or $n$ in different notations), and models that condition on this number. Conditionallikelihood models do not directly estimate abundance; abundance is estimated as a derived parameter (Schwarz and Arnason 1996). Fulllikelihood models include abundance as a parameter. The choice of formulation has virtually no effect on the parameter estimates[^footnote5]. The conditional likelihood form is somewhat faster and easier to fit (Schwarz and Arnason 1996), and it focuses on parameters that are estimated robustly (apparent survival, seniority, population growth rate).
The conditional models discussed by Pradel (1996), Link and Barker (2005), Schofield and Barker (2016) and others lack a distinguishing label to indicate their collective similarity. The label PradelLinkBarker PLB was suggested by Efford and Schofield (2020).
[^footnote5]: this may not be true for spatial models with spatially varying density, but these models are not considered in openCR.
Historically the CJS and JS likelihoods have been expressed in terms of 'sufficient statistics' that are timespecific counts of animals in different categories, such as the number caught, the number marked etc. This approach is used in the openCR function JS.direct
and with the Pradel model type in openCR.fit
. The likelihood may also be computed as a product over terms, one for each observed capture history[^footnote6]. Modelling of individual capture histories, is slower, but it is extremely flexible, allowing direct inclusion of censoring, learned responses, individual covariates, secondary sessions and other extensions. This is the approach used in MARK and openCR.fit
.
[^footnote6]: strictly, the product over observed histories is only one component of the likelihood
Most published formulations of CJS and JSSA models admit only one secondary session per primary session. Data collected according to a robust design with multiple secondary sessions must be collapsed to a single sample per primary session. However, it is simple to adapt the capturehistory models for multiple secondary occasions, and this makes better use of the data. MARK offers many specific robust design models. A robust design is assumed in openCR; data with a single secondary session per primary session are merely a special case.
Models may be spatially explicit or not. Nonspatial models ignore the spatial distribution of animals. Spatial models use the spatially explicit capturerecapture paradigm of Efford (2004), Borchers and Efford (2008) and Royle et al. (2014). Open population spatial models using MCMC were published by Gardner et al. (2010), Chandler and Clark (2014), Ergon and Gardner (2014), Whittington and Sawaya (2015) and others. Glennie et al. (2019) proposed a frequentist hidden Markov formulation. The spatial models in openCR are described by Efford and Schofield (2020) and provide very similar estimates to those of Glennie et al. (2019).
There are three major motivations for open spatial models
openCR fits spatial analogues of CJS and JSSA models by maximizing the likelihood. The abundance parameter is density $D$ (animals per hectare) rather than population size $N$.
Recruitment in spatial models may be modelled using parameterizations to those described above for nonspatial models, replacing 'number' by 'density'. The locations at which animals recruit are not modelled.
By definition, the interval between primary sessions is long enough for turnover due to births and deaths. It is also possible that resident animals shift their home ranges (i.e. disperse). Spatial models may either ignore such movement (Gardner et al. 2010, Chandler and Clark 2014, Whittington and Sawaya 2015) or attempt to model it (Ergon and Gardner 2014). There are good arguments for modelling movement:
Data should be provided to openCR.fit
as secr 'capthist' objects. The occasions of a singlesession[^footnote7] dataset are treated as openpopulation temporal samples. For spatial analyses, the capthist object should use a point detector type ('multi', 'proximity' or 'count').
[^footnote7]: The terms 'singlesession' and 'multisession' are here used in the secr sense (secrmultisession). openCR uses these data structures, but interprets them differently as explained here.
openCR mostly uses the terminology of primary and secondary sessions (Pollock 1982) rather than 'session' and 'occasions' as in secr. Where 'session' appears without qualifier it refers to a primary session composed of one or more secondary sessions.
The optional intervals
attribute of the capthist object defines the structure. If intervals are not specified then they default to 1.0 and each occasion is treated as a primary session. If intervals are specified then some may be zero; occasions separated by 'zero' intervals are treated as secondary sessions within the same primary session, as in MARK.
# load code from Appendix 1
# run code from Appendix 1
Fig. 1. Structure of data for openpopulation analysis in openCR. Primary sessions initially correspond to the sessions (components) of a multisession secr capthist object; each primary session may have one or more secondary sessions as numbered (top). For model fitting in openCR.fit
the multisession capthist is 'joined' to form a singlesession capthist with an 'intervals' attribute; nonzero intervals indicate breaks between primary sessions (bottom). The join
step is automatic when a multisession capthist is provided to openCR.fit
if stratified = FALSE
. (See Appendix 1 for code to make this figure).
To construct your own capthist objects 
secr::unRMarkInput
, orread.inp
.Examples of data input code also appear on the help pages for data objects FebpossumCH
, fieldvoleCH
, microtusCH
and dipperCH
.
A multisession capthist object will be converted automatically to a singlesession object using function secr::join
unless stratified = TRUE
(see below). An appropriate intervals attribute is constructed, using the intervals attribute of the multisession object for the intervals between primary sessions (1.0 if not specified), and setting other intervals to zero.
Table 1. Input formats for \textbf{openCR} 2.2
 Input  stratified
 Interpretation 
:::
 singlesession capthist  not used  secondary sessions split into primary sessions by 'intervals' 
 multisession capthist FALSE  single stratum (each 'session' is one primary session) 
 multisession capthist TRUE  multistratum (each 'session' is one stratum) 
From openCR 2.0 onwards any model may be stratified. For stratified models (stratified = TRUE
in openCR.fit
) each session of a multisession capthist object is interpreted as an independent stratum that contributes one component of the log likelihood. Each stratum (session) has its own detectors and capture data. This assumes that primary sessions within each stratum have previously been joined manually in a nominally 'singlesession' capthist. The function stratify
helps you construct stratified capthist objects from collections of singlesession objects.
# run code from Appendix 1
Fig. 2. Structure of data for stratified openpopulation analysis in openCR. Each stratum is a prejoined (singlesession) component of a multisession capthist object. The internal structure and detector may differ between strata.
Stratified models may use 'stratum' as a factorvalued predictor. Groups of strata may be contrasted using stratumlevel covariates as described later.
The various models available in openCR are named to encode the distinctions made in the 'Brief survey'. Names are formed by concatenating four components:
Thus 'JSSAsecrfCL' is a spatial JSSA model parameterized in terms of per capita recruitment $f$ and fitted by maximizing the conditional likelihood (a spatial version of Link and Barker (2005), minus parameter covariation). Any movement model is specified separately with the 'movementmodel' argument of openCR.fit
.
Models of the form 'JSSA...CL' are variations on the PradelLinkBarker models. openCR 2.2 recognises labels of the form 'PLB...' as an alias for each of these models. Thus 'PLBf' is synonymous with 'JSSAfCL', and 'PLBsecrl' is synonymous with 'JSSAsecrlCL'.
Parameters vary with the type of model, as listed below. Each of these primary parameters ('real' parameters in MARK) may also be modelled as a linear combination of predictors on a suitable link scale, allowing the inclusion of covariates and constraints. The coefficients of the parameterspecific linear combinations are called 'beta' parameters in MARK; the likelihood is maximized with respect to the concatenated list of beta parameters.
Table 2. Parameter definitions and default link functions (nonspatial models)
 Parameter  Symbol  Link  Description  ::::  p  $p$  logit  capture probability (recapture probability for CJS)   phi  $\phi$  logit  apparent survival   b  $b$  mlogit  entry probability cf PENT in MARK   f  $f$  log  per capita recruitment rate   gamma  $\gamma$  logit  seniority (Pradel 1996)   lambda  $\lambda$  log  population growth rate (finite rate of increase)  superN  $N$  log  superpopulation size  BN  $B_N$  log  number of entrants  N  $N_j$  log  timespecific population size
* parameters marked with an asterisk are scaled by the interval between primary sessions.
Table 3. Parameters of nonspatial openCR models
 Type  Alias  p  phi  b  f  gamma  lambda  superN  BN  N 
:
 CJS  +  +        
 JSSAbCL  PLBb  +  +  +       
 JSSAfCL  PLBf  +  +   +      
 JSSAgCL  PLBg  +  +    +     
 JSSAlCL  PLBl  +  +     +    
 JSSAb   +  +  +     +   
 JSSAf   +  +   +    +   
 JSSAg   +  +    +   +   
 JSSAl   +  +     +  +   
 JSSAB   +  +       +  
 JSSAN   +  +        + 
Models with type ending in CL are of the PradelLinkBarker type, with aliases as shown.
openCR mostly fits models by modelling capture histories onebyone. An alternative faster method is to evaluate the likelihood expressed in terms of sufficient statistics. Sufficient statistics vary among models, but they are typically counts such as provided by the function JS.counts
. The 'sufficient statistics' approach is not compatible with individual covariates. The nonspatial model types 'Pradel' and 'Pradelg' are implemented in openCR using sufficient statistics (Pradel 1996) and therefore fall outside the main framework (Table 3). They correspond to 'JSSAlCL' and 'JSSAgCL' respectively, and estimate the same parameters as those models. Estimates should coincide except when there are losses on capture. 'Pradel' is parameterized in terms of population growth rate (lambda) and 'Pradelg' is parameterized in terms of seniority (gamma).
Additionally, the function JS.direct
computes classic JollySeber estimates using the sufficient statistics.
Table 4. Parameter definitions and default link functions (spatial models)
 Parameter  Symbol  Link  Description  ::::  lambda0  $\lambda_0$  log  detection function intercept   sigma  $\sigma$  log  detection function scale (m)   z  $z$  log  detection function shape parameter (HHR, HAN, HCG, HVP)   phi  $\phi$  logit  apparent survival   b  $b$  mlogit  entry probability (beta)   f  $f$  log  per capita recruitment rate  gamma  $\gamma$  logit  seniority (Pradel 1996)   lambda  $\lambda$  log  population growth rate (finite rate of increase)  superD  $D$  log  superpopulation density  BD  $B_D$  log  entrants per hectare   D  $D_j$  log  timespecific population density
* parameters marked with an asterisk are scaled by the interval between primary sessions.
Table 5. Parameters of spatial openCR models
 Type  Alias  lambda0  sigma  z  phi  b  f  gamma  lambda  superD  BD  D 
:
 CJSsecr   +  +  +  +        
 JSSAsecrbCL  PLBsecrb  +  +  +  +  +       
 JSSAsecrfCL  PLBsecrf  +  +  +  +   +      
 JSSAsecrgCL  PLBsecrg  +  +  +  +    +     
 JSSAsecrlCL  PLBsecrl  +  +  +  +     +    
 JSSAsecrb   +  +  +  +  +     +   
 JSSAsecrf   +  +  +  +   +    +   
 JSSAsecrl   +  +  +  +     +  +   
 JSSAsecrB   +  +  +  +       +  
 JSSAsecrD   +  +  +  +        + 
 secrCL   +  +  +         
 secrD   +  +  +       +   
Spatial models with type ending in CL have features in common with the PradelLinkBarker models, hence the aliases as shown.
Formulae define a linear model for each 'real' parameter (p, phi, sigma etc.) on the link scale (logit, log etc.). Alternative link functions not shown in Tables 2 and 4 are 'loglog' and 'sin', both as defined in MARK.
The default linear combination for each parameter is a constant, null model (~1, parameter constant over time, unaffected by individual differences etc.). To include other effects build formulae using either predefined (builtin) predictors listed here, or the names of covariates.
Table 6. Builtin predictors ('sessions' refers to primary sessions)
 Predictor  Parameters  Description 
: 
 stratum  all  Factor, one level per stratum (stratified = TRUE
) 
 session  all except 'superN', 'superD'  Factor, one level per primary session 
 t  all except 'superN', 'superD'  synonym of 'session' 
 Session  all except 'superN', 'superD'  Continuous time 
 b  p, phi, lambda0, sigma  learned response (persists across sessions) 
 B  p, lambda0, sigma  transient (Markovian) response across sessions
 bk  p, phi, lambda0, sigma  detectorspecific learned response (persists across sessions) 
 bsession  p, lambda0, sigma  learned response within sessions 
 Bsession  p, lambda0, sigma  transient (Markovian) response within sessions
 bksession  p, lambda0, sigma  detectorspecific learned response within sessions 
 Bksession  p, lambda0, sigma  detectorspecific transient (Markovian) response within sessions
 h2  all except abundance  2class finite mixture 
 h3  all except abundance  3class finite mixture 
 age  all except abundance  age factor 
 Age  all except abundance  linear effect on age 
 Age2  all except abundance  linear effect on age$^2$ 
Differences among the various learned responses may be understood by examining their effect on the parameter index array (PIA). This table illustrates the PIA slice corresponding to an individual with the nonspatial detection history shown (4 primary sessions, each of 4 secondary sessions). The values '1' and '2' refer to different parameter combinations, most commonly to levels of lambda0.
 Detection history :  0100 0000 0000 0100   :   ~bsession  1122 1111 1111 1122  persistent within primary session   ~Bsession  1121 1111 1111 1121  transient within primary session   ~b  1122 2222 2222 2222  persistent   ~B  1122 2222 1111 1122  transient across primary sessions 
IMPORTANT NOTE: Learned response predictors ('b', 'bsession' etc.) were redefined in openCR 1.3.0. Models fitted with earlier versions should be refitted.
The rules for covariates largely follow secr (secroverview.pdf). Covariates may be at the level of stratum, primary session, secondary session (detection parameters only), individual (CL models only), or detector (spatial models only). Further complexity may be modelled by providing custom design data cutting across these categories (see below).
Individual and detector covariates are named columns in the 'covariates' attributes of the respective capthist and traps object. Covariate names should differ from the builtin predictors (Table 6).
Stratum covariates are provided to openCR.fit
in the argument 'stratumcov'. That should be a dataframe with one row per stratum; the name of any column may be used in a model formula.
Primary session covariates are provided to openCR.fit
in the argument 'sessioncov', rather than associated with a data object. If 'sessioncov' is a vector (length equal to number of primary sessions) rather than a dataframe then it may be referenced as 'scov' in model formulae. For stratified data, 'sessioncov' may be a list with one component per stratum (the lazy option of providing a single vector or dataframe works only if all strata have the same sessions).
Covariates for detection parameters in secondary sessions are provided in the 'timecov' argument. If 'timecov' is a vector (length equal to total number of secondary sessions) rather than a dataframe then it may be referenced as 'tcov' in model formulae. For stratified data, 'timecov' may be a list with one component per stratum (the lazy option of providing a single vector or dataframe works only if all strata have the same primary and secondary sessions).
The types 'secrD' and 'secrCL' cause openCR.fit
to treat the data as if from a closed population (no mortality, no recruitment, no movement); the intervals attribute is ignored.
msk < make.mask(traps(captdata), buffer = 100, type = 'trapbuffer') fit_secr < secr.fit(captdata, detectfn = 'HHN', mask = msk, trace = FALSE) fit_openCR < openCR.fit(captdata, detectfn = 'HHN', mask = msk, type = 'secrD') # massage the predict.openCR results to the same format as predict.secr pred_openCR < plyr::rbind.fill(predict(fit_openCR)) pred_openCR < pred_openCR[c(2,1,3), !(names(pred_openCR) %in% c('stratum','session'))] rownames(pred_openCR) < fit_secr$realnames # compare estimates predict(fit_secr)[,1] pred_openCR
# compare timings in seconds c(secr = fit_secr$proctime, openCR = fit_openCR$proctime)
The maximised log likelihoods differ because openCR does not include the multinomial constant. secr has function logmultinom
that lets us add it back:
# compare maximised log likelihoods c(secr.logLik = logLik(fit_secr), openCR.logLik = logLik(fit_openCR) + logmultinom(captdata))
Two and threeclass finite mixtures (h2, h3) allow for individual heterogeneity in detection and turnover parameters (Pledger et al. 2003, 2010). Using one of these predictors in a formula causes a further real parameter 'pmix' to be added. pmix is the proportion in latent mixture class 2 for h2, and the proportions in classes 2 and 3 for h3 (the proportion in class 1 is obtained by subtracting from 1). The implementation in openCR assumes that class membership applies across all parameters. The posterior probabilities of class membership for all detected individuals are returned as the 'posterior' component of the fitted model.
Finite mixture likelihoods are prone to multimodality. Misleading estimates result when the numerical maximization settles on a local maximum (see also [secrfinitemixtures.pdf].
If age is modelled as a factor then it is useful to group older animals in a maximum age class ('maximumage'). 'minimumage', 'maximumage' and 'initialage' are optional components of the 'details' argument of openCR.fit
. 'initialage' can name an individual covariate to avoid the assumption that all animals are the minimum age at first detection.
The agecov argument of openCR.fit
may be used to specify a recoding of numerical age, analogous to the sessioncov argument of secr.fit
. Thus numerical ages may be bracketed into 'young', 'middle' and 'old'. The length of agecov should match the number of numerical ages (maximumage  minimumage + 1).
For a quadratic relationship with age, specify an additive model with both Age and Age2 terms (e.g., model = phi ~ Age + Age2).
We have seen the role of the intervals attribute in defining primary and secondary sessions. Betweensession intervals need to be specified only if they vary, or if you would like rates (phi, gamma, lambda, f) to be reported in time units other than the (implicitly constant) sampling interval. Scaling from the standardised parameter $\theta_j$ to the intervalspecific value $\theta^\prime_j$ uses $\theta^\prime_j = \theta_j^{T_j}$ where $\theta_j$ is one of $\phi_j$ or $\lambda_j$, and $T_j$ is the duration of interval $j$.
Scaling $\gamma$ follows the same pattern except that the relevant duration for $\gamma_j$ is $T_{j1}$. Scaling per capita recruitment $f_j$ is more tricky. We use $f^\prime_j = (\phi_j + f_j)^{T_j}  \phi_j^{T_j}$.
Occasionally there is a need for covariates that do not relate specifically to individuals, sessions or detectors, and are not included as canned predictors. For this you must construct your own dataframe of design data and pass it as the 'dframe' argument of openCR.fit
. Design data are used as input to the model.matrix
function (the 'data' argument); model.matrix
generates the design matrix for each real parameter. Design data are usually constructed internally in openCR.fit
from named covariates and other predictors that appear in model formulae; if 'dframe' is provided then the internally constructed design data are added as extra columns, overwriting any custom columns of the same name. The same design dataframe is used for all parameters.
Constructing 'dframe' is fiddly. The dataframe should have one row for each combination of unique capture history, secondary session, detector and latent class (mixture). For nonspatial models without finite mixtures this collapses to one row for each capture history and secondary session. The order of rows follows that of the elements in an array with dimensions ($n$, $S$, $K$, $X$) for $n$ unique capture histories, $S$ secondary sessions, $K$ detectors and $X$ latent classes[^footnote8]. The secr function insertdim
can help to expand data into the correct row order.
[^footnote8]: This rectangular (or cuboidal) configuration includes cells that are redundant and unused for a particular model type (e.g., cells corresponding to sessions at or before first capture in CJS models). However, the full complement of rows is required in dframe.
A warning: by default openCR.fit
replaces the input capthist with a more compact version using only unique capture histories (the number of each is kept in the individual covariate 'freq'; see the function squeeze
). Design data are in terms of the 'squeezed' capture histories.
In this example we define a function to construct custom design data for a learned response.
makedf.b < function (ch, spatial = FALSE, nmix = 1, naive = FALSE) { R < 1 # assume single stratum ch < squeeze(ch) # Construct matrix of logical values TRUE iff caught before detected < apply(abs(ch),1:2,sum)>0 detected < t(apply(detected, 1, cumsum)>0) if (naive) b < rep(FALSE, prod(dim(ch)[1:2])) else b < t(apply(detected, 1, function(x) {x[which.max(x)] < FALSE; x})) # For a simple nonspatial case: data.frame(customb = as.vector(b)) # More generally: n < nrow(ch) S < ncol(ch) K < if (spatial) dim(ch)[3] else 1 data.frame(customb = insertdim(b, c(2,3,1), c(R,n,S,K,nmix))) }
Now compare the result with the canned predictor 'b' for a persistent learned response.
ovenj < join(ovenCH) fitb < openCR.fit(ovenj, model = p ~ b) fitbc < openCR.fit(ovenj, model = p ~ customb, dframe = makedf.b(ovenj)) AIC(fitb, fitbc)
Our custom model gives exactly the same result as the canned predictor 'b' when type = 'CJS' because the precise secondary session of first capture is irrelevant for CJS models (recaptures are modelled only for subsequent primary sessions unless details$CJSp1 == TRUE
).
Discrepancies can arise with nonCJS models because these account for animals never detected. The corresponding likelihood component uses a distinct design matrix for a 'naive' animal. To customize nonCJS models a separate dframe should be provided that applies to naive animals:
fitb2 < openCR.fit(ovenj, model = p ~ b, type = 'JSSAfCL', start = fitb) fitbc2 < openCR.fit(ovenj, model = p ~ customb, type = 'JSSAfCL', dframe = makedf.b(ovenj), dframe0 = makedf.b(ovenj, naive = TRUE)) AIC(fitb2, fitbc2)
An ad hoc adjustment for transience may be programmed as follows (cf Pradel et al. 1997).
makedf.resident < function (ch, spatial = FALSE, nmix = 1) { nstrata < 1 # assume single stratum ch < squeeze(ch) n < nrow(ch) S < ncol(ch) K < if (spatial) dim(ch)[3] else 1 primary < primarysessions(intervals(ch)) detected < apply(abs(ch),1:2,sum)>0 nprimary < apply(detected, 1, function(x) length(unique(primary[x]))) data.frame(resident = insertdim(nprimary>1, 1, c(nstrata, n, S, K, nmix))) }
A simpler approach is to code an individual covariate that scores whether an individual was detected in more than one primary session.
addresidentcov < function (ch) { primary < primarysessions(intervals(ch)) detected < apply(abs(ch), 1:2, sum)>0 nprimary < apply(detected, 1, function(x) length(unique(primary[x]))) covariates(ch) < data.frame(residentcov = nprimary>1) ch }
Results are identical:
ovenj < join(ovenCH) ovenj < addresidentcov(ovenj) fitnull < openCR.fit(ovenj, model = phi ~ 1) fitcov < openCR.fit(ovenj, model = phi ~ residentcov) fitdf < openCR.fit(ovenj, model = phi ~ resident, dframe = makedf.resident(ovenj)) fits < openCRlist(fitnull, fitcov, fitdf) AIC(fits) pred < predict(fits, newdata = data.frame(resident = TRUE, residentcov = TRUE)) do.call(rbind, lapply(pred, '[[', 'phi'))
Hines et al. (2003) suggested extending the definition of residence to include animals captured at least $d$ days apart within a primary session; either of the approaches here may be modified accordingly. Here is the code for two individual covariates:
addresidentcov2 < function (ch, d = 1) { primary < primarysessions(intervals(ch)) secondary < secondarysessions(intervals(ch)) detected < apply(abs(ch), 1:2, sum)>0 nprimary < apply(detected, 1, function(x) length(unique(primary[x]))) dsecondary < apply(detected, 1, function(x) max(by(secondary[x], primary[x], function(y) diff(range(y))))) covariates(ch) < data.frame(residentcov1 = nprimary>1, residentcov2 = nprimary>1  dsecondary>=d) ch }
Factor predictors take a number of discrete values (levels). These are usually represented by columns of 0's and 1's in the design matrix, where the number of columns (and coefficients) relates to the number of levels. The default in R is to use 'treatment contrasts'; one coefficient describes a reference class (level) and other coefficients represent the effect size (difference from the reference class on the link scale). By default the first level is used as the reference: for time effects (t, session) the first primary session is the reference level[^footnote9].
This may lead to trouble if the parameter is not identifiable in the reference class. One workaround is to specify a session covariate with differently ordered levels. Another is to switch to dummy variable coding in which each coefficient represents the magnitude of one real parameter on the link scale (useful in itself). Dummy variable coding is achieved by removing the intercept from the formula (1), assuming the default contrast function for factor coding (contr.treatment
; check with options()$contrasts
). The following model fits yield the same estimates of 'real' parameters and the same loglikelihood, but with different 'beta' parameters:
fit0 < openCR.fit(ovenCH, model = p~t) fitd < openCR.fit(ovenCH, model = p ~ 1+t) coef(fit0) coef(fitd)
Dummy variable coding has proved useful for avoiding some maximization problems. From openCR 2.1.0, dummy variable coding can be selected with the 'details' argument 'dummyvariablecoding'. This updates the model to remove the intercept, and assigns the default starting value across all levels of a factor (rather than zero for nonreference levels). The following fit is therefore equivalent to the preceding fitd
.
fitd2 < openCR.fit(ovenCH, model = p~t, details = list(dummyvariablecoding = 't')) coef(fitd2)
[^footnote9]: This does not apply for times when a parameter can never be estimated  for example, openCR understands that seniority (gamma) is not estimated for the first session and uses the second session for the reference level.
Suppose you wish to estimate the average of a parameter across levels of a factor such as time (session). Cooch and White (2019 Section 6.15) advocate modifying the design matrix so that one beta parameter (coefficient) relates directly to the mean. This is achieved very simply in openCR.fit
[^footnote10] by setting the contrast function for the factor to contr.sum
in the details
argument[^footnote11]. With the resulting factor coding the first coefficient corresponds to the mean. Applying this to estimate the average timespecific survival rate for the dippers assuming constant recapture probability:
fit < openCR.fit(dipperCH, model = phi~t, details = list(contrasts = list(t = contr.sum))) invlogit(coef(fit)['phi',c('beta','lcl','ucl')])
The mean is backtransformed from the link scale. This results in some bias owing to the nonlinearity of link functions other than the identity function. Cooch and White take the position that the bias is often ignorable.
[^footnote10]: This also works in secr.fit
.
[^footnote11]: Helmert contrasts (contr.helmert
) also yield the mean as the first coefficient, but the coding is more obscure.
Potential movement of home ranges between primary sessions (= dispersal) is a critical part of openpopulation models (Efford and Schofield 2022). The argument movementmodel
of openCR.fit
allows the possibilities in Table 7. Two of these do not model movement at all. The default 'static' is a null model in which each animal retains the same home range. The 'IND' option models the locations of an animal independently in each primary session; information is sacrificed and the implied movement depends on the size of the habitat mask.
The remaining options ('BVN', 'BVE', 'BVC', 'BVT', etc.) fit a dispersal kernel (Nathan et al. 2012) to represent movement between primary sessions. This usually requires at least one more parameter to represent the spatial scale of dispersal.
Table 7. Models for movement between primary sessions.
 Movement model  Parameters Description (aliases in parentheses) 

 static
 0  Centres constant across primary sessions 
 BVN
 1  Bivariate normal ('normal', Gaussian)
 BVE
 1  Bivariate Laplace kernel ('exponential') 
 BVC
 1  Bivariate Cauchy kernel 
 BVT
 2  Bivariate $t$ ('t2D', 2Dt)
 RDE
 1  Exponential distance moved 
 RDG
 2  Gamma distance moved 
 RDL
 2  Lognormal distance moved 
 UNI
 0  Uniform within arbitrary radius 
 IND
 0  Centres unconstrained within habitat mask 
 (user function)  0,1,2  Usersupplied function (ncores = 1
only)
Note: 'BVN', 'BVE' and 'BVT' kernels were previously designated 'normal', 'exponential' and 't2D' respectively, and these names are still recognised.
Note: The online vignette openCRkernel.pdf covers movement kernels in detail. This section is retained for historical reasons and will be removed in future.
All movement kernels in openCR are radially symmetrical. Relative probability of movement is specified in terms of radial distance $r$ from the point of origin (Table 8). Four of the builtin kernels (BVN, BVE, BVC, BVT) are defined directly as bivariate probability density functions $g(r)$. Three others (denoted RDE, RDG, RDL) are defined indirectly by the univariate distribution of distance moved $f(r)$ where a point on the kernel has polar coordinates $(r, \theta)$ assuming direction $\theta$ uniform on $(0,2\pi)$ (Cousens et al. 2008, Ergon and Gardner 2014).
The extent of the kernel is controlled by the argument 'kernelradius' that gives the radius in terms of mask cells. Cellspecific probabilities are normalised so that they sum to 1.0 across the kernel. Dispersal probability effectively falls to zero at the boundary of the kernel, so the kernel radius is a critical part of the model. The uniform 'UNI' kernel has no parameters but depends critically on the userspecified kernel radius.
Table 8. Kernel probability density functions. 'move.a' and 'move.b' are the names used in openCR for scale and shape parameters, as indicated in the table. Based in part on Nathan et al. (2012, Table 15.1) and Clark et al. (1999) with adjustment for parameterisation in openCR. $g(r) = f(r) / (2 \pi r)$.
 Kernel  move.a
 move.b
 $f(r)$  Mean $r$* 

 BVN  $\alpha$    $\frac{r}{ \alpha^2} e^{ \frac{r^2}{2\alpha^2} }$ $\alpha \sqrt{\frac{ \pi}{2}}$ 
 BVE  $\alpha$    $\frac{r}{\alpha^2} e^{ \frac{r}{\alpha}}$ $2\alpha$ 
 BVT  $\alpha$  $\beta$  $\frac{2 \beta r}{\alpha^2} \left( 1 + \frac{r^2}{\alpha^2} \right)^{(\beta+1)}$ $\alpha \frac{\sqrt \pi}{2} \frac{\Gamma \left( \beta0.5 \right)}{\Gamma (\beta)}$
 RDE  $\alpha$    $\frac{1}{\alpha} e^{\frac{r}{\alpha}}$  $\alpha$ 
 RDG  $\alpha$  $\beta$  $\frac{1}{\Gamma(\beta)\alpha^\beta} r^{\beta1} e^{ \frac{r}{\alpha} }$  $\alpha\beta$ 
 RDL  $\exp \mu$  $1/(e^{\sigma^2}1)$  $\frac{1}{r \sigma \sqrt{2 \pi}} e^{ \frac{(\ln(r)\mu)^2}{2\sigma^2} }$  $e^{\mu + \frac{\sigma^2}{2}}$ 
 UNI
      
 (user)  $a$  $b$    
* Continuous, untruncated, kernel. Expected values for the discretized and truncated kernel will be less (see summary.kernel
).
\vspace{12pt}
The 'BVT' kernel is the same as '2Dt' of Clark et al. (1999) and Nathan et al. (2012). The parameter $\alpha$ (move.a
) corresponds to $a$ in Nathan et al. (2012) and $\sqrt u$ in Clark et al. (1999); the parameter $\beta$ (move.b
) corresponds to $b1$ in Nathan et al. (2012) and $p$ in Clark et al. (1999). Defining move.b
as $\beta \equiv b1$ is handy because the default link for move.b
(log) then ensures $b>1$. The degrees of freedom of the corresponding $t$distribution are given by $\nu = 2\beta$.
The `BVT' kernel approaches bivariate normal as $\beta \to \infty$ and Cauchy as $\beta \to 0$ (e.g., Clark et al. 1999). Clark et al. (1999 p. 1485) found it hard to fit this kernel to seed dispersal data. The mean is undefined for $\beta \le 0.5$.
openCR 2.2 introduces zeroinflated versions of kernels otherwise defined with one parameter or none. Zeroinflated kernels use the suffix 'zi', hence 'BVNzi', 'BVEzi', 'RDEzi', 'UNIzi'. The kernelfree independent movement model 'IND' also has a zeroinflated form 'INDzi' that is not strictly independent or uncorrelated. Each of these models has an additional zeroinflation parameter (move.b or move.a depending on whether the base kernel does or does not already have a parameter). Zeroinflated kernels often fit well, but it is common for the fitted scale parameter 'move.a' of BVNzi, BVEzi and RDEzi models to become large and essentially unidentifiable as the kernel for $r>0$ flattens.
A kernel function may be specified by the user and passed in the argument movementmodel
. The function should have argument $r$, and optionally $a$, or $a$ and $b$ (the last two correspond to openCR parameters move.a
and move.b
) It should return a vector of values one for each element of $r$, although length(r) = 1 when the likelihood is evaluated in C++ (details$R = FALSE
, the default). The code should give a valid result when $r = 0$ that will be used for the origin cell. With the default link ('log' for both move.a
and move.b
) there is no risk of $a \leq 0$ or $b \leq 0$.
A big problem with standard kernels as defined in openCR <2.0.0 is that the number of cells increases with the square of the radius. Processing time is roughly proportional to the number of cells, and kernels with many cells fit slowly. openCR 2.0.0 introduces novel 'sparse' kernels that include only those grid cells that lie on 4 axes (NS, EW, NWSE, NESW) (Efford 2022b). The number of cells then increases only linearly with radius. Cellwise movement probabilities are adjusted so that the distribution of dispersal distances is almost unchanged (essentially multiplying by $2\pi r$ at radius $r$, and adjusting cells on the oblique axes by $\sqrt 2$).
Sparse kernels are obtained by setting sparsekernel = TRUE
when fitting a model with openCR.fit()
. Here is an example.
par(mar = c(3,1,4,5)) k < make.kernel(movementmodel = 'BVN', kernelradius = 10, spacing = 10, move.a = 40, sparse = TRUE, clip = TRUE) plot(k) symbols(0,0, add = TRUE, circles = 100, inches = FALSE)
Note that the maximum on each axis is no longer at the centre. In fact, the central cell is assigned zero weight because $r = 0$ (this is undesirable and may be corrected in future). Each oblique arm in the example has only 7 cells; these cells have higher weighting to avoid orientation bias.
Somewhat surprisingly, sparse kernels appear to work about the same as full kernels, only faster.
When a kernel is applied to cells near the edge of a habitat mask some projected movements will lie outside the mask. This creates a problem for the model. Kernel cell values are probabilities summing to one; the cell probabilities of a truncated kernel will no longer be true probabilities and results are prone to bias.
openCR.fit
offers two approaches to resolve this problem:
If the mask is rectangular, the truncated cells (and their probabilities) may be 'wrapped' to the opposing edge of the mask. This works fine if the kernel is not too large. Wrapping does not impose a computational burden. A rectangular mask is generated by make.mask
with the default type = 'traprect'
.
For any mask, the cell probabilities of a truncated kernel may be scaled (normalized) so that they sum to 1.0. This requires substantial additional computation.
The edge method is chosen by setting the argument 'edgemethod' in openCR.fit
; the options are 'truncate' (default in 1.5.0 and later), 'wrap', and 'none'. Wrapping is fast, but it will cause an error if the mask is not rectangular. If a rectangular mask does not make sense (e.g., because the habitat is patchy) then you must use edgemethod = 'truncate'
for unbiased estimates.
Prior to 1.5.0 there was no adjustment (movement truncated without normalization, equivalent to edgemethod = 'none'
in later versions) and estimates from movement models could be biased because the probability of a null (allzero) history was estimated incorrectly.
A kernel may be constructed with make.kernel
and visualised with the plot
method. Use the summary
method to obtain a terse description.
par (mar = c(3,3,4,6), cex = 0.9) k < make.kernel (movementmodel = 'BVN', spacing = 10, move.a = 40, clip = TRUE) plot(k, contour = TRUE) summary(k)
Use the secr function spotHeight(k)
to display cell values on the plot.
The independent option 'IND' (previously 'uncorrelated') is not recommended. It discards information on the continuity of home ranges between primary sessions, and estimates may vary with the (often arbitrary) extent of the habitat mask.
Kernelbased movement models require extreme care. Definitive advice cannot yet be given on the safe use of these models. Longdistance movements will usually be poorly sampled and poorly modelled.
Userdefined functions cannot be used with multithreaded C++, so they will be slow to fit; always set ncores = 1
.
When one of the kernel models is used (movementmodel
not 'static', 'IND' or 'INDzi') the movement model may be modified by weighting settlement according one or more local (mask) covariates. The weights are given by a new parameter 'settle'. By default the logarithm of 'settle' is a linear combination of mask covariates with no intercept. 'settle' may also vary by
The settlement model is invoked by setting details = list(settlemodel = TRUE)
.
This code compares a model with uniform settlement to one in which settlement varies with the 'y' coordinate of the mask. The data are the ovenbird mistnetting data compressed to a single occasion per year.
ovenCHb < reduce(ovenCHp, by = 'all', outputdetector = 'count') msk < make.mask(traps(ovenCHp[[1]]), buffer = 500, spacing = 40, type = 'trapbuffer') # uniform settlement fit0 < openCR.fit(ovenCHb, type = 'PLBsecrf', mask = msk, binomN = 1, movementmodel = 'BVN', details = list(settlemodel = FALSE)) # logarithmic NS gradient in settlement fit1 < openCR.fit(ovenCHb, type = 'PLBsecrf', mask = msk, binomN = 1, movementmodel = 'BVN', details = list(settlemodel = TRUE), model = settle~y)
load('settlement.RData')
In this example there is no evidence for a NS gradient in settlement. The intercept is removed automatically when the link function is the default ('log'), as shown in the resulting formula for 'settle'.
AIC(fit0, fit1)[,6] coef(fit1)
The cumulative effect of differential settlement may be visualised with function cumMove
by specifying a single mask covariate with values in the range 0 to 1.
NOTE: Settlement models are a novelty in openCR 2.2.2 and their value is limited because survival and recruitment are not allowed to vary across space. The implementation may change.
Various derived parameters may be computed from a fitted model. Specifically,
HTbysession = TRUE
) the HT estimate may be based on the number detected in each session and the corresponding sessionspecific estimates of $p$ or $a$.[^footnote12]: However, the effect of a constraint (e.g., parameter constant over sessions) will vary depending on the parameter to which it is applied.
Both goals are served by the derived
method for openCR
objects. Among other outputs, this generates a summary table with point estimates of all relevant parameters. We demonstrate this with a new dipper model, fitted using conditional likelihood:
dipperCL < openCR.fit(dipperCH, type = 'JSSAlCL', model = list(lambda~t, phi~t)) # only these parameters are in the model and estimated directly, names(predict(dipperCL)) # but we can derive b, f, gamma and N, as well as the superpopulation N d < derived(dipperCL) print(d, digits = 3, legend = TRUE)
The print
method for objects from derived
provides some control over formatting, as shown. Use the Dscale argument to change area units (spatial models only).
derived
does not yet provide deltamethod SE or confidence intervals for derived parameters. A reliable workaround for abundance paramaters (N, D)[^footnote13] is to (i) infer the point estimates with derived
[^footnote14], (ii) assemble a start vector on the link scale(s) for an equivalent fulllikelihood openCR.fit
model that includes the derived abundances, and (iii) run openCR.fit
with method = "none"
to compute the hessian at the MLE, and hence the full variancecovariance matrix.
[^footnote13]: This may sometimes be feasible for derived recruitment parameters, but given the doubts introduced by differing constraints (e.g. constant $f$ vs constant $\lambda$) it is better just to refit the model.
[^footnote14]: These are also the MLE when distribution = "poisson" (e.g., Schofield and Barker 2016).
The secr functions sim.popn
and sim.capthist
provide the means to generate spatial openpopulation data with known survival probability, population trend $\lambda$ and detection parameters. Open population data are generated by setting nsessions > 1 in sim.popn
and specifying a value for $\lambda$. Turnover settings are controlled by components of the 'details' argument of sim.popn
. The secr help page ?turnover should be consulted. sim.capthist
should be called with renumber = FALSE (otherwise individual capture histories cannot be matched across primary sessions).
Use the openCR function sim.nonspatial
to generate nonspatial openpopulation data. openCR also provides these functions to streamline simulation 
 Function  Purpose 

 runsim.nonspatial
 Generate data with sim.nonspatial
and fit models using openCR.fit

 runsim.spatial
 Generate data with sim.popn
and sim.capthist
, and fit models using openCR.fit

 sumsims
 Summarise list output from runsim.nonspatial
or runsim.spatial

runsim.nonspatial
and runsim.spatial
are essentially wrappers; the user must provide appropriate argument values for each of the nested functions.
It is common for some sessionspecific parameters of open capturerecapture models to be nonidentifiable, either for structural reasons or because the particular dataset is uninformative (e.g., Gimenez et al. 2004).
The main diagnostic is the rank of the Hessian matrix. If the rank is less than the number of parameters then the model is not fully identifiable and the estimates of some parameters will be confounded or unreliable. Matrix rank is determined numerically by counting nonzero eigenvalues. Computed eigenvalues of nonidentifiable parameters may appear as small positive numbers, so it is necessary to apply an arbitrary numerical threshold.
Exactly which parameter estimates are unreliable can usually be discerned from computed variances (SE and confidence intervals). Data cloning (Lele et al. 2010) is also helpful; function cloned.fit
implements the method for nonspatial models.
Sessionspecific turnover parameters may become nonidentifiable if home ranges are allowed to move freely between primary sessions (movementmodel = 'uncorrelated'
). Intuitively, this is because radical changes in individual detection probability (due to proximity to detectors) cannot be separated from mortality and recruitment.
Bad estimates (zero, very large, close to starting values or zero variance) may merely indicate a problem with the maximization algorithm rather than nonidentifiability.
Numerical maximization of the likelihood requires appropriate starting values for the parameters. If starting values are poor then initial evaluations of the likelihood may return an infinite value, or otherwise provide inadequate direction for the numerical algorithm.
openCR.fit
provides a mechanism for recycling earlier estimates as starting values: simply provide the name of a previously fitted model as the start
argument. Parameters shared between the models will be set to the old estimates, while unmatched parameters will be set to defaults. A list of two previous models may be provided; values from the first take precedence.
Variance estimation based on the Hessian matrix fails if the estimate lies on a boundary of the parameter space. Computed SE are then extreme, and confidence limits are implausible. This commonly happens when apparent survival (phi) approaches 1.0. Boundary estimates are more benign than other reasons for failure (the estimates themselves may be reliable). Alternative methods for variance estimation in this case have not been implemented.
Using the "sin" link for parameters bounded by 0 and 1 (the probability parameters p and phi) can be helpful.
The default method for maximizing the likelihood function is NewtonRaphson as implemented in the R function nlm
. This relies on numerical gradient estimates, which can cause trouble. Avoid gradient estimation entirely by using the somewhat slower 'NelderMead' method of function optim
e.g.,
fitnr < openCR.fit(ovenCH, type = 'JSSAlCL', model = list(phi ~ t, lambda~t)) fitnm < openCR.fit(ovenCH, type = 'JSSAlCL', model = list(phi ~ t, lambda~t), method = "NelderMead", details = list(control = list(maxit = 5000)))
The default maximum number of likelihood evaluations for the NelderMead algorithm (500) is often too small and results in a "probable maximization error" warning. Here we increase it to 2000 by setting the details argument "control" that is passed to optim
.
Somewhat alarmingly, the NM algorithm settles on a lower log likelihood and different estimates:
AIC(fitnm,fitnr)
We can fix that by feeding NelderMead the starting values from another model:
fitnm < openCR.fit(ovenCH, type = 'JSSAlCL', model = list(phi ~ t, lambda~t), method = "NelderMead", details = list(control = list(maxit = 2000)), start = fitnr) AIC(fitnm,fitnr)
In the longer term, better maximizers are needed.
Exceeding the iteration limit of nlm
results in code 4. The default number of iterations was increased to 300 in openCR 2.1.0 from the previous default of 100. Each iteration in nlm
completes one evaluation of the gradient of the likelihood function, and requires multiple likelihood evaluations. For many parameters 300 iterations may still not be enough: increase it with the details argument control
(e.g., details = list(control = list(iterlim = 500))
). See the help for optim
for the relevant control settings of other maximizers such as 'NelderMead'.
Boundary values of 'beta' coefficients may sometimes be avoided by changing the default factor coding, particularly for sessionspecific estimates (~t, ~session). This may be tackled manually, but from 2.1.0 there is a shortcut using the details
argument 'dummyvariablecoding'. See the section above on Factor coding.
Spatial models are slow to fit. Consider these options
derive
can give estimates of abundance (superN, N, superD, and D) from PLB models, as well as alternative measures of recruitment.openCR $\ge$ 1.2 uses multiple threads to run some calculations in parallel. Multithreading uses RcppParallel. A couple of tuning parameters are available. The number of threads is set with the 'ncores' argument of openCR.fit
, or with the setNumThreads
function of secr that sets the environment variable RCPP_PARALLEL_NUM_THREADS. By default openCR $\ge$ 1.5.0 uses only 2 cores, for compliance with CRAN rules. You can increase this up to the number of (virtual) cores available (i.e. 8 on a quadcore desktop with hyperthreading), or some lesser number if you want to multitask:
# RCPP_PARALLEL_NUM_THREADS # recommended for quadcore Windows PC setNumThreads(7)
The `grain size' tuning parameter (see RcppParallel) may be varied with details$grain, but it seems to have little effect.
Full models (not CL or Pradel) include superpopulation size $N$ as a variable. The default in openCR for both nonspatial and spatial models is to treat $N$ as a Poisson variable, from which it follows that the number of individuals detected at least once ($n$) is also Poisson. This is also the default in secr. However, estimates from POPAN models in MARK treat $N$ as fixed and $n$ as binomial. The assumption of fixed $N$ leads to narrower confidence intervals and estimates of detection and turnover parameters that differ slightly from conditional likelihood models (see e.g. Schofield and Barker 2016). To obtain JSSA estimates from openCR that match those from MARK it is necessary to set distribution = "binomial"
.
Several examples of analyses with openCR are given in the associated vignette openCRexamples.pdf. These use data already formatted as secr capthist objects in R; the objects are provided in one or other package. All are available immediately openCR is loaded with library
. Each has its own help page.
Table 9. Data objects in openCR. 'RD' indicates robust design with multiple secondary sessions. See openCRexamples.pdf for references.
 Data object  Spatial  RD  Species etc.  Source  :::::  microtusCH etc.  No  Yes  Meadow vole Microtus pennsylvanicus USA  Nichols et al. (1984), Williams et al. (2002)   FebpossumCH  No  Yes  Brushtail possum Trichosurus vulpecula New Zealand  M. Efford unpubl.   dipperCH  No  No  European dipper Cinclus cinclus France  Lebreton et al. (1992), MARK   gonodontisCH  No  No  Moth Gonodontis bidentata England  Bishop et al. (1978), Crosbie (1979)   fieldvoleCH  Yes  Yes  Field vole Microtus agrestis* Norway  Ergon and Lambin (2013) 
Table 10. Multisession data objects in secr.
 Data object  Spatial  RD  Species etc.  Source  :::::  OVpossumCH  Yes  Yes  Brushtail possum Trichosurus vulpecula New Zealand  M. Efford unpubl.   ovenCHp  Yes  Yes  Ovenbird Seiurus aurocapilla USA  D. Dawson and M. Efford unpubl. 
This is generally an undeveloped field for spatially explicit capturerecapture models. Demonstrating that assumptions were not satisfied may also be of no consequence: we would usually ignore such a finding if the estimator is reasonably robust.
For CormackJollySeber (nonspatial) models there is an established suite of tests following Burnham et al. (1987). The tests have been implemented in the UCARE software of Choquet et al. (2009), recently translated into R by Gimenez et al (2018) as package R2ucare. Program RELEASE (Burnham et al. 1987) also implements the core CJS tests and is available through MARK.
The openCR function ucare.cjs
is a wrapper for relevant functions in R2ucare, which should be installed. We briefly demonstrate it here for the dipper data of Marzolin (1988).
if (requireNamespace("R2ucare")) ucare.cjs(dipperCH, verbose = FALSE, by = 'sex')
This invocation of ucare.cjs
calls the R2ucare functions test3sr
, test3sm
, test2ct
, test2cl
and overall_CJS
for each sex and provides a condensed report. For interpretation see the original papers, the R2ucare vignette, and Chapter 5 of the MARK book (Cooch and White 2019). Lebreton et al. (1992: 86) indicate only Test 3SR is meaningful for these data (see also openCRexamples.pdf).
openCR does not do
ucare.cjs
(above).Parameter counting and overdispersion adjustment are probably the most critical omissions. See Cooch and White (2019) for detailed coverage in the context of MARK.
Defaults for some arguments differ between openCR.fit
and secr.fit
. For openCR.fit

trace = FALSE
details$multinom = FALSE
)distribution
has been elevated to a full argument rather than merely a component of details
. This argument describes the distribution of the number of individuals detected (default distribution = "poisson") (see here).
When details$LLonly = TRUE
, openCR.fit
returns a vector with the log likelihood in position 1, followed by the named starting values of the coefficients (beta parameters) (secr.fit
returns only the log likelihood).
In secr the argument CL
is used in secr.fit
to switch between full and conditionallikelihood models. In openCR conditionallikelihood models are given a separate type
with the suffix CL (or use PLB alias).
The predictor 't' is used in secr models to indicate a factor with one level for each secondary session. In openCR it is a synonym for 'session', i.e. a factor with one level for each primary session. This is consistent with the use of 't' in Lebreton et al. (1992) and makes for more compact model specification. In the unlikely event that you want to code a model with one level for each secondary session in openCR, use the 'timecov' argument.
Arguments to be passed to nlm() cannot merely be appended as in the \dots argument of secr.fit(), but must be passed as a named list in the details argument control
. See here for an example.
Parts of openCR are coded in C++, via the R package Rcpp. The Rcpp interface requires less copying of data, and enables the use of multiple threads via RcppParallel. openCR also duplicates some C++ functions in native R code, which is useful for debugging. Select the R version by setting details = list(R = TRUE)
in openCR.fit
. This currently works for most models except those with detector type 'multi' and some exotic movement models.
Strata (openCR >=2.0) are analogous to sessions in secr in that they are treated as independent with no redetections of animals between strata. The total loglikelihood in openCR is the sum of stratum log likelihoods, just as the total is the sum of session loglikelihoods in secr.
These features of secr are not available in openCR
make.table
may do the job)The nonspatial capability of openCR largely duplicates MARK and RMark. Several of the nonspatial model types have exact matches in MARK (Table 11).
Table 11. Relationship of nonspatial openCR models to MARK model types
 openCR type  MARK model  Reference  :  CJS  CJS  Seber (1982)   JSSAb  POPAN  Schwarz and Arnason (1996)   JSSAfCL  LinkBarker  Link and Barker (2005)   Pradel  Pradlambda  Pradel (1996)   Pradelg  Pradsen  Pradel (1996) 
The R package marked (Laake, Johnson and Conn 2013) also overlaps substantially with the nonspatial features of openCR. Its interface echoes RMark just as openCR echoes secr. marked has some fancy features for individual covariates and random effects, and promises fast processing of large datasets. marked 1.2.6 includes fulllikelihood JSSA (POPAN) models parameterized in terms of entry probabilities (type JSSAb)[^footnote15], but not the other JSSA options in Table 3.
[^footnote15]: dipper example in openCRexamples.pdf.
Borchers, D. L. and Efford, M. G. (2008) Spatially explicit maximum likelihood methods for capturerecapture studies. Biometrics 64, 377385.
Burnham, K. P., Anderson, D. R., White, G. C., Brownie, C. and Pollock, K. H. (1987) Design and Analysis Methods for Fish Survival Experiments Based on ReleaseRecapture. American Fisheries Society Monograph 5. Bethesda, Maryland, USA.
Chandler, R. B. and Clark, J. D. (2014) Spatially explicit integrated population models. Methods in Ecology and Evolution 5, 13511360.
Choquet, R., Lebreton, J.D., Gimenez, O., Reboulet, A.M. and Pradel, R. (2009) UCARE: Utilities for performing goodness of fit tests and manipulating CAptureREcapture data. Ecography 32, 10711074.
Clark, J. S., Silman, M., Kern, R., Macklin, E., and HilleRisLambers, J. (1999) Seed dispersal near and far: patterns across temperate and tropical forests. Ecology 80, 14751494.
Cooch, E. and White, G. (eds) (2019) Program MARK: A Gentle Introduction. 19th edition. Available online at http://www.phidot.org/software/mark/docs/book/.
Crosbie, S. F. and Manly, B. F. J. (1985) Parsimonious modelling of capturemarkrecapture studies. Biometrics 41, 385398.
Efford, M. G. (2004) Density estimation in livetrapping studies. Oikos 106, 598610.
Efford, M. G. (2022a) secr: Spatially explicit capturerecapture models. R package version 4.5.5. https://CRAN.Rproject.org/package=secr/
Efford, M. G. (2022b) Efficient discretization of movement kernels for spatiotemporal capturerecapture. Journal of Agricultural, Biological and Environmental Statistics In press. https://doi.org/10.1007/s13253022005034
Efford, M. G., Borchers D. L. and Byrom, A. E. (2009) Density estimation by spatially explicit capturerecapture: likelihoodbased methods. In: D. L. Thomson, E. G. Cooch, M. J. Conroy (eds) Modeling Demographic Processes in Marked Populations. Springer. Pp 255269.
Efford, M. G., Borchers D. L. and Mowat, G. (2013) Varying effort in capturerecapture studies. Methods in Ecology and Evolution 4, 629636.
Efford, M. G. and Fewster, R. M. (2013) Estimating population size by spatially explicit capturerecapture. Oikos 122, 918928.
Efford, M. G. and Schofield, M. R. (2020) A spatial openpopulation capturerecapture model. Biometrics 76, 392402.
Efford, M. G. and Schofield, M. R. (2022) A review of movement models in open population capturerecapture. Methods in Ecology and Evolution In press. https://doi.org/10.1111/2041210X.13947
Ergon, T. and Gardner, B. (2014) Separating mortality and emigration: modelling space use, dispersal and survival with robustdesign spatial capturerecapture data. Methods in Ecology and Evolution 5, 13271336.
Fletcher, D. J. (2012) Estimating overdispersion when fitting a generalized linear model to sparse data. Biometrika 99, 230237.
Gardner, B. J., Reppucci, J., Lucherini, M. and Royle, J. A. (2010) Spatiallyexplicit inference for open populations: estimating demographic parameters from cameratrap studies. Ecology 91, 33763383.
Gimenez, O., Viallefont, A., Catchpole, E. A., Choquet, R. and Morgan, B. J. T. (2004) Methods for investigating parameter redundancy. Animal Biodiversity and Conservation 27, 561572.
Gimenez, O., Lebreton, J.D., Choquet, R. and Pradel, R. (2018) R2ucare: An R package to perform goodnessoffit tests for capturerecapture models. Methods in Ecology and Evolution 9, 17491754.
Gimenez, O., Lebreton, J.D., Choquet, R. and Pradel, R. (2022) R2ucare: GoodnessofFit Tests for CaptureRecapture Models. R package version 1.0.2. https://github.com/oliviergimenez/R2ucare/
Glennie, R., Borchers, D. L., Murchie, M. Harmsen, B. J. and Foster, R. J. (2019) Open population maximum likelihood spatial capturerecapture. Biometrics 75, 13451355.
Hines, J. E., Kendall, W. L. and Nichols, J. D. (2003) On the use of the robust design with transient capturerecapture models. The Auk 120, 11511158.
Huggins, R. M. (1989) On the statistical analysis of capture experiments. Biometrika 76, 133140.
Laake, J.L., Johnson, D. S. and Conn, P.B. (2013) marked: An R package for maximumlikelihood and MCMC analysis of capturerecapture data. Methods in Ecology and Evolution 4, 885890.
Laake, J. and Rexstad E. (2014) Appendix C. RMark  an alternative approach to building linear models in MARK. In: Cooch, E. and White, G. (eds) Program MARK: A Gentle Introduction. 13th edition. http://www.phidot.org/software/mark/docs/book/.
Lebreton, J.D., Burnham, K. P., Clobert, J., and Anderson, D. R. (1992) Modeling survival and testing biological hypotheses using marked animals: a unified approach with case studies. Ecological Monographs 62, 67118.
Lele, S.R., Nadeem, K. and Schmuland, B. (2010) Estimability and likelihood inference for generalized linear mixed models using data cloning. Journal of the American Statistical Association 105, 16171625.
Link, W. A. and Barker, R. J. (2005) Modeling association among demographic parameters in analysis of openpopulation capturerecapture data. Biometrics 61, 4654.
Link, W. A. and Barker, R. J. (2010) Bayesian Inference with Ecological Applications. Academic Press, Amsterdam.
Marzolin, G. (1988) Polygynie du CineIe plongeur (Cinclus cinclus) dans les cotes de Lorraine. L'Oiseau et la Revue Francaise d'Ornithologie 58:277286.
Nathan , R., Klein, E., RobledoArnuncio, J. J. and Revilla, E. (2012) Dispersal kernels: a review. In: J Clobert et al. Dispersal Ecology and Evolution. Oxford University Press. Pp 187210.
Nichols, J. D., Pollock, K. H. and Hines, J. E (1984) The use of a robust capturerecapture design in small mammal population studies: a field example with Microtus pennsylvanicus. Acta Theriologica 29, 357365.
Nichols, J. D. (2016) And the first one now will later be last: timereversal in CormackJollySeber models. Statistical Science 31, 175190.
Pledger, S., Pollock, K. H. and Norris, J. L. (2003) Open capturerecapture models with heterogeneity: I. CormackJollySeber model. Biometrics 59, 786794.
Pledger, S., Pollock, K. H. and Norris, J. L. (2010) Open capturerecapture models with heterogeneity: II. JollySeber model. Biometrics 66, 883890.
Pollock, K. H. (1982) A capturerecapture design robust to unequal probability of capture. Journal of Wildlife Management 46, 752757.
Pradel, R. (1996) Utilization of capturemarkrecapture for the study of recruitment and population growth rate. Biometrics 52, 703709.
Pradel, R., Hines, J. E., Lebreton, J.D. and Nichols, J. D. (1997) Capturerecapture survival models taking account of transients. Biometrics 53, 6072.
Royle, J. A., Chandler, R. B., Sollmann, R. and Gardner, B. (2014) Spatial capturerecapture*. Academic Press.
Schwarz, C. J. (2001) The JollySeber model: more than just abundance. Journal of Agricultural, Biological, and Environmental Statistics 6, 195205.
Schwarz, C. J. and Arnason, A. N. (1996) A general methodology for the analysis of capturerecapture experiments in open populations. Biometrics 52, 860873.
Schofield, M. and Barker, R. (2016) 50yearold curiosities: ancillarity and inference in capturerecapture models. Statistical Science 31, 161174.
Seber, G. A. F. (1982) The estimation of animal abundance and related parameters. 2nd edition. Griffin.
White, G. C. and Burnham, K. P. (1999) Program MARK: survival estimation from populations of marked animals. Bird Study 46, Supplement S120S139.
Whittington, J. and Sawaya, M. A. (2015) A comparison of grizzly bear demographic parameters estimated from nonspatial and spatial open population capturerecapture models. PLoS ONE 10, e0134446. https://doi.org/10.1371/journal.pone.0134446
Williams, B. K., Nichols, J. D. and Conroy, M. J. (2002) Analysis and management of animal populations. Academic Press, San Diego.
\pagebreak
Code used to generate schematic diagrams of data structure.
onemulti < function(st = c(0,6,11,15), le = c(5,4,3,5), yb = 7, col=col1, outer = TRUE) { col < rep(col, le) xl < unlist(mapply(":",st,le+st1)) yb < rep(yb,length(xl)) xr < xl + width yt < yb + height rect(xl,yb,xr,yt,col=col) text(xl+width/2, yb+height/2, unlist(mapply(":", 1, le))) xl < st  margin yb < rep(yb[1], length(xl))  margin xr < st+le1+width+margin yt < yb+height+2*margin rect(xl,yb,xr,yt) text(st+le/2, rep(yb[1]+2*margin,length(st))+height+0.5, paste('session',1:length(st))) if (outer) { rect(st[1]3*margin, yb[1]2*margin, tail(st+le1,1)+width+3*margin, yb[1]+height+8*margin) } } onejoined < function(offset = 1.5, le = c(5,4,3,5), yb = 2.2, col=col1, intervals = TRUE, intlabel = 'intervals', leftlabel = '', outer = TRUE) { col < rep(col, le) xl < 0:(sum(le)1)+offset yb < rep(yb,length(xl)) xr < xl + width yt < yb + height rect(xl,yb,xr,yt,col=col) text(xl+width/2, yb+height/2, c(1:length(xl))) if (intervals) { xi < offset + (1:(length(xl)1))  (1width)/2 xip < cumsum(le)[length(le)] # intermediate between primary sessions intervals < rep(0,length(xi)) intervals[xip] < 1 text(xi, yb [1]0.8, intervals) text(0.2, yb[1]0.8, intlabel) segments(xi[xip], rep(yb[1]0.4,length(xip)), xi[xip], rep(yb[1]+0.4, length(xip))+height) } text (0.4, yb[1]+height/2, leftlabel, adj = c(1,0.5)) if (outer) { rect(offset2*margin, yb[1]2*margin, sum(le)1+offset+width+2*margin, yb[1]+height+2*margin) } }
# Fig. 1 Singlestratum data par(cex=1, xpd = TRUE, mfrow = c(1,1), mar=c(1,4,1,4)) width < 0.85 height < 1.1 margin < 0.15 col1 < c('salmon','pink','brown', 'red') col2 < c('green','lightgreen','darkgreen', 'lightblue') MASS::eqscplot(0,0,xlim=c(0,20), ylim=c(0,8), type='n', axes=F,xlab='',ylab='') onemulti(col = col1) text(9, 5.2, 'join()', cex=1.1) arrows (10.7,6.2,10.7,4.2) onejoined(leftlabel='')
# Fig. 2 Multistratum data par(cex = 0.9, xpd = TRUE, mfrow = c(1,1), mar = c(1,4,1,4)) MASS::eqscplot(0,0,xlim=c(3,20), ylim=c(2,8), type='n', axes=FALSE, xlab = '',ylab='') onejoined(leftlabel='stratum 1', yb = 6.5, intlabel='') onejoined(leftlabel='stratum 2', yb = 3, intlabel='') onejoined(leftlabel='stratum 3', yb = 0.5, le = c(4,3,4,4), intlabel='', col = col2) rect(3, 2, 19.3, 8.7)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.