A review of Bayesian state-space modelling of capture–recapture–recovery data

1.1. State-space models

State-space models separate an observed process into two components: a system process that models the underlying biological process over time and an observation process that accounts for imperfect detection of the system process, such as measurement error. State-space models have become increasingly popular within ecological problems, often separating the nuisance parameters (observational error) from the biological processes of interest.

In general, the state-space framework takes the following form. Let y₁, … , y_T denote the (possibly multivariate) observations of the process; x₁, … , x_T the true underlying states at times t = 1, … ,T and θ the set of model parameters. The state-space model can be described by,

When the underlying states of the system process, x₁, … , x_T, are discrete-valued the state-space model reduces to a traditional hidden Markov model (although we note that some authors have extended the term hidden Markov model to include continuous state variables, see Cappé et al. [32]). For such models with finitely many states, there is an explicit closed form likelihood function for the observed data, denoted by f(y|θ) (by summing over all possible states in the system process) which can be maximized to obtain parameter estimates. Alternatively, the EM algorithm [33] can be implemented to obtain parameter estimates (see also van Deusen [34] for an application of the EM algorithm to capture–recapture data). For a continuous system process, in general there is no such closed form expression for the likelihood function of the observed data (although it can be expressed in the form of an integral over the x values) with the Kalman filter [35] traditionally used to obtain parameter estimates, assuming a linear Gaussian state-space model. More recently, Bayesian model-fitting approaches have been used that permit a greater degree of flexibility including, for example, nonlinear and non-Gaussian processes. We focus on the Bayesian approach and briefly describe a data augmentation approach for fitting state-space models. See, Newman et al. [36] for further discussion.

1.2. Bayesian model fitting

In order to fit a state-space model, we use a Bayesian data augmentation approach [37], also referred to as an auxiliary or latent variable approach, and, more recently, a complete data likelihood approach [38,39]. We treat the underlying states x as parameters (or auxiliary variables) to be estimated in addition to the model parameters, θ. The joint posterior distribution of the model parameters and auxiliary variables can be written in the form,

However, the necessary integration to obtain the posterior distribution is analytically intractable. We will use Markov chain Monte Carlo (MCMC) to obtain a sample from the joint posterior distribution of x and θ. Considering only the sampled values of θ (and ignoring the x values) indirectly performs the necessary integration to obtain a set of sampled values from the marginal distribution π(θ|y) which can be used to obtain summary statistics of interest. We note that the posterior (marginal) distribution of x is often of interest itself, relating to the posterior distribution of the true underlying states of the system under study.

For some studies, some of the underlying states x may be observed without error. In such cases, we let x = {x_obs, x_mis}, where x_obs and x_mis denote the observed and unobserved states, respectively. x_obs is observed data and so regarded in the same manner as the observations y, while x_mis are auxiliary variables and we can write

1.1

At each iteration of the MCMC algorithm, the auxiliary variables x_mis (in addition to the model parameters θ) are updated.

The introduction of auxiliary variables is often used in two particular circumstances: (i) when the likelihood function of the observed data, f(y|θ), is analytically intractable or computationally expensive; and/or (ii) when the inclusion of the additional auxiliary variables simplifies or improves the MCMC algorithm. Thus, we note that an auxiliary variable approach may be used even if an explicit likelihood is available. For example, for hidden Markov models (i.e. discrete system processes), there is an explicit likelihood. However, using an auxiliary variable approach (where x_mis corresponds to the discrete underlying states) leads to a complete data likelihood, f(y|x,θ), that can be significantly computationally faster to calculate than the observed data likelihood, f(y|θ). The auxiliary variable approach may also lead to posterior conditional distributions of standard form for the model parameters (and auxiliary variables) and use of the Gibbs sampler compared with the Metropolis–Hastings algorithm when using the observed data likelihood due to non-standard posterior conditional distributions. See for example [27,36,39–42] for further discussion of Bayesian model fitting in the context of a range of ecological state-space models.

