Simulation-based assessment of the performance of hierarchical abundance estimators for camera trap surveys of unmarked species

Knowledge on animal abundances is essential in ecology, but is complicated by low detectability of many species. This has led to a widespread use of hierarchical models (HMs) for species abundance, which are also commonly applied in the context of nature areas studied by camera traps (CTs). However, the best choice among these models is unclear, particularly based on how they perform in the face of complicating features of realistic populations, including: movements relative to sites, multiple detections of unmarked individuals within a single survey, and low detectability. We conducted a simulation-based comparison of three HMs (Royle-Nichols, binomial N-mixture and Poisson N-mixture model) by generating groups of unmarked individuals moving according to a bivariate Ornstein–Uhlenbeck process, monitored by CTs. Under a range of simulated scenarios, none of the HMs consistently yielded accurate abundances. Yet, the Poisson N-mixture model performed well when animals did move across sites, despite accidental double counting of individuals. Absolute abundances were better captured by Royle-Nichols and Poisson N-mixture models, while a binomial N-mixture model better estimated the actual number of individuals that used a site. The best performance of all HMs was observed when estimating relative trends in abundance, which were captured with similar accuracy across these models.


Accuracy and precision of parameter estimates
We assessed estimator quality by inspecting the relative bias, RMSE and 95% CI coverage for both observation and state parameters for each model under various combinations of study design, population size N and HRA.Importantly, we found that increasing the number of grid cells that were sampled beyond 25% had a negligible impact on estimator quality (results not shown).

Detection parameters
All HMs consistently overestimated θ det , and the relative bias increased as the HRA increased (Fig. 1a; Table 1).Furthermore, none of the models reached the 95% CI coverage threshold for any of the simulation scenarios.The BP is the only model that reached coverages > 10 in most of the scenarios (Supplementary table S1).Nevertheless, the RMSE for BernP and BP is relatively small, with median values across all scenarios around 0.05 and 0.02 respectively.Overall, θ det obtained through the BP were the most accurate according to the median proportion of simulations with |Relative Bias|≤ 0.5 (Table 1).For BernP and PP, patterns in relative bias were similar for constant θ det , θ det with spatial variation analysed by naïve and covariate models (Fig. 2a).However, for BP the relative bias in θ det decreased when spatial variation was present.In addition, these scenarios led to lower variability in relative bias in θ det .

Abundance
Regardless of the simulation scenario, BP consistently overestimated abundance (Fig. 1b).Both BernP and PP slightly underestimated in most of the simulation scenarios.Overall, the relative bias in was nearly unaffected by the HRA.However, the relative bias in increased for BP, and decreased for BernP and PP when population sizes N become smaller.When HRA ≥ 10.38 and ρ = 0.70, 95% CI coverages were close to the threshold for BernP and PP (Supplementary table S2).None of the 95% CI coverages in BP reached the threshold.Furthermore, RMSEs were generally low for BernP and PP, but not for BP.According to the median, PP most accurately estimated , followed by the BernP and BP, with respectively 58%, 21% and 3% of simulations resulting in |Relative Figure 1.Mean relative bias (dots/ triangles), together with 2.5% and 97.5% quantiles (solid lines) in the estimated (a) detection parameters θ det , (b) abundances , (c) site-use frequencies use and (d) trends in .Results are displayed for all combinations of model used (BernP, BP, PP), population size N (60, 120, 240, 480), closure (closure, non-closure), movement parameters σ (150, 300, 600, 1200) and ρ (0.7: dots, 0.95: triangles), and the emerging home range area in km 2 (0.65, 2.61, 10.38).Line of equality (dotted line).Average relative bias for each HM (dashed line).For visual clarity, x-scales are different for the subpanels in (b) and (c).
Vol:.( 1234567890 1).Positive bias in from BP was stronger when spatial variation in θ det was present, especially when this information was used in a covariate model (Fig. 2b).For all models, the variability in relative bias in lacked a clear pattern.

