Hidden Markov models for extended batch data

Batch marking provides an important and efficient way to estimate the survival probabilities and population sizes of wild animals. It is particularly useful when dealing with animals that are difficult to mark individually. For the first time, we provide the likelihood for extended batch‐marking experiments. It is often the case that samples contain individuals that remain unmarked, due to time and other constraints, and this information has not previously been analyzed. We provide ways of modeling such information, including an open N‐mixture approach. We demonstrate that models for both marked and unmarked individuals are hidden Markov models; this provides a unified approach, and is the key to developing methods for fast likelihood computation and maximization. Likelihoods for marked and unmarked individuals can easily be combined using integrated population modeling. This allows the simultaneous estimation of population size and immigration, in addition to survival, as well as efficient estimation of standard errors and methods of model selection and evaluation, using standard likelihood techniques. Alternative methods for estimating population size are presented and compared. An illustration is provided by a weather‐loach data set, previously analyzed by means of a complex procedure of constructing a pseudo likelihood, the formation of estimating equations, the use of sandwich estimates of variance, and piecemeal estimation of population size. Simulation provides general validation of the hidden Markov model methods developed and demonstrates their excellent performance and efficiency. This is especially notable due to the large numbers of hidden states that may be typically required


Introduction
The standard protocol for capture-recapture studies of animals is to use individually numbered tags so that the capture history can be constructed, determining whether an individual was captured at each sampling occasion. Many sophisticated models have been built for this type of data (McCrea and Morgan, 2014). However, it may not be feasible to use individually numbered tags for some species because the small size of the animal makes it difficult to use a unique tag, or for reasons of cost and convenience. In these cases, batch marks can be used instead.
In extended batch marking, individuals are captured at several sampling occasions. Some of the unmarked individuals are given a common batch-mark, with different marks applied at different sampling occasions; without loss of generality, we shall talk in terms of different colored tags. At each subsequent capture occasion, marked individuals are counted and their tag colors are noted. A random sample of unmarked individuals are given a new colored mark and then all marked animals are released. Some unmarked individuals do not receive tags, and these animals may or may not be released; in the case of the weather-loach data that we analyze in this article, the unmarked individuals that are not tagged were not released.
There are many examples of batch marking in the literature, involving, for example, species of fish, insects, and amphibians. A wide variety of batch marks are available that are often species dependent, including tattoo, brand, fin-clip, o-rings, dyes, polymer, antibiotic, and radioisotope marks as possible fisheries batch marks, and dust, paint, dye, painted labels, self-marking baits, wire tags, genetic markers, and mutilation for marking insects.
Closed population models have been developed for extended batch-mark studies, aiming at estimating population size. For two-sample studies the Petersen estimator can be used to estimate population size and the Schnabel estimator (under model M t in Otis et al., 1978) is used with more than two sample occasions (Pine et al., 2003). Individual fish can be marked and released upriver, then caught downriver to enumerate salmon runs; these studies can be analyzed using a ratio estimator developed by Laplace (Rawson, 2009). For open populations, Skalski et al. (2009) review marking methods for small fish and analysis methods to estimate survival in both batch and uniquely tagged individuals. However, models for open populations and the use of standard likelihood methods are difficult because of the need to account for potential, non-identifiable, multiple detections of the same individuals over the different sampling occasions.
An important advance was made by Huggins et al. (2010), hereafter denoted HWK, who provide methods for an extended batch mark study. Conditional on the number of individuals released at each sample time, they developed a pseudo likelihood for the recaptured individuals. They then derived estimating equations to estimate survival and capture probabilities, with error estimation obtained from a sandwich estimator. They show how population size can be derived using a Horvitz-Thompson-like estimator also accompanied by a large-sample sandwich variance estimator. HWK do not model individuals that were captured but not marked; however, they did include these when estimating population size. During a week-long industrial modeling camp, Dang et al. (2009) derived a likelihood for unmarked individuals, but ignored the dependencies between the numbers of individuals captured at successive times. Cowen et al. (2014) developed a likelihood approach for the marked data, and compared the efficiency of the extended batch-mark study, analyzed using both likelihood and pseudo-likelihood, with a traditional capture-recapture study (using the Jolly-Seber model). Although the analysis only involved marked individuals, direct computation of the likelihood was a formidable computational exercise.
In this article, we present hidden Markov models (HMMs) for the extended batch-mark survey incorporating both marked and unmarked individuals captured and released at each sampling occasion. We illustrate our methods using data on the oriental weather-loach, Misgurnus anguillicaudatu, studied by HWK.
The article is structured as follows: the data motivating the work are described and presented in Section 2; Section 3 lists the notation used, including model specification, and explains the relevance of HMMs; Section 4 presents the two likelihood components, corresponding to marked and unmarked individuals. It is explained how the component likelihoods are efficiently computed using the methodology of HMMs. A framework for model selection and evaluation is provided in Section 5 and alternative procedures for estimating population size are presented in Section 6. The results are given in Section 7 and the article ends with discussion and avenues for future research in Section 8.

