Model types

openCR fits nonspatial open-population models of the Cormack-Jolly-Seber (CJS) and Jolly-Seber-Schwarz-Arnason (JSSA² 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. Conditional-likelihood JSSA models are also called Pradel–Link–Barker (PLB) models. Pradel analyses are also provided.

Spatial versions of the CJS and JSSA model types are also provided³. 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.

Data

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

Model specification and fitting

Models are specified using formula notation as in secr. Possible predictors include both pre-defined variables for learned responses, trend over time, etc., and user-provided 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⁴. Within-session variation in detection parameters may also be modelled (see field vole example in openCR-examples.pdf).

Parameterization

A selection of parameterizations is offered for recruitment in JSSA models. Models can also be parameterized in terms of the time-specific population size (non-spatial models) or density (spatial models), avoiding the super-population parameter.

Super-population 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 time-specific population sizes and densities.

Features and limitations

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.

CJS vs JS

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

Parameterization of recruitment in JSSA models

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 time-specific probability β_j (PENT in MARK), an idea from Crosbie and Manly (1985). Other parameterizations are

• number of new entrants at each time j
• per capita fecundity (new entrants at time j scaled by 1/number in population at j − 1)
• seniority (reverse-time survival Pradel 1996, Nichols 2016)
• population growth rate λ
• (relative) number in population at each time j

Estimates of recruitment or implied recruitment from any one of these six parameterizations can be used to infer the others⁵. 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).

Schwarz (2001) is illuminating (see also chapter on Jolly-Seber 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 Jolly-Seber models’ (‘Recruitment Parameters’ in the help index).

Conditional (PLB) vs full likelihood JSSA

For each JSSA recruitment parameterization there is a choice between models that include the total number of detected individuals (u_⋅ or n in different notations), and models that condition on this number. Conditional-likelihood models do not directly estimate abundance; abundance is estimated as a derived parameter (Schwarz and Arnason 1996). Full-likelihood models include abundance as a parameter. The choice of formulation has virtually no effect on the parameter estimates⁶. 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 Pradel–Link–Barker PLB was suggested by Efford and Schofield (2020).

Sufficient statistics vs capture histories

Historically the CJS and JS likelihoods have been expressed in terms of ‘sufficient statistics’ that are time-specific 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⁷. 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.

Robust design

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 capture-history 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.

Spatial vs nonspatial

Models may be spatially explicit or not. Nonspatial models ignore the spatial distribution of animals. Spatial models use the spatially explicit capture–recapture 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

• allowance for varying extent of sampling area
• modelling of individual heterogeneity due to differential access to detectors
• separation of emigration and mortality

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 non-spatial models, replacing ‘number’ by ‘density’. The locations at which animals recruit are not modelled.

Home-range shifts between primary sessions

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:

• Movement that is ignored inflates estimates of the within-session scale of detection σ, with flow-on effects on demographic parameters.
• If the distribution of dispersal distances can be inferred from the detection histories of residents then it is possible in principle to separate actual mortality from losses due emigration (Ergon and Gardner 2014). However, the robustness and data requirements of movement models have yet to be fully understood.

Stratification

From openCR 2.0 onwards any model may be stratified. For stratified models (stratified = TRUE in openCR.fit) each session of a multi-session 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 ‘single-session’ capthist. The function stratify helps you construct stratified capthist objects from collections of single-session objects.

Fig. 2. Structure of data for stratified open-population analysis in openCR. Each stratum is a pre-joined (single-session) component of a multi-session capthist object. The internal structure and detector may differ between strata.

Stratified models may use ‘stratum’ as a factor-valued predictor. Groups of strata may be contrasted using stratum-level covariates as described later.

Non-spatial openCR models

Spatial openCR models

Table 4. Parameter definitions and default link functions (spatial models)

* parameters marked with an asterisk are scaled by the interval between primary sessions.

Table 5. Parameters of spatial openCR models

Spatial models with type ending in CL have features in common with the Pradel–Link–Barker models, hence the aliases as shown.

Built-in predictors

Table 6. Built-in predictors (‘sessions’ refers to primary sessions)

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 non-spatial 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.

IMPORTANT NOTE: Learned response predictors (‘b’, ‘bsession’ etc.) were re-defined in openCR 1.3.0. Models fitted with earlier versions should be re-fitted.

User-provided covariates

The rules for covariates largely follow secr (secr-overview.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 built-in 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).

Closed populations

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]

