Incorporating capture heterogeneity in the estimation of autoregressive coefficients of animal population dynamics using capture–recapture data

Abstract Population dynamic models combine density dependence and environmental effects. Ignoring sampling uncertainty might lead to biased estimation of the strength of density dependence. This is typically addressed using state‐space model approaches, which integrate sampling error and population process estimates. Such models seldom include an explicit link between the sampling procedures and the true abundance, which is common in capture–recapture settings. However, many of the models proposed to estimate abundance in the presence of capture heterogeneity lead to incomplete likelihood functions and cannot be straightforwardly included in state‐space models. We assessed the importance of estimating sampling error explicitly by taking an intermediate approach between ignoring uncertainty in abundance estimates and fully specified state‐space models for density‐dependence estimation based on autoregressive processes. First, we estimated individual capture probabilities based on a heterogeneity model for a closed population, using a conditional multinomial likelihood, followed by a Horvitz–Thompson estimate for abundance. Second, we estimated coefficients of autoregressive models for the log abundance. Inference was performed using the methodology of integrated nested Laplace approximation (INLA). We performed an extensive simulation study to compare our approach with estimates disregarding capture history information, and using R‐package VGAM, for different parameter specifications. The methods were then applied to a real data set of gray‐sided voles Myodes rufocanus from Northern Norway. We found that density‐dependence estimation was improved when explicitly modeling sampling error in scenarios with low process variances, in which differences in coverage reached up to 8% in estimating the coefficients of the autoregressive processes. In this case, the bias also increased assuming a Poisson distribution in the observational model. For high process variances, the differences between methods were small and it appeared less important to model heterogeneity.

models seldom include an explicit link between the sampling procedures and the true abundance, which is common in capture-recapture settings. However, many of the models proposed to estimate abundance in the presence of capture heterogeneity lead to incomplete likelihood functions and cannot be straightforwardly included in state-space models. We assessed the importance of estimating sampling error explicitly by taking an intermediate approach between ignoring uncertainty in abundance estimates and fully specified state-space models for density-dependence estimation based on autoregressive processes. First, we estimated individual capture probabilities based on a heterogeneity model for a closed population, using a conditional multinomial likelihood, followed by a Horvitz-Thompson estimate for abundance.
Second, we estimated coefficients of autoregressive models for the log abundance.
Inference was performed using the methodology of integrated nested Laplace approximation (INLA). We performed an extensive simulation study to compare our approach with estimates disregarding capture history information, and using R-package VGAM, for different parameter specifications. The methods were then applied to a real data set of gray-sided voles Myodes rufocanus from Northern Norway. We found that density-dependence estimation was improved when explicitly modeling sampling error in scenarios with low process variances, in which differences in coverage reached up to 8% in estimating the coefficients of the autoregressive processes. In this case, the bias also increased assuming a Poisson distribution in the observational model. For high process variances, the differences between methods were small and it appeared less important to model heterogeneity.