Sampling Design and Data
Data from marked individuals can be summarized in a similar way to the m-arrays used in band-recovery experiments (Brownie et al., 1985), but individuals can be captured on multiple occasions, and significantly the use of multinomial distributions is no longer appropriate.
The weather-loach study is described in detail by HWK. Here, different colored batch tags were given to a random sample of unmarked individuals at each occasion. The data are provided in Table 1 for this 11-sample occasion study. We do not use the information on the numbers lost ( t ), but it is presented to illuminate the data as a whole. The analysis of Table 1 Weather-loach batch mark data array taken from HWK. The number unmarked is the number of individuals caught without a mark, before marking takes place. Thus, trivially this number is 306 at the first sampling occasion. Of this number, 280 individuals are given a batch mark, and 32 of them are recaptured at sample occasion 2, 22 at occasion 3, etc. In order to understand the table structure, note, for example, that at sample occasion 4, 207 = 23 + 17 + 20 + 147, etc. To illustrate the notation, T = 11, G = 10, S 2 = 219, R 2 = 139, u 2 = 187, r 23 = 28, r 24 = 17, etc. At each sample occasion a number of the fish sampled are not marked, and they are not returned to the water, becoming lost to the study. Thus, for example, at recapture occasion 2, 48 = 187 − 139, etc.

Recapture occasion Sample
Number marked individuals is conditional upon the numbers of marked individuals released, and the analysis of unmarked individuals is based upon the numbers of unmarked individuals sampled, and so the lost individuals do not bias the analyses. Note that the number of sampling occasions (T = 11) and the number of individuals marked at the first occasion (R 1 = 280) are sufficiently large for direct computation of the likelihood, as in Cowen et al. (2014), to be computationally prohibitive, as discussed below.

Models and Assumptions
We allow for immigration and emigration/death between sample times, but all emigration is assumed permanent. Other assumptions are similar to those of typical capture-recapture experiments namely: all individuals have the same probability of survival between sample times, all individuals have the same probability of capture at a sample time, tags are not lost, the color of the tags is identifiable, sampling is instantaneous, and individuals are independent with respect to capture and survival.
As ages of individuals are unknown in the case study, we only consider possible time-dependence in model parameters. We specify models using a standard notation so that the four models for marked individuals are (φ, p), (φ t , p), (φ, p t ), and (φ t , p t ), where φ and p denote, respectively, survival and recapture probabilities, which may be constant or time dependent. The (φ t , p t ) model is the batch-marking version of the Cormack-Jolly-Seber model, see Buckland and Morgan (2016) and McCrea and Morgan (2014, p. 70). There are also no covariates in the case study, but if relevant time-varying covariates are available then these might be accommodated by appropriate logistic regressions.
Hidden Markov models (HMMs) are a particular type of state-space model where the state space is discrete. Markrecapture clearly fits under the HMM framework; see King (2012). Here, a capture-recapture observation is generated by a distribution that is dependent on the state of an unobserved Markov process (Zucchini et al., 2016). For typical capture-recapture models there are two hidden states for each animal, and the true sequences of states are often only partially observed; when an animal is captured then it is known to be alive, whereas when it is not then its true state may be unknown. This is explained in detail in Laake (2013). A more complex example is provided by Chapter 24 of Zucchini et al. (2016), in which there are three hidden states at each time, according to whether individuals are alive, have died since the previous sampling time, or have died prior to that. In that application observed capture histories are then described in terms of survival, recapture, and recovery probabilities. For the extended batchmarking experiment, the hidden states are at the batch, rather than the individual level, and the true sequence of states is entirely unobserved. This gives rise to a large number of states.

Primary Notation
The following notation is used in specifying the likelihoods:

Constants and Statistics
T the number of capture-recapture occasions. G the number of batch-marked release groups; G ≤ T . S g the number of individuals sampled at sampling occasion g; g = 1, 2, . . . , G. R g the number of individuals marked and released at sampling occasion g from batch group g; g = 1, 2, . . . , G. We condition on the {R g } when we form the likelihood for marked individuals. r gt the number of individuals from batch group g recaptured at recapture occasion t; g = 1, 2, . . . , G, t = g + 1, . . . , T . u t the number of individuals captured at sampling occasion t that were not marked; t = 1, . . . , T .
t the number of individuals lost at sampling occasion t;

Latent Variables
X gt the number of individuals present at occasion t from marked group g; g = 1, 2, . . . , G, t = g + 1, . . . , T . d gt the number of individuals from marked group g that die at sampling occasion t.
Parameters gt the (R g + 1) × (R g + 1) state transition probability matrix of the Markov chain for the marked individuals of group g, describing transitions between sample occasions t and t + 1. u t the (U max + 1) × (U max + 1) state transition probability matrix of the Markov chain for the unmarked individuals describing transitions between sample occasions t and t + 1. P gt (m) the (R g + 1) × (R g + 1) diagonal matrix containing the state-dependent probabilities of observing m recaptures at occasion t for group g.
. δ g the initial distribution of the Markov chain for marked group g. δ u the initial distribution of the Markov chain for the unmarked population. φ t the probability of surviving and remaining in the population between occasions t and t + 1, given an individual was alive and in the population at occasion t. We use φ to indicate the set of survival parameters in a model. p t the probability of capture at occasion t. We use p to indicate the full set of recapture probabilities in a model, and p 2:T when p 1 is omitted. λ the initial mean abundance (at occasion 1) for the unmarked population.
η the recruitment rate into the unmarked population. Note that we shall consider constant and densitydependent versions of recruitment. U t the total number of unmarked individuals in the population available for capture at occasion t. We may use U to denote the full set of the numbers of unmarked individuals in the population at times t = 1, . . . , T. In the Student model of Section 4.2.1, these numbers are parameters which are estimated directly from the likelihood, whereas in the open N-mixture model the U are obtained as derived variables, which are functions of the other model parameters. U max the maximum number of unmarked individuals available for capture on any occasion. N t the total population size at time t.

Likelihood Constructions
As we can see from Table 1, in examples of extended batch-mark studies, not all individuals that are captured are marked. This may occur for reasons such as handling time constraints where animals must be released after being contained for a set period of time, or practical constraints, such as having a limited number of marks available. We shall provide the likelihood for marked individuals for extended batch-marking experiments. We shall also incorporate information on unmarked individuals into the likelihood by developing separate models for such individuals. This was not done by HWK or Cowen et al. (2014); it will allow us to include additional information on the capture probabilities in the analysis and also to estimate simultaneously the numbers of unmarked individuals in the population.

