Model construction

The manner in which models are constructed in the marked package is different than in MARK. We start by describing how MARK specifies and constructs models and how the marked package differs. The likelihood for CJS/JS is a product multinomial and the multinomial cell probabilities are functions of the parameters. We will focus on CJS which has two types of parameters: \(p\) (capture/recapture probability) and \(\phi\) (apparent survival). Each cell is associated with a release cohort and recapture time (occasion). The cohorts can be further split into groups based on categorical variables (e.g., sex). In MARK, a parameter index matrix (PIM) is used to specify the parameters. Each cell can have a unique index for each type of parameter type (e.g., \(p\) and \(\phi\)). As a very simple example, we will consider a CJS model for a single group with k=4 occasions. For this structure the PIMS are triangular with a row for each cohort released at each occasion (time) and the column representing the occasions (times) following the release. Using an unique index for each potential parameter, the PIMS would be:

The index 5 represents survival of the second release cohort between occasions 3 and 4 and index 8 represents recapture probability of the first release cohort on the third occasion. A separate set of PIMS is created for each group defined in the data. For this example, if there were \(g\) groups there would be \(g\times 12\) possible unique indices.

Some reduced models with constraints on parameters can be constructed by modifying the PIMS. For example, a model with constant survival and capture probability (Phi(.)p(.)) could be specified by setting all of the PIM values to 1 for the \(\phi\) PIM and 2 for the \(p\) PIM. Not all reduced models can be specified by modifying the PIMS and more generally reduced models are constructed with a design matrix which is a set of linear constraints. Each parameter index specifies a row in the design matrix which can be used to apply constraints on the parameters and to relate covariates to the parameters. The PIM approach minimizes the work of manually creating a design matrix by reducing the set of potential parameters and thus the number of rows in the design matrix. For our example, there would be 12 rows in the design matrix and each column in the design matrix (\(X\)) represents one of the parameters in the vector\(\beta\). If a logit link is used, the real parameters \(\theta\) (e.g., \(\phi\) and \(p\) in CJS) are:

\[\begin{equation} \theta=\frac{1}{1+\exp(-X\beta)} \tag{1} \end{equation}\]

To specify the Phi(.)p(.) model, the design matrix (\(X\)) would have 12 rows and 2 columns. The first column in \(X\) would have a 1 in the first six rows (indices 1-6) and a 0 in the last six rows (indices 7-12). The second column would have 0 in the first six rows and a 1 in the last 6 rows. The resulting model would have two parameters with the first representing \(\phi\) and the second \(p\).

To automate the creation of design matrices with formula, the RMark package creates “design data” which are data about the parameters like cohort, time, group, and any related variables. However, having a single \(X\) becomes problematic with individual covariates (variables that differ for animals in the same group/cohort). In MARK, these individual covariates are entered into the design matrix as a covariate name. When the real parameters are computed, the name of the covariate is replaced with the value for an animal. If there are only a few individual covariates and small number of animals, the time required for replacement and computation is minor, but when the individual covariates differ in time, there is a different covariate name for each time and most rows in the design matrix must be recomputed which results in longer execution times for large models.

Even if computation time is not limiting, we agree with McDonald et al. (2005) that the PIMs used in MARK can be quite confusing to the novice user because they are an additional layer of abstraction from the data. PIMS are useful for manual model creation, but they become an unnecessary nuisance when models and design matrices are created with formula. Explicit PIMS can be avoided by having an underlying real parameter for each animal for each occasion regardless of the release cohort. It can be viewed as a rectangular matrix with a column for each animal and a row for each occasion. That was the solution of McDonald et al. (2005) who in their mra package use a rectangular covariate matrices as predictor for the real parameters as in standard regression.

We use a similar approach that we believe is simpler for the user. For a specified model structure (e.g., CJS), the marked package creates from the raw data, a list of dataframes. Each dataframe contains the covariate data for a type of parameter in the model (e.g., \(\phi\) and \(p\)). The dataframes contain a record for each animal-occasion and the columns are the covariates. Each record contains the covariate values for that animal-occasion. If the covariate is constant across time then the value is the same value in each record for an animal and the value may be different for each occasion if the covariate is time-varying. Time-varying covariates must be named in the raw data in a certain manner and specified in the processing step. Each column in the dataframe is a vector that is equivalent to one of the covariate matrices in mra; however, rather than specifying the model as a set of covariate matrices, the standard R model formula (e.g., ~time+sex) can be used with the dataframe for each parameter which is even closer to a typical regression analysis.