Finally, we comment on the issue of parameter redundancy, or local identifiability, within the Bayesian framework. A parameter redundant model has a ridge in the corresponding likelihood surface leading to non-unique maximum-likelihood estimates. Traditionally, to remove this problem, restrictions are specified on the parameter space to allow estimation of all identifiable parameters. In general, parameter redundancy can be difficult to identify, particularly for models with a large number of parameters and/or processes acting, although some advances have been made using symbolic algebraic approaches [43–47]. However, within the Bayesian framework, the likelihood is combined with the prior distribution to form the posterior distribution of the parameters. Thus, parameter redundancy still permits the estimation of the posterior distribution (assuming a proper prior distribution), but due to the identifiability issues, parameters (or functions of parameters) suffer from strong prior sensitivity. For example, suppose that a given parameter is unidentifiable (so that the likelihood is flat with respect to this parameter), then the posterior distribution for this parameter is simply the corresponding prior for the parameter. The idea extends to when functions of parameters are unidentifiable, and the individual parameters within this function will typically exhibit strong prior dependence. Gimenez et al. [48] used this idea to assess potential identifiability issues, by considering the overlap between the posterior and prior distributions. A large overlap may provide an indication that the parameters are not being well estimated and potential parameter identifiability issues. See also King et al. [41] and Kéry & Schaub [42] for further discussion and applications.

1.2.1. Model discrimination

The Bayesian approach can be extended to incorporate model uncertainty. For example, suppose that there are a total of M possible models. We introduce a discrete model indicator parameter, denoted by m ∈ {1, … , M} with corresponding parameters θ_m (since the set of parameters is model dependent). The joint posterior distribution over the combined parameter and model space is given by

where p(m) denotes the prior specified on model m. Marginal posterior model probabilities are often of particular interest, providing a quantitative comparison of competing models, with the posterior probability for model m given by

Bayes factors, defined to be the ratio of posterior model probabilities to prior model probabilities, are also frequently used to compare competing models [49]. Posterior model probabilities (or Bayes factors) are a quantitative comparison of competing models and hence naturally permit posterior model-averaged estimates of parameters of interest. For example, suppose that the parameter θ₁ is common to all models. The posterior model averaged distribution for θ₁ is given by

MCMC methods can again be used to obtain a sample from the posterior distribution of interest to obtain Monte Carlo estimates of interest. For example the posterior model probabilities are estimated as the proportion of time the chain is in each model. However, an extension to the standard MCMC algorithm is necessary in order for the constructed Markov chain to move between models with a different number of parameters. The most common such algorithm is the reversible jump MCMC algorithm [50] but also see, Godsill [51] for additional discussion. See King et al. [41] and references therein for further discussion of ecological applications.

2.1. State-space model formulation

The state-space model formulation of the CJS model has been proposed independently by Royle [18], Gimenez et al. [61] and Schofield & Barker [62] (although this modelling approach dates back to (at least) Dupuis [63]—see §2.3.1). The state-space approach separates the survival process from the observation processes of individuals observed within the study (either alive or dead). We note that (as pointed out by an anonymous referee) this underlying concept can be seen within the original seminal papers, although not formally expressed in a (conditional) state-space representation. For example, both Jolly [3] and Seber [4] describe (in words) the two separate processes acting on each individual member of the population of surviving from one capture event to the next and subsequently being captured (or not).

Formally, within the state-space formulation, capture–recapture–recovery data can be viewed as the combination of two different processes: (i) individuals surviving between successive time points and (ii) individuals being observed (either alive or dead) with some probability. Thus, we can regard each encounter history as a combination of the processes of survival between each capture event and observations of whether they are recaptured (if alive) or recovered (if dead), and thus dissect w into these distinct (conditional) components. Mathematically, to represent the survival process, we define x to be the n × T matrix such that