Site-use frequency
For open populations, relative bias in use differs considerably with that observed for for all models (Fig. 1c).Both BernP and PP strongly underestimated use in all simulation scenarios, while BP slightly underestimated or overestimated use depending on the simulation scenario (Fig. 1c).Negative bias in use most frequently occurs when ρ = 0.70, while positive bias dominated when ρ = 0.95.Negative bias in use was observed for BernP and PP for all simulation scenarios.None of the HMs was able to reach the 95% CI coverage threshold in any of the simulation scenarios, and even for the most accurate model, i.e., BP, the median RMSE across all simulations was substantial (Supplementary table S3).According to the median proportion of simulations resulting in |Relative Bias|≤ 0.5, BP accurately estimated use in half of the cases, while both PP and BP failed to produce accurate estimates of use (Table 1).The patterns in relative bias in use were affected by spatial variation in θ det in a similar way as those in , that is, stronger positive relative bias for BP (Fig. 2c).

Relative abundance trends
When relative trends in abundance were of prime interest, BernP and PP outperform BP in nearly all scenarios (Table 2).Despite some variability in accuracy, Fig. 1d revealed that the mean estimated trends in were close to the true trends for all models and for all simulation scenarios.However, the |Relative Bias|≤ 0.5 revealed a tendency for poorer trend estimates with stronger departures from closure.According to the median, the models yielded respectively 52% and 69% (BernP), 44% and 56% (BP), and 50% and 65% (PP) of trend estimates (10% and 20% decline in population size N ) with a |Relative Bias|≤ 0.5 (Table 2).Average relative bias in trend estimates were nearly unaffected by spatial variation in θ det (Fig. 2d).However, the variability in trend estimates from BernP and PP, but not BP, was substantially larger when spatial variation in θ det was present.Finally, trend estimates from this scenario had a lower accuracy across all these models (Supplementary table S5).

