Improved methods for estimating abundance and related demographic parameters from mark‐resight data

Abstract Over the past decade, there has been much methodological development for the estimation of abundance and related demographic parameters using mark‐resight data. Often viewed as a less‐invasive and less‐expensive alternative to conventional mark recapture, mark‐resight methods jointly model marked individual encounters and counts of unmarked individuals, and recent extensions accommodate common challenges associated with imperfect detection. When these challenges include both individual detection heterogeneity and an unknown marked sample size, we demonstrate several deficiencies associated with the most widely used mark‐resight models currently implemented in the popular capture‐recapture freeware Program MARK. We propose a composite likelihood solution based on a zero‐inflated Poisson log‐normal model and find the performance of this new estimator to be superior in terms of bias and confidence interval coverage. Under Pollock's robust design, we also extend the models to accommodate individual‐level random effects across sampling occasions as a potentially more realistic alternative to models that assume independence. As a motivating example, we revisit a previous analysis of mark‐resight data for the New Zealand Robin (Petroica australis) and compare inferences from the proposed estimators. For the all‐too‐common situation where encounter rates are low, individual detection heterogeneity is non‐negligible, and the number of marked individuals is unknown, we recommend practitioners use the zero‐inflated Poisson log‐normal mark‐resight estimator as now implemented in Program MARK.

aerial (eg, Keech et al., 2011), and camera (eg, Alonso et al., 2015) surveys. Because mark-resight does not require repeated physical capture of individuals, the main advantage of this approach is that costs associated with marking and recapturing can be minimized. In many circumstances, it can therefore be a less-invasive and lessexpensive alternative to traditional mark recapture for population monitoring. Relative to capture recapture, mark-resight abundance estimators also tend to be more robust to individual heterogeneity in detection probability because marked individuals are selected independent of the sighting process (ie, "all-zero" detection histories are observed).
Mark-resight estimators have seen much development over the past decade (McClintock and White, 2012;Matechou et al., 2013;Alonso et al., 2015;Lyons et al., 2016), and the models that are currently implemented in the freeware Program MARK (White and Burnham, 1999) are widely used by practitioners. Because it does not require knowledge of the exact number of marked individuals in the population, one of the most popular mark-resight models for estimating abundance and survival parameters is the zero-truncated Poisson lognormal estimator (ztPNE; ). Here, we identify several potential deficiencies of the ztPNE currently implemented in Program MARK. Among them are its failure to utilize valuable information about the number of marks ever released into the population and, under Pollock's robust design, its conditioning on first sighting instead of first encounter (eg., capture). The ztPNE in Program MARK also relies on an ad-hoc adjustment term to account for marks that are sighted but not identified to individual identity (see McClintock et al., 2014). Lastly, while the ztPNE implemented in Program MARK can include individual-level random effects in detection probability, these individual heterogeneity effects are (perhaps unrealistically) assumed to be independent across primary sampling periods.
As a solution to these potential deficiencies of ztPNE, we propose the zero-inflated Poisson lognormal estimator (ziPNE). By allowing for zero inflation in the marked individual sighting data, ziPNE can incorporate all marked individuals ever released into the population and thus does not need to condition on first sighting. We also extend the likelihoods for these models to accommodate unidentified marks and individual-level random effects across primary sampling occasions. After describing ziPNE in the next section, we compare the performance of ziPNE relative to ztPNE under a broad range of simulated sampling scenarios. We then revisit the New Zealand Robin (Petroica australis) analysis of  and compare inferences from the various estimators.