##         estimate SE.estimate    lcl     ucl
## D          5.485     0.64703  4.356  6.9059
## lambda0    0.307     0.03413  0.247  0.3815
## sigma     28.764     1.30055 26.326 31.4283

pred_openCR

##         estimate SE.estimate    lcl     ucl
## D          5.485     0.64480  4.356  6.9059
## lambda0    0.307     0.03403  0.247  0.3815
## sigma     28.764     1.29988 26.326 31.4283

# compare timings in seconds
c(secr = fit_secr$proctime, openCR = fit_openCR$proctime)

##   secr.elapsed openCR.elapsed 
##          1.655          1.322

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

##   secr.logLik openCR.logLik 
##        -758.9        -758.9

Finite mixtures

Two- and three-class 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 [secr-finitemixtures.pdf].

Age

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.

Specify the details argument ‘agebreaks’ to group numeric ages into age classes. Breaks are used with the cut function to generate a factor from the numeric ages; the cut argument ‘right’ is set to FALSE to include the lower limit in each age class. Extreme ages are shrunk to the interval [‘minimumage’, ‘maximumage’] before grouping, so these arguments must be compatible with ‘agebreaks’ (e.g., ‘maximumage’ >= lower bound of oldest group). Check the grouping by applying it to the matrix of numeric ages in your data. For example,

agebrk <- c(0, 2, Inf)

# construct matrix of numeric ages (animal x secondary session)
age <- age.matrix(join(ovenCH), maximumage = 2, unborn = NA)

# tabulate the grouped ages
aclass <- cut(age, breaks = agebrk, right = FALSE)
table(age, aclass)

##    aclass
## age [0,2) [2,Inf)
##   0   680       0
##   1   620       0
##   2     0    1000

The notation [0,2) indicates ages in the interval 0 ≤ age < 2. The numbers in older groups include animals never seen again and possibly dead.

This example illustrates how to use grouped ages in a model. The data are an undocumented non-spatial selection from the same brushtail possum study as OVpossumCH (trapping in February, June and September, 1980–1988, known-age females only; interval in years). The individual covariate ‘age’ records the age of each possum at first capture, in years.

# retrieve capthist object
datadir <- system.file('extdata', package = 'openCR')
CH <- readRDS(paste0(datadir,'/poss8088F.RDS'))

# model with grouped ages; any maximumage>=6 OK
fit <- openCR.fit(CH, model = list(phi ~ age), details = list(
  agebreaks = c(0,2,4,6,Inf), initialage = 'age', maximumage = 6))

# show results for first session only, as no time effect fitted
# levels of age grouping factor are stored in 'design' object
newdat <- data.frame(age = fit$design$agelevels)
predict(fit, newdata =  newdat)$phi

##   session     age estimate SE.estimate    lcl    ucl
## 1 Feb1980   [0,2)   0.5622     0.06095 0.4415 0.6760
## 2 Feb1980   [2,4)   0.9152     0.03530 0.8157 0.9634
## 3 Feb1980   [4,6)   0.8775     0.04529 0.7583 0.9424
## 4 Feb1980 [6,Inf)   0.8220     0.02079 0.7776 0.8592

If ‘agebreaks’ is omitted then the default uses bins of width 1 time unit from ‘minimumage’ up to, but not including, ‘maximumage’, with an extra bin for ‘maximumage’ and above (7 bins [0,1),…,[6,Inf) in this example).

The older ‘agecov’ mechanism is limited to sampling sessions separated by one time unit and is deprecated from 2.2.6.

For a quadratic relationship with age, specify an additive model with both Age and Age2 terms (e.g., model = phi ~ Age + Age2).

Sampling intervals

We have seen the role of the intervals attribute in defining primary and secondary sessions. Between-session 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 θ_j to the interval-specific value θ_j^′ uses θ_j^′ = θ_j^T_j where θ_j is one of ϕ_j or λ_j, and T_j is the duration of interval j.

Scaling γ follows the same pattern except that the relevant duration for γ_j is T_j − 1. Scaling per capita recruitment f_j is more tricky. We use f_j^′ = (ϕ_j + f_j)^T_j − ϕ_j^T_j.