Discussion
In this study, we have evaluated three hierarchical abundance estimators (BernP, BP and PP) in the context of complicating features of realistic animal populations simulated through random walks.
Increasing trends in per-capita detectability and detection rate with population size N , as indicated in Table 3, are within the line of expectations, since denser populations are expected to result in more frequent detection  Results are displayed for all combinations of model used (BernP, BP, PP), the submodel and covariate structure for θ det (constant, spatial variation (naïve model), spatial variation (covariate model)), population size N (120, 240), ρ (0.7: dots, 0.95: triangles) for σ = 300 (emerging home range area of 2.61 km 2 ).Line of equality (dotted line).Average relative bias for each HM (dashed line).For visual clarity, x-scales are different for the subpanels in (b) and (c).www.nature.com/scientificreports/events (i.e., an individual crossing the viewshed).We follow Neilson, et al. 18 , in assuming that populations of individuals occupying large HRAs, such as apex predators, are generally less dense (Table 3).This assumed confounding of HRA and population size, explains the decreasing trend of θ det with HRA.Evidently, the per- capita detection probability is consistently lower than the per-capita detection rate, as some individuals may be detected more than once on a survey day.Regardless of geographical closure and the HRA, the abundance equals N/N site , resulting in the densities indicated in Table 3 ("Abun.").Site-use frequency, use , is a product of both population size N and HRA.In theory, use increases with both increasing N and HRA.However, in our simulations these parameters are confounded to simulate realistic animal populations 18 .Hence, use is seemingly unaffected by N and HRA, as their effects cancel out.Only, the effect of closure violations on use is clearly visible (Table 3; "Site-use freq.").
In our simulations, true θ det are extremely small (i.e., always θ det < 0.05).Hence, it is perhaps not surprising that in these scenarios positive bias in θ det persisted across most of the HMs considered (Figs.1a and 2a).It is notoriously difficult for these HMs to accurately estimate parameters near the boundary of their parameter space 20,21 .In cases where θ det is accurately estimated, abundance estimates appear to suffer from increased positive bias (Fig. 2a-c BP-covariate model) in order to provide an explanation for the observed counts (i.e.some incidences of high counts occur because of accidental double counting).The use of a beta-binomial distribution for p in the BP case, which can alleviate bias resulting from correlated detections within samples, may alleviate some of this bias in abundance, while conserving accuracy in θ det 12 .Overall, we would like to warn practitioners that very small detection probabilities/ detection rates, which are frequent in camera trapping studies, often do not contain sufficient information to confidently analyse within complex modelling frameworks.However, modelling detections of multiple species jointly may hold the potential to overcome issues related to low detectability 22 .
Furthermore, our results reveal that none of the HMs considered in this study were able to estimate accurately under all scenarios (Figs.1b and 2b).The BernP model slightly underestimates in most scenarios as opposed to simulation results of Royle and Nichols 6 .Nonetheless, BernP estimates more accurately than BP, despite only using binary input data.A possible explanation for the binary input model to outperform a count model lies in the occurrence of double counting of individuals.Indeed, we assumed that animals were unmarked, and hence could not be individually identified by CTs, which led to badly inflated BP , a well-known issue for BP 10,13 .However, this problem does not apply to binary models, since double counts do not affect a binary encoding.PP, which alleviates most of the bias in observed under BP, produces the most accurate estimates of from any of the models considered here.In line with our expectations, a PP even yields accurate estimates in over 80% of the cases when populations are closed (Table 1).This highlights the importance of accounting for double counts when absolute abundances are desired.Nonetheless, PP are still somewhat negatively biased in open populations.Specifying a distribution for the latent N i 's that allows for overdispersion and/or underdisper- sion can be a solution to improve inference when extra-variation in counts is present among sites.We have used zero-inflated Poisson distribution for N i 's 19 , which did not result in more accurate 's (results not shown).On the contrary, both BZIP and PZIP display some of the strongest biases that we observed in our simulations.This is likely a consequence of these models being very complex, i.e., they require the estimation of two latent parameters in addition to a detection probability/ detection rate, for the sparse data simulated in our study.Despite their inaccurate estimates BZIP and PZIP models were preferred by LOO-ELPDs, which is possibly because they provide a good description of the many zero-observations present in the data 3,19 .Negative binomial models have also been suggested to deal with extra-variation in counts, but suffer from a similar "good-fit-bad-prediction dilemma" 4,23 .Hence, there seems no reason to assume that HMs with negative binomial state models perform any better than the zero-inflated Poisson HMs evaluated here.
Given that is not well estimated for open populations, the use of HMs for inferring the site-use frequency use could be an alternative that provides reliable answers to many wildlife related questions.Site-use frequency is effectively a product of population size and HRA, hence reflects both the true number ( N ) of animals present, as well as their movement pattern 11,15,24 .We expected that this quantity ( use ) would be better estimated than abundance ( ) when populations are open, yet our results suggest that this is not the case (Figs.1b-c and 2b-c).The BernP, which uses only binary input, is too conservative to produce estimates in the order of magnitude of use .PP on the other hand can attribute additional counts resulting from groups moving across sites to a higher detection rate ( µ ).This seems to be supported by estimates of detection rates in the PP that steadily increase with increasing HRA (Fig. 1a).Site-use frequency estimated from BP outperform those of BernP and BP, but, on average, they still only yield accurate results in 40% of the simulations.Contrasting with our results, Nakashima 11 found that BernP and PP, and not BP, were able to accurately estimate the site-use frequency.Fogarty and Fleishman 15 , on the other hand, also find that a BP, which accounts for permanent movements among sites (i.e., immigration and emigration) yields accurate estimates of site-use frequency ("superpopulation abundance").However, they point out that the ecological interpretation of site-use frequency is difficult when the movement characteristics of animals are unknown.The problem is that there is no way of discriminating between sites that are frequently used because they are on a travelling route and between sites that are frequently used because they offer important resources.Thus, site-use frequency does not appear to provide a reliable solution for the interpretation of abundances of population that are open relative to the sampling frame.Abundance (or density) of open populations may be estimated more accurately by spatially-explicit models for count data, which leverage observed spatial correlation in detections 25 .However, for entirely unmarked populations, which we simulated in this study, the precision of abundance/density estimates based on this method is typically low 26,27 .Therefore, future developments in spatially-explicit models for count data should focus on improving precision of their estimates.
Regardless of the capability of HMs to yield unbiased estimates of (or use ), it is often more interesting to explore relative trends in .Moreover, while there are concerns regarding identifiability within HMs for abundance, their ability to estimate relative abundances has not been questioned 8,17 .Here, we explored trend estimates obtained through comparing abundances N / 0.9N and N / 0.8N to induced declines of respectively, 10 and 20%.Overall, simulated declines in population size N of 10 and 20%, are accurately estimated more frequently than standalone (Table 2).However, when a relative bias up to 50% is tolerated, 48%, 66%, 50% and 31%, 44%, 35% of the trends will be misleading for respectively 10% and 20% declines estimated by BernP, BP and PP.Thus, HMs for abundance are better suited to capture relative trends in abundance, but obtaining them post-hoc, as is often done by practitioners, may still be misleading in many cases.Estimates of relative trends in abundance may be biased particularly when they are confounded with between-year variation in detection probability that is unaccounted for.Because of our focus on a single sampling season within a single year, we did not evaluate the impact of temporal variation in the probability of detecting an individual that crosses a camera viewshed, i.e., P(detected|in viewshed), which is a study limitation.However, the true detection probability is a product of P(detected|in viewshed) and the probability that an individual crosses a camera viewshed if it is present, i.e., P(in viewshed).Inherent variation in the simulated movement trajectories of our study induce small-scale variation in P(in viewshed), and from our results, it appears that relative trends are robust against this form of variation.Our results suggest that when spatial variation in P(detected|in viewshed) exists, the accuracy of relative trends from BernP, BP and PP decrease owing to a larger variability (Fig. 2d; supplementary table S5).Future research should determine if relative trends in abundance suffer from similar losses in accuracy when stronger temporal variation in detection probability exist.
While there are potentially many other candidate models that are useful to estimate animal abundance (or density) in ecological surveys, we explicitly chose not to include them in our comparison.We did not include distance sampling 28 or random encounter (and staying time) models 29,30 as they are hampered by a need for additional information, either distance data or movement data that are often absent or of low quality in camera trapping studies.We did not consider capture-recapture methods 31 since identifying individuals through photographs is difficult when individuals are unmarked, which was one of the assumptions of our study.For the same reason, we disregarded methods based on removal sampling 32,33 , where individuals need to be identified to be removed from the study after their first detection.Finally, we did not include models in Moeller, et al. 34 : the timeto-event model was disregarded because it assumes perfect detectability, the space-to-event and instantaneous sampling models require time-lapse photos which may lead to many missed encounters.
In conclusion, we have shown that under realistic settings, and with 24-h as the time unit of temporal replicates, the probabilities of detecting an unmarked species like wild boar using camera traps are extremely low.In addition, when individuals cannot be individually identified, and thus double counting cannot be excluded, we find strong biases in estimates of abundances for BP.Through accommodating false-positives, the PP model was able to estimate more accurately.Furthermore, we reveal that shifting the interpretation to frequency of site-use does not improve the inference in general, but leads to better estimates in BP.Finally, we report that relative trends in abundance are estimated with greater precision than absolute abundances, but that they can still suffer from bias especially when spatial variation in detection probability exist.Practitioners should thus avoid using absolute abundance estimates, regardless of their interpretation, from static HMs and turn to relative trend estimates, but also these should be used with caution.Depending on the context, model-based approaches, taking into account temporal trends, spatial trends or a combination thereof, might further improve the accuracy of inference from HMs.From this perspective, it would be valuable to assess the estimator quality of temporally explicit HMs in a simulation study similar to ours.