We note that x_it = 0 corresponds to the cases where an individual is either dead or has not yet entered the study (i.e. before it is first observed). We let f(i) denote the first time individual i is observed in the study and ϕ_it the probability individual i survives until time t + 1, given they are alive at time t, with ϕ = {ϕ_it, i = 1, … , n; t = 1, … , T − 1}. Submodels are obtained by setting restrictions on these parameters. For example, a constant survival model is obtained by setting ϕ_it = ϕ for all i = 1, …, n and t = 1, … , T − 1. However, we retain the more general notation as we will consider a range of models allowing for temporal and/or individual heterogeneity. The underlying system process can be easily described as follows:

2.1

For capture–recapture–recovery data, we partition the observation process into two distinct processes corresponding to the (re-)capture (of live animals) and recovery (of dead animals). We initially consider the recapture process and let y denote the n × T matrix such that

We note that for an individual to be recaptured, it is conditional on the individual being alive in the study at that time (i.e. if x_it = 0 then y_it = 0). Let p_it denote the probability individual i is recaptured at time t, given they are alive in the study area at this time, for i = 1, … , n and t = 2, … , T. The observation process can be written as

2.2

The recovery process is defined similarly. For simplicity, we assume that the recapture and recovery pro-cesses are observed over the same period of time, though this can be easily relaxed. We let z denote the n × T matrix with elements

2.3

The observed encounter histories w can be separated into the data components x_obs, y and z. Note that a glossary of these latent variables, or auxiliary variables (and others introduced later in the manuscript), are listed in table 4 in appendix A. In addition, we note that this state-space formulation is simply a (partial) hidden Markov model, with the system process composed of the two discrete states ‘alive’ and ‘dead’, with some of the states known. The corresponding underlying state-space (or hidden Markov model) representation is provided graphically in figure 1.

• Download figure
• Open in new tab
• Download PowerPoint

We fit the models using a Bayesian data augmentation approach to form the joint posterior distribution (cf. equation (1.1)),

Using a state-space framework, when considering deviations from the above model and/or additional complexities, it is not necessary to explicitly derive a new observed likelihood expression, but simply consider each individual process in turn. For example, the assumption of recovery only in the capture event following the death of an individual is often reasonable, due to the decay of marks used to uniquely identify individuals. However, in some cases, this assumption is not valid, as for abalone (see Catchpole et al. [64], who derive an explicit observed data likelihood for the model where individuals can be recovered at any point following death). The model proposed by Catchpole et al. [64] can be easily fitted within the state-space framework changing only the observation process for recoveries such that

2.2. Estimating abundance

Removing the conditioning of initial capture from the CJS model permits the estimation of birth rates and abundance via the JS model [3,4]. An observed data likelihood can be obtained [58–60] by considering a superpopulation corresponding to all possible individuals who may be observed within the study (i.e. individuals present and available for capture for at least one capture event) and summing over all possible combinations of their entry and exit times (i.e. ‘birth’ and ‘death’ times). Entry to the system could be via birth or immigration and exit from the system, death or permanent emigration. However, for simplicity, we use the terms ‘birth’ and ‘death’ to include such migration. Thus, for the JS model, there is an additional entry process to be included within the system process. We let the number of individuals in the superpopulation be denoted by N, which is a parameter to be estimated. We note that (N − n) individuals are available for capture at least once within the study but are not observed and that x is of dimension N × T, and hence a function of N. We let b_it denote the probability that individual i enters the study population in the interval (t − 1,t] and so is available for capture for the first time at capture event t, for t = 1, … , T. This b_i1 is interpreted as the probability individual i is in the population and available for capture at the first capture event. Since we condition on all N individuals being available for capture at least once within the study period, we have that ∑_t=1^Tb_it = 1 for all i = 1, … , N.

Within the state-space model, the survival process is replaced by a ‘presence’ process allowing for individuals to both join and leave the study (assuming individuals can only join and leave the study a maximum of one time), corresponding to ‘birth’ and ‘death’. In particular, we have that for t = 1, … , T and i = 1, … , N,

2.4

The term ensures that an individual can only be ‘born’ once within the study and cannot be ‘reborn’ (i.e. re-enter the study once it has died). The parameter γ_it can be interpreted as the probability individual i enters the population within the interval (t − 1,t], conditional on not entering the study population before or at time t − 1. The recapture and recovery processes essentially remain the same as for the CJS model, given in equations (2.2) and (2.3), but are now specified for t = 1, … , T, since the JS model removes the conditioning on first capture. The posterior distribution for the JS model, assuming no recoveries, is given by