Likelihood for Marked Individuals
HWK regarded the likelihood for marked individuals as "intractable," and suggested that a possible EM approach would be "complicated and computationally intensive." Instead, they present the pseudo likelihood shown below, i=g φ i denotes the probability that an individual released from group g is recaptured at occasion t. In fact, we can obtain an expression for the likelihood for the data from marked individuals by conditioning upon the unknown numbers of dead individuals, {d gt }, to obtain Evaluating the conditional probabilities in equation (1) then results in the explicit expression for the likelihood, where π gk = (1 − φ k ) k−1 j=1 φ j , and by convention, 0 j=1 = 1. Direct computation of the likelihood L m can be slow, due to the evaluation of the multiple summations involved, which is O((R 1 + 1) T ), equating to O(281 11 ) for the case study. Consequently, we have only been able to maximize the likelihood in this form for the weather-loach data from the first 7 samples only. However, the results provide a useful check of the HMM computations which follow.
The model for each release group is a HMM, with the {d gt } being the hidden information, described by the conditional multinomial distributions in equation (1). It is the fact that the {d gt } are unknown/hidden that results in the expensive summations in equation (2). We can therefore make use of the efficiency of the standard forward probability approach for HMMs to compute the likelihood, in addition to other benefits; see Zucchini et al. (2016, p. 37). Therefore, the HMM likelihood component for release group g of the marked individuals, L m,g , can be written in the usual way as a product of the initial distribution vector δ g , the appropriate transition probability matrices { gt }, and the state-dependent probability matrices {P gt } for each sample occasion t. Thus, from Zucchini et al. (2016, p. 37) we can write, where 1 denotes the unit row vector, with P(r gt |X gt = i) = 0 for i < r gt , and otherwise we have the binomial forms: For each sample time t, the state transition probability matrix gt has elements P(X g,t+1 = j|X gt = i), and has the form, .
Each row of gt corresponds to a binomial distribution, and the rows sum to unity, as required. As the number of animals alive at t = 0 for group g is known to be R g for each g, the initial state distribution for group g is given by = (0, 0, . . . , 0, 1).
If we assume independence between release groups, the likelihood for the marked individuals is then given as This is exactly the same expression as equation (2) This likelihood has the same structural form as the first part of the Jolly-Seber likelihood, (McCrea and Morgan, 2014, p. 150), and we shall refer to the underlying model as the Student model. The model includes U, the elements of which denote the numbers of unmarked individuals in the population at the different times, as a set of parameters to be estimated. By itself of course this likelihood is over parameterized, as there are two unknowns for each degree of freedom. There are various ways of dealing with this, and a referee has suggested using Bayesian modeling with a joint regularization prior on the {U t } with a positive correlation structure to constrain the {U t }, for example, a random walk of order 2 or a Gaussian Markov random field; see Schmidt et al. (2015). Here, we combine the likelihood with the likelihood for the marked individuals. The likelihoods for the marked and unmarked individuals can be multiplied to get the overall joint likelihood, as there are no individuals in common: This is an illustration of integrated population modeling; see Besbeas et al. (2002). It is possible to maximize L j (φ, p, U; {r gt }, {R g }, {u t }), after some minor adjustments described later. However, if the {U t } have no structure such a model is still too general, and in particular there is no involvement of φ in L u (p, U; {u t }). An alternative structural approach is described below.

Open N-mixture approach.
Here, we introduce structure, by relating the {U t } over time in a stochastic manner, using an open N-mixture model for the unmarked individuals. This model adopts a first-order Markov dependence; see Dail and Madsen (2011). The likelihood for the unmarked individuals is then given by where P(U 1 ) is the probability function for U 1 and P(U t |U t−1 ) is the probability of U t conditional upon the value of U t−1 . The total number of unmarked individuals in the population at time t is unknown, and following Dail and Madsen (2011), after experimentation we assume this to be less than some large number, U max , for all t. Thus, U max can be taken as the upper limit of all of the summations in equation (6). For related discussion, see Dennis et al. (2015). One might similarly introduce an appropriate U min term, which could be beneficial in some cases. We can see how the binomial expression of the Student model, in equation (4), is the basis for the binomial term in equation (6), although the interpretation of {U t } is different between the two models. Once again we have a HMM, due to the Markov structure, and the computationally expensive summations in the likelihood expression can be seen to be a consequence of the unknown numbers of uncaptured individuals at each occasion; cf equation (1). In order to obtain P(U t |U t−1 ), we write U t = A t + C t where A t is the number of individuals that have survived from occasion t − 1 to t, and C t is the number of individuals recruited into the population at occasion t. We model A t |U t−1 ∼ Binomial(U t−1 − S t−1 , φ t−1 ) and C t |U t−1 ∼ Poisson(U t−1 η), although other distributions are possible, for example to account for overdispersion. The conditional distribution of C t provides a form of density dependence of recruitment, and an alternative is simply to use Poisson(η); we have used both in the case study of the next section. The initial number of unmarked individuals in the population is modeled as U 1 ∼ Poisson(λ), and P(U t |U t−1 ) is given by the convolution sum, (7) and zero otherwise, where Bin(x; n, p) and Pois(x; μ) denote the probability functions of the binomial and Poisson distributions, respectively. In view of the typically very large number of terms in the summations of equation (7) the Fast Fourier Transform was used to evaluate the summations efficiently.
As the open N-mixture model is also a HMM, we can write the likelihood of equation (6) in the form, and P(u t |U t = i) = 0 for i < u t , and otherwise we have the binomial forms: P(u t |U t = i) = i ut p ut t (1 − p t ) i−ut , for i ≥ u t . This is clearly the same approach as that adopted in Section 4.1.
The state transition probability matrix for the unmarked population u t , describing both the survival and immigration of individuals, has elements given by equation (7), and the initial state distribution, δ u , is Poisson(λ), as stated above.
The joint likelihood in this case has the form, We thank a referee for pointing out similar work in Schmidt et al. (2015). In that article data are collected on known-fate individuals which are radio marked, and also unmarked individuals. Integrated modeling is employed: an open N-mixture model is used to model the information on unmarked individuals, and that likelihood component is also formed using a HMM.