A separate dataframe is created for each parameter to allow different data (e.g., time values for \(\phi\) and \(p\)) and number of occasions (e.g., \(\phi\) and \(p\) in CJS) for each parameter. For \(\phi\) and \(p\) in the CJS model, there are \(n(k-1)\) rows for \(n\) animals and \(k-1\) occasions for each parameter but time is labeled by 1 to \(k-1\) for \(\phi\) and 2 to \(k\) for \(p\). For the JS model, there are \(n(k-1)\) records for survival probability and entry probability, but for capture probability there are \(nk\) records because the initial capture event is modeled.

The parameter-specific animal-occasion dataframes are automatically created from the user’s data which contain a single record per animal containing the capture history and any covariates. Some covariates like cohort, age, and time are generated automatically for each record. Other static and time-varying covariates specified in the data are added. Static variables (e.g., sex) are repeated for each occasion for each animal. Time-varying covariates are specified in the data using a naming convention of vt where v is the covariate name and t is a numeric value for the time (occasion) (e.g., td19 contains the value of covariate td at time 19). The time-varying covariates are collapsed in the animal-occasion dataframes to a single column with the name identified by v and each record contains the appropriate value for the occasion. All of this is transparent to the user who only needs to specify a dataframe, the type of model (e.g., CJS), the variables that should be treated as time-varying, and the formula for each parameter.

With the R function model.matrix, the formula is applied to the dataframe to create the design matrix for each parameter (e.g, \(X_{\phi}\),\(X_{p}\)). They are each equivalent to a portion of the design matrix in MARK for an animal which are then “stacked on top of one another” to make a matrix for all animals. For maximum likelihood estimation (MLE), equation (1) (or similar inverse link function) is applied and the resulting vector is converted to a matrix with \(n\) rows and a column for each required occasion (\(k-1\) or \(k\)). For Bayesian MCMC inference in marked, only\(X\beta\) are needed for updating.

Model fitting

Currently there are only three types of models implemented in the marked package: 1) CJS, 2) JS, and 3) probitCJS. The first two are based on MLE and the third is an MCMC implementation. The likelihoods for CJS and JS were developed hierarchically as described by Pledger, Pollock, and Norris (2003) for CJS. For simplicity we only consider a single animal to avoid the additional subscript. Let \(\omega\) be a capture history vector having value \(\omega_{j}\)= 1 when the animal was initially captured and released or recaptured at occasion \(j\) and \(\omega_{j}\)= 0 if the animal was not captured on occasion \(j\), for occasions \(j=1,\dots,k\). The probability of observing a particular capture history \(Pr(\omega)\) for CJS can be divided into two pieces (\(Pr(\omega)=Pr(\omega_{1})Pr(\omega_{2})\): 1) \(\omega_{1}\) is the portion of the capture history from the initial release (i.e., first 1) to the last occasion it was sighted (i.e., last 1), and 2) \(\omega_{2}\) is from the last 1 to the last occasion. \(Pr(\omega_{1})\) is easily computed because the time period when the animal was available for capture is known but \(Pr(\omega_{2})\) is more difficult to compute because it is unknown whether the animal survived until the last occasion, or if it died, when it died. Typically, \(Pr(\omega_{2})\) has been computed recursively (Nichols 2005). A hierarchical construction is more direct and understandable. As described in Pledger, Pollock, and Norris (2003), we let \(f\) be the occasion an animal was released and let \(d\) be the occasion after which the animal is no longer available to be recaptured due to death or termination of the study. Then \(Pr(\omega|\, f)=\sum_{d}Pr(\omega|\, f,\, d)Pr(d\,|\, f)\) where the conditional probability \(Pr(\omega|\, f,\, d)\) is: \[ Pr(\omega|\, f,\, d)=\prod_{j=f+1}^{d}p_{j}^{\omega_{j}}(1-p_{j})^{(1-\omega_{j})} \] and the departure probability is \[ Pr(d\,|\, f)=\left(\prod_{j=f}^{d-1}\phi_{j}\right)(1-\phi_{d}) \] Viewed in this way, it is possible to construct the likelihood values for all of the observations with a set of matrices as we describe in the Appendix.