We note that when specifying the model within the set of BUGS (and associated) computer packages, it is not possible to specify x as a function of a stochastic node (i.e. N). The solution proposed is to specify a fixed upper bound for N, denoted by N*, and let x be of fixed dimension (N* × T) allowing for (N* − N) individuals not available for capture within the study period, or, in other words, do not enter into the study population. See earlier studies [40,42,65–68] for further discussion of this approach, alternative parametrizations and associated computer codes. Finally, we note the link of the JS model to two other common models of occupancy models and closed population models.

Occupancy models relate to the occupancy of a species (rather than individuals) within a given area. Species can colonize an area (enter the population), die out from the area (exit the population), as for capture–recapture data, but also recolonize an area (i.e. be ‘reborn’). The state-space formulation of the JS model can be easily extended to an occupancy model, by adding the term ) (1 − X_{it − 1}) to ν_it in (2.4) for capture times t = 3, … , T to allow re-entry into the population (see Royle and Dorazio [40] for further discussion and references therein in relation to occupancy models). Finally, we note that removing the presence process (assuming that all individuals are present for the full duration of the study, so that x_it = 1 for all i = 1, … , N and t = 1, … , T) and considering only a recapture process, leads to the special case of a closed population where the estimation of the total population, N, is of particular interest. See earlier research [39–42,67–69] for a range of different closed population models.

2.3. Individual heterogeneity

In many biological systems, the model parameters may be individual specific, in that they differ between individuals despite identical environmental conditions. For example, individuals typically differ in body condition (which itself may be a function of feeding or breeding behaviour, parasite load, etc.) which in turn has a direct impact on their corresponding survival (animals in better condition are more able to compete for resources and/or survive harsh environmental conditions). Similarly, animals may disperse through an area, with individuals more/less likely to be recaptured/recovered in more remote or dense vegetation areas. To incorporate such heterogeneity individual level factors (or covariates) and/or random effects are often used to model such differences.

For simplicity, we consider capture–recapture data and the CJS model, extensions allowing for further system and/or observation processes follow by including additional separate processes. We consider the model parameters as a function of individual level time-varying covariates. The matrix of observed encounter histories, w, has elements given by

The survival and recapture processes for x and y given in equations (2.1) and (2.2) have the same definitions as before but where we now have the survival and recapture probabilities expressed as a function of the covariate values. However, we need an additional system process corresponding to the transition process for the underlying covariates. For the transition process, we let u be an n × T matrix, with elements

2.5

2.3.1. Multi-state data

Discrete-valued time-varying covariates, such as breeding state or discrete location, are most commonly referred to as multi-state (or multi-site) data. These models are an example of a (partially) hidden Markov model, as the system processes have a finite number of states. Notationally, we let the set of possible states for the covariates be denoted by 𝒦 = {1, … , K} (this does not include the state of death). The traditional Arnason–Schwarz (AS) model extends the CJS model by considering a first-order Markov transition process between the different states, with an observed data likelihood available [9,10,70,71], allowing parameter estimation. Dupuis [63] proposed an alternative Bayesian state-space modelling approach, although did not phrase the approach as a state-space model but using an auxiliary variable approach (see also Kéry & Schaub [42, ch. 9]; Marin & Robert [72, ch. 5.5]; Clark et al. [73] and references therein and King & Brooks [54] for an extension to allow for model uncertainty in terms of time and covariate dependence on the model parameters).

We let the transition parameters ψ_it(j,k) denote the probability that individual i is in state k at time t + 1, given they are alive and in state j at time t for t = f(i), … , T − 1 and set ψ_it = {ψ_it(j,k): j,k ∈ 𝒦}. We initially consider a first-order Markovian structure for the transition process, but discuss relaxations to this assumption below. In particular, we can write

where

