7.3 Mark-recapture models
A widely applied field method in ecology is to capture (or sight) animals, mark them (e.g., by tagging), then release them. This process is then repeated one or more times, and is often done for populations on an ongoing basis. The resulting data may be used to estimate population size.
The first subsection describes a simple mark-recapture model that does not involve any latent discrete parameters. The following subsections describes the Cormack-Jolly-Seber model, which involves latent discrete parameters for animal death.
Simple mark-recapture model
In the simplest case, a one-stage mark-recapture study produces the following data
- \(M\) : number of animals marked in first capture,
- \(C\) : number animals in second capture, and
- \(R\) : number of marked animals in second capture.
The estimand of interest is
- \(N\) : number of animals in the population.
Despite the notation, the model will take \(N\) to be a continuous parameter; just because the population must be finite doesn’t mean the parameter representing it must be. The parameter will be used to produce a real-valued estimate of the population size.
The Lincoln-Petersen (Lincoln 1930; Petersen 1896) method for estimating population size is \[ \hat{N} = \frac{M C}{R}. \]
This population estimate would arise from a probabilistic model in which the number of recaptured animals is distributed binomially, \[ R \sim \textsf{binomial}(C, M / N) \] given the total number of animals captured in the second round (\(C\)) with a recapture probability of \(M/N\), the fraction of the total population \(N\) marked in the first round.
data {
int<lower=0> M;
int<lower=0> C;
int<lower=0, upper=min(M, C)> R;
}
parameters {
real<lower=(C - R + M)> N;
}
model {
R ~ binomial(C, M / N);
}A probabilistic formulation of the Lincoln-Petersen estimator for population size based on data from a one-step mark-recapture study. The lower bound on \(N\) is necessary to efficiently eliminate impossible values.
The probabilistic variant of the Lincoln-Petersen estimator can be directly coded in Stan as shown in the Lincon-Petersen model figure. The Lincoln-Petersen estimate is the maximum likelihood estimate (MLE) for this model.
To ensure the MLE is the Lincoln-Petersen estimate, an improper uniform prior for \(N\) is used; this could (and should) be replaced with a more informative prior if possible, based on knowledge of the population under study.
The one tricky part of the model is the lower bound \(C - R + M\) placed on the population size \(N\). Values below this bound are impossible because it is otherwise not possible to draw \(R\) samples out of the \(C\) animals recaptured. Implementing this lower bound is necessary to ensure sampling and optimization can be carried out in an unconstrained manner with unbounded support for parameters on the transformed (unconstrained) space. The lower bound in the declaration for \(C\) implies a variable transform \(f : (C-R+M,\infty) \rightarrow (-\infty,+\infty)\) defined by \(f(N) = \log(N - (C - R + M))\); the reference manual contains full details of all constrained parameter transforms.
Cormack-Jolly-Seber with discrete parameter
The Cormack-Jolly-Seber (CJS) model (Cormack 1964; Jolly 1965; Seber 1965) is an open-population model in which the population may change over time due to death; the presentation here draws heavily on Schofield (2007).
The basic data are
- \(I\): number of individuals,
- \(T\): number of capture periods, and
- \(y_{i,t}\): Boolean indicating if individual \(i\) was captured at time \(t\).
Each individual is assumed to have been captured at least once because an individual only contributes information conditionally after they have been captured the first time.
There are two Bernoulli parameters in the model,
- \(\phi_t\) : probability that animal alive at time \(t\) survives until \(t + 1\) and
- \(p_t\) : probability that animal alive at time \(t\) is captured at time \(t\).
These parameters will both be given uniform priors, but information should be used to tighten these priors in practice.
The CJS model also employs a latent discrete parameter \(z_{i,t}\) indicating for each individual \(i\) whether it is alive at time \(t\), distributed as \[ z_{i,t} \sim \mathsf{Bernoulli}(z_{i,t-1} \; ? \; 0 \: : \: \phi_{t-1}). \]
The conditional prevents the model positing zombies; once an animal is dead, it stays dead. The data distribution is then simple to express conditional on \(z\) as \[ y_{i,t} \sim \mathsf{Bernoulli}(z_{i,t} \; ? \; 0 \: : \: p_t). \]
The conditional enforces the constraint that dead animals cannot be captured.
Collective Cormack-Jolly-Seber model
This subsection presents an implementation of the model in terms of counts for different history profiles for individuals over three capture times. It assumes exchangeability of the animals in that each is assigned the same capture and survival probabilities.
In order to ease the marginalization of the latent discrete parameter \(z_{i,t}\), the Stan models rely on a derived quantity \(\chi_t\) for the probability that an individual is never captured again if it is alive at time \(t\) (if it is dead, the recapture probability is zero). this quantity is defined recursively by \[ \chi_t = \begin{cases} 1 & \quad\text{if } t = T \\ (1 - \phi_t) + \phi_t (1 - p_{t+1}) \chi_{t+1} & \quad\text{if } t < T \end{cases} \]
The base case arises because if an animal was captured in the last time period, the probability it is never captured again is 1 because there are no more capture periods. The recursive case defining \(\chi_{t}\) in terms of \(\chi_{t+1}\) involves two possibilities: (1) not surviving to the next time period, with probability \((1 - \phi_t)\), or (2) surviving to the next time period with probability \(\phi_t\), not being captured in the next time period with probability \((1 - p_{t+1})\), and not being captured again after being alive in period \(t+1\) with probability \(\chi_{t+1}\).
With three capture times, there are eight captured/not-captured profiles an individual may have. These may be naturally coded as binary numbers as follows.
| profile | 1 | 2 | 3 | probability |
|---|---|---|---|---|
| 0 | \(-\) | \(-\) | \(-\) | n/a |
| 1 | \(-\) | \(-\) | \(+\) | n/a |
| 2 | \(-\) | \(+\) | \(-\) | \(\chi_2\) |
| 3 | \(-\) | \(+\) | \(+\) | \(\phi_2 \, p_3\) |
| 4 | \(+\) | \(-\) | \(-\) | \(\chi_1\) |
| 5 | \(+\) | \(-\) | \(+\) | \(\phi_1 \, (1 - p_2) \, \phi_2 \, p_3\) |
| 6 | \(+\) | \(+\) | \(-\) | \(\phi_1 \, p_2 \, \chi_2\) |
| 7 | \(+\) | \(+\) | \(+\) | \(\phi_1 \, p_2 \, \phi_2 \, p_3\) |
History 0, for animals that are never captured, is unobservable because only animals that are captured are observed. History 1, for animals that are only captured in the last round, provides no information for the CJS model, because capture/non-capture status is only informative when conditioned on earlier captures. For the remaining cases, the contribution to the likelihood is provided in the final column.
By defining these probabilities in terms of \(\chi\) directly, there is no need for a latent binary parameter indicating whether an animal is alive at time \(t\) or not. The definition of \(\chi\) is typically used to define the likelihood (i.e., marginalize out the latent discrete parameter) for the CJS model (Schofield 2007).
The Stan model defines \(\chi\) as a transformed parameter based on
parameters \(\phi\) and \(p\). In the model block, the log probability is
incremented for each history based on its count. This second step is
similar to collecting Bernoulli observations into a binomial or
categorical observations into a multinomial, only it is coded directly
in the Stan program using target += rather than
being part of a built-in probability function.
The following is the Stan program for the Cormack-Jolly-Seber mark-recapture model that considers counts of individuals with observation histories of being observed or not in three capture periods
data {
array[7] int<lower=0> history;
}
parameters {
array[2] real<lower=0, upper=1> phi;
array[3] real<lower=0, upper=1> p;
}
transformed parameters {
array[2] real<lower=0, upper=1> chi;
chi[2] = (1 - phi[2]) + phi[2] * (1 - p[3]);
chi[1] = (1 - phi[1]) + phi[1] * (1 - p[2]) * chi[2];
}
model {
target += history[2] * log(chi[2]);
target += history[3] * (log(phi[2]) + log(p[3]));
target += history[4] * (log(chi[1]));
target += history[5] * (log(phi[1]) + log1m(p[2])
+ log(phi[2]) + log(p[3]));
target += history[6] * (log(phi[1]) + log(p[2])
+ log(chi[2]));
target += history[7] * (log(phi[1]) + log(p[2])
+ log(phi[2]) + log(p[3]));
}
generated quantities {
real<lower=0, upper=1> beta3;
beta3 = phi[2] * p[3];
}Identifiability
The parameters \(\phi_2\) and \(p_3\), the probability of death at time 2
and probability of capture at time 3 are not identifiable, because both
may be used to account for lack of capture at time 3. Their product,
\(\beta_3 = \phi_2 \, p_3\), is identified. The Stan model defines
beta3 as a generated quantity. Unidentified parameters pose a
problem for Stan’s samplers’ adaptation. Although the problem posed
for adaptation is mild here because the parameters are bounded and
thus have proper uniform priors, it would be better to formulate an
identified parameterization. One way to do this would be to formulate
a hierarchical model for the \(p\) and \(\phi\) parameters.
Individual Cormack-Jolly-Seber model
This section presents a version of the Cormack-Jolly-Seber (CJS) model cast at the individual level rather than collectively as in the previous subsection. It also extends the model to allow an arbitrary number of time periods. The data will consist of the number \(T\) of capture events, the number \(I\) of individuals, and a boolean flag \(y_{i,t}\) indicating if individual \(i\) was observed at time \(t\). In Stan,
data {
int<lower=2> T;
int<lower=0> I;
array[I, T] int<lower=0, upper=1> y;
}The advantages to the individual-level model is that it becomes possible to add individual “random effects” that affect survival or capture probability, as well as to avoid the combinatorics involved in unfolding \(2^T\) observation histories for \(T\) capture times.
Utility functions
The individual CJS model is written involves several function definitions. The first two are used in the transformed data block to compute the first and last time period in which an animal was captured.18
functions {
int first_capture(array[] int y_i) {
for (k in 1:size(y_i)) {
if (y_i[k]) {
return k;
}
}
return 0;
}
int last_capture(array[] int y_i) {
for (k_rev in 0:(size(y_i) - 1)) {
int k;
k = size(y_i) - k_rev;
if (y_i[k]) {
return k;
}
}
return 0;
}
// ...
}These two functions are used to define the first and last capture time for each individual in the transformed data block.19
transformed data {
array[I] int<lower=0, upper=T> first;
array[I] int<lower=0, upper=T> last;
vector<lower=0, upper=I>[T] n_captured;
for (i in 1:I) {
first[i] = first_capture(y[i]);
}
for (i in 1:I) {
last[i] = last_capture(y[i]);
}
n_captured = rep_vector(0, T);
for (t in 1:T) {
for (i in 1:I) {
if (y[i, t]) {
n_captured[t] = n_captured[t] + 1;
}
}
}
}The transformed data block also defines n_captured[t], which is
the total number of captures at time t. The variable
n_captured is defined as a vector instead of an integer array
so that it can be used in an elementwise vector operation in the generated
quantities block to model the population estimates at each time point.
The parameters and transformed parameters are as before, but now there
is a function definition for computing the entire vector chi, the
probability that if an individual is alive at t that it will
never be captured again.
parameters {
vector<lower=0, upper=1>[T - 1] phi;
vector<lower=0, upper=1>[T] p;
}
transformed parameters {
vector<lower=0, upper=1>[T] chi;
chi = prob_uncaptured(T, p, phi);
}The definition of prob_uncaptured, from the functions block,
is
functions {
// ...
vector prob_uncaptured(int T, vector p, vector phi) {
vector[T] chi;
chi[T] = 1.0;
for (t in 1:(T - 1)) {
int t_curr;
int t_next;
t_curr = T - t;
t_next = t_curr + 1;
chi[t_curr] = (1 - phi[t_curr])
+ phi[t_curr]
* (1 - p[t_next])
* chi[t_next];
}
return chi;
}
}The function definition directly follows the mathematical definition
of \(\chi_t\), unrolling the recursion into an iteration and
defining the elements of chi from T down to 1.
The model
Given the precomputed quantities, the model block directly encodes the
CJS model’s log likelihood function. All parameters are left with
their default uniform priors and the model simply encodes the log
probability of the observations q given the parameters p
and phi as well as the transformed parameter chi defined
in terms of p and phi.
model {
for (i in 1:I) {
if (first[i] > 0) {
for (t in (first[i]+1):last[i]) {
1 ~ bernoulli(phi[t - 1]);
y[i, t] ~ bernoulli(p[t]);
}
1 ~ bernoulli(chi[last[i]]);
}
}
}The outer loop is over individuals, conditional skipping individuals
i which are never captured. The never-captured check depends
on the convention of the first-capture and last-capture functions
returning 0 for first if an individual is never captured.
The inner loop for individual i first increments the log
probability based on the survival of the individual with probability
phi[t - 1]. The outcome of 1 is fixed because the individual
must survive between the first and last capture (i.e., no zombies).
The loop starts after the first capture, because all
information in the CJS model is conditional on the first capture.
In the inner loop, the observed capture status y[i, t] for
individual i at time t has a Bernoulli distribution
based on the capture probability p[t] at time t.
After the inner loop, the probability of an animal never being seen
again after being observed at time last[i] is included, because
last[i] was defined to be the last time period in which animal
i was observed.
Identified parameters
As with the collective model described in the previous subsection,
this model does not identify phi[T - 1] and p[T], but
does identify their product, beta. Thus beta is defined
as a generated quantity to monitor convergence and report.
generated quantities {
real beta;
// ...
beta = phi[T - 1] * p[T];
// ...
}The parameter p[1] is also not modeled and will just be uniform
between 0 and 1. A more finely articulated model might have a
hierarchical or time-series component, in which case p[1] would
be an unknown initial condition and both phi[T - 1] and
p[T] could be identified.
Population size estimates
The generated quantities also calculates an estimate of the population
mean at each time t in the same way as in the simple
mark-recapture model as the number of individuals captured at time
t divided by the probability of capture at time t. This
is done with the elementwise division operation for vectors
(./) in the generated quantities block.
generated quantities {
// ...
vector<lower=0>[T] pop;
// ...
pop = n_captured ./ p;
pop[1] = -1;
}Generalizing to individual effects
All individuals are modeled as having the same capture probability, but this model could be easily generalized to use a logistic regression here based on individual-level inputs to be used as predictors.
References
An alternative would be to compute this on the outside and feed it into the Stan model as preprocessed data. Yet another alternative encoding would be a sparse one recording only the capture events along with their time and identifying the individual captured.↩︎
Both functions return 0 if the individual represented by the input array was never captured. Individuals with no captures are not relevant for estimating the model because all probability statements are conditional on earlier captures. Typically they would be removed from the data, but the program allows them to be included even though they make not contribution to the log probability function.↩︎