The hierarchical approach is easily extended to JS. Now the capture history has an additional component with \(Pr(\omega)=Pr(\omega_{0})Pr(\omega_{1})Pr(\omega_{2})\) where \(\omega_{0}\) is the portion of the capture history vector from the first occasion (\(j\)=1) to the occasion the animal was first seen (\(j=f\)). The \(Pr(\omega_{0})\) is similar to \(Pr(\omega_{2})\) but now it is \(e\), the occasion at which the animal was first available for capture (i.e., entered in the prior interval) that is unknown. The CJS likelihood provides \(Pr(\omega|\, f)\) so we only need to compute \(Pr(\omega_{0})=\sum_{e=0}^{f-1}Pr(\omega_{0}|\, e)\pi_{e}\) where \(\pi_{e}\) is the probability of an animal enters the population in the interval between occasion \(e\) and \(e+1\) (\(\sum_{e=1}^{k-1}\pi_{e}=1-\pi_{0}\); pent parameters in MARK and specified as \(\beta\) by ) and \[ Pr(\omega_{0}|\, e)=\left(\prod_{j=e+1}^{f-1}(1-p_{j})\right)p_{f} \]

The JS likelihood also contains a component for animals that entered but were never seen. The details for the JS likelihood construction are provided in the Appendix.

The MLEs are obtained numerically by finding the minimum of the negative log-likelihood using optimization methods provided through the optimx R package (Nash and ). The default optimization method is BFGS but you can specify alternate methods and several methods used independently or in sequence as described by Nash and (). Initial values for parameters can either be provided as a constant (e.g., 0) or as a vector from the results of a previously fitted similar model. If initial values are not specified then they are computed using general linear models (GLM) that provide approximations. Using the underlying idea in Manly and Parr (1968) we compute initial estimates for capture-probability \(p\) for occasions 2 to \(k-1\) using a binomial GLM with the formula for \(p\) which is fitted to a sequence of Bernoulli random variables that are a subset of the capture history values \(y_{ij}\) \(i=1,...,n\) and \(j=f_{i}+1,...,l_{i}-1\) where \(f_{i}\) and \(l_{i}\) are the first and last occasions the \(i^{th}\) animal was seen. A similar but more ad-hoc idea is used for \(\phi\). We know the animal is alive between \(f_{i}\) and \(l_{i}\) and assume that the animal dies at occasion \(l_{i}+1\leq k\). We use a binomial GLM with the formula for \(\phi\) fitted to the \(y_{ij}^{*}\) \(i=1,...,n\) and \(j=f_{i}+1,...,l_{i}+1\leq k\) where \(y_{ij}^{*}\)=1 for \(j=f_{i}+1,...,l_{i}\) and \(y_{ij}^{*}\)=0 for \(j=l_{i}+1\leq k\). For \(\phi\) and \(p\), a logit link is used for MLE and a probit link for MCMC. For JS, a log link is used for \(f^{0}\), the number never captured and the estimate of super-population size is the total number of individual animals plus \(f^{0}\). Also, for JS a multinomial logit link is used for \(\pi_{e}\). The initial value is set to 0 for any \(\beta\) that is not estimated (e.g., \(p_{K}\) in CJS or \(p_{1}\),\(\pi_{e}\), \(N\) for JS).

Likelihood details

Here we provide some details for the likelihoods that are computed in cjs.lnl and cjs.f for the CJS model and in js.lnl for the JS model. Let \(\phi\),\(M\) \(p\) be \(n\times k\) matrices which are functions of the parameters where\(\phi_{ij}\) is the survival probability for animal \(i\) from occasion \(j-1\) to \(j\) (\(\phi_{i1}=1\) ), \(m_{ij}\) is the probability of dying in the interval \(j\) to \(j+1\) (\(m_{ik}=1)\) and \(p_{ij}\) is the capture probability for animal \(i\) on occasion \(j\) (\(p_{i1}=0\)). In addition we define a series of matrices and vectors computed from the data in the function process.ch. Let \(C\) be an \(n\times k\) matrix of the capture history values and \(f_{i}\)and \(l_{i}\) are the first and last occasions the \(i^{th}\) animal was seen. Derived from those values are \(n\times k\) matrices \(L\), \(F\) and \(F^{+}\) which contain 0 except that \(L_{ij}=1\; j\geq l_{i}\),\(F_{ij}=1\; j\geq f_{i}\) and \(F_{ij}^{+}=1\; j>f_{i}\). The likelihood calculation has the following steps where the matrix multiplication is element-wise and \(1\) and represents an \(n\times k\) matrix where each element is 1:

1. Construct the \(n\times k\) matrix \(\phi'=1-F^{+}+\phi F^{+}\)which is the \(\phi\) matrix modified such that \(\phi_{ij}=1\) for \(j\leq f_{i}\).
2. Construct the \(n\times k\) matrix \(\phi*\)from \(\phi'\) where \(\phi_{ij}*=\prod_{i=1}^{j}\phi_{ij}'\).
3. Construct the \(n\times k\) matrix \(p'=1-F^{+}+F^{+}(Cp+(1-C)(1-p))\).
4. Construct the \(n\times k\) matrix \(p*\)from \(p'\) where \(p_{ij}*=\prod_{i=1}^{j}p_{ij}'\).
5. Compute \(n\times1\) vector of probabilities of observed capture histories \(Pr(\omega_{i})\; i=1,...,n\) which are the sums of the rows of \(LM\phi^{*}p^{*}\).
6. Compute the log-likelihood \(\sum_{i=1}^{n}\ln(Pr(\omega_{i}))\)

The R code translates from the mathematics quite literally as:

 Phiprime = 1-Fplus + Phi*Fplus
 Phistar = t(apply(Phiprime,1,cumprod))
 pprime = (1-Fplus)+Fplus*(C*p+(1-C)*(1-p))
 pstar = t(apply(pprime,1,cumprod))
 pomega = rowSums(L*M*Phistar*pstar)
 lnl = sum(log(pomega))

While this code is simple it is faster if the apply is replaced with a loop because there are far more rows than occasions or replaced with compiled code which we did with a FORTRAN subroutine (cjs.f called from cjs.lnl).

The Jolly-Seber likelihood can be partitioned into 3 components: 1) CJS likelihood for \(Pr(\omega_{1})Pr(\omega_{2})\) treating the first ``1’’ as a release, 2) a likelihood component for \(Pr(\omega_{0})\); entry and first observation, 3) a component for those that entered before each occasion but were never seen. We define an additional \(n\times k\) matrix \(\pi\) which are the entry probabilities into the population (specified as \(\beta\) by Schwarz and Arnason (1996)) with the obvious constraint that \(\sum_{j=0}^{k-1}\pi_{ij}=1\) and \(N_{g}\: g=1,...,G\) is the abundance of animals in each of the defined \(G\) groups (e.g., male/female) that were in the population at some time (super-population size). \(N_{g}=n_{g}+f_{g}^{0}\) where \(n_{g}\) is the number observed in the group and \(f_{g}^{0}\) are the estimated number of animals in the group that were never seen.

The same hierarchical approach used for CJS can be used for the second component. Construct the capture history probability for a given entry time and then sum over all possible entry times. . The second component is constructed with the following steps:

1. Construct the \(n\times k\) matrix \(E=(1-p)\phi(1-F)+F\),
2. Construct the \(n\times k\) matrix \(E*\)where \(E_{ij}*=\prod_{l=j}^{k}E_{il}\),
3. Compute \(n\times1\) vector of probabilities \(Pr(\omega')\) which are the sums of the rows of \(E*(1-F^{+})\pi\) multiplied by the vector \(p_{if_{i}}\: i=1,...,n\),
4. Compute the log-likelihood \(\sum_{i=1}^{n}\ln(Pr(\omega_{i}'))\).

For the third likelihood component for missed animals, we constructed \(G\times k\) dummy capture histories of all 0’s except for a ``1’’ at the occasion the animals entered. From the CJS portion of the code, we obtained \(p_{gj}^{0}\) the probability that an animal released in group \(g\) on occasion \(j\) would never be captured. The final log-likelihood component is: \[ \sum_{g=1}^{G}f_{g}^{0}\ln\left[\sum_{j=1}^{k}\pi_{g,j-1}(1-p_{gj})p_{gj}^{0}\right]+\ln(N_{g}!)-\ln(f_{g}^{0}!) \]

The JS log-likelihood is the total of the 3 components plus \(-\sum_{g=1}^{G}\sum_{j=1}^{k}\ln(n_{gj}!)\) which do not depend on the parameters but is added after optimization in function js, to be consistent with the output from MARK.

With the structure we have used, the design matrices can become quite large and available memory could become limiting. For the design matrices, we use sparse matrices with the R package Matrix (citation). In addition, for MLE analysis, we use the following to reduce the required memory:

1. We reduce the data to \(n^{*}\leq n\) individuals by aggregating records with identical data and using the frequencies (\(f_{1},...,f_{n^{*}};n=\sum_{i=1\:}^{n^{*}}f_{i}\)) in the likelihood.
2. We construct the design matrix \(X\) incrementally with a user-specified size of data chunk that is processed at one time.
3. We retain only the rows of \(X\) which are unique for the design and including any fixed parameters which can be animal-specific.

For Bayesian MCMC inference, we cannot aggregate records but we reduce the required memory by:

1. Eliminating the unused animal-occasion data for occasions prior to and including the release occasion for the animal, and
2. Storing only the unique values of the real parameters which are unique rows of \(X_{\phi}\)and \(X_{p}\).

Appendix References