Material and methods
Our simulation study is designed with a specific unmarked target species (wild boar) in mind, in a specific park (Hoge Kempen National Park (NPHK) situated in Belgium).In that park, the species commonly investigated throughout a relatively short survey period within a single year, a strategy we adopt as well.Yet, our findings, based on this design, are applicable to a wider range of camera trapping studies since most of them take place in areas of comparable size, with similar CT spacing and with a focus on mammal species, many of which live in groups and display home range tendencies 1,35 .In fact, ungulates comprised the second most studied group of species according to Burton, et al 1 .Moreover, we do not simulate climatic conditions or landscape features that are specific to the NPHK.

Simulated space
We simulated a space with a total surface area of 116.64 km 2 (10.8 km * 10.8 km; Fig. 3), which roughly equals that of the NPHK 36 .Next, we simulated a grid layer consisting of 144 grid cells ( i = 1, 2, . . ., 144 ) of 0.9 km × 0.9 km, for the placement of CTs.Our grid cells are nine times coarser than the grid cells in Wevers, et al. 36 , such that the area (0.81 km 2 ) of a single grid is closer to the empirical home range sizes of wild boar 37,38 .Next, we placed a CT in 25% (36) of the grid cells, according to a randomised regularly spaced sampling design.Within sampled grid cells, a CT "detector" was then simulated by placing a single CT at its centroid with a viewshed radius r = 15m and an angle of view θ = 42 • , based on camera specifications of the Reconyx Hyperfire HC600 39 .All CT views- heds, in which passing animals may be detected, were simulated facing North.

