# Robust abundance estimation of harvested populations from low quality time series data: A red deer case study

Stig W Omholt, Marlène Gamelon, Erling L Meisingset

► To cite this version:

Stig W Omholt, Marlène Gamelon, Erling L Meisingset. Robust abundance estimation of harvested populations from low quality time series data: A red deer case study. Ecological Indicators, 2025, 173 (113398), 10.1016/j.ecolind.2025.113398. hal-05022701

HAL Id: hal-05022701

https://cnrs.hal.science/hal-05022701v1

Submitted on 7 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Original Articles

# Robust abundance estimation of harvested populations from low quality time series data: A red deer case study

Stig W. Omholt<sup>a</sup>, Marlène Gamelon<sup>b</sup>, Erling L. Meisingset<sup>c, *, e</sup>

ARTICLE INFO

Keywords:

Abundance estimation

Citizen data

Hunter observations

Male emigration

Population model

Red deer

Time series data

ABSTRACT

Reliable estimates of the size and composition of harvested populations over time are key to designing adequate population management plans, regardless of management objectives. In Norway, a national system for collecting and analysing hunter-reported data on red deer (*Cervus elaphus*) has been operational for about 20 years. The system was expected to provide population metrics that would substantially improve deer population management routines at the municipal level. This has proven to be challenging when using existing state-of-the-art estimation methodology. The main reasons are that the variation in the observation data is generally much larger than population abundance variability, and that one does not have a clear understanding of the stochastic process generating the observation data. Here, using hunter-reported observation data and harvest data from six Norwegian municipalities collected in the period 2007–2023, we show that a straightforward estimation methodology based on population modelling can produce robust abundance estimates despite frequent low quality of the observation data. Its major assets are that it does not involve strong assumptions about the stochastic processes underlying the observation process and that it does not involve assumptions about initial population size and structure in terms of prior statistical distributions. We anticipate that the method can be applied in several other population management contexts, and we think that the results offer fresh perspectives on to what extent noisy citizen-collected time series data can be used to inform management decisions.

# 1. Introduction

Ungulates are important members of most ecosystems in the north- ern hemisphere, and questions related to their conservation and man- agement receive considerable scientific and public attention. The availability of reliable methods to estimate their abundance and popu- lation structure is a key to understanding their population dynamics and how they affect ecosystem processes (Apollonio et al., 2017; Forsyth et al., 2022; Langvatn & Loison, 1999; Meriggi et al., 2009; Morellet et al., 2010). Such methods are also crucial for efficient management through hunting, regardless of management objectives (Forsyth et al., 2022).

In several European countries, hunters gather non-survey game species population metrics to inform wildlife management (Forsyth et al., 2022; Hušek et al., 2021; Morellet et al., 2010). In 2002, Nor- wegian authorities implemented a national system for annual collection

and compilation of moose (*Alces alces*) and red deer (*Cervus elaphus*) data provided for free by hunters (seen deer data), the so-called Cervid Register (CR) (Solberg et al., 2022). The CR was expected to provide crucial population metrics for the formulation of efficient deer popula- tion management routines at the municipal level. However, due to a large, and partly biologically implausible, variation in the reported CR population indices from year to year, it has proven challenging to extract reliable estimates of population size and structure. Therefore, the administrative and scientific value of continuing to collect these data has been questioned (Hušek et al., 2021).

State-space models (SSMs) are increasingly used in ecology to model population dynamics time series data (Auger-Méthé et al., 2021). SSMs are hierarchical models that decompose observed time series data into an observation error component and a process variance (Kéry & Schaub, 2012). This is a very strong concept, and both a frequentist approach and a Bayesian approach can now be used to fit complex SSMs (Auger-Méthé

<footnote><sup>a</sup> Department of Biology, Centre for Biodiversity Dynamics (CBD), Norwegian University of Science and Technology (NTNU), NO-7491, Trondheim, Norway</footnote>

<footnote><sup>b</sup> Laboratoire de Biométrie et Biologie Evolutive, UMR 5558, CNRS, Université Lyon 1, Campus de la Doua, Bâtiment Gregor Mendel, 43 Boulevard du 11 novembre 1918, 69622 Villeurbanne, France</footnote>

<footnote><sup>c</sup> Department of Forest and Forest Resources, Norwegian Institute of Bioeconomy Research (NIBIO), Tingvoll gard, 6630 Tingvoll, Norway</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: erling.meisingset@nibio.no (E.L. Meisingset).</footnote>
et al., 2021). SSMs have been successfully applied in a variety of contexts across a wide range of species from the 1980s in fisheries (Maunder & Punt, 2013) and the last 20 years in terrestrial systems (Schaub & Abadi, 2011; Schaub & Kéry, 2022). They are also widely acknowledged as useful tools in conservation and management (Arnold et al., 2018; Gamelon et al., 2021; Lee et al., 2015; Lieury et al., 2015; Tenan et al., 2016; Zipkin & Saunders, 2018).

Most SSMs combine capture-recapture data from marked individuals with other data types such as population counts or hunting bags (but see Zipkin et al. (2017) for an integration of counts and detection-nondetection data from unmarked individuals). The use of capture-recapture data is arguably the most informative approach to tracking populations (Lebreton et al., 1992; Zipkin et al., 2017). However, in exploited populations, the use of unmarked individuals to generate population indices is the rule, rather than the exception, potentially limiting the applicability of SSMs for such populations. A further issue is that the use of an SSM approach on the Norwegian red deer time series data would require a good prior understanding of the observation process underlying the annual aggregated population abundance indices. However, these indices reflect a very complex process, and the quality of the data, in terms of representing unbiased counting processes, can vary substantially for a number of reasons that cannot be controlled for. Moreover, in the Norwegian red deer time series data, the process variability is much less than the observation data variability, and complex SSMs are in general prone to suffer from parameter- and state-estimation problems when the measurement error is larger than the biological stochasticity (Schaub & Kéry, 2022). Finally, even if the above issues turned out to be less severe than anticipated, the use of SSMs in a local management context would be quite demanding due to the need for substantial computational know-how and the need for locally adapted model specification (Frost et al., 2022; Maunder & Punt, 2013). Taken together, this motivated us to probe whether it was possible to develop an alternative estimation approach less exposed to the above concerns attached to SSMs.

Using a deterministic one-equation population model that projects the development of a red deer population from one year to the next and the L2 norm (or Euclidian distance) as a metric to quantify the degree of compatibility between population projection and observation record, we found that the range of population sizes at time point T compatible with the observation data (seen-deer-per-hour) can, when $T = 17$ years, be estimated with reasonable accuracy. Encouraged by this result, we developed a stochastic sex- and stage-structured discrete dynamic population model informed by harvest data and generated a repository of hypothetical long-term population projections based on a large full factorial set of initial condition parameters (population size and age/sex structure). By letting the upper and lower bounds of the initial condition set to be well outside of what is biologically feasible, we avoided the use of prior statistical distributions regarding initial condition parameters. We then fitted these projections to the census data using the L2 norm as a metric. Using six municipalities as case studies, we show that robust abundance estimates for the period 2007–2023 can be extracted from the current CR data by using the population projections that are most compatible with the observation data.

We anticipate that our estimation method is simple enough to be put into practical use at a local management level. Considering that uncertainty about the impacts of interventions is a major challenge facing adaptive management, the managerial value of having access to reasonably reliable annual population estimates of harvested populations based on long-term time series data is likely to be substantial, as it opens for large-scale use of what may be called *mathematically informed population management*. By this we mean the formulation and execution of short-term management plans targeting a specific set of objectives that are based on locally adapted population models projecting the development of a population from one year to the next as a function of the chosen harvesting plan. Since our estimation method allows for independently obtained population estimates to be fed into

projection models on an annual basis and thus serve as a corrective negative feedback signal, the targeted objectives are likely to be met more efficiently and reliably than what can be achieved by current practices.

# 2. Material and methods

## 2.1. Study areas and data sources

We used reported hunter-data from six municipalities in Norway during the period 2007–2023, five from the county Møre and Romsdal (Averøy, Tingvoll, Surnadal, Sunndal and Vestnes) and one from the neighbouring county Vestland (Lærdal), to assess the reliability of our estimation method (Fig. 1(a)). Although these municipalities are located mainly close to each other, they vary substantially in terms of landscape, topography, and area (Supplementary material, Table S1). The six municipalities are located in the main areas of red deer distribution in Norway and red deer populations are heavily managed by hunting (Loe et al., 2016).