Custom design data

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⁹. The secr function insertdim can help to expand data into the correct row order.

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 non-spatial 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)

##                 model npar rank logLik   AIC  AICc dAIC AICwt
## fitb        p~b phi~1    3    2 -254.6 515.2 515.6    0   0.5
## fitbc p~customb phi~1    3    2 -254.6 515.2 515.6    0   0.5

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 non-CJS models because these account for animals never detected. The corresponding likelihood component uses a distinct design matrix for a ‘naive’ animal. To customize non-CJS 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)

##                      model npar rank logLik  AIC AICc dAIC AICwt
## fitb2        p~b phi~1 f~1    4    4 -660.9 1330 1330    0   0.5
## fitbc2 p~customb phi~1 f~1    4    4 -660.9 1330 1330    0   0.5

Transience

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)

##                       model npar rank logLik   AIC  AICc  dAIC AICwt
## fitcov  p~1 phi~residentcov    3    2 -225.8 457.6 458.0  0.00     1
## fitnull           p~1 phi~1    2    2 -254.6 513.2 513.4 55.56     0
## fitdf      p~1 phi~resident    3    2 -254.6 515.2 515.6 57.56     0

pred <- predict(fits, newdata = data.frame(resident = TRUE, residentcov = TRUE))
do.call(rbind, lapply(pred, '[[', 'phi'))

##         session resident residentcov estimate SE.estimate    lcl    ucl
## fitnull    2005     TRUE        TRUE   0.4630     0.05473 0.3590 0.5703
## fitcov     2005     TRUE        TRUE   0.7387     0.08767 0.5372 0.8732
## fitdf      2005     TRUE        TRUE   0.4630     0.05473 0.3590 0.5703

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 coding

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¹⁰.

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 log-likelihood, but with different ‘beta’ parameters:

fit0 <- openCR.fit(ovenCH, model = p~t)
fitd <- openCR.fit(ovenCH, model = p ~ -1+t)
coef(fit0)

##          beta SE.beta     lcl     ucl
## p    -1.54953  0.2459 -2.0315 -1.0675
## p.t3  0.32963  0.3280 -0.3133  0.9725
## p.t4 -1.42728  0.5259 -2.4581 -0.3965
## p.t5 -0.14375  0.4489 -1.0236  0.7361
## phi  -0.03141  0.2399 -0.5016  0.4388

coef(fitd)

##          beta SE.beta     lcl     ucl
## p.t2 -1.54955  0.2459 -2.0316 -1.0675
## p.t3 -1.21990  0.2188 -1.6487 -0.7911
## p.t4 -2.97677  0.4663 -3.8907 -2.0628
## p.t5 -1.69325  0.3783 -2.4347 -0.9518
## phi  -0.03142  0.2399 -0.5016  0.4387

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 non-reference levels). The following fit is therefore equivalent to the preceding fitd.

fitd2 <- openCR.fit(ovenCH, model = p~t, details = list(dummyvariablecoding = 't'))
coef(fitd2)

##          beta SE.beta     lcl     ucl
## p.t2 -1.54955  0.2459 -2.0316 -1.0675
## p.t3 -1.21990  0.2188 -1.6487 -0.7911
## p.t4 -2.97677  0.4663 -3.8907 -2.0628
## p.t5 -1.69325  0.3783 -2.4347 -0.9518
## phi  -0.03142  0.2399 -0.5016  0.4387

Mean of a parameter across levels of a factor

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¹¹ by setting the contrast function for the factor to contr.sum in the details argument¹². With the resulting factor coding the first coefficient corresponds to the mean. Applying this to estimate the average time-specific 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')])

##       beta   lcl    ucl
## phi 0.5633 0.505 0.6199

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.

Movement kernels

Note: The online vignette openCR-kernel.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 built-in 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, θ) assuming direction θ uniform on (0, 2π) (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. Cell-specific 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 user-specified 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πr).

* Continuous, untruncated, kernel. Expected values for the discretized and truncated kernel will be less (see summary.kernel).

The ‘BVT’ kernel is the same as ‘2Dt’ of Clark et al. (1999) and Nathan et al. (2012). The parameter α (move.a) corresponds to a in Nathan et al. (2012) and $\sqrt u$ in Clark et al. (1999); the parameter β (move.b) corresponds to b − 1 in Nathan et al. (2012) and p in Clark et al. (1999). Defining move.b as β ≡ b − 1 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 ν = 2β.