Movement trajectories
We assumed that N individuals of a species of interest (with N fixed a priori according to Table 3) live in the region of interest, say A .Moreover, we assumed that these animals moved in groups, which we monitored dur- ing 25 consecutive days within a single season in which environmental factors potentially affecting probabilities of detecting an individual were assumed to be more or less constant over time.We started by simulating forty independent, group-specific and activity-adjusted random walks according to a bivariate Ornstein-Uhlenbeck process (OU).To generate realistic animal trajectories for a species with home range tendencies, we varied three parameters: the population size ( N ) and parameters that control serial autocorrelation ( ρ ) and variance in move- ments ( σ ) (Table 3).Note that the size of the HRA emerges from the choice of σ.Let u represent a point location, defined by a set of coordinates {x, y} .For g = 1, . . ., G groups and t = 1, . . ., 600 one-hour steps (25 days of 24 h), the probability for a given set of group-specific positions u g(t+1) in the next hour t + 1 according to a discrete-time OU can be expressed by: where s g are the group-specific home range centers,ρ controls the tendency to move in the direction of the home range center, σ determines the diffusion of the movement and I represents the 2 × 2 identity matrix.However, when the species was resting (here we assumed a nocturnal species, but the results would be similar for a diurnal species with similar proportions of active time), we assumed that all groups stayed in place, hence u g(t+1) = u gt if (t + 1)%24 ∈ {0 : 5, 18 : 23} , where '%' represents the modulo.
We generated group sizes N g for each group g by sampling empirical group sizes of the target population that satisfy N = G g=1 N g .In addition, we simulated 10% and 20% declines in population size for all populations (i.e., N 0.9 = 0.9N and N 0.8 = 0.8N ) to explore the ability of the models to estimate these declines by comparing abundances post-hoc ( N / 0.9N and N / 0.8N ).Group-specific home range centers s g were (i) sampled from the fixed set of 0.9 km × 0.9 km grid centroids (resulting in closure) or (ii) sampled randomly from the entire study space (resulting in closure violations).When simulating hourly displacements, information on finer scale movements is lost.However, these may be important when they result in the passing or not passing through the camera viewshed 40 .To ensure that we generated more realistic, finer scale movements, we simulated Brownian motion between consecutive positions x = u gt and y = u g(t+1) using a Brownian Bridge (Supplementary section 1).

Simulating count and detection-non detection data
After simulating animal movements, we generated sightings for each hour h = 1, 2, ..., 24 within each day j = 1, 2, . . ., 25 at locations where a group of individuals crossed the CT viewshed, for i = 1, 2, . . ., 36 CT loca- tions.We assumed that all individuals in a group were detected when they pass this viewshed (i.e., P(detected|in viewshed) = 1).Note that the probability to detect an individual also depends on the probability that an individual that is present enters a CT viewshed (i.e., P(in viewshed), which will be < < 1).Additionally, we assessed the impact of imperfect P(detected|in viewshed) ≈ 2/3, which also display spatial variation according to a unique Matérn process 41 per simulation for a subset of the simulation scenarios (i.e., σ = 300 resulting in HRA = 2.61 km 2 ) that best reflect empirical home range areas of the target species.In both cases, we obtained hourly counts as: y ijh ∼ binomial(n ijh , P[det|inview] i ), where n ijh is the number of individuals crossing a CT viewshed located at grid cell i on day j and hour h .Hourly counts y ijh were then aggregated across days to yield daily summaries of counts y ij = 24 h=1 y ijh .Since we have assumed that individuals are unmarked, the same group could potentially contribute to y ij on each hour if they re-appeared in front of the CT.When only binary inputs were required (i.e. for BernP), we reduced count data to detection/non-detection data.

Assessment of criteria for the model inferences
For each simulation replicate, we calculated the true abundance ( ) as the average number of individuals per grid cell, assuming that all individuals are bounded to the grid cells that contain their activity centers (i.e., the average number of activity centers per grid).This quantity was used as ground truth to assess the ability of the HMs to make accurate inferences on N .In simulation scenarios with closure violations, individuals move across their initial grid cells such that the number of individuals per grid cell changes over time.Hence, we also calculated site-use frequency ( use ) as the average number of individuals that have used a grid cell at least once during the simulated study period, given that all individuals completely covered their HRA (i.e., the average number of HRAs overlapping a grid cell) 11,15 .Finally, we assessed the ability of HMs to accurately estimate relative abundances by comparing how accurately N / 0.9N and N / 0.8N captured true induced population declines of respectively 10 and 20% for different choices of N.