for i = 1, … ,n, t = f(i) + 1, … , T, ψ_it(k,·) = {ψ_it(k, j): j = 1, … , K} and s_if(i) = u_if(i). In other words, if x_it = 1 and u_it−1 = s_it−1 = j, the probability that individual i is in state k at time t is equal to ψ_it−1 (j,k); if x_it = 0 (individual is dead), then u_it = 0 (although s_it ∈ 𝒦). The survival and recapture parameters can be extended to be dependent on the additional state parameter, with the analogous Bernoulli processes occurring for the observation and system survival processes.

The joint probability mass function of the survival and transition processes, f(x, u|ϕ,ψ) given in equation (2.5) can be specified in conditional form

The first term corresponds to the Bernoulli survival process and the second to the multinomial transition process. In the standard AS model, if an individual is recaptured (i.e. y_it = 1), the corresponding state of the individual, u_it, is also known. However, if this is not the case and the state may not be observed when an individual is observed then within the state-space formulation the model specification remains the same. It is simply that the set of u_mis and u_obs values change to reflect the set of missing and observed states, and that any unobserved u_it values are auxiliary variables and updated within the MCMC algorithm.The state-space model can be extended to include additional system or observation processes. For example, within the traditional AS model, we assume that there are no recoveries within this process, but these can be easily included by considering an additional observation process corresponding to the recovery process [42,61,74]. Similarly, a birth process can be included (and a parameter for total population size) leading to the multi-state JS model [55,65].

Using a state-space framework, it is straightforward to fit more complex model structures, including a second-order transition process, such that the state of an individual at time t is dependent not only on their state at time t − 1 but also t − 2. An explicit observed data likelihood for a second-order Markov model is available [10,75] but is computationally expensive to calculate compared with the complete data likelihood, conditional on the unobserved states within the transition process. In addition, within the state-space framework the survival and observation processes remain the same, with only the transition system process needing to be modified. Let Ψ_it(j,k,l) denote the probability that individual i is in state l ∈ 𝒦 at time t + 2, given that they are alive and in states j ∈ 𝒦 and k ∈ 𝒦 at times t and t + 1, respectively, for t = f(i), … , T − 2. We let Ψ_it = {Ψ_it(j,k,l): j, k, l ∈ 𝒦}. The transition process for t = f(i) + 1, … , T can be described in the form

where

for t = f(i) + 2, … ,T and Ψ_it−2(s_it−2, s_it−1,·) = {Ψ_it−2 (s_it−2, s_it−1,l),l ∈ 𝒦}; and for t = f(i) + 1, typically, a first-order process is considered, so that

using the notation as for the first-order Markovian process and where s_if(i) = u_if(i).

Possible error within the state observation process can also be incorporated, leading to multi-event models proposed by Pradel [76]. In such models, the state of an individual is not observed directly, but an event relating to the state is observed. For example, states may refer to ‘breeding’ and ‘non-breeding’ and events ‘seen near nest’ and ‘not seen near nest’. In this case, there are no observed states (u_obs = ∅︀). The state-space framework extends the AS model specification by including an additional event observation process, conditional on the underlying states. We note that the event process can, equivalently, be incorporated into the state process, corresponding to the event an individual is in at each capture event (with the introduction of auxiliary variables for the unobserved events). However, for simplicity, we retain the event process within the observation process (but note that if the model parameters such as the recapture probabilities are dependent on the event process, rather than the true underlying states, the event process needs to be incorporated into the system process and the unobserved events treated as auxiliary variables).

Let ℰ = {1, … , E} denote the set of possible observed events, and let v_it denote the event observation process for individual i at time t, conditional on the individual being observed at this time (i.e. y_it = 1). For statistical convenience, we let v_it = 0 if y_it = 0 (i.e. the individual is not observed). Thus, we have

We define the additional model parameters β_it(k_it, e_it) to be the probability individual i is in event e_it ∈ ℰ at time t, given they are in state k_it ∈ 𝒦 and observed at this time and let β_it = {β_it(k,e): k ∈ 𝒦,e ∈ ℰ}. Then, for i = 1, … , n and t = f(i), … , T, we have