For each municipality, we used publicly available CR data on seen-deer-per-hunting-hour, the aggregated ratio of number of seen hinds per seen stag during the hunting season, and the maximum annual spring count for the period 2007–2023 (Fig. 1(b) gives an outline of the data collection process). In the cases where data within this time window were missing in the CR, these data were kindly provided by the municipal wildlife authorities. Annual harvest figures for the period 2007–2023 for each municipality were obtained from Statistics Norway (https://www.ssb.no/en).

In Norway, municipalities are divided into separate local management units (LMUs) that are further divided into several hunting fields (Meisingset et al., 2018). The seen-deer-per-hour data and the hinds per stag ratio data used in this study come from daily recordings by hunters in each hunting field throughout the hunting season (September 1 – December 23). For each day, hunters record the number of deer seen (divided into age and sex classes, including days with no observations), the number of hours hunted, the number of hunters that participated, the hunting habitat (farmland or forest), and the sex, age, and slaughter weight of deer shot. Since 2017, the recorded data from each hunting field have been reported online to the CR using a specific web reporting scheme. The data of a municipality are not publicly released before the chairman of each LMU has approved the reports of the hunting fields and the municipal authority has approved the reports of the LMUs.

Several municipalities also report spring census data annually (Fig. 1(b)), although this is not yet mandatory. These data are obtained by direct vehicle counts of deer along fixed routes according to a fixed procedure on more or less all farmland known to be used by red deer during spring (Godvik et al., 2009). The counting is organized by municipalities (or LMUs) and is done during the same time window (a few hours) one to three times from March to May (Mysterud et al., 2007). Of the six municipalities studied, only Averøy had not recorded spring census data annually in the period 2007–2023. The spring census data are distinctly different from the seen-deer data, since they are based on a direct count of a substantial, though unknown, proportion of the deer population.

## 2.2. Description of the data

Regarding the number of females, number of males and total number of red deer harvested per year in the period 2007–2023 (Fig. 2, left panel column), the municipalities Averøy and Tingvoll stand out by having systematically harvested more females than males throughout the whole period. The resolution of the observation data available for each municipality varied somewhat (Fig. 2, right column).

Slaughter weight data for the entire period were available only for the municipalities Surnadal and Tingvoll (Supplementary material, Fig. S1), but they are probably representative for all six. Although there
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** (a). Locations of the six Norwegian municipalities studied (Averøy, Tingvoll, Surnadal, Sunndal, Vestnes, Lærdal) (areas highlighted, see Table S1 for detailed geographical information). (b). Study outline. Showing types of data collected, when they are collected and the time points when the population model records the population status (see Table S2).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Parts of the data underlying the study. **Left panel column:** Reported sex-specific harvest figures and total harvest over the period 2007–2023 for all six municipalities studied. **Right panel column:** Reported seen deer data combined and divided into farmland and forest observations and spring census data over the period 2007–2023 for the six municipalities. Times series data divided into farmland and forest observations were not available for the municipalities Averøy and Sunndal, and sufficient spring census data were not available for Averøy.</figcaption></figure>
have been some minor changes in mean body weight during the study period, the data do not support the presence of temporal variation in food supply that would have markedly affected reproduction and survival rates.

The relationship between annual harvest quota and annual total harvest (Fig. S2) indicates that the population size has been decreasing the last few years in four of the six municipalities (Averøy, Tingvoll, Surnadal and Vestnes). Lærdal municipality increased the quota and the harvest substantially in 2017.

Each seen-deer-per-hour time series can be viewed as a realization of a discrete time stochastic process {$X_t: t \in \Psi$}, where X is a real random variable (i.e. the population index), $\Psi$ is an index set {$1, ..., T$} of the stochastic process, and where each collection of random variables {$X_t: t \in \Psi$} is defined on a common probability space ($\Omega$, A, P), where $\Omega$ is a sample space, A is a $\sigma$-algebra to be interpreted as the collection of events which can be assigned probabilities, and P is a probability measure (satisfying the Kolmogorov axioms) that can, somewhat informally, be viewed as a probability distribution on $\Omega$, the collection of events A usually being understood. Each of the six time series is thus one trajectory out of a very large sample of trajectories that could have been realized by the process {$X_t: t \in \Psi$} specific to each municipality. By viewing the seen-deer-per-hour time series as realizations of stochastic processes, we can apply well-established tools that provide indications about the structure of the underlying probability spaces ($\Omega$, A, P).

An Augmented Dickey-Fuller (ADF) test for stationarity (Seabold & Perktold, 2010) indicates we can reject $H_0$ (non-stationarity) only for the series coming from the municipalities Averøy and Sunndal (p < 0.002). A Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test for stationarity (Seabold & Perktold, 2010), however, gives consistent results since we cannot reject $H_0$ (stationarity) in any of the six cases (p > 0.1). These results suggest that the discrete time stochastic processes having generated the data from Averøy and Sunndal can be considered to have been stationary (Seabold & Perktold, 2010), causing the mean and variance to have stayed fairly constant throughout the time series period. If a KPSS test indicates stationarity and an ADF test on the same data indicates non-stationarity, the time series is trend stationary (Seabold & Perktold, 2010). Thus, the discrete time stochastic processes underlying the observation data from the municipalities Tingvoll, Surnadal, Vestnes and Lærdal can be considered to have been trend stationary, causing the mean to have been time-dependent, while the variance has remained fairly constant. However, since a Mann-Kendall trend test ($H_0$: no monotonic trend) did not pick up a trend in any of the six municipalities (p > 0.09), the time-dependency is arguably quite weak. Thus, it appears fair to conclude that the statistical properties of the six stochastic processes that have generated the six seen-deer-per-hour time series have not changed much over time. According to a Shapiro-Wilk test ($H_0$: normally distributed variable) (Virtanen et al., 2020), we cannot reject that the data points of each time series are drawn from a normal distribution (p > 0.06 for Vestnes, and p > 0.15 for the rest). In contrast, a Kolmogorov-Smirnov test ($H_0$: uniformly distributed variable) (Virtanen et al., 2020) strongly suggests that the underlying distributions are not uniform (p < 0.0004).

### 3. Theoretical foundation for the estimation method used

Under a highly quality-assured observation regime, there is a linear relationship between deer seen per hour and population size (Van Deelen & Etter, 2003). However, the observation data at hand are generated by observation routines that are neither well scrutinized nor highly standardized. Thus, based on the variation in reported harvest numbers and knowledge about annual variation in fertility and natural mortality, the impact from several uncontrollable factors causes the observations from one year to the next to frequently vary much more than what the population abundance plausibly can. This makes the relationship between an actual population size and its associated population index sometimes quite dubious, which warrants an estimation

method that does not depend on a linear relationship. As will be shown in the following, an estimation method based on the L2 norm seems fit for our purpose.

#### 3.1. The resolution capacity of an L2 norm metric

Assume a population develops according to:

$$N_t = (1+r)N_{t-1}hN_0 \quad (1)$$

where $t \in \{1, 2, ..., T\}$, $N_t$ is the number of individuals at end of year $t$, $t \in \{0, 1, 2, ..., T\}$, r is the population growth rate (accounting for natural mortality), $N_0$ is the population size at $t = 0$, h = constant, and $hN_0$ is the number of individuals harvested each year. The population development can then be expressed in closed form as:

$$N_t = \left[ (1 + r)^t - h \sum_{i=0}^{t-1} (1 + r)^i \right] N_0. \quad (2)$$

Associate a discrete time stochastic process {$X_t: t \in \Psi$} with this population development, where $X_t$ is a real-valued random variable, for example in units of seen-deer-per-hour. The L2 norm measuring the distance between the population vector and the observation vector can then be expressed as:

$$L2 = \|\lambda N - X\| = \left( \sum_{t=1}^{T} (\lambda N_t - X_t)^2 \right)^{1/2}, \quad (3)$$

where $\lambda$ is a scaling factor that makes the $N_t$ values comparable with the $X_t$ values. To obtain a value of $\lambda$ that minimizes the distance between the population vector and the observation vector, one can differentiate the above expression for L2 with respect to $\lambda$ and solve the equation dL2/d$\lambda$ = 0 for $\lambda$. This gives:

$$\lambda = \frac{\sum_{t=1}^{T} N_t X_t}{\sum_{t=1}^{T} N_t^2} = \frac{\sum_{t=1}^{T} [(1 + \beta)(1 + r)^t - h \sum_{i=0}^{t-1} (1 + r)^i] N_0 X_t}{\sum_{t=1}^{T} [(1 + \beta)(1 + r)^t - h \sum_{i=0}^{t-1} (1 + r)^i] N_0} \quad (4)$$

Since $N_0$ cancels out in the product $\lambda N_0$, this L2 norm metric is independent of the size of $N_0$. The value of L2 will depend on the observation series {$X_t: t \in \Psi$}, that is, the realized trajectory of the underlying stochastic process. This means that given a large number of hypothetical trajectories determined by the probability space ($\Omega$, A, P) associated with a given population development, we will obtain a distribution of L2 norms. If the development of the population had been based on a slightly different value of $N_0$ ($N_0'$), we would have obtained a new L2 norm distribution (L2') that would overlap with the original one. The resolution capacity of the L2 norm metric in the context of a large number of possible observation series will therefore be defined by the $N_0'$ values below and above $N_0$ having L2 distributions clearly distinguishable from that of L2.

To quantify these upper and lower bound values, we made use of the equation

$$N_t = \left[ (1 + \beta)(1 + r)^t - h \sum_{i=0}^{t-1} (1 + r)^i \right] N_0, \quad (5)$$

which projects the population development as a function of $\beta$, i.e. the fraction of $N_0$ by which the initial population size is increased or decreased. We first assumed that the population size remained constant between $t = 0$ and $t = T$ ($N_t = N_0$), and made 1000 hypothetical observation series {$X_t: t \in \Psi$} based on the uniform distribution $U(N_0 - \alpha N_0, N_0 + \alpha N_0)$ and 1000 series based on the normal distribution $N(N_0, \tau N_0)$. We let $\alpha = 0.5$ and $\tau = 0.3$. These values are, respectively, well above the mean maximum distance relative to the mean (0.34) and the mean coefficient of variation (0.19) of the six reported seen-deer-per-hour observation series (Fig. 2). We included the about 5 % of cases where the sampling gave an observation time series with a trend
(assessed by a Mann-Kendall monotonic trend test). Each hypothetical
observation series was scaled by dividing each sampled value with N₀.
The net (potential) growth rate *r* of an unharvested Norwegian red deer
population is about 0.2 (Langvatn & Loison, 1999). In accordance with
the time series available we let T = 17. Although the time series at hand
do not seem to have been sampled from an uniform distribution, we
included it because we wanted to assess whether our estimation
approach was dependent on the data being normally distributed.

With the above parametrization, we found that the resolution capacity of the L2 norm metric is quite good irrespective of the distribution used to produce the observation series (Fig. S3). Panels (a) (increasing N₀) and (c) (decreasing N₀) are based on the uniform distribution, and panels (b) and (d) are based on the normal distribution. The red line shows the 95 % quantile of the L2 norm distribution when β = 0, and the blue line marks the β value where the 5 % quantile of its associated L2 norm distribution is clearly higher.

When we let the population grow (N₁七年 = 1.5 N₀) or decline (N₁七年 = 0.5 N₀) by adjusting the parameter h, and using U(Nᵢ - αN₀, Nᵢ + αN₀), with α = 0.5, to generate the observation series, we obtained quite similar results (not shown).

3.2. *Estimation error as a function of time series length*

We then probed how well the L2 norm metric would perform in an
estimation setting as a function of time series length and population
development pattern. To this end, we tested three population develop-
ment scenarios (N<sub>T</sub> = N<sub>0</sub>, N<sub>T</sub> = 1.5 N<sub>0</sub>, N<sub>T</sub> = 0.5 N<sub>0</sub>) with T ∈ {10, 20,
30}. For each of the three scenarios and for each T we defined a true
projection vector N<sub>t</sub> t ∈ {0, 1, 2, ..., T}, and generated a large ensemble of
population projections by use of Eqn. (2) and by use of all hypothetical
initial population sizes (N<sub>0</sub>) in the set {N&#771;<sub>0</sub> ∈ Z: 0.25 N<sub>0</sub> ≤ N&#771;<sub>0</sub> ≤ 2 N<sub>0</sub> and
N&#771;<sub>0</sub>/5 ∈ Z}. For Scenario 1 (N<sub>T</sub> = N<sub>0</sub>, T ∈ {10, 20, 30}), we generated
1000 observation series from both uniform and normal distributions as
described above. For Scenarios 2 and 3, we used the uniform distribution
U(N<sub>t</sub> - αN<sub>0</sub>, N<sub>t</sub> + αN<sub>0</sub>) to show that even in these cases we do not need to
assume a normally distributed variable. In concordance with the
observation data (Fig. 2), we now let α = 0.35 and τ = 0.2, and we
ensured that no observation series had a trend according to a Mann-
Kendall trend test. We then calculated the L2 norm for each projection
in each scenario for each T using Eqns. (3) and (4) and found the pro-
jection with the lowest value. The boxplots of the distributions of the
initial population size and the end population size associated with each
scenario (Fig. S4) show that the performance of this estimation method
is, as expected, heavily dependent on the length of the time series. The
Scenario 1 results (Fig. S4, panels (a) and (b)) clearly show that the
performance of the estimation method is insensitive to the type of dis-
tribution used to generate the observation series. We see that the esti-
mation method is unbiased for all scenarios.

The variation in the estimates for T = 10 and T = 20 is due to the
variation in the observation series since Eqn. (2) is deterministic. This
means that the accuracy of the estimation method depends on the
observation series pattern up to a certain T. Since T in our case is 17
years, the population estimates will be exposed to uncertainty even
before we take into account that stochastic processes will influence the
population development. However, even when T = 17, there will be a
wide range of observation series patterns where the L2 norm metric will
work very well. To obtain an indication of whether the reported
observation series were within this range, we used the seen-deer-per-
hour series (combination of farmland and forest observations) of the
Sunndal municipality to estimate N₀ and N_T. The harvest data, the quota
data, the seen-deer-per-hour data, and the results of the KPSS and ADF
tests indicate that the population size has remained fairly constant in the
period 2007–2023, suggesting that r ≅ h. Testing for N₀ values in the
range 1500–2500 when using Eqn. (3) (i.e. N_t = N_0 for all t), and
selecting the population projection with the lowest L2 norm among the
700 projections tested, the estimated N₀ is almost identical to the

assumed N₀, and the estimated N_T is in the range of 5–8 % less than the
assumed N_T, depending on the value of N₀ assumed to be true. Using the
seen-deer-per-hour data from Vestnes municipality (combination of
farmland and forest observations) gave the same result.

Assuming that the observation process over the T years has been
sufficiently standardized to prevent a systematic trend in the observa-
tion series that does not reflect the development of the population (false
trend) and that the time series is of some length, the above analysis
suggests that the method will work sufficiently well as long as the
relationship between X_t and the population size P_t is strong enough to
reveal a trend in {X_t: t ∈ Ψ} if there is a clear overall trend in {P_t: t ∈ Ψ}.

It is beyond the scope of this paper to do an exhaustive analysis of the
types of stochastic process that allow for sufficiently good estimates
when using the L2 norm metric. However, it is worth noting that even in
the case where we let N<sub>T</sub> = 1.5 N<sub>0</sub> and allowed heteroscedasticity, by
letting the uniform distribution be U(N<sub>t</sub> - αN<sub>t</sub>, N<sub>t</sub> + αN<sub>t</sub>), we obtained
unbiased estimates. In this case, the variance of the estimates increased
at T = 10 and T = 20 compared to when we have homoscedasticity
(Fig. S4(c)), but at T = 30 there is negligible variance (not shown). This
suggests that even when there is a monotonic increase in the degree of
heteroscedasticity, the estimation accuracy can be maintained by
extending the length of the time series.

The results suggest that given a low quality population index time
series of sufficient length that fulfils the requirements specified above,
one can estimate the red deer population abundance with reasonable
precision by confronting the population index time series with a large
family of model-generated hypothetical population projections by use of
the L2 norm metric. Most notably, the estimation can be done without
making any assumption about initial population size beyond providing
upper and lower bounds that confine the actual value with a wide
margin. The population projection model will, of course, have to be
biologically sound, but the design of such a model is arguably within
reach when one has access to high-quality harvest data and empirical
data on variation in fertility rates and natural mortality. Since these
three premises were met in our red deer case, we developed a sex- and
stage-structured population model that allowed the generation of a large
family of hypothetical population projections in each of the six munic-
ipalities in the period 2007–2023.

3.3. *Design of the sex and stage-structured population model*

3.3.1. *Outline of the model*

We built a sex and stage-structured red deer population model that
projected the number of individuals in each population category (female
calves, young hinds, adult hinds, male calves, young stags and adult
stags) in each municipality at three time points per year over a period of
17 years. The resolution of the model was deliberately tuned to the
resolution of the data at hand. The three time points each year were ≈
May 15 (spring before the birth peak), August 31 (just before the start of
the hunting season), and December 31 (just after the end of the hunting
season) (Fig. 1(b)). We assumed a 1:1 sex ratio among born calves. A calf
became a young hind or young stag when it turned one year old, and
young hinds and stags became adult hinds and stags when they turned
two years old. The projections were made as a function of the initial
population size, the number of individuals in each category in the initial
population summing up to the initial population size, a stags per hind
threshold value defining stag emigration (see below), the number of
harvested individuals in each category for each year reported by the
municipality under study, the mortality risk during winter for each
category, and the fertility rates of 2-year and ≥ 3-year old hinds. A
formal description of the full projection model containing 6 population
categories x 3 equations per population category respectively covering
the periods t₃ → t₁, t₁ → t₂, t₂ → t₃ = 18 projection equations, is given in
Supplementary material (Table S2), where the number of female calves,
young hinds, adult hinds, male calves, young stags and adult stags at the
three time points each year are given by the equations in rows (1–3),
(4–6), (7–9), (10–12), (13–15) and (16–18), respectively.

In the model, no upper limit was set for how old the hinds and stags
could be. To assess the effect of this, we calculated a survival graph
based on the risk of winter mortality and the probability of being shot
during the hunting season each year. For animals that survived the first
year of life, the analysis showed that the probability of living to 15 years
and 20 years was about 2 % and 0.7 %, respectively. Since these figures
are in agreement with age registration data from Norwegian deer pop-
ulations (Solberg et al., 2017), it suggests that there is no accumulation
of unrealistically old animals in the model.

3.3.2. The set of hypothetical initial population states

We defined a large putative initial total population size range (900–4000 individuals) for each municipality and made a set of initial population sizes in increments of 10 individuals within this range ({ x ∈ Z: 900 ≤ x ≤ 4000 and x /10 ∈ Z}). Based on the reported harvest data and the reported seen-deer-per-hour and spring census data, this population span clearly covered the actual initial population size in all municipalities studied (Fig. 2). The set of possible initial hinds per stag ratios was defined to be { x ∈ R: 1.3 ≤ x ≤ 2.2 and x /0.1 ∈ Z}. This range of values most likely covers the actual initial hinds per stag ratio in all municipalities studied (Fig. 2).

The sets of possible initial proportion of calves in the total popula-
tion, the proportion of young hinds of the total hind population, and the
proportion of young stags in the total stag population were defined to be
{ x ∈ R: 0.2 ≤ x ≤ 0.35 and x /0.01 ∈ Z}, { x ∈ R: 0.09 ≤ x ≤ 0.27 and x
/0.01 ∈ Z} and { x ∈ R: 0.18 ≤ x ≤ 0.54 and x /0.02 ∈ Z}, respectively.
For every combination of elements from these three sets, elements from
the set of initial total population sizes, and the set of hinds per stag ra-
tios, the initial number of individuals in all six categories was uniquely
determined. All possible combinations of the elements from the above
five sets (that is, a full factorial design) defined the set of hypothetical
initial population states for each municipality. It should be emphasized
that beyond assuming that the range of each initial condition parameter
is wide enough to cover the true initial value in each of the six municipal-
ities, we make no assumption about the location of the true value in
this range in terms of a prior statistical distribution.

3.3.3. Allowing for net emigration of stags

The harvest regimes in two municipalities (Averøy and Tingvoll) have been consistently female biased throughout the 17-year period (Fig. 2, left panel column). This bias is not reflected in the reported hinds per stag ratio data (both farmland and forest registrations combined). In fact, instead of being lower, as one would expect, the ratios in these two municipalities are consistently higher than the ratios in the other four municipalities (Fig. 2, right panel column) (see also Discussion, Section 5.3). Considering published data on the winter mortality of hinds and stags (Loison & Langvatn, 1998) and on the sex ratio of calves born under normal conditions (Kruuk et al., 1999; Mysterud et al., 2000), this led us to suspect that there could have been a significantly enhanced net emigration of stags from these two municipalities compared to the other four throughout the observation period. To account for this possibility, we included a stags per hind threshold parameter in the model, expressing that when the ratio between the number of adult stags and the number of adult hinds exceeded a specific value, the excess proportion of stags emigrated during summer. We let the set of possible stags per hind threshold values be { x ∈ R: 0.5 ≤ x ≤ 0.8 and x /0.05 ∈ Z} for all six municipalities and included in addition a scenario with no net emigration of stags (denoted by the set {0}), no matter how large the ratio was between number of stags and hinds in early summer. In this way, we preserved an unbiased approach, while allowing the data to tell whether there was reason to assume that female-biased harvest causes an enhanced net emigration of stags. As many areas experience both an influx and an efflux of stags in connection with regular seasonal migrations (Meisingset et al., 2018), it should be noted that the above mechanism addresses only the annual net emigration. Each of the

hypothetical initial population states was combined with all elements in
the two sets of possible stags per hind threshold values.

3.3.4. Annual winter mortality and fertility parameters

To account for annual variation in winter survival and fertility rates, we let the winter survival probability (Langvatn & Loison, 1999; Loison & Langvatn, 1998) and the fertility rates of hinds (Langvatn et al., 2004; Veiberg et al., 2010) become normally distributed variables. We used the following mean and standard deviation values for winter survival: female calf (0.8 ± 0.04) (Table S2, Eqn. (1); young hind (0.92 ± 0.03) (Eqn. (4); adult hind (0.93 ± 0.03) (Eqn. 7); male calf (0.8 ± 0.04) (Eqn. 10); young stag (0.87 ± 0.04) (Eqn. 13); adult stag (0.93 ± 0.03) (Eqn. 16) (in the very few cases where the obtained values were > 1.0, the values were set to 1.0). The mean and standard deviation values are based on a critical assessment of reported data (Loison & Langvatn, 1998).

In Norway, the proportions of 1-year-old, 2-year-old and > 3-year-
old hinds ovulating are estimated to be 0.5–0.7, 0.95 and close to 1.0,
respectively (Langvatn et al., 2004). Due to lack of data, we did not
explicitly include embryonic mortality and postnatal mortality (during
summer before the start of the hunting season). This could possibly lead
to an overestimation of the number of calves reaching age 1. Never-
theless, we assumed that it would be captured in large part by letting the
mean fertility be somewhat lower than the above ovulation figures and
by allowing a fairly large standard deviation, particularly for the
youngest hinds (2 years old when giving birth). We therefore let the
mean and standard deviation values for calves born per hind be
respectively 0.4 ± 0.05 for 2-year-old hinds and 0.975 ± 0.008 for
hinds ≥ 3 years (Table S2, Eqns. 2 and 11). We used the same mean and
standard deviation values for the six municipalities.

3.3.5. Evaluation of population projections by use of the L2 norm

For each initial condition, we let the model project population
development over the 17-year period (2007–2023), incorporating the
harvest figures reported annually for each population category and
letting the probabilities of winter survival and fertility rates for each
year be sampled from the eight normal distributions described above.
We calculated the L2 norm (Eqn. (3) of each hypothetical population
projection using the seen-deer-per-hour data and the spring census data.
The seen-deer-per-hour data were compared with the predicted total
population size before and after hunting.

3.3.6. Computer implementation of the estimation approach

Data and computer code needed to produce all figures and supple-
mentary figures can be downloaded from the Zenodo repository
(https://zenodo.org/records/15113932). The code was developed in
Python (v. 3.8.5) in a Jupyter lab (v. 3.6.1) environment, a general
framework for interactive computing and data science. It makes exten-
sive use of the data manipulation package Pandas (v. 1.2.1). The repo
contains four scripts that produce hypothetical population projections
from the population model in the period 2007–2023 (based on varying
premises), using a full factorial design of initial parameter values con-
tained in a predefined parameter value hypercube. The scripts select the
hypothetical population projections with the lowest L2 norm values and
store the results in Pandas Dataframes. In addition, the repo contains
five scripts that produce Figs. 2-6 in the main paper, and four scripts that
produce Figs. S1-S5 in the Supplementary material. To reduce
computing time, all four scripts doing the population projection and L2
norm evaluation make use of IPython Parallel (v. 7.20.0) to allow trivial
parallelization.

4. Results

4.1. Population size estimations based on the seen-deer data

Because the predicted total population size before hunting (time T<sub>2</sub>,
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Estimation results based on seen-deer-per-hour data from the six municipalities. **Left panel column:** Reported seen-deer-per-hour data (red dots) and the scaled pre-hunt population sizes (blue dots) of the hypothetical population trajectory with the lowest L2 norm. **Middle panel column:** Boxplots showing the distributions of the predicted stags per hind threshold value when we use the 20, 100 and 1000 population projections with the lowest L2 norm (stags_per_hind_t_20, stags_per_hind_t_100, stags_per_hind_t_1000), and the predicted total population sizes (divided by 1000) at the end of December 2006 (pop_2006_20/1000, pop_2006_100/1000, pop_2006_1000/1000) and the end of December 2023 (pop_2023_20/1000, pop_2023_100/1000, pop_2023_1000/1000). **Right panel column:** Boxplots showing the distribution of the predicted annual percentages of stags emigrating (s_emig.%) from the six municipalities for the 20 population trajectories with the lowest L2 norms. In all municipalities except Surnadal, we used the aggregated seen-deer-per-hour figures where the farmland and forest registrations are combined (see Fig. 2). In the Surnadal case, we used the farmland data only as use of the aggregated farmland and forest data gave highly unrealistic population estimates based on the annual harvest numbers (see Section 5.2).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The effect on population structure by imposing no stag emigration. **Left panel column:** The predicted population sizes (mean ± SEM) of calves, young hinds, adult hinds, young stags and adult stags) at the end of December (post-hunt) each year based on the 20 hypothetical population trajectories for each municipality that gave the lowest L2 norms when compared with the seen-deer-per-hour data and when we allowed for net emigration of stags. **Right panel column:** The predicted population sizes (mean ± SEM) for the same groups at the end of December (post-hunt) each year based on the 20 best hypothetical population trajectories when we did not allow for any net emigration of stags.</figcaption></figure>

Fig. 1(b)) overall gave a lower L2 norm than the total population size after hunting (time T3, Fig. 1(b)), we chose to use the pre-hunt value (time T2) as a predictor in the seen-deer-per-hour case. Comparison of the seen-deer-per-hour data and the scaled pre-hunt population trajectory with the minimum L2 norm for each municipality (Fig. 3, left panel column) shows that in several cases there is a significant distance between the predicted and the observed annual value. Even though the depicted process variability is tuned by the estimation process, it will by no means be comparable to the variation in the observation data which frequently indicate a 50–100 % drop or increase in the population size from one year to the next without any marked change in harvest numbers (Fig. 2). Considering the numerous and uncontrollable factors influencing the magnitude of the annual seen-deer-per-hour figures due to the way they were collected, such a discrepancy may be expected.

In each municipality, the 25 %, 50 % and 75 % quantiles of the predicted initial and end population sizes are practically speaking insensitive to whether we make use of the 20, 100 or 1000 population projections with the lowest L2 norm (Fig. 3, middle panel column).

Considering the range of the initial population sizes (900–4000 individuals) used to produce the hypothetical population trajectories, the spread in the population estimates in terms of distance between the 25 % and 75 % quantiles at each of the two time points are encouragingly low for all municipalities. And the results resonate very well with those obtained by using the one-equation population model (Eqn. (1)).

4.2. Female-biased harvesting induces enhanced stag emigration

The predicted median stags per hind threshold values are remarkably
similar across all six municipalities irrespective of whether we used the
20, 100 or 1000 population projections with the lowest L2 norm (Fig. 3,
middle panel column), suggesting that all six have experienced a net
emigration of stags. However, the predicted percentages of emigrating
stags (Fig. 3, right panel column, boxplots showing the distribution of
the top 20 population trajectories) were clearly associated with the
degree of female-biased harvesting (Fig. 2, left panel column). This
might explain why a systematic female-biased harvesting regime is not
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Assessment of the initial population size estimate. Population projections (dotted cyan and blue lines) based on a deterministic projection model vs the seen-deer-per-hour data (red dots). The projections were scaled such that the start value corresponded with the mean of the seen-deer-per-hour series. In each case, the initial population fractions were identical with the median of the distributions underlying *pop_2006_20/1000* (Fig. 3, middle panel column), the *stags_per_hind* threshold value was identical with the mean of the distribution *stags_per_hind_t_20* (Fig. 3, middle panel column), and the initial population size was either 10% lower than the 25% quantile of the distribution *pop_2006_20/1000* (Fig. 3, middle panel column) (blue dotted line) or identical with the 75% quantile (cyan dotted line).</figcaption></figure>

reflected in the observed hinds per stag ratios (Fig. 2).

When we assumed no net emigration of stags, we obtained highly improbable population development patterns over the 17-year period compared to when we allowed the data to determine the stags per hind threshold value, and particularly so for the two municipalities with the most pronounced female-biased harvest regimes, Averøy and Tingvoll (Fig. 4). It appears challenging to reconcile these patterns with the observation and the harvest data from these two municipalities (Fig. 2), providing further support for the notion that a female-biased harvest regime might cause an enhanced net emigration of stags.

Furthermore, since the predicted mean stags per hind threshold value for each municipality varied slightly around 0.58, it follows that the predicted hinds per stag ratio is approximately 1.7. This estimate is consistent with the reported hinds per stag ratios for four of the six municipalities (Fig. 2, right panel column). However, it is puzzling that the hinds per stag ratios for Averøy and Tingvoll are so consistently high (see Discussion, Section 5.3, for possible explanations).
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Predicted population development 2007–2023 in the six municipalities. The boxplots show the distribution of predicted total pre-hunt population size in each municipality each year when using the 20 hypothetical population projections with the lowest L2 norm values obtained by use of the stochastic population model (Table S2, Fig. 3).</figcaption></figure>

### 4.3. Assessing the effect of stochasticity

Based on the very consistent result presented above, we now fixed the hinds per stag and stags per hind values at 1.7 and 0.58, respectively, and repeated the estimation of initial and end population sizes by employing either a single winter survival and fertility record (as in Fig. 3) or 50 different records, for each initial condition. Use of 50 different records instead of one had a marginal impact on the distributions when using the top 20 projections (Fig. S5). Fixing the hinds per stag and stags per hind values had also overall a negligible effect (Fig. S5).

To assess the effect of using stochastic winter survival and fertility figures versus fixed ones, we estimated the initial and end population sizes by using a deterministic version of the population model where the winter survival and fertility parameter values were identical with the mean values of the original normal distributions (Fig. S5, red-coloured boxplots). Otherwise, the initial parameter value ranges were identical to the original stochastic model. The predicted mean values from the deterministic model are overall quite similar to those obtained by the stochastic model, but the spread of the predicted distributions are expectedly much less. Hence, with the observation data at hand, the main effect of stochasticity appears to be just an increase in the range of initial and end population sizes among the projections that are most compatible with the data.

### 4.4. Estimations based on the spring census data

Among the five municipalities for which we had spring census data, the spring counts for Tingvoll and Sunndal show an increasing trend according to a Mann-Kendall test, while the seen-deer-per-hour data show no trend (Fig. 2). Although the four other spring censuses show no trend, a semi-average analysis (ratio (R) between the mean of the eight last observations and the eight first observations) shows that they are also clearly different from the associated seen-deer-per hour records since the {$R_{seen_deer}$, $R_{spring_census}$} values for the five municipalities are: Tingvoll: {0.97, 1.45}; Surnadal: {0.90, 1.14}; Sunndal: {1.0, 1.72}; Vestnes: {0.98, 1.16}; Lærdal: {1.02, 1.23}. We therefore expected that using the spring census data would give end-of-2023 population estimates that would be higher than the estimates obtained by using the seen-deer-per-hour data. In particular, we expected estimates for Tingvoll and Sunndal to be substantially higher. This turned out to be the case, as the median end population estimates for Tingvoll and Sunndal were around 4000, while the estimates for the three other municipalities were much more similar to those based on the seen-deer-per-hour data. As the end-of-2023 population estimates for Tingvoll and Sunndal were clearly unrealistic considering the quota data and the harvest data (Fig. S2), the current spring census data appear to be less reliable than the seen-deer data.
### 4.5. Assessment of the uncertainty tied to the initial population size estimate

The terms that make up the L2 norm (Eqn. (3)) are power functions, where the contribution of a term to the total sum increases nonlinearly with the distance between the predicted and the observed value. If the hypothetical initial population size deviates much from the actual one, the population model's inherent logic will force a larger and larger deviation for each year between the actual population size and the predicted one. As long as there is a reasonable overall correspondence between the actual population size time series and the reported seen-deer-per-hour time series (see Section 3.2), this dramatically narrows down the set of initial population-size values capable of producing population trajectories that stay close to the observed values across the whole time series.

To illustrate the above reasoning, we first found the median six initial population fraction values and the median stags per hind threshold value among the 20 projections with the lowest L2 norm obtained by use of the original stochastic model (Fig. 3). For each municipality, we then made use of these values in the deterministic projection model described in Section 4.3 to generate two contrasting population trajectories. In one case we let the initial population size be 10 % lower than the 25 % quantile of the distribution pop.2006_20/1000 (Fig. 3). In the other case we let the initial population size be identical with the 75 % quantile of this distribution. We scaled both population projections so that the 2007 prediction corresponded to the mean of the seen-deer-per-hour observation series (red dots) (Fig. 5). In all six cases, the model predicts that the population will either die out (blue) or expand to a size far beyond what is sensible (cyan) (Fig. 5).

In full correspondence with the analytical results (Fig. S4), the results show very clearly the inverse relationship between the length of the time series and the range of initial population sizes that can give rise to population projections compatible with the observation data. The results support the notion that there is quite a narrow range of initial population sizes that can give rise to the reported seen-deer-per-hour data.

### 4.6. Assessment of the sensibility of the predicted 2007–2023 population development

The relationships between the quota and harvest data (Fig. S2) indicate that in the municipalities that are geographically close to each other (Averøy, Tingvoll, Surnadal, Sunndal, and Vestnes), the populations expanded up to 2010 before they declined somewhat. Furthermore, the data indicate that in the last few years the populations in the municipalities Averøy, Tingvoll, Surnadal and Vestnes have declined markedly, while the Lærdal population has remained fairly constant after the 2017 culling.

To test whether our estimation method was capable of recapitulating these patterns, we used the top 20 projections obtained by using the original stochastic model (Table S2, Fig. 3) to display the annual distribution of predicted total pre-hunt population size across the entire time series (Fig. 6). Considering the low quality of the observation series, and that we used a very wide range of initial population structures, the correspondence is quite remarkable. Only the predicted decline in the Sunndal population since 2020 is not well reflected in the quota and harvest data (Fig. S2), although these data suggest that there has been a decline the last couple of years.

## 5. Discussion

### 5.1. Practical applicability of the estimation method

Without access to reasonably good population abundance estimates, it is challenging to formulate harvest regimes that will reliably lead to the fulfilment of a set of management objectives for many game species.

We have shown that noisy but long-term hunter-reported data can be successfully used to estimate the abundances of a game species by simple means without access to detailed individual data. Our approach should be quite easy to exploit in a practical management context for any management body possessing long-term hunter-reported time series data of reasonable quality.

With access to repeatedly updated survey and harvest data, our estimation approach used together with short-term projection models may become a handy device in the toolbox of adaptive management (Nichols, 2021; Walters & Hilborn, 1978). A sex- and stage-structured population model can easily be revised to serve as a tool for projecting short-term population development as a function of the harvest regime, and thus be used to identify a harvest regime that will meet a given set of management objectives. And it can be adapted to local conditions and populations: for instance, if specific information on sex- and stage-specific vital rates (e.g. survival) is available, then this information can be used to feed the population model. Since independently obtained population abundance estimates based on time series data can be fed into such short-term projection models on an annual basis and thus serve as a corrective negative feedback signal, our estimation approach may facilitate more widespread use of mathematically informed population management. Finally, our study suggests that even noisy hunter-reported data, when combined with low-resolution dynamic models, can generate fresh perspectives on the behavioural ecology of a species.

Our results suggest that in a heavily managed red deer population, the end-of-time-series population size estimate is quite robust to annual variation in winter mortality and fertility within the empirically observed bounds. However, it should be noted that the current model does not incorporate density-dependent effects on winter mortality and fertility. Such effects can easily be incorporated if this is deemed necessary. Since the observation data contain as much noise as they do, there is clearly a limit to how well the predicted 17-year population trajectories should track the reported time series data without becoming biologically unrealistic. When we allowed all the winter survival and fertility parameters to have much larger standard deviations than supported by the available data, we obtained a slight reduction in the L2 norms of the 20 best trajectories. However, despite that the results did not differ much from those presented above (Fig. 3), the predicted population trajectories started to show clear signs of overfitting (data not shown). It is therefore clearly a trade-off between obtaining as low L2 norms as possible and biological realism regarding demographic rates when using noisy empirical time series data.

### 5.2. The estimation method is sensitive to systematically unrepresentative data

We used the combined farmland and forest seen-deer data for all municipalities, except for Surnadal, where we used the farmland registrations only. The reason for this was that the combined seen-deer data gave a spuriously large end-population estimate compared to Tingvoll, where the harvest numbers have been substantially higher over the whole time period. In Surnadal, the combined data and the forest data (Fig. 2, second right panel) show a systematic upward trend in the period 2017–2019, while the farmland data do not. The upward trend coincides with the introduction of individual hunter registrations, which introduced the possibility of multiple registrations of the same animal by individual hunters in the same hunting team (Solberg et al., 2022). When the probability of registering deer during forest hunting is as low as it is for Surnadal and the invested effort is so high that it dominates the combined data, this may partly explain the discrepancy between recorded farmland and forest data. Therefore, we chose to use the farmland data for Surnadal. This gave an end-population estimate on par with that of Sunndal (having comparable harvest figures) but substantially lower than that of Tingvoll. This reasoning suggests that the results obtained by our estimation approach for a given municipality may benefit from being evaluated in the context of population abundance
estimates and reported data from other comparable municipalities.

To further substantiate that the introduction of new reporting routines from 2017 may have introduced a systematic change in the seen-deer-per-hour data not reflecting a corresponding change in the actual population numbers, we also pre-analysed a data set (2007–2021) from the municipality Stranda, belonging to the county Møre and Romsdal. Here, the seen-deer-per-hour figures (farmland, forest, and combined) have increased substantially from 2019 onwards. Despite the fact that the harvest figures in 2016 and 2017 were very high and that the harvest numbers in the period 2018–2021 were on par with those before 2016, the combined seen-deer-per-hour figure in 2021 was almost twice as high as in 2018. It seems therefore justified to claim that in some cases the introduction of new reporting routines has caused an unfortunate disruption of the seen-deer time series data (Forsyth et al., 2022).

Another concern is the rapidly increasing use of thermal binoculars that allow a larger proportion of the deer population to be counted in the spring census. This may explain the recent increase in spring census counts for Tingvoll and Sunndal, which is not reflected in their seen-deer data records. The same concern applies to the increased use of thermal binoculars during dusk and dawn hunting, as it will affect the seen-deer-per-hour index. In case, the introduction of thermal binoculars has introduced another unfortunate disruption of the CR time series data.

### 5.3. Considerations on the proposed stag emigration mechanism

Red deer have spread steadily in Norway over the past several decades (Meisingset, 2015). This means that there has been an emigration of both stags and hinds from established populations, even though it is a clear sex-biased emigration toward males (Loe et al., 2009). What triggers this emigration is likely a complicated function of genetics, environment, and population density and –structure (Clutton-Brock et al., 2002; Loe et al., 2009) that we do not fully understand. However, the constancy of the predicted stags per hind threshold figure suggests that we may have captured a key mechanism regulating stag emigration. Assuming this to be the case, it suggests that stags are capable of sensing stag and hind population structure relationships that are at least partly reflected by the relative abundance of stags versus hinds, and that they use this population structure information as a cue for making behavioural decisions tied to prospects for future reproductive success. It follows that if a population is exposed to a highly female-biased harvest regime, the emigration of stags will become accentuated and lead to emigration of young stags that under normal conditions would have stayed.

Other factors influencing stag emigration, related to population density (Clutton-Brock et al., 2002) and grazing opportunities in high-altitude areas during summer (Bischof et al., 2012; Mysterud et al., 2017), could have systematically differed between the municipalities studied. Such differences may have accentuated the large difference in predicted proportions of stags emigrating from municipalities with a female-biased harvest regime compared to those with a 50:50 regime. This issue can be resolved by analysing a much larger set of populations that have sported a female-biased harvest regime over many years and where the municipalities vary substantially in terms of topographical features and population density. Unfortunately, we have not been able to establish such a data set from the CR data. Although we acknowledge this as a shortcoming, we still find it pertinent to point out that in a municipality with a low population density and excellent grazing opportunities for stags at higher altitudes during summer, all stags will benefit from this in terms of growth rate. This implies that none of them will get a competitive edge from this enhanced growth rate as long as they compete for the same hinds. Thus, since stags are most likely able to assess their status relative to others, the decision mechanism driving emigration alluded to above may also in this case cause enhanced stag emigration rates under a female-biased harvesting regime comparable to those areas where the stags are generally smaller due to less favourable grazing conditions. We did not include the possibility of emigration of

hinds in our model, as we think it has been too small to impact our key findings.

As already alluded to, despite that the harvest regimes in the two municipalities Averøy and Tingvoll have been consistently female biased throughout the 17-year period (Fig. 2, left panel column), the reported hinds per stag ratio values (farmland and forest registrations combined) are consistently higher than the ratios in the other municipalities (Fig. 2, right panel column). It may be due to the low quality of the data or that many farmland registrations early in the hunting season dominate (before the rut) more than they do in the other municipalities. Another possibility is that female-biased harvesting triggers the emigration of a higher proportion of stags. This would imply that the female biased harvesting has induced processes that are not captured by the simple mechanism we embedded in the population model, which for example does not account for the effect of genetic variation or nonlinear stimulus-response relationships. Thus, it is conceivable that the stag emigration has been even higher than we predict, which would in case affect the population estimates. However, the fact that the predicted stags per hind threshold value is practically identical in the six municipalities speaks against this interpretation.

### 5.4. Quality assurance of seen-deer data and spring census data

Even though we have established that the current CR data can indeed be used to get narrowly distributed population size estimates, there are limits to how well our estimation method can handle nonrepresentative data. We think that there is ample space for improving data collection procedures in most Norwegian municipalities reporting to the CR and, in particular, the managerial value of high-quality spring census data appears to be underestimated in many cases. Our results will hopefully incentivize authorities to implement routines that will improve the representativeness of the reported data. The extra amount of work associated with this refinement is very modest compared to the resulting gain in the managerial usefulness of the data. To avoid multiple interpretations of what constitutes best practice routines regarding representative collection of seen-deer data and spring census data, it would be preferable that the same set of routines had to be adhered to by all municipalities harbouring a red deer population.

### Credit authorship contribution statement

**Stig W. Omholt:** Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Conceptualization. **Marlène Gamelon:** Writing – review & editing, Writing – original draft, Conceptualization. **Erling L. Meisingset:** Writing – review & editing, Writing – original draft, Data curation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

We thank two anonymous referees for useful comments on an earlier draft. Vebjørn Veiberg, Norwegian Institute for Nature Research (NINA) provided valuable advice in an early phase of the model development, and Atle Mysterud, University of Oslo, provided valuable advice concerning important aspects of the deer population in the Lærdal municipality. We would also like to thank the game managers from the different municipalities that provided crucial data not present in the CR. EM was partly founded by Research Council of Norway project AgriDeer (318575) and MG by the French National Research Agency ANR PURE project (ANR-23-CE02-0028).
# Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ecolind.2025.113398.

## Data availability

Data and Python scripts that produce the results are archived in Zenodo (https://zenodo.org/records/15113932).

## References

Apollonio, M., Belkin, V.V., Borkowski, J., Borodin, O.I., Borowik, T., Cagnacci, F., Danilkin, A.A., Danilov, P.I., Faybich, A., Ferretti, F., Gaillard, J.M., Hayward, M., Heshtaut, P., Heurich, M., Hurynovich, A., Kashtalyan, A., Kerley, G.I.H., Kjellander, P., Kowalczyk, R., Yanuta, G., 2017. Challenges and science-based implications for modern management and conservation of European ungulate populations. *Mammal Research* 62 (3), 209–217. https://doi.org/10.1007/s13364-017-0321-5.

Arnold, T.W., Clark, R.G., Koons, D.N., Schaub, M., 2018. Integrated population models facilitate ecological understanding and improved management decisions. *The Journal of Wildlife Management* 82 (2), 266–274. https://doi.org/10.1002/jwmg.21404.

Auger-Méthé, M., Newman, K., Cole, D., Empacher, F., Gryba, R., King, A.A., Leos-Barajas, V., Mills Flemming, J., Nielsen, A., Petris, G., Thomas, L., 2021. A guide to state-space modeling of ecological time series. *Ecological Monographs* 91 (4), e01470. https://doi.org/10.1002/ecm.1470.

Bischof, R., Loe, L.E., Meisingset, E.L., Zimmermann, B., Van Moorter, B., Mysterud, A., 2012. A Migratory Northern Ungulate in the Pursuit of Spring: Jumping or Surfing the Green Wave? *The American Naturalist* 180 (4), 407–424. https://doi.org/10.1086/667590.

Clutton-Brock, T.H., Coulson, T.N., Milner-Gulland, E.J., Thomson, D., Armstrong, H.M., 2002. Sex differences in emigration and mortality affect optimal management of deer populations. *Nature* 415 (6872), 6872. https://doi.org/10.1038/415633a.

Forsyth, D.M., Comte, S., Davis, N.E., Bengsen, A.J., Côté, S.D., Hewitt, D.G., Morellet, N., Mysterud, A., 2022. Methodology matters when estimating deer abundance: A global systematic review and recommendations for improvements. *The Journal of Wildlife Management* 86 (4), e22207. https://doi.org/10.1002/jwmg.22207.

Frost, F., McCrea, R., King, R., Gimenez, O., Zipkin, E., 2022. Integrated Population Models: Achieving Their Potential. *Journal of Statistical Theory and Practice* 17 (1), 6. https://doi.org/10.1007/s42519-022-00302-7.

Gamelon, M., Nater, C.R., Baubet, É., Besnard, A., Touzot, L., Gaillard, J.-M., Lebreton, J.-D., Gimenez, O., 2021. Efficient use of harvest data: A size-class-structured integrated population model for exploited populations. *Ecography* 44 (9), 1296–1310. https://doi.org/10.1111/ecog.05738.

Godvik, I.M.R., Loe, L.E., Vik, J.O., Veiberg, V., Langvatn, R., Mysterud, A., 2009. Temporal scales, trade-offs, and functional responses in red deer habitat selection. *Ecology* 90 (3), 699–710. https://doi.org/10.1890/08-0576.1.

Husek, J., Boudreau, M.R., Panek, M., 2021. Hunter estimates of game density as a simple and efficient source of information for population monitoring: A comparison to targeted survey methods. *PLOS ONE* 16 (8), e0256580. https://doi.org/10.1371/journal.pone.0256580.

Kéry, M., & Schaub, M. (Eds.). (2012). Front Matter. In *Bayesian Population Analysis using WinBUGS* (pp. i–ii). Academic Press. doi: 10.1016/B978-0-12-387020-9/000024-9.

Kruuk, L.E.B., Clutton-Brock, T.H., Albon, S.D., Pemberton, J.M., Guinness, F.E., 1999. Population density affects sex ratio variation in red deer. *Nature* 399 (6735), 459–461.

Langvatn, R., Loison, A., 1999. Consequences of harvesting on age structure, sex ratio and population dynamics of red deer Cervus elaphus in central Norway. *Wildlife Biology* 5 (4), 213–223. https://doi.org/10.2981/wlb.1999.026.

Langvatn, R., Mysterud, A., Stenseth, N.C., Yoccoz, N.G., 2004. Timing and synchrony of ovulation in red deer constrained by short northern summers. *American Naturalist* 163 (5), 763–772.

Lebreton, J.-D., Burnham, K.P., Clobert, J., Anderson, D.R., 1992. Modeling Survival and Testing Biological Hypotheses Using Marked Animals: A Unified Approach with Case Studies. *Ecological Monographs* 62 (1), 67–118. https://doi.org/10.2307/2937171.

Lee, A.M., Bjørkvoll, E.M., Hansen, B.B., Albon, S.D., Stien, A., Sæther, B.-E., Engen, S., Veiberg, V., Loe, L.E., Grøtan, V., 2015. An integrated population model for a long-lived ungulate: More efficient data use with Bayesian methods. *Oikos* 124 (6), 806–816. https://doi.org/10.1111/oik.01924.

Lieury, N., Gallardo, M., Ponchon, C., Besnard, A., Millon, A., 2015. Relative contribution of local demography and immigration in the recovery of a geographically-isolated population of the endangered Egyptian vulture. *Biological Conservation* 191, 349–356. https://doi.org/10.1016/j.biocon.2015.07.008.

Loe, L.E., Mysterud, A., Veiberg, V., Langvatn, R., 2009. Negative density-dependent emigration of males in an increasing red deer population. *Proceedings of the Royal Society b: Biological Sciences* 276 (1667), 2581–2587. https://doi.org/10.1098/rspb.2009.0224.

Loe, L.E., Rivrud, I.M., Meisingset, E.L., Bøe, S., Hamnes, M., Veiberg, V., Mysterud, A., 2016. Timing of the hunting season as a tool to redistribute harvest of migratory deer across the landscape. *European Journal of Wildlife Research* 62 (3), 315–323. https://doi.org/10.1007/s10344-016-1004-2.

Loison, A., Langvatn, R., 1998. Short- and long-term effects of winter and spring weather on growth and survival of red deer in Norway. *Oecologia* 116 (4), 489–500.

Maunder, M.N., Punt, A.E., 2013. A review of integrated analysis in fisheries stock assessment. *Fisheries Research* 142, 61–74. https://doi.org/10.1016/j.fishres.2012.07.025.

Meisingset, E.L., 2015. Space use of red deer and its implications for management, Doctoral dissertation, University of Oslo. https://www.duo.uio.no/bitstream/handle/10852/48260/1/PhD-Meisingset-DUO.pdf.

Meisingset, E.L., Loe, L.E., Brekkum, Ø., Bischof, R., Rivrud, I.M., Lande, U.S., Zimmermann, B., Veiberg, V., Mysterud, A., 2018. Spatial mismatch between management units and movement ecology of a partially migratory ungulate. *Journal of Applied Ecology* 55 (2), 745–753. https://doi.org/10.1111/1365-2664.13003.

Meriggi, A., Sotti, F., Lamberti, P., Gilio, N., 2009. A review of the methods for monitoring roe deer European populations with particular reference to Italy. *Hystrix, the Italian Journal of Mammalogy* 19 (2). https://doi.org/10.4404/hystrix-19.2-4421.

Morellet, N., Klein, F., Solberg, E., & Andersen, R. (2010). The census and management of populations of ungulates in Europe. *Ungulate Management in Europe: Problems and Practices*, 106–143.

Mysterud, A., Meisingset, E.L., Veiberg, V., Langvatn, R., Solberg, E.J., Loe, L.E., Stenseth, N.C., 2007. Monitoring Population Size of Red Deer Cervus Elaphus: An Evaluation of Two Types of Census Data from Norway. *Wildlife Biology* 13 (3), 285–298. https://doi.org/10.2981/0909-6396(2007)13[285:MPORD]2.0:CO:2.

Mysterud, A., Vike, B.K., Meisingset, E.L., Rivrud, I.M., 2017. The role of landscape characteristics for forage maturation and nutritional benefits of migration in red deer. *Ecology and Evolution* 7 (12), 4448–4455. https://doi.org/10.1002/ece3.3006.

Mysterud, A., Yoccoz, N.G., Stenseth, N.C.H.R., Langvatn, R., 2000. Relationships between sex ratio, climate and density in red deer: The importance of spatial scale. *Journal of Animal Ecology* 69 (6), 959–974. https://doi.org/10.1111/j.1365-2656.2000.00454.x.

Nichols, J. D. (2021). Adaptive management: Making recurrent decisions in the face of uncertainty. In R. Salguero-Gomez & M. Gamelon (Eds.), *Demographic Methods across the Tree of Life* (pp. 313–328). Oxford University Press.

Schaub, M., Abadi, F., 2011. Integrated population models: A novel analysis framework for deeper insights into population dynamics. *Journal of Ornithology* 152 (1), 227–237. https://doi.org/10.1007/s10336-010-0632-7.

Schaub, M., & Kéry, M. M. M. (2022). *Integrated population models: Theory and ecological applications with R and JAGS*. Elsevier, Academic Press.

Sebold, S., Perktold, J., 2010. Statsmodels: Econometric and statistical modeling with python. 9th Python in Science Conference, Austin, 28 June-3 July, 2010, 57-61. https://doi.org/10.25080/Majora-92bf1922-011ScPy.

Solberg, E.J., Strand, O., Veiberg, V., Andersen, R., Heim, M., Rolandsen, C.M., Solem, M. I., Holmstrøm, F., Jordhøy, P., Nilsen, E.B., Granhus, A., Eriksen, R., 2017. Hjortevilt 1991–2016—Oppsummeringsrapport fra Overvåkningsprogrammet for hjortevilt. NINA report 1388. http://hdl.handle.net/11250/2453679.

Solberg, E.J., Veiberg, V., Strand, O., Hansen, B.B., Rolandsen, C.M., Andersen, R., Heim, M., Solem, M.I., Holmstrøm, F., Granhus, A., Eriksen, R., Bøthun, S.W.. 2022. Oppsummeringsrapport fra Overvåkningsprogrammet for hjortevilt 1991-2021. NINA report 2141. https://brage.nina.no/nina-xmlui/handle/11250/3020056.

Tenan, S., Iemma, A., Bragalanti, N., Pedrini, P., De Barba, M., Randi, E., Groff, C., Genovart, M., 2016. Evaluating mortality rates with a novel integrated framework for nonmonogamous species. *Conservation Biology* 30 (6), 1307–1319. https://doi.org/10.1111/cobi.12736.

Van Deelen, T., Etter, D., 2003. Effort and the Functional Response of Deer Hunters. *Human Dimensions of Wildlife* 8 (2), 097–108. https://doi.org/10.1080/10871200304306.

Veiberg, V., Nilsen, E.B., Ueno, M., 2010. Fruamtid forvaltning av norske hjortebestandar—Utfordringar knytt til bestandstettleik og demografi. NINA report 571, 1–40.

Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., van Mulbregt, P., 2020. SciPy 1.0: Fundamental algorithms for scientific computing in Python. *Nature Methods* 17 (3), 261–272. https://doi.org/10.1038/s41592-019-0686-2.

Walters, C.J., Hilborn, R., 1978. Ecological Optimization and Adaptive Management. *Annual Review of Ecology and Systematics* 9 (1), 157–188. https://doi.org/10.1146/annurev.es.09.110178.001105.

Zipkin, E.F., Rossman, S., Yackulic, C.B., Wiens, J.D., Thorson, J.T., Davis, R.J., Grant, E. H.C., 2017. Integrating count and detection—nondetection data to model population dynamics. *Ecology* 98 (6), 1640–1650. https://doi.org/10.1002/ecy.1831.

Zipkin, E.F., Saunders, S.P., 2018. Synthesizing multiple data types for biological conservation using integrated population models. *Biological Conservation* 217, 240–250. https://doi.org/10.1016/j.biocon.2017.10.017.