Statistical models
Within the current study, we assessed the performance of HMs, that are commonly used to make inference on abundance.Depending on the data at hand, BernP 6 (for replicated detection/ non-detection data) and BP/PP 7,11 (for replicated count data) models are commonly used.In this study, we adopt a Bayesian estimation framework, as it allows the flexible modelling of HMs.Model fitting was performed using Stan (via the R package cmdstanr), a probabilistic programming language that enables Bayesian estimation through a dynamic Hamiltonian Monte Carlo sampler 42 .
All three statistical models start from the series of observations y ij .Observations y ij can either consist of binary detection (y ij = 1) / non-detection (y ij = 0) data or count (y ij ≥ 0) data.The mathematical structure and distributional assumptions for y ij and N i , the number of individuals present at site i , are given in Table 4.Note that in practice, N i 's are obtained by taking the finite sum over K possible latent abundances to make the maximisation of the likelihood numerically tractable 7 .Here we chose an upper bound K , such that K = 100 for detection/ non-detection data and K = max y ij + 100, ∀{i, j} for count data.Moreover, we construct the likelihood by marginalising over N i 's with upper bound K given that Stan cannot sample discrete latent variables.Importantly, p -the per capita detection probability -represents the probability that an individual is detected (assuming that there are no false-positives), while µ -the per capita detection rate -corresponds to the rate at individuals are detected each day (when false-positives occur).For notational simplicity, we will refer to these detection parameters jointly as θ det .For simulation runs with spatial variation in P(detected|in viewshed), the p u g(t+1) |u gt , s g = Norm s g + e ρI u gt − s g , σ 2 I − (e ρI • σ 2 I • e ρI ) 'naïve' models as in Table 4 were fitted alongside 'covariate' models that used the values from the Matérn process, which generated P(detected|in viewshed), as a covariate x : logit(θ det ) = β 0 + β 1 x.
To assess goodness of fit for competing HMs, we calculated Bayesian P-values 10,43 .Additionally, we computed and compared the LOO-ELPD for each model 44 .For more information on the prior specification, goodnessof-fit evaluation and MCMC convergence of HMs in Stan we refer to Supplementary sections 2-3.To give a representative view on the results, we did not exclude simulation runs where convergence and/or identifiability issues occurred from further analysis.

Relative bias, RMSE and coverage probability
In this study, we assessed estimator performance from three Bayesian HMs of interest by exploring the bias, relative bias, RMSE and 95% CI coverage.We calculated the bias as Bias(θ) = ( θ − θ) ; the relative bias as Rel.Bias(θ) = ( θ − θ)/θ ; the RMSE(θ) = n i=1 ( θ i − θ i ) 2 /n ); and the 95% CI coverage as the proportion of simulations where the true parameter value is enclosed by the 95% CI.We used |relative bias|< 0.5 as the threshold for simulations that yielded accurate estimates.Moreover, we regard 0.95 of simulations that reached the 95% CI coverage as an indicator for good uncertainty quantification.

Figure 3 .
Figure 3. Graphical overview of the simulation study.The dimensions of the state space represented by the grid layers are 10.8 km × 10.8 km.OU: Ornstein-Uhlenbeck, BernP: Bernoulli-Poisson mixture model, BP: Binomial-Poisson mixture model, PP: Poisson-Poisson mixture model.

Table 1 .
Summary table for estimator quality of detection parameters θ det , abundances and site use frequencies use obtained from three Bayesian hierarchical models (BernP, BP, PP).Cells display the proportion of simulation replicates that satisfy |Relative Bias|≤ 0.5.

Table 3 .
Design of the simulation study, based on three ecological/ behavioural parameters, home range area, population size and mean hourly displacement of (a group of) individual animals.The range of true detection probabilities p and rates µ , resulting from each scenario, are also indicated.Vol.:(0123456789) Scientific Reports | (2023) 13:16169 | https://doi.org/10.1038/s41598-023-43184-w

Table 4 .
Mathematical structure of the hierarchical models used in this study, their main reference and applications in wildlife research.− 1 − p Ni y ij |N i ∼ Bernoulli(P) N i ∼ Poisson( ) 6 45,46 BP Binomial/Poisson Counts y ij |N i ∼ Binomial(N i , p) Poisson y ij |N i ∼ Poisson(N i • µ)