| METHODS
There are several different study designs for the markresight models in Program MARK (McClintock and White, 2012). The mark-resight scenario under consideration consists of T sampling occasions where sightings of both marked and unmarked individuals are recorded. The observable population size N ( ) t is defined as the number of individuals that are available to be observed (eg, alive and within the study area) during sighting survey t T = 1,…, . The population is assumed to be geographically and demographically closed within each survey, but not between surveys. To account for temporary emigration off the study area (or any other process that causes a segment of the population to be temporarily unobservable), we define the unobservable population size as the number of individuals that were alive during survey t but had temporarily emigrated off the study area between times t − 1 and t. During is the total number of marked individuals ever released into the population prior to survey t, and R t is the number of (unmarked) individuals that were newly marked between times t − 1 and t. The mark-resight data therefore consist of the number of resightings for each marked individual ) and the total number of unmarked individual sightings ( )  is the set of M t marked individuals ever released into the population prior to survey t. In the common circumstance where the identity of an individual detected as marked cannot always be determined (McClintock et al., 2014), then the observed data for the marked segment of the population consist of y , for s t  ∈ and the total number of marked sightings that The number of marks alive and within the study area can sometimes be reasonably assumed to be known exactly (eg, if all n t marks are introduced immediately prior to sighting survey t), but more commonly n t is unknown due to attrition of marks over time from mortality and permanent emigration.
When n t is unknown for at least one sampling occasion, the zero-truncated Poisson log-normal estimator (ztPNE, McClintock et al., 2009; can be used to estimate abundance N ( ) t , apparent survival from time t to time t ϕ + 1 ( ) t , the transition probability from an observable state (eg. in the study area) at time t to an unobservable state (eg, temporarily off the study area) at time t γ + 1 ( ″) t , and the transition probability from an unobservable state at time t to an unobservable state at time t γ + 1 ( ′) t . However, the ztPNE does not capitalize on all of the available information afforded by this mark-resight study design. Notably, the model does not incorporate M t into estimation of detection probability parameters, and, under Pollock's robust design (Kendall et al., 1997;, it conditions on the first sighting of each marked individual (instead of first encounter; eg, capture) for estimating the open population parameters (ϕ γ , ′, and γ′′). To make better use of the available information, we propose the zero-inflated Poisson log-normal estimator (ziPNE). The assumptions of ziPNE are identical to the standard mark-resight assumptions of ztPNE ). Unmodeled individual sighting heterogeneity is well known to induce bias in abundance estimators, and we develop two ziPNE formulations based on whether or not individual-level effects for sighting rates are assumed independent ("within" heterogeneity sensu McClintock et al., 2009) or not ("across" heterogeneity sensu Gimenez and Choquet, 2010) and r t is the conditional probability that a mark is individually identified (given detected as marked). The individual sighting rate λ ( ) . When the latent individual-level effects include time dependence, we refer to this as the "within" heterogeneity model. For an "across" heterogeneity model, we simply remove time dependence from the individual effects: Under this scenario, likelihood-based estimation of N t and related demographic parameters is challenging for several reasons. The individual-level effects clearly induce dependencies between y Y = ( ) and e e e = ( ,…, ) T 1 , and, for the "across" heterogeneity models, within Y e , , and u u u = ( ,…, ) T 1 . A further complication for the "across" heterogeneity models is that when R > 0 t for t T {2,…, } ∈ , then y e , s t T t T , : : , and u t 1: −1 may not be independent for s t−1  ∉ . One possible solution is to formulate a complete data likelihood and use data augmentation (eg, McClintock et al., 2014;Rankin et al., 2016), but such approaches can be time-consuming and difficult to fit (see Web Appendix A). We instead adopt a practical perspective and focus our efforts on a composite likelihood (eg, Varin et al., 2011) that is amenable to maximum likelihood estimation and implementation in Program MARK. Web Table B1 includes definitions for all parameters, latent variables, and random variables in our composite likelihood.
To make the problem tractable, we propose to ignore any covariance between Y e , , and u and simply integrate over the latent Z for each of three components of a composite observed data likelihood as if they were independent: where bracket notation data parameter [ ] | denotes a probability density (or mass) function. We can account for zero inflation in y s by assuming marked individuals are in one of three states at time t (1 = dead or permanent emigrant; 2 = alive, not a permanent emigrant, off the study area; 3 = alive, not a permanent emigrant, on the study area). Let t t be a row vector of initial state probabilities, where w t is the probability that a randomly selected individual from the R t new marks introduced between sighting occasion t − 1 and t was alive and not a permanent emigrant during sighting occasion t, and g t is the conditional probability that this individual was within the study area (given alive and not a permanent emigrant). Let P s t , be a 3 × 3 matrix where the diagonal elements denote the probability mass function of y s t , conditional on the state of individual s at time t, and let G t be a 3 × 3 state transition probability matrix with rows corresponding to the state at time t and the columns corresponding to the state at time t + 1. We then have α σ ϕ γ γ y z r w g w P G P 1 [ , , , , , , , is the first sighting occasion for which individual s was marked (and thus potentially available for resighting), I () is the indicator function, 1 = [1, 1, 1] 3 T , and , ,

|
. By allowing y s t , for the R t newly marked individuals to be zero inflated (with probability w g 1 − t t ), we can include M T and condition on first encounter (instead of first sighting) in our composite likelihood.
For both the "within" and "across" models, we wish to integrate over the latent variables z s t , and z s , respectively, in Equation (1). While this is straightforward for pose difficulties because sums of Poisson log-normal random variables are not of standard form (eg, Poisson) and they involve multidimensional integrals that are not amenable to numerical approximation. Similar to ztPNE, we rely on a left-truncated normal distribution for of a Poisson lognormal random variable after integrating over the individual-level effects: represents a normal distribution truncated at a and b. Note that for the "across" heterogeneity model, this composite likelihood also ignores the covariance induced in u by the latent individual-level effects.
We employ a similar left-truncated normal distribution for α σ e r [ , , ] | : where e α σ r nλ r n κ κ λ r * , ,~( , and e e = * t t is a non-negative continuous random variable. While n 1 can sometimes be reasonably assumed to be known exactly (ie, n M = 1 1 and w g = = 1 1 1 ), for all occasions where the number of marks alive and within the study area is not known exactly, we replace n t in Equation (2) with a Horvitz-Thompson-type estimator: is the probability of a marked individual being resighted at least once.
We therefore have the following composite likelihood for the "within" heterogeneity ziPNE: , , , , , , , , . For the "across" heterogeneity ziPNE: , , , , , , , , , . For maximum likelihood analysis, the integrals in Equations (3) and (4) can be easily approximated using Gauss-Hermite quadrature: , , | , and f v ( , ) l l are the corresponding weights and nodes, respectively, for quadrature points l L = 1,…, (Abramowitz and Stegun, 1964). When n t is unknown, we have The size of the observable population (eg, alive and within the study area) during sighting survey can be derived from the estimated parameters after fitting Equations (3) or (4) using maximum likelihood methods. As with all derived parameters in Program MARK, standard errors for N t are estimated using the delta method and finitedifference approximations of the first derivative. Program MARK also facilitates modeling any of the parameters as a function of covariates using link functions (see Section 4). We note that in a fully time-dependent model for w t and g t T ( = 1,…, ) t , only the product w g T T is estimable. When T = 1, only the product w g 1 1 is estimable. We also note that because the "within" and "across" heterogeneity models use the same data, they can be compared using model selection criteria such as Akaike's Information Criterion (AIC; Burnham and Anderson, 2002).

| SIMULATION STUDY
Using data generated from the exact model, we performed a series of simulation experiments to evaluate the performance of our composite ziPNE relative to the ztPNE as , some of the U t−1 unmarked individuals became newly marked (with probability 0.02) in order to maintain a somewhat consistent marked population size over time. However, no additional unmarked individuals were introduced to the population after the first occasion, so N t is expected to steadily decline through time due to mortality and permanent emigration. Scenarios with λ = 100 t represent the "best case" scenario where sighting rates are very high and individuals that are alive and within the study area are virtually guaranteed to be sighted (i.e. x > 0 s t , ). Under this scenario, we would expect ziPNE and ztPNE to perform similarly. Scenarios with both λ = 1.5 t and τ = 5 represent the pathological case where mean sighting rates are low and individual heterogeneity is high, whereby the vast majority of individuals are never expected to be sighted (x = 0 s t , ) but some individuals will be sighted many times (x 0 s t , ≫ ). In terms of α σ , , and N estimation, we would generally expect all estimators to perform less than nominally under this extreme scenario (particularly when r = 0.5 t ). While we would not expect ϕ estimators to be biased by extreme heterogeneity coupled with low mean sighting rates, we suspect it could also be deleterious to reliable estimation of γ′′ and γ′.
For each combination, we simulated data from both the "within" and "across" heterogeneity models, thus yielding 2 × 2 × 2 × 2 × 3 × 4 = 192 simulated scenarios. For each scenario, we fit the two models corresponding to the data-generating heterogeneity model: 1) ziPNE within heterogeneity Equation (3) and ztPNE within heterogeneity; or 2) ziPNE across heterogeneity Equation (4) and ztPNE across heterogeneity. We investigated heterogeneity model mis-mispecification in an additional series of simulation experiments (see Web Appendix E). We assumed that n 1 was known, n t T ( = 2,…, ) t was unknown, and M n = 1 1 . For scenarios with τ = 0, we fixed σ = 0 t for each of the fitted models. For simplicity, we also fixed w g = = 1 1 1 (because n 1 is known and n M = 1 1 ) and otherwise assumed no time dependence in parameters (eg, ). All simulations used L = 101 quadrature points for numerical integration. We evaluated estimator performance based on bias, 95% confidence interval coverage of the true parameter values, and root mean squared error (RMSE). Simulation results are presented using the nested loop plot of Rücker and Schwarzer (2014). Similar to a time series plot, nested loop plots serve to present a large number of simulation results by putting all scenarios into a lexicographical order and arranging them consecutively along the horizontal axis. The quantity of interest (eg, bias, coverage, and RMSE) is then plotted on the vertical axis.
We found the "within" heterogeneity ziPNE to be far superior to ztPNE for the estimation of the closed population parameters (N α , , and σ), particularly as λ t decreased, τ increased, and r decreased. In nearly all of the scenarios, ziPNE tended to be less biased with higher 95% confidence interval coverage and lower RMSE (Figure 1, Web Figures C1-C4). Although it performed much better than ztPNE, ziPNE still exhibited an increase in bias and decrease in coverage under the pathological scenarios with low mean resighting rates and extreme heterogeneity. While little difference was found between ziPNE and ztPNE in terms of apparent survival ϕ ( ) estimation (Figure 2), ziPNE performed much better than ztPNE in the estimation of both γ′′ and γ′ (Web Figures C5-C8).
We did not find such drastic differences in the performance of the "across" heterogeneity ziPNE relative to ztPNE, but ziPNE did tend to perform better in the estimation of the closed population parameters (Figure 1, Web Figures C2-C4). While the "across" ztPNE tended to exhibit slightly less bias and lower RMSE for N in the most extreme scenarios, the "across" ziPNE tended to yield higher 95% confidence interval coverage (particularly when r = 0.5 t ; Figure 1, Web Figure C1). Nevertheless, coverage MCCLINTOCK ET AL.

| 803
of N still tended to be less than nominal when τ > 0 and r = 0.5 t for all λ t . For the open parameters ϕ γ γ ( , ′′, ′), we found that both ziPNE and ztPNE performed equally well in terms of bias and confidence interval coverage (Figure 2), but RMSE for ϕ and γ′ in particular tended to be larger for ztPNE (Web Figures C5-C8).

| EXAMPLE: NEW ZEALAND ROBIN
As an example application, we return to the New Zealand Robin (P. australis) data originally analyzed using ztPNE by . This mark-resight study was conducted by the New Zealand Department of t estimators based on 192 simulated scenarios. Scenarios are ordered from outer to inner loops by λ {100, 5, 1.5} t ∈ ("lambda"; three levels, decreasing), τ {0, 0.5, 1, 5} ∈ ("tau"; four levels, increasing), r {0.9, 0.5} ∈ (two levels, decreasing), R E( ) {120, 60} 1 ∈ ("Number marked"; two levels, decreasing), and U E( ) {600, 300} 1 ∈ ("Number unmarked"; two levels, decreasing) such that the left-and right-most scenarios are the best-and worst-case sampling scenarios, respectively. Estimators include the zero-inflated Poisson log-normal estmator (dark orange = ziPNE within, light orange = ziPNE across) and the zero-truncated Poisson log-normal estimator (dark purple = ztPNE within, light purple = ztPNE across). Results are based on 400 simulated data sets for each scenario. This figure appears in color in the electronic version of this article, and color refers to that version Conservation on three 100 ha study areas ("Knobs Flat", estimators based on 192 simulated scenarios. Scenarios are ordered from outer to inner loops by λ {100, 5, 1.5} t ∈ ("lambda"; three levels, decreasing), τ {0, 0.5, 1, 5} ∈ ("tau"; four levels, increasing), r {0.9, 0.5} ∈ (two levels, decreasing), R E( ) {120, 60} 1 ∈ ("Number marked"; two levels, decreasing), and U E( ) {600, 300} 1 ∈ ("Number unmarked"; two levels, decreasing) such that the left-and right-most scenarios are the best-and worst-case sampling scenarios, respectively. Estimators include the zero-inflated Poisson log-normal estmator (dark orange = ziPNE within, light orange = ziPNE across) and the zero-truncated Poisson log-normal estimator (dark purple = ztPNE within, light purple = ztPNE across). Results are based on 400 simulated data sets for each scenario. This figure appears in color in the electronic version of this article, and color refers to that version within a study area. Sighting surveys did not begin in Smithy until August 2005, so there were T = 6 sighting occasions for both Knobs Flat and Walker Creek, but only T = 5 occasions for Smithy. Notably, ship rat (Rattus rattus) control was conducted on the Walker Creek study area from June 2006 to April 2007. Complete details of the study design and data can be found in .
We performed a ziPNE analysis for comparison to the "within" ztPNE analysis of McClintock and White (2009)  Covariates examined for sighting rate and open population parameters included age class (age), sex, age class at time of capture (agecap), study area (area), time, breeding season (season), and rat control on Walker Creek. As in McClintock and White (2009), we also included an interaction between age, sex, and breeding season (age*sex*season) to investigate potential effects of reproduction on adult survival. For model selection, we used AIC adjusted for small sample sizes (AICc; Burnham and Anderson, 2002). Three covariate models for rat control on Walker Creek were examined: (1) an immediate and constant effect during the entire program (rat); (2) an immediate effect with no effect thereafter (rat1); and (3) no immediate effect with a delayed effect thereafter (rat2). We also investigated models under completely random emigration, γ γ ′′ = ′ (Kendall et al., 1997), and models ignoring individual heterogeneity σ ( = 0) or temporary emigration γ γ ( ′′ = ′ = 0). Similar to McClintock and White (2009), we first investigated various models for the sighting rate parameters (α σ , , and r) under the most general (but identifiable) time-and group-dependent structure for U and the open parameters (ϕ γ , ′′, and γ′). We proceeded by selecting the model structure for α σ , , and r best supported by AICc and then using only this structure for investigating more parsimonious models for ϕ γ , ′′, and γ′ (for an evaluation of a similar model selection procedure, see Doherty et al., 2012). We then fit all combinations of the three best supported structures for each parameter based on each stage of this selection process. We assumed the same structure for w and ϕ and the same structure for g γ , ′′, and γ′. However, because they pertain to the entire capture period beginning in September 2003, we also investigated w 1 and g 1 effects (cap1) for the Knobs Flat and Walker Creek study areas, and w 2 and g 2 effects (cap2) for the Smithy study area. Unlike the "within" ztPNE analysis of McClintock and White (2009), we investigated both "within" and "across" heterogeneity models, and we again note that, unlike the ziPNE proposed here, the ztPNE of McClintock and White (2009) does not include the r parameter and instead relies on an ad-hoc adjustment term for unidentified marks. When fitting the models, we found the Markovian temporary emigration models including additive age, agecap, or area effects to have suspect convergence and thus excluded these from our model set. We based our comparisons on model-averaged parameter estimates and unconditional variances calculated from AICc weights.
AICc overwhelmingly supported "across" heterogeneity models for ziPNE, with the minimum AICc model including: agecap, sex, season, and area effects on ϕ and w; cap1 and cap2 effects on w; agecap and season effects on random emigration (ie, γ γ ′′ = ′) and g; age, sex, and rat effects on α; sex effects on σ; and area by time interaction effects on r. Model selection was somewhat similar for the "within" ztPNE model of , but the ziPNE analysis exhibited considerably more model selection uncertainty (Web Table D1) and weak evidence of rat control effects on ϕ.
For ziPNE, model-averaged mean sighting rates λ ( ) t were the lowest in August 2007 (2.2, 95% CI: 1.9 − 2.5) and the highest on Walker Creek during rat control in March 2007 (3.2, 95% CI: 2.6 − 4). For the "across" ziPNE models, model-averaged σ = 0.52 (SE = 0.18) for females, and σ = 0.28 (SE = 0.07) for males. Individual mark identification probabilities r ( ) varied by study area and time but were generally high, with modelaveraged estimates ranging from 0.81 (SE = 0.07) to 1.00 (SE = 0.00). Point and variance estimates for N t tended to be larger for ziPNE, and, except for Walker Creek in August 2006, all estimated 95% confidence intervals for N t were greater than the minimum number known alive based on territory mapping that was independent of the mark-resight methodology ( Figure 3). As with the original ztPNE analysis, the ziPNE analysis found evidence of higher apparent survival probabilities for individuals that were first captured as adults (Figure 4) when compared to individuals that were first captured as juveniles (Web Figure D1), but ziPNE exhibited less seasonal variation in adult female apparent survival relative to ztPNE. Temporary emigration parameter (γ′′ and γ′) point estimates were similar between ziPNE and ztPNE, but ziPNE was generally more precise (Web Figures D2-D4). Estimates for w after longer capture intervals tended to be the lowest, particularly for individuals that were first captured as juveniles. Estimates for g were generally high and tended to be somewhat higher for juveniles (Web Figures D5-D6).
We have proposed a new zero-inflated version of the Poisson log-normal mark-resight estimator in Program MARK. Using simulation, we demonstrated ziPNE has superior performance relative to the zero-truncated model of . By accounting for zero inflation in the marked individual sighting data, ziPNE is able to utilize valuable information about detection rates provided by marked individuals that were released but never subsequently sighted. This not only makes ziPNE more robust to individual heterogeneity for reliable estimation of population size but also allows the likelihood for the open population parameters (ϕ γ γ , ′′, ′) to condition on first encounter (instead of first sighting). ziPNE thus makes better use of available data and, unlike ztPNE, allows γ′ 1 to be estimated under Pollock's robust design.
In terms of abundance estimation, previous simulation studies also found the "within" heterogenity zerotruncated model of  to perform increasingly poorly with low sighting rates, high levels of individual heterogenity, and low-marked individual identification probabilities McClintock et al., 2014). Our extensive simulations included a broader range of more extreme scenarios, and the new ziPNE performed near nominally and better than ztPNE in almost every scenario examined. However, the average coverage of N could still be as low as 81% for the "across" ziPNE when high levels of individual heterogeneity were combined with low mark identification probabilities. While we have not encountered actual mark-resight data as information poor as some of the scenarios included in our simulations, these pathological sampling scenarios provide valuable insights about the limitations of our composite likelihoods. Practitioners should carefully consider study designs that increase sighting rates, decrease individual heterogeneity, and increase marked individual identification whenever possible.
Unless changes in individual behavior (or sighting survey design) effectively eliminates correlations in the marked individual sighting data, the "across" model may often be a more natural choice than the "within" model. Indeed, our New Zealand Robin example demonstrated overwhelming support for the "across" models and thus strong correlations in the individual sighting data across sighting surveys. It deserves note that both the "across" and "within" models tended to perform poorly under heterogeneity model mis-specification (see Web Appendix E Creek study area. This figure appears in color in the electronic version of this article, and color refers to that version therefore important to examine both "within" and "across" heterogeneity models. Fortunately, these models can be compared using standard model selection criteria such as AIC and Bayesian information criterion (BIC). In some cases, it may be desirable to simultaneously model both "within" and "across" heterogeneity; this extension remains an avenue for future research. Our ziPNE is a composite likelihood that ignores some of the dependencies induced by individual heterogeneity among Y e , , and u in order to facilitate maximum likelihood estimation in Program MARK. While this composite likelihood proved robust under a broad range of simulation scenarios, an alternative (exact) solution can be implemented using a data-augmented complete data likelihood in a Bayesian framework (see Web Appendix A). While relatively easy to code using standard software, poor mixing can cause such implementations to be time consuming and challenging to fit. Thus, while conceptually appealing, they can often be impractical for simulation studies or real-world application (particularly for large T or N ). To aid practitioners in determining whether or not the performance of ziPNE can be expected to meet their specific objectives under sampling scenarios not covered in our simulations, we provide R code for simulating data from the exact model and fitting ziPNE in the Supporting Information.
We recommend practitioners use ziPNE for markresight analyses in Program MARK whenever the number of marks in the population is unknown due to attrition arising from mortality and permanent emigration. As demonstrated in our New Zealand Robin example, we also recommend previous ztPNE analyses be revisited with ziPNE where appropriate. Program MARK includes both individual heterogeneity model formulations ("within" and "across"), and has many additional analysis features, such as covariate modeling via link functions, model selection using AIC or BIC, multimodel inference, and plotting. MARK also provides the ability to simulate data under the mark-resight models, a useful feature for designing studies. For advanced users, these mark-resight models can also be implemented through the R package RMark (Laake, 2013) interface for Program MARK.

ACKNOWLEDGMENTS
The findings and conclusions in the manuscript are those of the author(s) and do not necessarily represent the views of the National Marine Fisheries Service, NOAA. Any use of trade, product, or firm names does not imply an endorsement by the US Government.