The `BVT’ kernel approaches bivariate normal as β → ∞ and Cauchy as β → 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 β ≤ 0.5.

Zero-inflated kernels

openCR 2.2 introduces zero-inflated versions of kernels otherwise defined with one parameter or none. Zero-inflated kernels use the suffix ‘zi’, hence ‘BVNzi’, ‘BVEzi’, ‘RDEzi’, ‘UNIzi’. The kernel-free independent movement model ‘IND’ also has a zero-inflated form ‘INDzi’ that is not strictly independent or uncorrelated. Each of these models has an additional zero-inflation parameter (move.b or move.a depending on whether the base kernel does or does not already have a parameter). Zero-inflated 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.

User-defined kernel

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 ≤ 0 or b ≤ 0.

Sparse kernels

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 (N-S, E-W, NW-SE, NE-SW) (Efford 2022b). The number of cells then increases only linearly with radius. Cell-wise movement probabilities are adjusted so that the distribution of dispersal distances is almost unchanged (essentially multiplying by 2π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.

Edge effects

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:

1. 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'.

2. 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 (all-zero) history was estimated incorrectly.

Plotting and summary

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)

## Kernel radius (cells)     :  10 
## Spacing (side of cell)    :  10  (m)
## Number of cells           :  349 
## Movement model            :  BVN 
## Parameter(s)              :  move.a = 40 
## Proportion truncated      :  0.03189 
## Movement as truncated at edge of kernel
## Empirical mean distance   :  47.91  (m)
## Expected distance         :  47.19  (m)
## 50th percentile (median)  :  45.61  (m)
## 90th percentile           :  79.39  (m)
## Movement, untruncated kernel
## Expected distance         :  50.13  (m)
## 50th percentile (median)  :  47.1  (m)
## 90th percentile           :  85.84  (m)

Use the secr function spotHeight(k) to display cell values on the plot.

Warnings

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

2. Kernel-based movement models require extreme care. Definitive advice cannot yet be given on the safe use of these models. Long-distance movements will usually be poorly sampled and poorly modelled.

3. User-defined functions cannot be used with multithreaded C++, so they will be slow to fit; always set ncores = 1.

Nonidentifiability

It is common for some session-specific parameters of open capture–recapture 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 non-zero eigenvalues. Computed eigenvalues of non-identifiable 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.

Session-specific 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.

Failure of numerical maximization

Bad estimates (zero, very large, close to starting values or zero variance) may merely indicate a problem with the maximization algorithm rather than nonidentifiability.

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 = "Nelder-Mead", details = list(control = list(maxit = 5000)))

AIC(fitnm,fitnr)

##                    model npar rank logLik  AIC AICc  dAIC  AICwt
## fitnr p~1 phi~t lambda~t    9    9 -656.7 1331 1334 0.000 0.9898
## fitnm p~1 phi~t lambda~t    9    9 -661.3 1341 1344 9.154 0.0102

fitnm <- openCR.fit(ovenCH, type = 'JSSAlCL', model = list(phi ~ t, lambda~t),
                    method = "Nelder-Mead", details = list(control = list(maxit = 2000)),
                    start = fitnr)
AIC(fitnm,fitnr)

##                    model npar rank logLik  AIC AICc dAIC AICwt
## fitnm p~1 phi~t lambda~t    9    9 -656.7 1331 1334    0   0.5
## fitnr p~1 phi~t lambda~t    9    9 -656.7 1331 1334    0   0.5

fitnr2 <- openCR.fit(ovenCH, type = 'JSSAlCL', model = list(phi ~ t, lambda~t),
    details = list(stepmax = 2))

Speed

Spatial models are slow to fit. Consider these options

• Use no more mask points than necessary. Typically about 1000 will do (may not apply for kernel movement models).
• Data with many occasions (secondary sessions) should be collapsed.
• Use the conditional likelihood (PLB) models: estimates of phi and lambda are often all you need, and derive can give estimates of abundance (superN, N, superD, and D) from PLB models, as well as alternative measures of recruitment.
• Avoid individual covariates with many levels. This applies especially to continuous individual covariates: normally these should be discretized (coded as a few ordered categories, but not converted to factor).
• First fit the most simple model and then add complexity; use a simpler related model for ‘start’.
• ‘secr’ data with detector type ‘multi’ fit much faster than ‘proximity’ data; use this option if it makes sense (and even maybe when it doesn’t).
• For problems with many parameters, ‘cyclic fixing’ may be useful (Schwarz and Arnason 1996; Pledger et al. 2003).

openCR ≥ 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 ≥ 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 quad-core desktop with hyperthreading), or some lesser number if you want to multitask:

# RCPP_PARALLEL_NUM_THREADS 
# recommended for quad-core Windows PC
setNumThreads(7)

The grain size' tuning parameter (see [RcppParallel](https://rcppcore.github.io/RcppParallel/)) may be varied withdetails$grain`, but it seems to have little effect.