where

and β_it(u_it,·) = {β_it(u_it,r_it): r_it = 1, … , E}, such that ∑_e∈ℰ β_it(u_it, e) = 1 for all u_it ∈ 𝒦 (assuming that on initial capture the observed event is independent of the fact that it is their first encounter).

The underlying system process is the same as for the AS model, with the event process incorporated in the observation process. However, we note that within the AS model, the likelihood conditions on the state an individual is in when first observed. Within multi-event models, this conditioning is removed (since we do not observe the true states), and an initial state distribution needs to be specified. Let κ_i(k) denote the probability individual i is in state k when it is first encountered. The initial state distribution is given by

where κ_i = {κ_i(k): k = 1, … , K}. More generally, the multi-event framework allows for both observed known state and partially observed state where an observed event restricts the set of possible true states (see Smout et al. [77] who consider this framework in relation to tag-loss models and King & McCrea [78] who derive an efficient closed form observed data likelihood).

Finally, we consider the special case of mixture models proposed by Pledger et al. [79] to allow for unobservable heterogeneity within the population. Within these models, the population is composed of K homogeneous sub-populations, with each individual belonging to one sub-population and there are no transitions between sub-populations. Typically, the sub-populations are unobservable, so that there is no directly observable covariate information (u_obs = ∅︀). The covariate values are time-invariant, so that u_it = u_i for all, t = 1, … , T, and

where κ = {κ(k): k = 1, … , K} and κ(k) denotes the probability an individual belongs to sub-population k. In general, K may be known or be a parameter to be estimated [79].

2.3.2. Mixed/random effect models

We now consider continuous individual covariates. In these cases, the demographic parameters are typically specified as a parametric function of covariate values [11,12,15–20,41,42,65,80] or recently, more flexible semi-parametric spline models have also been proposed to model the relationship between the demographic parameters and covariates [14,81]. For example, let u_it = {u_it(j): j = 1, … , J} denote the set J covariate values for individual i at time t. The survival probability of individual i at time t may be (parametrically) logistically regressed on the set of covariates, such that

We note that if an individual is unobserved at any capture event (or following capture for the CJS model), their corresponding time-varying covariate(s) are also generally unobserved. This typically leads to an analytically intractable observed data likelihood, specified only in integral form, where the integration is over the unknown covariate values. Catchpole et al. [15] proposed an alternative conditional likelihood approach, conditioning on only the observed covariate values while Bonner et al. [16] compared this approach with an alternative Bayesian approach. The comparative study conducted suggested that, in general, the conditional likelihood approach suffers from poorer precision (as a result of discarding a potentially large proportion of the data) and can have boundary value problems for the parameter estimates, but has the potential advantage that no underlying model needs to be specified on the underlying covariate values, as in the Bayesian approach. In particular, the transition process for the covariate values are written in the form of the probability density function associated with the covariates (potentially allowing for dependence between the values). For example, Bonner & Schwarz [13] considered only one time-varying covariate and set

for t = f(i) + 1, … , T, dropping the notation on covariate j and where δ_t denotes an additive year effect. A probability density function is specified for u_if(i) if there are unobserved values at initial capture. The missing covariate values are treated as auxiliary variables and updated within the MCMC algorithm.

Including individual random effects to the model allows for additional unexplained individual heterogeneity, such that

where ε_i ∼ N(0, σ_ε²). Gimenez & Choquet [20] use a numerical integration procedure to approximate the observed data likelihood for random effect models on survival probabilities. Such an approach can be computationally very expensive for more than two individual random effects. Alternatively, within the Bayesian model-fitting approach, the additional random effect terms ε = {ε_i: i = 1, … , n} are treated as auxiliary variables that are updated within the Markov chain. The joint posterior distribution over the model parameters, unobserved system state values and random effect terms is given by

where denotes the probability density function for the random effects (independent Gaussian for the above random effects model). Model selection can be performed on the covariates present in the model as a variable selection problem [12,19,41].

2.4. Efficiency