Estimation, Model Selection, and
Goodness-of-Fit Hidden Markov modeling is challenging in the work of this article in view of the very large numbers of states that might be involved. We discuss this issue in Section 7.2. However, the computer programs written were effective and efficient. Parameters are estimated using numerical maximum likelihood methods, and we use the observed information matrix and the bootstrap to provide estimates of uncertainty. Confounded parameters (e.g., the product φ T −1 p T ) can be dealt with by introducing single parameters denoting such a product. To obtain estimated standard errors in general, we need to be able to deal with the Hessian being singular due to boundary estimates (e.g.,φ 7 = 1). We obtain the singularvalue decomposition (SVD) of the Hessian, H, which in standard notation can be written as H = WDV T ; see Searle (1982, p. 316). We then determine which values of the diagonal matrix D are close to zero, and remove the corresponding rows and columns of the SVD, resulting in H * . We then obtain standard errors from the inverse of H * using H * −1 = V * D * −1 W * T (Searle 1982, p. 318), and set standard errors equal to zero for the estimates corresponding to the deleted rows.
Model selection can be achieved using standard likelihood methods, see Besbeas et al. (2015), and we illustrate the model selection process for the case study using the Akaike information criterion (AIC) in Section 7. The only model fitted by HWK, to data on marked individuals, was that where survival probability varies by time and capture probability is constant, (φ t , p), which they justify in terms of how the sampling was carried out.
To study goodness-of-fit for marked individuals, observed and expected cell counts of the marked individuals can be compared straightforwardly, where expected cell counts are calculated as E(r gt ) = R g Q gt (φ, p), as done by HWK, who showed that for the weather-loach data the Pearson goodnessof-fit statistic for model (φ t , p) has the value X 2 = 21.62, which is insignificant (p-value > 10%) when referred to χ 2 tables on 16 degrees of freedom, after appropriate pooling of cells. For the model for unmarked individuals, we plot normal pseudo-residual segments. These arise as we are dealing with discrete random variables, so that the standard pseudo residuals, which are points defined for continuous random variables, are replaced by intervals (Zucchini et al., 2016, p. 101). Ordinary normal pseudo residuals (Zucchini et al., 2016, p. 102) could also be used for the model for marked individuals, but cohort lengths are typically short, and so instead we plot observed versus expected values in that case.

Estimating Population Size
Population size was estimated by HWK using a Horvitz-Thompson-like approach, resulting in wherep was estimated from the data on marked individuals. However, when we also model the number of unmarked individuals, we can form alternative estimates of population size. This is because we can estimate the numbers of both marked and unmarked individuals in different ways. For the open N-mixture model the estimated expected numbers of unmarked individuals are derived variables, functions of the estimated model parameters, and given by recursion; see Dail and Madsen (2011, p. 4). For instance, for the case of constant survival and density-dependent recruitment, we havê We shall refer back to this equation later in the article. We could use the multivariate delta method to estimate the variances of the estimated expectations in equation (11), but in the following we shall instead use the bootstrap. An alternative approach for estimating the numbers of unmarked individuals arises because of the HMM nature of the model for unmarked individuals, and results from application of the Viterbi algorithm, a dynamic programming algorithm which produces the most likely set of {Û t }; see Zucchini et al. (2016, p. 89). In the following, we shall refer to {Û t }, but when the Dail/Madsen approach is used, rather than the Viterbi algorithm, this is to be interpreted as {Ê(U t )}. In addition, whichever estimate of unmarked individuals is used, by adopting a Horvitz-Thompson-like approach, the number of individuals alive at occasion t can be obtained from,N Cf equation (10). Alternatively, we can estimate the numbers of marked individuals alive directly from the numbers of individuals marked in appropriate samples: We shall consider the relative merits of alternative methods for estimating population sizes in the analyses of the next section. However, we note here that our experience has been that the expression of equation (13) results in smaller standard errors than those of equation (12), and we only illustrate the former in the case study which follows. A further approach, suggested by a referee, would result from posterior sampling of numbers of marked and unmarked individuals following a Bayesian analysis, such as that outlined in Schmidt et al. (2015).