Sampling variance warning

Full models (not CL or Pradel) include superpopulation size N as a variable. The default in openCR for both non-spatial 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".

Example datasets

Several examples of analyses with openCR are given in the associated vignette openCR-examples.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 openCR-examples.pdf for references.

Table 10. Multi-session data objects in secr.

Testing assumptions

This is generally an undeveloped field for spatially explicit capture–recapture 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 Cormack-Jolly-Seber (nonspatial) models there is an established suite of tests following Burnham et al. (1987). The tests have been implemented in the U-CARE 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')

## Loading required namespace: R2ucare

## $Male
## $Male$components
##          stat df p_val sign_test
## test3sr 6.778  5 0.238    -1.530
## test3sm 0.000  2 1.000        NA
## test2ct 4.284  2 0.117    -1.035
## test2cl 0.000  0 1.000        NA
## 
## $Male$overall_CJS
##                          chi2 degree_of_freedom p_value
## Gof test for CJS model: 11.06                 9   0.271
## 
## 
## $Female
## $Female$components
##          stat df p_val sign_test
## test3sr 4.985  5 0.418     1.428
## test3sm 2.041  3 0.564        NA
## test2ct 3.250  4 0.517    -0.901
## test2cl 0.000  0 1.000        NA
## 
## $Female$overall_CJS
##                          chi2 degree_of_freedom p_value
## Gof test for CJS model: 10.28                12   0.592

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 openCR-examples.pdf).

Limitations of openCR

openCR does not do

1. Continuous random effects (consider finite mixtures as an alternative)
2. Parameter counting to adjust AIC
3. Overdispersion adjustment (chat, QAIC) or goodness-of-fit tests, except for ucare.cjs (above).
4. MCMC
5. Bootstrap confidence intervals
6. Temporary emigration parameterizations of non-spatial robust-design models
7. Age-specific survival curves (Weibull etc.)
8. Mark-resight
9. SE for derived parameters and estimates with mlogit link (to be fixed)

Parameter counting and overdispersion adjustment are probably the most critical omissions. See Cooch and White (2019) for detailed coverage in the context of MARK.

Differences from secr

Defaults for some arguments differ between openCR.fit and secr.fit. For openCR.fit –

1. trace = FALSE
2. By default the reported log likelihood and AIC do not include the multinomial constant (details$multinom = FALSE)
3. The default criterion for AIC() and modelAverage() is ‘AIC’, not ‘AICc’ as in secr.
4. The maximum number of iterations used by openCR.fit() to maximize the likelihood defaults to 300 rather than 100.

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 conditional-likelihood models. In openCR conditional-likelihood 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 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 re-detections of animals between strata. The total log-likelihood 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

1. Hybrid mixture models (hcov in secr)
2. Groups (use strata, or CL and individual covariates, or see marked)
3. Regression splines from mgcv
4. Density surfaces and other spatial density models
5. Post-hoc probability density of activity centres (fxi in secr)
6. Non-point detectors (polygon, polygonX etc. in secr; discretize instead)
7. ‘collate’ function (make.table may do the job)
8. Variable effort for nonspatial models (cf Efford, Borchers and Mowat 2013) (The ‘usage’ attribute of traps objects is applied in spatial openCR models).
9. Negative binomial counts (binomN<0)

Relationship to other software

The non-spatial 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 non-spatial openCR models to MARK model types

The R package marked (Laake, Johnson and Conn 2013) also overlaps substantially with the non-spatial 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 full-likelihood JSSA (POPAN) models parameterized in terms of entry probabilities (type JSSAb)¹⁶, but not the other JSSA options in Table 3.