In general, when there are competing model-fitting algorithms that can be used, their relative performance will be both model and data dependent. In order to investigate the computational efficiency of a state-space modelling approach, we consider capture–recapture data relating to dippers (see Lebreton et al. [82] for further details) where there are T = 7 capture occasions and compare different model-fitting algorithms. Brooks et al. [53] compare different time-dependent models for these data, identifying the model with constant recapture probability, p, and two survival probabilities, corresponding a flood effect year (years 2 and 3) with survival probability, ϕ_f, and non-flood effect year (years 1, 4, 5 and 6) with survival probability, ϕ_n, as having largest posterior support (assuming uniform priors on the recapture and survival probabilities). Using this temporal model, we specify additional individual random effects, so that the model parameters are of the form

2.6

We consider three model-fitting algorithms:

Of the three models fitted including individual heterogeneity, model {p(h), ϕ(f)}, corresponding to only individual random effects on the recapture probabilities, is preferred in terms of posterior model probabilities. Bayes factors of 3.9 and 32 are obtained for model {p(h), ϕ(f)} versus {p, ϕ(f + h)} and {p(h), ϕ(f + h)}, respectively. (However, we do note that the model excluding any random effects, {p, ϕ(f)}, is preferred with a Bayes factor of 5.2 versus model {p(h), ϕ(f)}.) We use model {p(h), ϕ(f)} to compare the performance of the different model-fitting algorithms including some level of individual heterogeneity. Note that for this model, we set q = logit⁻¹b, corresponding to the underlying mean recapture probabilities of individuals in the presence of individual random effects, and ϕ = {ϕ_n, ϕ_f} the survival probabilities for non-flood and flood years, respectively (no individual random effects present).

We consider the convergence of the Markov chains to the stationary distribution, using the Brooks–Gelman–Rubin (BGR) statistic [85]. The BGR statistic for parameters ϕ, b and σ_p² consistently tend to one within 10 000 iterations for repeated simulations starting from different initial values, though ϕ_n and ϕ_f always tend to one much quicker than the other parameters. This appears to be a result of the poorer mixing of the Markov chain for σ_p² (and b). We return to this issue below in relation to autocorrelation. To be conservative, we specify a burn-in of 20 000 iterations for the MCMC simulations, to ensure the chain has sufficiently converged.

To compare the relative performances of the MCMC algorithm for the different model-fitting approaches, accounting for autocorrelation present within the Markov chain, we consider the effective sample size (ESS) [86] of each of the parameters. The ESS is essentially the number of independent samples within the MCMC sampled values, correcting for autocorrelation within the Markov chain and are presented in table 3. There are clear differences in the ESS across the different parameters with the consistent signal that all the MCMC algorithms perform poorly in the updating of the random effect variance term σ_p, with very high autocorrelation. Generally, the L1 model-fitting approach appears to have the largest ESS, suggesting a better mixing Markov chain; however, we recall that this was also the slowest model-fitting algorithm. Table 3 provides the corresponding ESS per unit time (100 s), which calculates the effective number of independent samples per unit of computational time. Using this ESS per unit time, the state-space approaches SS1 and SS2 each perform better than the L1 approach, for all parameters in this example. However, the SS2 state-space parametrization (auxiliary variables correspond to the time of death of each individual) significantly outperforms the SS1 parametrization (auxiliary variables correspond to alive/dead state of each individual at each capture time following final capture). However, we note the very poor performance of the MCMC algorithm for the random effect variance across all algorithms, and subsequently q (or b), due to the high posterior correlation between b and σ_p.

We emphasize that, in general, the performance of different model-fitting approaches will be both problem and model dependent. For example, King [87] compares a state-space approach (using the Gibbs sampler) with the observed data likelihood approach (using a Metropolis–Hastings random walk algorithm) for the AS model using data analysed by Dupuis [63] consisting of 96 lizards, six capture events and three sites. In this case (using bespoke written code), the state-space approach was computationally faster, exhibited lower autocorrelation and did not require pilot-tuning compared with the observed data likelihood approach.