Case Study
We use the data of Table 1 to illustrate the methods of the article.

Analysis of Marked Data
To compare the HMM likelihood method with the pseudolikelihood (PL) method of HWK we present in Table 2 parameter estimates for the model (φ t , p) for data on marked individuals only, the sole case for which results were presented by HWK. We found that parameter estimates were in good agreement between the two methods for this case study; however, standard errors for the likelihood method were usually smaller, as independently suggested by the simulation study of Cowen et al. (2014). We note that our estimates of error for PL differ substantially from what HWK reported, partly due to errors both in their R code for obtaining standard error estimates and in the delta method that HWK used, and partly due to the likelihood method being more efficient. Note also that the sandwich standard errors are harder to compute than those resulting from inverting a Hessian in the normal way. This is due to the need to construct manually a modeldependent derivative matrix. Additionally, although most of the HMM standard errors are smaller than those from the PL approach, we have found the difference to be far greater, with HMM resulting in smaller errors, when smaller data sets are analyzed, as found, for example, in Haynes and Robinson (2011).
Several models were fitted using the HMM and their AIC values were compared (Table 3a). Interestingly, we found model (φ t , p) to be the best fitting model, the only model fitted by HWK, and the model (φ, p t ) performs well. Note that model (φ t , p t ) is parameter redundant, as only the product φ T −1 p T can be estimated, rather than either of the component terms.

Table 2
Parameter estimates and estimated standard errors (given below in parentheses) produced by the pseudo-likelihood method (PL) and the likelihood method (HMM) for model (φ t , p) for the marked weather-loach data. In the PL case we also give, in square parentheses, the incorrect standard errors from HWK.

Parameter estimates
Method

. Combining Marked and Unmarked Data
A complication with the analysis of the unmarked data is that a large value of U max was required. We experimented with several values and found U max = 1800 to work well for these data. However, this requires operations with matrices of dimension 1801 × 1801 for the computation of L u , which is memory (RAM) intensive for a personal computer. To deal with this, we allocate the states to bins of equal size w: we let ζ 1 , ζ 2 , . . . , ζ n be a partition of the state space such that each ζ i is of length w; the midpoint of the interval is taken to represent the state; see Zucchini et al. (2016, p. 158). For this case study, we experimented with values of w = 4, 10, 20, and 50, and found that taking w = 20 appeared to be satisfactory, though there was little difference between taking w = 4 and w = 20. This resulted in 90 hidden states, whereas, for example, taking w = 4 increased the number of states to 450. Binning in this way introduces an element of approximation to the formation of the likelihood component for unmarked individuals, the approximation increasing with w. For illustration, we present results for the density-dependent model for new individuals, as in equation (7), as these resulted in slightly better AIC values than the alternative, of constant augmentation. However, which alternative is more realistic might also depend on which agrees better with the biology of the study. When the unmarked individuals are incorporated into the likelihood, model selection may be done in the same way as before, using AIC; see Table 3b. Here, we find that model (φ t , p t ) fits the data best in both the Student and N-mixture approaches. However, uncritical use of AIC can result in overly complex models, and the models (φ, p t , U t ) and (φ, p t , λ, η) perform well in terms of AIC. Therefore, we shall use these for ease of illustration. We see in particular that models with constant recapture probability are not supported, in contrast to the findings of Table 3a, and the belief of HWK that this would be an acceptable assumption for the weather-loach data, due to the use of constant-effort electrofishing. Table 4 presents parameter estimates and estimated standard errors for the models (φ, p t , U t ) and (φ, p t , λ, η). We note that the average of the estimates of recapture probability from the N-mixture model is 0.19, in comparison with the value of 0.18 in Table 2, and conversely, from Table 2 the average estimate of φ is given by 0.65, in comparison with the value of 0.63 from Table 4. Standard errors are presented from using an estimated Hessian in the usual way, and are very similar to those obtained from a parametric bootstrap approach (not shown). The estimates of U from the Viterbi analysis of the N-mixture model are rounded due to the binning of states. They are in agreement with the predictions from equation (11), as they should be. As observed earlier, the Student estimates of U differ in interpretation compared with those obtained from the N-mixture model. The differences are generally small in comparison with the standard errors. The Student models generally are hard to fit, being subject to confounding, boundary estimates and convergence to local optima. For example, it is not possible to estimate p 1 and U 1 separately, nor φ 10 and p 11 in the (φ t , p t , U t ) model, which can only be estimated as a product, as observed also by HWK. However, they are estimable separately for model (φ t , p t , λ, η). We note in Table 4 the much smaller standard errors for the components of {U} from using the N-mixture approach, compared with the Student model, as expected, due to the chaining of the components of U in the open N-mixture model; standard errors are formed by using a parametric bootstrap, based on simulating HMMs. We can see that what is driving the need for the recapture probabilities to change with time is the estimate for the 7th sample occasion. It is interesting to note from Table 4 that the Dail/Madsen estimates of error are slightly larger than those from the Viterbi algorithm. In addition the Viterbi approach is more general. The very good agreement between the Viterbi and Dail/Madsen estimates, obtained by quite different methods, provides a validation for the methods used.
We show in Table 5 how we can estimate {N t }, in this case making use of the Viterbi estimates of the numbers of unmarked individuals for illustration. We see that the standard errors for N t and the population size estimates using "HT with p-vals," resulting from using the Horvitz-Thompson-like approach of equation (10) with estimated time-varying recapture probabilities, are similar in size, and generally different from those resulting from HWK. A remarkable feature of Table 5 is the very close agreements of both the estimates of N t and their corresponding estimated standard errors, whether one forms the estimates by evaluating the marked and unmarked components of N t separately, as in this article, or simply uses a HT approach. We note also the curious feature that the estimated standard errors are close to the sums Figure 1 illustrates goodness-of-fit plots for model (φ, p t , λ, η), with w = 20 and density dependence, separately for each likelihood component, as recommended by Besbeas and Morgan (2014). Overall there does not appear to be a serious lack of fit for the models fitted. The patterns in Figure 1a are a simple consequence of the repeated values of observed counts.
We note finally that the expected number of new unmarked individuals entering the population at occasion t can be derived usingB as in the Jolly-Seber model; see McCrea and Morgan (2014, p. 150). We would expectB t ≈ E[C t |U t−1 ], whereB t is estimated in equation (14), and this relationship is found to hold.