| INTRODUC TI ON
Models used to analyze population dynamics include a combination of density dependence and environmental effects. Ignoring the uncertainty in abundance estimates biases estimates of the strength of density dependence, and different approaches exist to achieve better accuracy (see Lebreton & Gimenez, 2012 for a review). In particular, state-space models combining an observation model-linking the observations such as counts to the true abundance-and a process model-describing the processes driving population dynamics-have become a standard approach in many analyses (Dennis & Taper, 1994). However, these models rarely include an explicit model of the link between how counts were obtained and true abundance, often relying on a nonspecific observation model such as log normal or Poisson distribution (for instance, Ono, Langangen, & Chr. Stenseth, 2019, but see below).
Capture-recapture methods have been extensively used to estimate abundance and density dependence, and many methods have been developed to incorporate different sources of variability into capture probability estimation, such as environmental information, survival, or trophic interactions (Barker, Fletcher, & Scofield, 2002;Lebreton & Gimenez, 2012;Schofield & Barker, 2008;Yackulic, Korman, Yard, & Dzul, 2018). Estimating abundance is a challenging statistical problem (Link, 2003), and heterogeneity in capture probabilities can lead to large biases in abundance estimates when using models assuming no heterogeneity (Carothers, 1973;Otis, Burnham, White, & Anderson, 1978). However, many of the models that have been proposed to estimate abundance in the presence of heterogeneity do not lead to observation models that can be included in state-space models as they do not lead to likelihood functions in a closed form (Chao & Huggins, 2006;Huggins & Hwang, 2011).
Many studies investigating density dependence have used simple process models such as the Gompertz model-that is, a model which is a first-order autoregressive model on a log scale (Ono et al., 2019;Thibaut & Connolly, 2019). However, ecological processes such as trophic interactions (Bjørnstad, Falck, & Stenseth, 1995) or intrinsic ecological properties such as age structure (Lande, Engen, & Saether, 2002) may lead to more complex process models such as a second-order autoregressive model (AR (2)). An important case is the population cycles observed in many small mammal populations, particularly in northern environments (Bjørnstad & Chr. Stenseth, 1999;Elton, 1924;Stenseth, 1999). These quasi-periodic fluctuations are quite well approximated by AR(2) models on a logarithmic scale (Bjørnstad et al., 1995). Whereas most analyses have ignored the uncertainty in abundance estimates (Bjørnstad et al., 1995), some have used state-space models (Cornulier et al., 2013;Ims, Yoccoz, & Killengreen, 2011;Kleiven, Henden, Ims, & Yoccoz, 2018;Stenseth et al., 2003). However, approaches using a capture-recapture framework and including capture heterogeneity have relied on integrating out random effects describing capture heterogeneity (King, Brooks, & Coulson, 2008;Schofield & Barker, 2013) and using superpopulation data augmentation (Royle, 2008); these approaches did not consider the conditional likelihood approach to estimating population size, which can easily handle, for example, individual covariates (Huggins & Hwang, 2011). Moreover, a fully MCMCbased Bayesian approach is computationally intensive on large data sets and requires that careful considerations are given to choices of priors and superpopulation sizes for data augmentation (Royle, Dorazio, & Link, 2007).
Here, we investigated the performance of an intermediate approach between ignoring uncertainty in abundance estimates (i.e., using the raw population counts) and fully specified state-space models. Specifically, we first used a multinomial observation model to estimate capture probabilities followed by estimating abundance at each time point using the Horvitz-Thompson estimator (Horvitz & Thompson, 1952). Second, we fitted an AR(2) process model to the log abundance to estimate direct and delayed density dependence given by the first and second coefficients of the AR(2) model, respectively. Both estimation steps were performed in a unified way, incorporating the models within the general class of latent Gaussian models (Rue, Martino, & Chopin, 2009). Full Bayesian inference was then obtained using the methodology of integrated nested Laplace approximation (INLA) (Rue et al., 2009. We based our analyses on a large-scale study of population dynamics of the dominant small mammal species in northern Fennoscandia, the gray-sided vole Myodes rufocanus (Ims et al., 2011). This species shows large fluctuations with a 4-to 5-year periodicity (Ims et al., 2011;Marolla et al., 2019). We monitored populations of gray-sided voles along a 155-km gradient from coast to inland, using live capture-recapture methods, starting in 2000.
Previous analyses have shown that there was large heterogeneity in capture probabilities (Yoccoz & Ims, 2004). In this rodent study, the goal was to understand spatial patterns of population dynamics, assessing potential seasonal effects on the density-dependence estimates. For this, we first needed to assess the robustness of using an approach based on estimated abundances but without implementing a full state-space model. In this paper, we therefore use a simulation study built around the case study (adaptable to other situations from the code provided) to assess the estimation accuracy of the density dependence, both including and excluding capture history information.
The structure of this paper is as follows. Section 2 provides our methodological background to analyze capture-recapture data and describes the Bayesian framework to perform parameter estimation. This includes using INLA to estimate individual capture probabilities and the direct and delayed density dependence given by the K E Y W O R D S abundance, capture probability, closed population models, density dependence, INLA, process variance coefficients of AR(2) models. Section 3 contains an extensive simulation study, investigating how density-dependence estimates are influenced when individual capture probabilities are taken into account.
In Section 4, we study the population cycles of gray-sided voles. We first compare different observation models in estimating individual capture probabilities and then assess whether incorporation of individual capture probabilities influences density-dependence estimates. A summary and concluding remarks are given in Section 5.

| ME THODOLOGY
Capture-recapture experiments are important to assess heterogeneity in individual capture probabilities. This section describes our approach to incorporate capture-recapture information in the estimation of density dependence. First, we define an observation model in which capture probabilities are modeled in terms of individual features and then used to estimate abundance. Second, we fit an AR(2) process model to the estimated log abundance to assess density dependence. When using state-space approaches, the parameters of the observation and process model are estimated simultaneously. This is not possible in our case as the capture probabilities are estimated based on a conditional multinomial likelihood, due to individuals that were not observed. Instead, we apply a sequential approach, first estimating the capture probabilities and then the AR(2) coefficients. This allows us to use an explicit sampling model to estimate capture probabilities, instead of assuming that the observed counts have a Poisson or log normal distribution. The given sequential approach is computationally efficient using the R-INLA package which is freely available at www.r-inla.org.

| Statistical background on capturerecapture data
Assume a closed population with a total of N individuals and a capture-recapture experiment with τ capture sessions. Let denote the capture history for the ith individual. If w ij = 1, the individual was captured at the jth capture session, while w ij = 0 otherwise, that is, w ij ~ Bernoulli(p ij ), j = 1, …, τ. For each individual, the probability of a given capture history is then Assuming that all individuals are captured independently, the complete likelihood becomes where both N and the set of probabilities {p ij } are unknown. Due to the unknown number of noncaptured individuals, computation of the likelihood is unfeasible. This is a well-known problem (Huggins & Hwang, 2011) and requires alternative strategies to perform parameter estimation.
A commonly applied approach is to maximize the conditional likelihood for the n individuals that were captured at least once. Let c ik , k = 0, …, 2 τ − 1, denotes the probability that the capture history of individual i is equal to category k. The different categories are defined by all possible permutations of the capture session vector, giving a total of m = 2 τ − 1 categories for the captured individuals.
From here onwards, we will refer to data sets with only two capture events, in which mortality and emigration are disregarded considering capture events on adjacent days. The event that an individual is never captured is then defined as category 0, while the categories 1, 2, and 3 refer to the capture histories (1,0), (0,1), and (1,1), respectively. To perform parameter estimation, we need to make realistic assumptions on the capture probabilities for different capture sessions. Otis et al. (1978) propose a total of eight different models characterizing capture probabilities for different sessions depending on time, behavior, and homogeneity of the individuals, also including combinations of these three factors. Here, we consider a heterogeneity model including a temporal effect, Mth. This implies that the capture probabilities depend on different features of the individuals.
Further, we assume that the capture probability on the first and second capture sessions is independent. The probabilities for the different categories are then specified as To estimate abundance based on individuals that were captured, we use the Horvitz-Thompson estimator (Horvitz & Thompson, 1952) where ĉ i0 denotes the estimated probability that individual i was not captured. This probability is estimated using a regression model as explained in the next section.

| A multinomial capture-recapture regression model including a Poisson transformation
An important question in analyzing population processes from capture-recapture data is whether features of the captured individuals give valuable information in further analysis of density dependence.
To estimate the probabilities in (2), it is natural to assume a multinomial regression model for the captured individuals, incorporating covariate information which helps to separate different capture categories. We define the vector Y ′ i = (Y i1 , …, Y im ), where Y ik = 1 for an individual classified to category k, while the remaining elements of Y i are 0. Each of the vectors Y 1 , …, Y n has a multinomial distribution.
Based on (1), probabilities for the m = 3 observed categories are defined by c = c ∕ 1 − c i0 , k = 1, …, m, ensuring that the probabilities sum to 1. These probabilities can then be modeled in terms of observed individual features such as weight, sex, and age.
We denote the individual features or covariates by z � r = z 1r , …, z . Further, we define the linear predictor where the coefficient γ kr is specific for category k and covariate r, while v is the number of covariates. The scaled probabilities for the captured individuals are then expressed as The resulting multinomial likelihood is Notice that in maximizing (5), the denominator of c does not simplify using the ordinary logarithmic transformation. It is therefore common to apply the well-known multinomial Poisson transformation (Baker, 1994) in which the likelihood is rewritten as Here, = e V + i represents the rate of a Poisson distributed random variable Y ik . The given transformation from a multinomial likelihood to the Poisson likelihood introduces auxiliary parameters β′ = (β 1 , . This is just a technical detail to make the approximation work correctly. The likelihood L P (.) is proportional to L M (.) and gives the same maximum-likelihood estimates for the coefficient vectors γ r . The resulting regression model is then summarized in terms of linking the expectation of the Poisson variables to the linear predictor using the log transform, that is, where i ∼ N 0, −1 denotes small independent random error terms.
In fitting the given model to a data set, the vectors { r } v r=1 will not be identifiable. However, in our case we only need estimates of the differences in these coefficients as these represent ratios of log probabilities between the different categories. For categories k and l, we notice that In estimating the parameters of the model, this implies that the auxiliary parameters and error terms disappear, but these are still included in fitting (6) to a data set. In the case of assuming (1), the estimated individual probabilities are then given by or equivalently These probabilities are then used to estimate ĉ i0 in (2).

| Implementation using a Bayesian framework
To fit (6) to a data set and estimate the capture probabilities, we choose to apply a Bayesian approach. This implies that all parameters in (6) are viewed as random variables. Specifically, the resulting regression model can be incorporated within the computational framework of latent Gaussian models. This is a flexible class of threestage hierarchical models, which can be analyzed in a unified way using INLA. Subsequently, the model in (6) is reformulated in terms of having conditionally independent observations, given a latent field and hyperparameters.
The three stages of a latent Gaussian model are expressed as follows, where π(.) is generic notation for probability densities: 1. The first stage specifies the likelihood where the observations are assumed conditionally independent given a latent field x and hyperparameters θ. In our case, let y � = (y � 1 , …, y � n ) denote the stacked vector of the m categories for the n individuals.
The likelihood is then expressed as 2. The latent field x collects all random variables of the linear where we could also include the predictor itself. The latent field models the dependency structure of the observations and is assigned a multivariate Gaussian prior The precision (inverse covariance) matrix Q is typically sparse such that x has Markov properties and is then referred to as a Gaussian Markov random field.
3. The hyperparameters θ of a latent Gaussian model are usually assigned non-Gaussian priors. Here, we only have one hyperparameter being the precision parameter of the random error terms, θ = κ. This parameter is assigned a penalized complexity prior (Simpson, Rue, Riebler, Martins, & Sørbye, 2017), implying that κ −1/2 has an exponential density.
The joint posterior for all elements of the latent field and the additional hyperparameter is then described as The main interest is to calculate the marginal posteriors for each of the latent field components and each of the hyperparameters.

| Estimating density dependence
Our final step is to fit a process model to study population dynamics of a species. Specifically, we focus on estimating density dependence by fitting an AR(2) model to a given time series, reflecting the population cycle for the relevant species. Let ln(N t ) denote the true log abundance at time t. The AR(2) model is then defined by where ln(η) denotes an offset, while the noise terms are independent Gaussian variables, t ∼ N(0, 2 ). T denotes the length of the time series, while the coefficients ϕ 1 and ϕ 2 characterize the direct and delayed density dependence of the series. The given process is stationary when −1 ≤ 2 ≤ 1 − | | 1 | | < 1 and has pseudoperiodic behavior when 2 1 + 4 2 ≤ 0. Estimation of the coefficients of AR (2) is not influenced by the offset ln(η). This implies that if the number of captured individuals at different time points is proportional to the underlying true abundance, we would get identical parameter estimates.
The AR(2) model is fitted within the framework of latent Gaussian models using INLA. In this case, the model has three hyperparameters, including = −2 and the coefficients ϕ 1 and ϕ 2 . These parameters are all assigned PC priors . Of main interest is to study how the estimates of ϕ 1 and ϕ 2 vary when capture heterogeneity is accounted for using the multinomial observational model.
Often, simplifying assumptions regarding the data generating process are made, for example, by assuming a Poisson process (Stenseth et al., 2003)  Here, β 0 denotes an intercept, while e 1 , …, e T denotes independent and identically distributed random variables, e i ∼ N(0, −1 e ). These error terms are included to model random variation as a function of time. As detailed in the next section, the AR(2) model will be fitted either to the estimated log abundance ln(N 1 ), …, ln(N T ) or to the posterior means of the log rates of the corresponding Poisson process, denoted ̂ 1 , …,̂ T .

| S IMUL ATION S TUDY COMPARING ME THODS TO E S TIMATE DEN S IT Y DEPENDEN CE
This section provides an extensive simulation study to assess how the inclusion of capture history information influences estimation of density dependence. We start by simulating data to approximate a realistic capture-recapture sampling scenario. The underlying log population of the sampled species is generated as an AR(2) process in time, using different fixed combinations of the coefficients (ϕ 1 , ϕ 2 ) and the innovation variance 2 , from here onwards referred to as the (population) process variance. Each resulting individual is then assigned a random weight, and a two-day capture history according to a multinomial model with probabilities defined by (1). We then fit an AR(2) process model to the estimates of log abundance or log rates obtained by different methods. These different methods are described in Section 3.1, while Section 3.2 specifies the simulation procedure and the method performance criteria used. Finally, Section 3.3 provides simulation results and an evaluation of the different methods.

| Estimation methods
An overview of the different estimation methods used in the simulation study is given in Figure 1. The left-hand side of the figure shows the additional steps needed to implement the observation model, incorporating sampling error in terms of capture history information. We employ two methods of estimating individual capture probabilities. The first is described in Sections 2.2 using INLA (method: CR-INLA) and corresponds to our suggested approach. The second, for comparison, estimates individual capture probabilities using the R-package VGAM (Yee, 2019). Among other utilities, the VGAM (vector generalized additive model) framework can be used to analyze closed population capture-recapture data, allowing the incorporation of individual covariates while using the conditional likelihood (Yee, Stoklosa, & Huggins, 2015). This application of VGAM allows for a flexible and efficient estimation of capture probabilities for all of the eight heterogeneity models given by Otis et al. (1978) (method: CR-VGAM). From the estimated capture probabilities from either of the two methods, we proceed to estimate the true log abundance using the Horvitz-Thompson estimator in (2). At this point, we have two possible variants in estimating density dependence: We either fit the AR(2) model to the times series of estimated log abundance {ln(N t )} T t=1 (A variant), or we fit the AR(2) model to the corresponding estimated log rate of a Poisson process, {̂ t } T t=1 (P variant). The right-hand side of Figure 1 illustrates the approach disregarding capture history, fitting the AR(2) model directly to the observed log counts, or to the log rate of the corresponding Poisson process (method: ObsCount).
Finally, the performance of the different estimation methods is compared with the results fitting the AR(2) model to the true generated log abundance or estimated log rate (method: Baseline).

| Simulation procedure
For each combination of AR (2) (12) where T = 20, using different fixed combinations of (ϕ 1 , ϕ 2 ). To remove the effect of sample size on the estimation of capture probability, we assumed that E(N t ) = 20 by using an offset ln( ) = ln(20) − 1 2 Var(ln(N t )). The series was rounded to give integer values for {N t } T t=1 , representing the abundance of an animal population. The total number of individuals generated for one simulated AR(2) process was then Ñ = ∑ T t=1 N t . 2. For each of the Ñ individuals, we generated a random weight where σ w = 1.2, while μ t ~ Log normal(ln(30), ln (5)). The weight was then scaled by the sample standard deviation of the generated weights to make it dimensionless. The resulting variable was used as an individual-specific covariate in (3). In this context, weight is a proxy for detectability. We varied the expected value of weight with time to model varying detectability, reflecting changes in the composition of the population at different time points. Thus, the varying mean reflects biological variation which we considered more realistic than assuming constant capture probabilities for different time points. The parameters relating to the weight distribution were here chosen to illustrate this biological variation.

AR(2) Population
3. Assume a temporal effect Mth for the capture-recapture process with τ = 2. To assign a capture history to each individual, we first assumed that the capture probabilities for days 1 and 2 were p i1 ≡ p 1 = 0.55 and p i2 ≡ p 2 = 0.75 for the total generated population. These probabilities were used to find reasonable values for the specific coefficients for the observed categories in terms of The final individual capture probabilities were then computed according to (9)-(10) including the generated random weight as a covariate, implying v = 1.
4. Remove individuals with capture history according to category 0 (undetected). 5. Estimate abundance using each of the methods described in Section 3.1 and fit an AR(2) model to the resulting time series including both the A and P variants.
The choices made in this simulation study intended to approximate the characteristics of a real ecological data set. Specifically, z | t ∼ Log normal ln t , ln w , 31 − 21 = ln p 1 1 − p 1 and 31 − 11 = ln p 2 1 − p 2 .

TA B L E 1
The estimated average coverage and RMSE for all combinations of (ϕ 1 , ϕ 2 ) in the four methods, using five levels of 2 . The AR(2) process was either fitted to the log abundance (A) or the log rate of the corresponding Poisson process (P) we have chosen to simulate rather short time series, having similar length as the real data set used in Section 4. Also, the initial capture probabilities for day 1 and day 2 were close to the proportions of captured individuals in the real data set (being 0.55 and 0.77, respectively).
Our next step was to apply INLA and fit the AR (2) Table 1 displays the average performance in terms of coverage and RMSE for the different methods used to estimate density dependence, including the two variants A and P. The averages were computed across all the given combinations of (ϕ 1 , ϕ 2 ) and for each of the five fixed values of 2 . Due to the short time-series length, coverage using the Baseline method will not achieve the nominal level of 0.95 (only nominal for the A variant). It is well known that estimators for the coefficients of AR processes are biased for small sample sizes (Shaman & Stine, 1988). Furthermore, the Baseline method for the P variant is not optimal considering it is based on the true log counts rather than alternatively generating true log rates of a Poisson process. provided generally higher coverage than using the A variants. When comparing the different methods using RMSE, which considers both bias and variance, we see that CR-INLA had the smallest error in all cases, while the method ObsCount had the largest error. Again, the differences between the methods were very small except for the lowest levels of the process variance. In general, RMSE was reduced for all methods as the process variance increased. Moreover, RMSE was higher for the P variants compared with the A variants at the two lowest process variance levels, using all methods. This was due to both an increased variance and bias, which explains why the P variants had higher coverage.

| Simulation results
The estimation bias of the different methods can be assessed explicitly in Figure Figures S1-S9). For the two lowest levels of process variance, the estimation bias using CR-INLA was slightly lower than using the other methods for all combinations of (ϕ 1 , ϕ 2 ). When the process variance was increased, the different methods gave approximately the same point estimates. The bias was slightly larger using the P variants compared with the A variants. This was in correspondence with the higher averages of the RMSE values for the P variants, as already observed.
To further study coverage and RMSE for each of the 15 combinations, we computed a joint coverage being the proportion of times both of the estimated 95% credible intervals contained ϕ 1 and ϕ 2 , respectively. We also computed a joint RMSE for both parameters, defined by The results for coverage and RMSE are shown in Figures 3 and 4, respectively. Our results showed that the coverage was smallest and RMSE is largest when |ϕ 1 | = 1 in most combinations of the AR coefficients. CR-INLA was seen to give the highest coverage and lowest RMSE for most of the combinations when 2 = 0.08, at least for the A variants. When 2 = 0.32, the results were very similar for all methods.
In summary, we can conclude that including capture history information improved the estimation of density dependence in process models having low process variance. Out of the tested methods, our suggested approach CR-INLA performed best, followed by CR-VGAM. For the given simulated data, the additional step of estimating log rates of a Poisson process resulted in larger RMSE.
Finally, we notice that both of the two AR coefficients were underestimated, and this bias increased with the absolute values of the coefficients.
The given simulation study was based on certain choices to illustrate a capture-recapture scenario using an AR(2) process model.
Here, we have assumed independent capture probabilities for the two capture sessions. The given approach could have easily been adapted to other models given by Otis et al. (1978), such as to also include a behavioral effect. Longer time series would have improved the estimation results using all of the suggested methods, albeit being less realistic from an ecological point of view.

| E S TIMATING DEN S IT Y DEPENDEN CE US ING A RE AL DATA S E T
In this section, we estimated density dependence for a real capturerecapture data set of small mammals, collected at 20 different spatial locations over a period of 18 years. Our main focus was to assess density-dependence estimates, studying how inclusion of capture history influenced the estimation. Using the CR-INLA approach, we estimated capture probabilities by the regression model in (6), including individual-specific covariate information and random effects. We proceeded to estimate the true abundances at each time point for each spatial location using (2). Finally, we fitted the AR(2) model to estimate density dependence and compared the results with using the methods CR-VGAM and ObsCount. For all three methods, we assessed both the A and P variants.

| Data description
The data included a total of 3,090 gray-sided voles, captured alive in

| Observation model selection, estimating capture probabilities
To estimate individual capture probabilities, we used the whole data set across time points and stations. Our first step was to select a reasonable observation model. Fitting the regression model in (6) To select which variables should be included, we evaluated different models using various information criteria. When applying CR-INLA, we used the estimates for the deviance information criterion (DIC) (Spiegelhalter, Best, Carlin, & van der Linde, 2002) and Watanabe-Akaike's information criterion (WAIC) (Watanabe, 2010).
When using CR-VGAM, we used the estimates of Akaike's F I G U R E 2 Posterior mean estimates of ϕ 1 and ϕ 2 , for the A variants on the left (panels a and c) and P variants on the right (panels b and d). The points of intersection of the dotted gray lines correspond to the true parameter values. The intersections, at which each set of dots lean to, correspond to the true value of that given set. Panels a and b show results when 2 = 0.08, whereas c and d correspond to 2 = 0.32 information criterion (AIC) (Akaike, 1973) and the Bayesian information criterion (BIC) (Schwarz, 1978).
An overview of the different models and the estimated information criteria is shown in  (6) is here given by F I G U R E 3 Joint coverage for different combinations of (ϕ 1 , ϕ 2 ) for 2 = 0.08 (panels a and b) and 2 = 0.32 (panels c and d). A variants are represented on the left (panels a and c) and P variants on the right (panels b and d). The results were split into 3 sets (ϕ 2 ∈ (−0.8, −0.5, −0.2), where each set includes the coverage results for ϕ 1 ∈ (−1, −0.5, 0, 0.5, 1)    (2)  using both the A and P variants. Station 9 did not have enough observations for the parameters to be estimated and was thus not included in the results.
The main results are displayed in Figure 7, showing the posterior mean estimates of the AR coefficients for the two seasons, for variants A and P. The estimates of both direct and delayed density dependence were very similar using all the given methods and were thus lumped together (see Figures S11 and S12 for detailed values).
Interestingly, the differences seen in the capture probability estimates between CR-INLA and CR-VGAM do not seem to have influenced the density-dependence estimates. This is in correspondence with the simulation study in Section 3, as the process variance 2 for all of the stations was quite high, with the overall estimated average being ̂ 2 = 0.9. This value likely corresponds to an overestimation of 2 , with the true average being likely significantly lower and closer to 0.6 (see Figure S10 for bias estimates at different variance levels).
In both spring and fall, the estimates of ϕ 1 varied from around −0.25 to 0.6, whereas the estimates of ϕ 2 ranged from around 0 to −0.8.
For all stations, except 3 and 13, the estimated time series showed a semi-periodic behavior. We also notice that the AR(2) coefficients varied with season for the same station, which suggests a seasonal effect in the density dependence. Additionally, during both seasons, the results indicate a decreasing trend in the value of ϕ 2 along the given transect (from coast to inland).

| D ISCUSS I ON
The main goal of this paper was to assess the importance of includ- with ObsCount (A) for the lowest tested process variance; see Table   S1). However, in scenarios with a large process variance, the methods that estimated capture probability did not stand out, producing extremely similar results compared with the observed counts approach. Furthermore, parameter estimates for both AR coefficients were generally biased toward 0, using all the methods, increasingly underestimating the absolute values of the parameters. In the context of quasi-periodic dynamics described by an AR(2) process, this means underestimating the strength of direct ϕ 1 and delayed ϕ 2 density dependence, and overestimating the process variance of the AR(2) model (see Figure S10).
The data collected in Porsanger showcased vole populations with very large fluctuations in abundance, as is typical of such systems (Cornulier et al., 2013;Henttonen & Hanski, 2000 Extending our approach to other observation process models (e.g., spatial capture-recapture models (Royle, Fuller, & Sutherland, 2017), including individual heterogeneity (Efford & Mowat, 2014), would provide a general approach to reducing biases in population dynamic models. One disadvantage of the CR-INLA method is that it would be cumbersome to apply in CR data sets with more than 3 days, given the data expansion necessary to fit multinomial likelihoods in INLA, where all the category combinations, observed and not, must be present. This could potentially be automatized as in Bayesian fitting of capturemark-recapture models (McCrea, 2014).
Two limitations of our simulation study were pointed out during the revision process of this manuscript. Given the complex nature of the simulation setup, involving 4 methods in 2 variants, and 75 combinations of parameters, the running time proved to be lengthy.
Even when running in parallel, we needed roughly 150 hr to obtain the full results, from running 200 simulations per unique combination of parameters. Time thus became a constraint and prevented us from running a higher number of simulations, such as 500 or even 1,000. Nonetheless, we believe our results show true patterns as we noticed early convergence in both coverage and RMSE, from around 50 simulations. Moreover, because we decided to split the results by process variance level, our aggregated averages, displayed in Table 1, combine 3,000 simulations for each method. In addition, it was noted that we did not propagate the uncertainty from our first modeling step into the second. We recognize that this would be a definite advantage, working with a Bayesian framework, and we have now investigated this. We first generated 200 posterior samples from the fitted multinomial model, and then, we fitted the AR(2) model to these samples. The resulting variance in estimating the AR coefficients was very small. This can be seen with an example displayed in In summary, we have found that capture-recapture information given by two capture events contributes to improve density-dependence estimates of AR(2) models with low process variance. In such cases, we recommended that capture heterogeneity is accounted for in the observation model, as it can constitute an important part of the total sampling error. Further analyses are required to assess how more capture events could impact process estimation.

ACK N OWLED G M ENTS
This study was supported by the COAT Tools project. We thank Prof.
Rolf Ims for taking part in the collection of the field data used to test the studied methods. We thank the anonymous reviewers for their comments that improved the quality of this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

O PE N R E S E A RCH BA D G E S
This article has been awarded Open Materials, Open Data Badges.
All materials and data are publicly accessible via the Open Science Framework at https://doi.org/10.5061/dryad.fj6q5 73rr.

DATA AVA I L A B I L I T Y S TAT E M E N T
The simulation code used to obtain the results in this paper and the real data set used to test the methods have been submitted to Dryad