Figure 1.
Goodness-of-fit plots for the weather-loach data, for model (φ, p t , λ, η): (a) we plot the difference, observed minus expected recapture counts for the marked individuals, divided by the square root of the expected counts, against the logarithm of the expected counts, and (b) the QQ plot for the normal pseudo-residual segments of the unmarked counts.

Discussion
Batch marking is fundamentally important in ecology, providing an important tool for studying population dynamics. However, there is a pressing need for effective new methods of analysis. We have developed a comprehensive model for the extended batch-marking experiment, using an efficient hidden Markov formulation that supersedes previous analyses. The approach incorporates information on unmarked individuals, leading to increased precision of parameter estimates, and improved model-selection. It has been interesting to see, in the case study, that the selected model changes as a result of the joint analysis. In addition one can devise and compare different models for the unmarked individuals, and useful descriptions of the population result. For instance, we have examined whether constraining the population sizes for the unmarked individuals in the Student model to be constant over time results in a realistic model for the weather-loach data, and it was not found to be competitive. An added advantage of using HMMs is the ability to examine pseudo residuals to check goodness-of-fit. In addition the Viterbi algorithm provides a convenient, general method for estimating numbers of unmarked individuals, which can be bootstrapped to provide estimates of standard errors. As the model is fitted using standard likelihood optimization methods, it is possible to easily undertake model selection and compute estimates of parameter uncertainty. The models presented provide a flexible platform for additional development, depending on what data are available. For instance, one can incorporate suitable covariates, which would reduce the number of parameters to be estimated and potentially increase biological understanding. Extensions to cope with age and size information, as discussed by HWK, would also be straightforward. Further, we know from discussions with batch marking users that tag loss can be non-negligible, though apparently this was not an issue with the weatherloach study. For some study designs, it might be possible to implement double marks on individuals and implement methods similar to Cowen and Schwarz (2006). Alternatively, one might be able to obtain auxiliary information on tag-loss rates through a holding study and adjust parameter estimates using methods of Seber and Felton (1981).

Supplementary Materials
Matlab code for the analyses of the article is available with this article at the Biometrics website on Wiley Online Library. It is written for use in the Parallel Computing toolbox, but is readily translated into R if necessary.