Research Article Aoife Hurley, Department of Mathematics and Statistics,
University of Limerick, V94 T9PX, Ireland<EMAIL_ADDRESS>Science
Foundation Ireland, Grant/Award Number: 18/CRT/6049

# Spatial Joint-Species N-mixture Models for Multi-Source Observational Data
with Application to Wild Deer Population Abundance in Ireland

Aoife K. Hurley  Department of Mathematics and Statistics, University of
Limerick, V94 T9PX, Ireland  Ruth F. Carden  School of Archaeology, University
College Dublin, D04 F6X4, Ireland  Sally Cook  School of Geography and
Environmental Science, Ulster University, BT52 1SA, United Kingdom  Irish Deer
Commission  Irish Deer Commission, Ireland  Ferdia Marnell  National Parks and
Wildlife Services (NPWS), Dublin, D07 N7CV, Ireland  Pieter A. J. Brama
School of Veterinary Medicine, University College Dublin, D04 F6X4, Ireland
Daniel J. Buckley  National Parks and Wildlife Services (NPWS), Dublin, D07
N7CV, Ireland  James Sweeney  Department of Mathematics and Statistics,
University of Limerick, V94 T9PX, Ireland

###### Abstract

Accurate predictions of the populations and spatial distributions of wild
animal species is critical from a management and conservation perspective.
When overpopulation is observed or suspected, culling of animals is an often
used and frequently controversial management strategy. Consequently, accurate
estimates of animal populations are essential for specifying, monitoring, and
justifying the impact of culling programmes. However, population estimates are
generally collated from disparate data sources and at differing spatial
resolutions. Citizen science projects typically provide point referenced data,
whereas site surveys, hunter reports, and official government data may be
aggregated and released at a small area or regional level. In this article we
present a multi-species N-mixture modelling framework for joint modelling of
species populations in the presence of spatially misaligned multi-source data.
Compared to traditional models for species population estimates which are
carried out for each species separately, our approach provides inferential
advantages from the perspective of making best use of the data available to
hand by modelling correlation in the underlying population distributions. We
illustrate the methodology with application to sparse, spatially misaligned
data for three deer species in Ireland.

###### keywords:

Bayesian hierarchical models, spatially misaligned count data, N-mixture
models, correlated species distribution models, deer abundance models

##  1 Introduction

Abundance data on wild animal species can be collected in numerous ways
including through motion cameras, tracking devices, pellet identification,
capture/mark release and site surveys. Site survey data is generally collected
at a fine spatial resolution, such as within a predetermined area or along a
predefined path, since collection is limited by human sense (vision and
hearing) and the costs of collecting data itself. Citizen science is an
emerging cost-effective method for collecting data on wild animal species
[59], though some projects have poor data practices including a lack of
accuracy, poor spatial or temporal representation, and insufficient sample
size [4]. Multi-species data have many complexities, including imperfect
detection, spatial autocorrelation, between-species correlation, and
preferential sampling [18]. Other potential complexities include data being
collected at different spatial scales, or having spatially misaligned data.

With spatially misaligned data, there is an additional challenge in combining
the data sources. Statistical downscaling approaches are often used in climate
models, where information obtained at lower spatial resolutions, say via
satellites, is extrapolated to much finer spatial scales. There are multiple
methods used for statistical downscaling including regression based and neural
network based methods [7]. [47] relate the issue of misaligned data to
integrated species distribution models, and provide a general solution for the
change-of-support problem. They make predictions at an ecologically relevant
scale of inference by leveraging spatial correlation and repeat observations
at multiple scales. In disease mapping, spatial disaggregation approaches are
used to make fine-scale predictions of disease risk from aggregated response
data by using high-resolution covariate data. [3] use simulation studies to
investigate the performance of disaggregation regression at making fine scale
predictions from various sized aggregated areas with varying levels of data.
These simulation studies, while helpful, possibly do not capture the nuances
of real-life response data. However, disaggregation involves downscaling areal
data to point level data, which does not always align with the available data
for a problem. [25] investigate the local prevalence of opioid misuse in Ohio
using a Bayesian hierarchical spatio-temporal framework. They incorporate
county-level data on opioid outcomes with state-level estimates on prevalence
of opioid misuse to estimate the latent county-level prevalence and counts of
people who misuse opioids.

Accurate estimates of wild animal populations are important for the management
and conservation of species, as well as investigating both the positive and
negative impacts these species have on the environment and ecosystems, and
their ability to harbour and spread diseases. [54] introduced N-mixture models
to model animal populations with point referenced data and imperfect
detection. N-mixture models account for detection and abundance
simultaneously, but require repeated observations at a given site to estimate
detection [28, 5]. However, some drawbacks to using N-mixture models include
issues of identifiability of parameters in models and computational
feasibility when populations are large [34]. [24] expresses the N-mixture
likelihood in closed form using a hypergeometric function, which is
algebraically tractable and computable to a high degree of accuracy, and does
not require the infinite sums seen in [54]. [22] develop a multi-species
N-mixture model framework for the estimation of abundance of rare species with
a focus on improvement in the detection probability. They do this by borrowing
strength from other species in the community and assuming the detection
probabilities are drawn at random from a beta distribution. [40] use a multi-
species expansion of the N-mixture model to estimate both positive and
negative correlations between 10 bird species in Alaska, with data taken from
the North American Breeding Bird Survey. By using their models, they find over
the decade they investigated, that bald eagle abundances increased,
particularly in the area of the Alexander Archipelago, albeit their
observational data consists solely of point referenced data as opposed to data
comprised from several spatial scales. A comprehensive review of N-mixture
models and it’s extensions is provided by [34].

Population estimates may also be used to evaluate and potentially influence
policies on culling programmes. Culling programmes are typically used to
manage species numbers and their negative impacts on habitats, when the
species has no natural predator in the area, or to halt the spread of
diseases. [60] investigated the effectiveness of the preemptive culling of
cattle herds that took place during the 2001 foot-and-mouth outbreak in the
UK. Although some aspects of this policy were controversial, their study
indicates that the adopted strategy was effective. Marine mammal culling
programmes were reviewed by [10], and they evaluated the extent of their
efficacy as a fisheries management measure. In the case of white-tailed deer,
[35] found that localised culling as a disease management strategy, can
maintain low disease prevalence while minimising impacts on recreational
hunting. [37] investigated the impact of declining recreational hunting has on
the population of wild boar across Europe between 1982 and 2012. They note
that the decreasing number of recreational hunters is a cause for concern, and
that reliance on recreational hunting alone is not sufficient to limit the
population growth of wild boar. The contribution of recreational hunters on
the population and abundance of invasive species in Australia was reviewed by
[8]. While they find insufficient evidence to resolve arguments that
contemporary recreational hunting programmes are effective at controlling
introduced mammal populations on public lands, they propose a framework with
various criteria for assessing the effects of recreational hunting on
introduced mammal populations.

In this article we focus on the development of models for estimating the total
population of wild animal species given information presented at two distinct
spatial levels, with spatially misaligned data at both levels. In the
following section we outline the motivating application which is deer
abundance mapping in the Republic of Ireland.

###  1.1 Motivating example

The proposed model is motivated by data collected on three wild deer species
in the Republic of Ireland: fallow (Dama dama), red (Cervus elaphus), and sika
(Cervus nippon). Fallow and sika deer were introduced to the Republic of
Ireland in the 1200’s and in 1860 respectively, while the Co. Kerry population
of red deer are descended from a 5,000 year old introduction other red deer
populations in Ireland are descended from introductions from UK and elsewhere
in the 19th century to modern times [11, 12, 6]. The wild deer population in
the Republic of Ireland has gained national media coverage, primarily due to
the potential for deer to transmit bovine tuberculosis to livestock [45, 30],
the threat posed to drivers [32], as well as their negative impact, both
economically and on biodiversity, on various habitats including conservation
ecosystems, agricultural lands, commercial forestry plantations and
semi/native woodlands [52, 11, 42]. More recently, the record number of culls
from the 2021/2022 season was highlighted [33], which lead to the conclusions
“the overall size of the deer population in Ireland is unknown because no
census of numbers has ever been conducted" [39].

Currently in the Republic of Ireland there is no national deer management
policy and no coordinated system to monitor population or density measurements
[52], though in 2022 the Government of Ireland established a Deer Management
Forum to try to address issues and provide recommendations to the national
management of wild deer. There is an ongoing debate about the reintroduction
of predators, such as wolves, as an alternative to controlling wild deer
numbers, though [26] believes the introduction of wolves in Ireland to be
unlikely. [36] note that a key driver in the debate to restore grey wolves to
the Highlands of Scotland is to control expanding numbers of the native red
deer in the area. Though they do not advocate for either side, they also
examine the ecological and landscape impact the reintroduction of grey wolves
would have. The introduction of the Eurasian lynx has also been muted in
regard to controlling wild deer numbers, though the analysis from [23], there
is a lack of suitable connected forestry and woodlands available in Ireland
currently and thus such an introduction would be unlikely to succeed in the
long term.

Disease transmission between wild deer and livestock in the Republic of
Ireland is an area of growing concern, though has long been speculated in the
transmission of bovine tuberculosis [30]. [20] examined the possibility of
deer being a transmitter of foot-and-mouth disease to livestock in the UK.
While their results show deer can carry the virus, it is unlikely to be an
important factor in the maintenance and transmission of the virus. A
literature review by [9] examined the different transmission routes posed by
deer, in particular the diseases present in the UK. They note that disease
transmission models will improve with better data collection, and the accurate
monitoring for all deer species population sizes. [13] conducted a literature
review of diseases impacting deer and livestock, with a focus on Australia,
with bovine tuberculosis and foot-and-mouth diseases ranking highly in their
risk assessment. Other diseases that ranked highly in their assessment, but
are not present in the Republic of Ireland at present, include malignant
catarrhal fever, surra, and screw-worm fly infestation.

Wild deer can cause a variety of damage, through bark stripping, browsing of
lateral and leader shoots, and adult males thrashing younger aged trees with
their antlers, to growing coniferous and deciduous trees. Such damage may
affect the individual tree’s overall growth, health (pathways for diseases and
pests), and timber quality and yields. There are limited studies on the
economic and biodiversity negative impacts, or damage, to plantation forestry
and semi/native woodlands in Ireland. These have determined estimations of
costs of deer damage to Sitka/other commercial plantations and deciduous
woodlands (commercial and semi/native) annually to the Irish economy at a
national level in the range of c.€1.3 million [42]. However, data was only
available from a few locations across the country and there is an absence of
quantitative baseline data, thus this annual cost estimate may be lower or
higher.

Multiple studies have been conducted regarding wild deer in Ireland including
modelling the distribution and range expansion of wild deer in Ireland [11];
presence-only mapping of wild deer in Ireland [41]; degree of hybridisation
between sika and red deer [38, 56], origins [38, 12, 6]; and the overlap
between wild deer species and endangered or vulnerable plant species [46].
[30] estimate the population of sika deer using hunter returned bag numbers,
but failed to account for uncertainty in the cull percentages themselves.
However, currently there is no modelling approach that accounts for multiple
sources of data to estimate the population of these three species in the
Republic of Ireland. To the best of our knowledge, no such comparable model
exists in any other country either. Accurate estimates of the populations of
the three deer species are vital for investigating forest health, conservation
of habitats, diseases, and disease transfer between wildlife and livestock, as
well as informing policies on culling programmes.

For this problem, data regarding the number of hunting licences granted, the
number of deer culled by recreational hunters and the number of deer observed
was made available (sources: Irish Deer Commission; National Parks and
Wildlife Services). These datasets are spatially misaligned, with the number
of licences granted and the recorded number of deer culled on a county or
regional level, and the observed data at a 10km $\times$ 10km grid resolution.
We also incorporate proportion data on 14 different land types from the Corine
land cover dataset at the 10km $\times$ 10km resolution [19]. Further
information on the available data is presented in Section 3.1.

The main aim of this study is to estimate the populations of the three species
of wild deer in the Republic of Ireland. Other important information such as
the spatial distribution of each species and the impact of the differing land
covers on each species is also of interest. The paper is structured as
follows: in Section 2 we outline the structure of a general model for joint
species population estimation given spatial misalignment in the observational
scales, as well as strategies for parameter inference. Section 3 shows the
model applied to our motivating example, and we examine the sensitivity of the
model to assumed cull percentages. Finally, in Section 4 we discuss our
findings and possible avenues for future work.

##  2 A Joint N-mixture model for species abundance with spatial misalignment
in observational data

###  2.1 Refined spatial scale

Let
$\bm{\mathcal{S}_{1}}=\left\\{\mathcal{S}_{1,1},\dots,\mathcal{S}_{1,n}\right\\}$
represent $n$ areal units at which we wish to estimate species abundance (e.g.
at the finest spatial resolution of available covariates). At this level, the
model uses the N-mixture framework proposed by [54]. Let
$\mathbf{N}_{a}=\left\\{N_{a,1},\dots,N_{a,n}\right\\}$ represent the true
population of species $a$ at each of the locations in $\mathbf{S}_{1}$, and
$\mathbf{y}_{a}=\left\\{y_{a,1},\dots,y_{a,n_{a}}\right\\}$ represent the
$n_{a}$ observations collected for species $a$. Our model at this spatial
level is then represented by

$\begin{split}y_{a,k}&\sim D_{w}(N_{a,k\in n},p_{a,k\in n},[\theta_{w}])\\\
N_{a,i}&\sim D_{b}(\lambda_{a,i},[\theta_{b}])\\\
\text{logit}(p_{a,i})&=\delta_{a,0}+G_{i}^{T}\bm{\delta}_{a}\\\
\log(\lambda_{a,i})&=\beta_{a,0}+X_{i}^{T}\bm{\beta}_{a}+\phi_{a,i}\\\
\end{split}$ (1)

where $D_{b}$ and $D_{w}$ are the probability distributions representing
between- and within-site variation respectively, $k\in n$ identifies the
location where the observation $k$ occurred, $p_{a,i}$ is the detection
probability for each site $i$, $\lambda_{a,i}$ is the mean abundance at site
$i$, $\mathbf{G}$ and $\mathbf{X}$ are the covariate matrices for detection
and abundance, respectively and $\bm{\phi}$ is an areal spatial process.
Similar to the N-mixture model shown in [21], each probability distribution
model ($D_{b}$ and $D_{w}$) can include an overdispersion parameter
($\theta_{b}$ and $\theta_{w}$) allowing the flexibility for various N-mixture
models to be applied. Equation (1) can easily be extended to account for
multiple site visits, which is a requirement in order to estimate site-level
detection probabilities.

###  2.2 Regional spatial scale

The spatial misalignment occurs with the inclusion of data at a second,
coarser spatial scale with $m$ areal units,
$\bm{\mathcal{S}_{2}}=\left\\{\mathcal{S}_{2,1},\dots,\mathcal{S}_{2,m}\right\\}$.
At this level, incompletely observed data
$\mathbf{z}_{a}=\left\\{z_{a,1},\dots,z_{a,m}\right\\}$, is available on the
population of species $a$ within each region. Using elements seen in Section
2.1, our model at this lower resolution spatial level is given by

$\begin{split}R_{a,j}&=\sum_{i\in j}\mathbf{N}_{a,i}\\\ z_{a,j}&\sim
D_{R}(R_{a,j},\kappa_{a,j},[\theta_{R}])\end{split}$ (2)

such that $R_{a,j}$ is the population of species $a$ in region $j$,
$\sum_{i\in j}$ sums the populations at the finer scale that are contained
within region $j$, $D_{R}$ is the probability distribution representing
within-region variation, and $\kappa_{a,j}$ represents the regional detection
level. Similar to Section 2.1, it is possible for the probability distribution
$D_{R}$ to include an overdispersion parameter $\theta_{R}$.

###  2.3 Model fitting and inference

We adopt a Bayesian paradigm, where model specification is complete after
assigning a prior distribution $p(\bm{\Theta})$ for the parameter vector
$\bm{\Theta}$, which represents all unknown model parameters. All models
outlined in this article are fit using the statistical software package `R`
[53] \- we have have good experience in using the `nimble` package [61] for
ecological applications in particular. In our experience, `RStan` is limited
in its capabilities for N-mixture models [58]. This is due to issues
surrounding the binomial $N$ parameter being an integer constrained random
variable, which is a function of an underlying Poisson count generating
process. As noted by [34], the Integrated Nested Laplace Approximation
approach -`INLA` [55] is also limited in its capabilities for N-mixture models
due to the requirement of averaging survey-level covariates to the site or
site-sampling level. Furthermore `INLA` is not a compatible inference approach
for applications involving data collected at an aggregated areal level as is
the case in our deer study. `JAGS` [49] is an alternative inference approach,
however, we encountered issues with parameter mixing in the joint 3 species
model due to issues in specifying Wishart priors. Therefore, we recommend
using `nimble`, which relies on Markov Chain Monte Carlo (MCMC) algorithms to
obtain samples from the posterior distribution. `nimble` adopts a similar
model language to `BUGS` [57] and `JAGS` [49], making it easy to use, and
allows for the customisation of sampling methods. [50] showcase the
customisable nature of `nimble` for a suite of models, while [31] provides a
tutorial on how to use `nimble` for disease mapping, including various spatial
models for both univariate and multivariate cases.

##  3 Case study of wild deer population in the Republic of Ireland

This section begins with a detailed description of the data used for model
fitting, following the introduction given in Section 1.1. We detail the model
for the case of wild deer in the Republic of Ireland in Section 3.2. Following
this, we show the results from fitting this model with 200 different sets of
cull percentages in Section 3.3. We also investigate the sensitivity of our
models to the cull percentages by comparing these results to those generated
from using just the mean cull percentages and changing the cull percentage in
one aggregated area in $\bm{\mathcal{S}_{2}}$ for one species. The main goals
of this analysis are: (a) to estimate, with quantifiable uncertainty, the
population of the three deer species in Ireland, (b) to assess the spatial
distribution of each species and (c) to investigate the impact of land types
on each species’ presence and abundance.

###  3.1 Data

As briefly outlined in Section 1.1, we have data on two distinct spatial
scales: the 26 counties in the Republic of Ireland and 732 10km $\times$ 10km
grid squares across the Republic of Ireland. At the county level,
$\bm{\mathcal{S}_{2}}$, we have data on the bagged returns or number of deer
culled. At our more granular level, $\bm{\mathcal{S}_{1}}$, we have count data
on the number of deer seen in a 10km $\times$ 10km grid during the 2017/2018
season, as well as land cover covariates from the Corine land cover data set
[19].

Private licensed hunters are required to report the number of deer they hunted
in a season to the National Parks and Wildlife Services (NPWS), who collate
this data to a county level noting both the overall number of culled deer and
the breakdown per species. Although these bagged numbers do not include all
deer deaths (e.g. natural causes, poaching, road deaths), using these numbers
has been shown to be a good predictor for deer populations, [1] use the
density of bagged numbers to predict the residual density of white-tailed deer
in Tennessee. Figure 1 visualises the returned bagged numbers for each deer
species in each of the 26 counties in the Republic of Ireland in the 2017/2018
hunting season.

Figure 1: The number of deer bagged by hunters in each of the 26 counties of
the Republic of Ireland during the 2017/2018 hunting season. While red and
sika deer culls were not recorded in every county (due to not being shot by
hunters, this does not preclude the absence of these deer), at least 3 fallow
deer were culled in every county.

Additionally, the number of licences granted within a hunting season at a
county level was also made available for this analysis (source: NPWS). The
number of culled or bagged deer per licence may serve as a proxy measure for
hunting effort, and alongside expert opinion, helps inform the ranges for our
cull percentages. The bagged number of deer per licence for the hunting season
2017/2018 for each of the 26 counties in the Republic of Ireland is shown in
Figure 2.

Figure 2: Bagged number of wild deer per licence for each species in the
hunting season 2017/2018. Similar pattern seen across each of the three
species as in Figure 1.

Data on sightings of these three main wild deer species has been collected at
a 10km $\times$ 10km grid level in the Republic of Ireland (source: Irish Deer
Commission). The sightings information includes repeated site visits where no
deer were spotted, resulting in presence-absence data and the ability to
estimate detection probabilities. The recorded number of visits differs for
each species, with fallow deer having the most number of squares accounted for
(96) and the most number of visits (184), followed by red (a total of 91
squares visited 163 times) and lastly sika deer having the smallest
observational data (61 squares visited a total number of 105 times). Figure 3
highlights the squares visited for each species and the number of times each
square was visited. In total, under 3% of grid squares have observations.

Figure 3: Squares coloured by the frequency of recorded visits.

Completing the data required for the analysis is the Corine Land Cover (CLC)
aggregated to a 10km $\times$ 10km resolution across 732 grid squares. For
this analysis, we used the most recent update to the CLC which was in 2018 and
is freely available to download [19]. The CLC contains three class levels with
each level having a different number of categories, of which we selected 14
different land types across the three levels. We account for the proportion of
each different land type within each of the 732 grid squares. Table 1 shows
the 14 different land types used in the matrix $\mathbf{X}$, as well as the
four groupings used for the matrix $\mathbf{G}$ to model the detection
probabilities. The groupings were decided based on the experience of deer
ecology experts and their informed knowledge from two decades of observations
of wild deer and their natural ecological habitat preferences in Ireland.

Table 1: CLC classification of land cover variables used and the associated groupings of these variables used in modelling the detection probabilities. CLC Classification |  Grouping
---|---
Burnt Areas | Bad
Urban Areas | Bad
Agricultural Land | Agriculture
Orchards | Good
Pastures | Agriculture
Broad Leaved Forests | Good
Coniferous Forests | Coniferous and Peat
Mixed Forests | Good
Natural Grasslands | Agriculture
Peat Bogs | Coniferous and Peat
Transitionary Woodlands with Scrub | Good
Beach | Bad
Marsh | Bad
Sea and Ocean | Bad

###  3.2 Model

In this section, we outline prior choice and other modelling specifications
that were not provided in the general outline shown in Section 2. For species
$a$, let $\mathbf{y}_{a}=\\{y_{a,1},y_{a,2}\dots,y_{a,n_{a}}\\}$ represent the
$n_{a}$ data points observed and
$\mathbf{N}_{a}=\\{N_{a,1},N_{a,2},\dots,N_{a,n}\\}$ represent the unknown
species population at the desired grid scale with $n$ locations, then

$\begin{split}y_{a,k}&\sim\text{Binom}(N_{a,k\in n},p_{a,k\in n})\\\
N_{a,i}&\sim\text{Neg Binom}(d_{a,i},r_{a,i})\text{.}\end{split}$

We assume a negative binomial distribution for the population counts of each
species, allowing for overdispersion in the observed noisy counts. We also
parameterise the negative binomial such that the mean abundance is
$\bm{\lambda}$. The rest of the grid-level model is given by

$\begin{split}\text{logit}(p_{a,i})&=\delta_{a,0}+G_{i}^{T}\bm{\delta}_{a}\\\
\log(\lambda_{a,i})&=\beta_{a,0}+X_{i}^{T}\bm{\beta}_{a}+\phi_{a,i}\\\
\bm{\phi}&=\Sigma\otimes\mathbf{u}\\\
\mathbf{u}_{a}&\sim\text{MVN}(0,[\tau_{a}(\mathbf{D}-\mathbf{A})]^{-1})\text{,}\end{split}$
(3)

where the covariate matrices, $\mathbf{X}$ and $\mathbf{G}$ are outlined in
Table 1 and $\bm{\phi}$ is a multivariate correlated spatial surface across
the 10km $\times$ 10km grid for the three deer species. We assign to
$\bm{\phi}$ a multivariate intrinsic conditional autoregressive (MICAR) prior
as seen in [27]. We select an ICAR rather than a proper conditional
autoregressive (PCAR) prior, as it allows for computational conveniences
including the use of Kronecker product to reduce computation time.

For both covariate matrices, $\mathbf{X}$ and $\mathbf{G}$, there is a unit
sum constraint on the rows as they consist of proportions data. Although this
is a natural and obvious property, it leads to difficulties for modelling
purposes due to these matrices not having full rank. As our data has
meaningful zeros, the log-ratio transformations originally proposed by [2] are
not applicable. We impose the Bayesian lasso penalty [48] on $\bm{\beta}$, the
coefficients of $\mathbf{X}$. An alternative to using the Bayesian lasso would
be to use principal component analysis, but this lacks the interpretability we
desire for this model.

At our regional level,

$\begin{split}R_{a,j}&=\sum_{i\in j}\mathbf{N}_{a,i}\\\
z_{a,j}&\sim\text{Poisson}(\kappa_{a,j}R_{a,j})\text{, }\end{split}$ (4)

where $z_{a,j}$ represents the number of culled deer for species $a$ in county
$j$, and $\kappa_{a,j}$ indicates the cull percentage of species $a$ in county
$j$. We use a Poisson approximation to the binomial distribution to evaluate
the likelihood for the observed cull numbers to overcome mixing issues in the
sampling chains observed in other model parameters. Grid squares are assigned
to counties based on the location of the centroid of the square

Due to confounding and identifiability issues between detection probabilities
at the 10km $\times$ 10km grid level and cull percentages at the county level,
we sample cull values for each species in each county based on expert
ecological opinion. Using the trend seen in the number of culled deer per
licence, as seen in Figure 4, we group the 26 counties into three distinct
bands per species. The appropriate intervals for each of the three bands were
suggested by deer ecologists, and are classed as Low (between 5% and 15%), Mid
(between 10% and 20%), and High (between 20% and 30%).

Figure 4: The number of culled deer per licence by species over 5 years.

When the culling percentage is believed to be below that of the reproductive
rate of the herd and the number of bagged deer per licence has increased, the
cull percentage is sampled from $\text{Low}\sim\text{Normal}(10,2.55)$ \-
roughly corresponding to a cull rate between 5% and 15%. An example of a
county classified as having a low culling percentage is County Donegal for red
deer, as over the course of 5 years, the culled number per licence has not
decreased at any stage. Counties where the reproductive rate and cull
percentage are similar, thus maintaining the size of the herd, are grouped
into the mid band, where $\text{Mid}\sim\text{Normal}(15,2.55)$ \- roughly
corresponding to a cull rate between 10% and 20%. Kerry is classified as
having mid culling percentages for sika deer, as there is little variation
between the start and end of the 5 year period. Finally, counties where
indications of reducing herd size are shown, are classified as counties with
high cull percentages, such that $\text{High}\sim\text{Normal}(25,2.55)$ \-
corresponding to a cull rate between 20% and 30%. A high culling percentage is
only assigned for fallow deer in counties Laois, Longford and Sligo, due to
the decreasing trend seen across the 5 year period. Table 2 gives the assigned
band by species to each of the 26 counties in the Republic of Ireland.

Table 2: Classification for each cull percentage by county associated with the three wild deer species in the Republic of Ireland. County |  Fallow |  Red |  Sika
---|---|---|---
Carlow | Mid | Mid | Low
Cavan | Low | Mid | Mid
Clare | Low | Low | Mid
Cork | Mid | Mid | Low
Donegal | Mid | Low | Mid
Dublin | Low | Mid | Low
Galway | Mid | Mid | Low
Kerry | Mid | Mid | Mid
Kildare | Mid | Low | Low
Kilkenny | Low | Mid | Mid
Laois | High | Mid | Mid
Leitrim | Low | Low | Mid
Limerick | Mid | Mid | Mid
Longford | High | Mid | Mid
Louth | Mid | Mid | Mid
Mayo | Mid | Low | Mid
Meath | Mid | Low | Mid
Monaghan | Low | Mid | Low
Offaly | Low | Mid | Low
Roscommon | Mid | Low | Mid
Sligo | High | Mid | Low
Tipperary | Low | Mid | Low
Waterford | Low | Mid | Low
Westmeath | Low | Mid | Mid
Wexford | Mid | Mid | Low
Wicklow | Low | Mid | Low

Lastly, to ensure the mean abundance of $N_{a,i}$ is $\lambda_{a,i}$, the
parameters $d_{a,i}$ and$r_{a,i}$ are modelled as

$\begin{split}d_{a,i}&=\frac{\nu_{a}}{\nu_{a}+\lambda_{a,i}}\\\
r_{a,i}&=\frac{\lambda_{a,i}d_{a,i}}{1-d_{a,i}}\text{, }\end{split}$

such that $\bm{\nu}$ controls overdispersion and the mean abundance of
$\mathbf{N}$ is $\bm{\lambda}$. Our remaining parameters for the fitted model
are contained in the vector
$\bm{\Theta}=\\{\bm{\beta_{0}},\bm{\beta},\bm{\tau},\bm{\delta},\bm{\theta}\\}$.
For the coefficients $\bm{\beta}$, we impose a Laplace(1,1) prior following
that of the Bayesian Lasso penalty outlined by [48] to alleviate the effects
of confounding between variables. The spatial precision parameter $\bm{\tau}$
is assigned a Gamma(shape = 10, scale = 0.5) prior, and $\bm{\nu}$, has an
Exponential(0.5) prior assigned. We assign a zero-mean Gaussian prior with a
standard deviation of 5 to the coefficients of the detection probability
$\bm{\delta}$ and species-specific intercept term $\bm{\beta_{0}}$. The values
for the components of the prior distributions were chosen after conducting
simulation studies to assess their impact on parameter estimates.

###  3.3 Results

This section interprets the results of running the model with 200 different
cull scenarios, with culling percentages sampled from the intervals previously
outlined in Section 3.1 and Table 2. Guided by simulation studies, we generate
200 cull scenarios as the credible intervals for the population estimates do
not change to any significant degree with the consideration of additional cull
scenarios. Each scenario contains a cull percentage for each of the 26
counties and each of the three species, resulting in 78 (26 $\times$ 3)
different cull percentages per scenario. Unlike the work conducted by [30], in
which only one fixed cull percentage is assumed at the county level, sampling
the cull percentages adds additional uncertainty in population estimates,
which reflects the uncertainty in the cull percentages themselves. For each
scenario, we generate 1,000 independent samples from the posterior
distribution for each model parameter, resulting in our analysis being
conducted on 200,000 samples from the posterior distribution in total. These
samples were generated using 10 chains run in parallel, each generating 100
independent samples post burn-in, taking approximately 80 minutes to run.
Posterior checks were carried out to ensure independence of the samples and
sufficient burn-in. The the total run-time for the analysis was approximately
11 days.

Table 3 shows the estimated species populations by county for the Republic of
Ireland, and estimated total populations. Across the Republic of Ireland, sika
deer are estimated to be the most populous species of deer (median estimate of
185,056), with the median estimated population of fallow deer lower at
155,964. Unsurprisingly, the native red deer have the smallest population with
a median estimate of 39,211. Red deer almost disappeared in the
20${}^{\text{th}}$ century, but due to rigorous protection and management,
their population has increased [12]. Figure 5 illustrates the estimated median
population for each species by county. Deer population estimates are also
available, but not presented here, at the 10 km square grid level.

Figure 5: Estimated median populations for each deer species shown on a map of
the Republic of Ireland.

Overall, the estimated population numbers from the model are aligned to what
we know about presence distribution of wild deer and anecdotal observations
from deer hunters, farmers, and conservation personnel, and reports on damage
caused by over localised populations of wild deer in many areas to various
land uses and conservation habitats and vehicular accidents involving deer on
many primary roads and secondary roads.

Fallow deer are the most widely distributed species in the Republic of
Ireland, as observed in Figure 5, and there has been little change since the
2008 survey by [11], though they are expected to increase their total range
distribution at relatively low rates (3% compound annual rate of expansion).
However, the herd may build up to high densities but remain in the area of
their original release, thus creating large populations beyond the carrying
capacity of the habitats. The landscape of the largely limestone-based Irish
lowlands is preferred and well suited to fallow deer. The highest
concentration of estimated numbers the model produced was in County Tipperary,
which is a county based on limestone and of enriched agricultural crops,
including grasslands for domestic livestock grazing and broadleaf
forestry/woodlands. The estimates presented here follow similar patterns of
presence as per [11], therefore indicating preferred habitats for foraging and
reproduction has not changed much in the previous 10 years. However, they have
expanded their range to the eastern midlands and this may be due to higher
densities being reached in areas.

Although red deer can expand their range at high rates (7% annually over a
30-year period [11]), red deer are the least widely distributed species, being
limited mainly to County Kerry (Killarney National Park and immediate
surrounds), counties Galway and Mayo, and Glenveagh National Park in County
Donegal. There are smaller outlier populations of red deer to the known
distribution hot spots across Ireland. However, there are no red deer left in
County Wicklow due to hybridisation with sika since 2008. Red deer have a
natural preference for deciduous broadleaf woodlands, but have been pushed out
into more open upland habitats by humans, thus limiting their population
through this ecological adaptation across the mountain regions. In the upland
peat areas in the west and northwest of the Republic of Ireland, there are
large commercial plantations of Sitka spruce, where red deer find shelter and
forage. They may negatively impact on the trees through browsing of lateral
and leader shoots of saplings and immature trees and bark stripping of semi-
mature trees.

It has been shown in a previous 30-year study that sika in Ireland expand
their total range by 5% per annum [11]. However, by far the most human-deer
conflict issues with sika occur in the east of the country in County Wicklow,
where anecdotes and evidence show they are causing damage to conservation
habitats, commercial plantation coniferous and broadleaf forestry and to
agricultural crops. It’s unsurprising that the model has identified County
Wicklow as it is a hot spot of high numbers of sika. There are large tracts of
uplands, which sika have adapted to, of peat bog and mountains with large
scale Sitka/other commercial forestry plantations in which they use for
shelter and foraging, as well foraging on agricultural crop lands and
grasslands adjacent to the plantations. Through much of County Wicklow, due to
a large issue of illegal shooting at night under a lamp, or even in daylight
hours, (poaching), the deer have adapted a more nocturnal lifestyle, emerging
from dense forest to feed in dark hours, and thus managing these herds by the
recreational hunters during legal shooting daylight hours has become quite
difficult, if impossible in many areas. A combination of dense habitats,
inaccessible uplands, the presence of a national park where recreational
hunters cannot shoot in, and the opportunistic ecological behaviour of sika
may have led to the overabundance of sika numbers in County Wicklow relative
to other counties where sika are found in much smaller estimated population
numbers.

Table 3: Population for each species of deer given by county and across the Republic of Ireland. Included in the table are the estimated median and the associated 95% credible intervals. County |  Fallow |  Red |  Sika
---|---|---|---
Carlow | 175 [106 - 297] | 11 [1 - 36] | 4,051 [2,591 - 9,100]
Cavan | 1,862 [1,197 - 4,331] | 61 [29 - 119] | 40 [17 - 84]
Clare | 20,349 [13,900 - 34,128] | 3,289 [2,171 - 7,063] | 322 [218 - 522]
Cork | 4,552 [3,445 - 6,992] | 340 [225 - 532] | 8,215 [5,565 - 15,314]
Donegal | 268 [176 - 428] | 11,037 [7,806 - 22,851] | 425 [289 - 705]
Dublin | 2,614 [1,650 - 4,947] | 21 [5 - 55] | 7,495 [5,156 - 16,306]
Galway | 15,873 [12,059 - 23,085] | 1,490 [1,072 - 2,275] | 872 [543 - 1,565]
Kerry | 414 [218 - 640] | 428 [365 - 524] | 13,036 [9,49970 - 18,737]
Kildare | 451 [330 - 641] | 1,191 [726 - 2,275] | 1,302 [855 - 2,786]
Kilkenny | 2,080 [1,362 - 3,959] | 305 [202 - 487] | 398 [278 - 637]
Laois | 4,974 [4,173 - 6,157] | 225 [147 - 356] | 425 [291 - 637]
Leitrim | 959 [608 - 1,833] | 642 [410 - 1,429] | 333 [219 - 534]
Limerick | 1,302 [950 - 2,100] | 620 [438 - 978] | 251 [165 - 404]
Longford | 225 [169 - 306] | 14 [3 - 43] | 9 [1 - 29]
Louth | 99 [68 - 147] | 521 [362 - 780] | 10 [2 - 33]
Mayo | 5,272 [3,976 - 7,649] | 11,665 [9,917 - 15,681] | 390 [260 - 623]
Meath | 188 [119 - 319] | 2,328 [1,494 - 4,709] | 88 [47 - 162]
Monaghan | 2,804 [1,924 - 6,369] | 68 [33 - 131] | 129 [61 - 312]
Offaly | 12,405 [8,103 - 23,123] | 234 [148 - 372] | 537 [328 - 1,239]
Roscommon | 2,562 [1,902 - 3,788] | 454 [267 - 1,007] | 194 [122 - 328]
Sligo | 2,348 [1,882 - 2,927] | 1,238 [865 - 1,759] | 990 [641 - 1,881]
Tipperary | 37,358 [24,705 - 78,751] | 1,083 [779, 1,548] | 2,037 [1,328 -4,348]
Waterford | 20,692 [14,011 - 38,657] | 107 [59, 193] | 1,806 [1,131 - 3,602]
Westmeath | 1,427 [921 - 3,393] | 149 [90 - 250] | 80 [42 - 145]
Wexford | 109 [71 - 174] | 91 [58 - 150] | 2,537 [1,542 - 5,302]
Wicklow | 8,791 [5,557 - 16,335] | 541 [371 - 875] | 136,130 [90,023 - 252,436]
Total | 155,964 [133,669 - 213,053] | 39,211 [34,033 - 52,264] | 185,056 [135,793 - 298,827]

Figure 6 illustrates the spatial surface for each species from the MICAR
component $\bm{\phi}$. The top row of Figure 6 shows the mean estimate for
each surface, and the lower row is the associated standard deviation. The
additional markings on the standard deviation maps indicate squares that have
had visits.

Figure 6: Estimated mean spatial surface for each wild species of deer and the
associated standard deviations at each grid square. The yellow triangles on
the standard deviation map indicates grid squares where we have at least one
recording from our observational data.

The pattern seen across all three species is similar to the presence-only maps
produced by [41] and the binary species distribution maps seen in [11]. The
hot-spots seen for all three species also align with expert opinion, in
particular the known hot spots for sika in the east, southwest and northwest,
with smaller population outliers elsewhere.

In Figure 7 the estimates for each coefficient of the 14 chosen CLC land
coverings, along with the associated 95% credible intervals are presented. We
observe that of the 42 estimated $\beta$’s associated with the CLC land types,
only two credible intervals do not contain zero. Unsurprisingly, Pastures is
positive for fallow deer. Fallow deer are a lowland species, preferring a
mosaic of open woodland and pasture to shelter and forage within. Due to the
forced migration of red deer to the uplands and bog lands across the Republic
of Ireland, ecologically the positive association between red deer and Peat
Bogs is tenable. Due to the sparse nature of the data, there is a large degree
of uncertainty in the parameter estimates for covariates. Given additional
observational data for each species, we would expect the credible intervals in
Figure 7 to increase in precision. We observe large, positive, median
estimates for the regression parameters corresponding to Agricultural Land,
Pastures and Coniferous Forests for red and sika populations. Based on expert
ecological opinion, these correspond to the preferred feeding habitats of
these species as well as reflecting the migration of these species to areas of
lower human habitation and uplands as represented by the percentage of
coniferous forest plantations. The large positive coefficient for Sea and
Ocean observed for fallow is an artefact of a small population based on Lambay
island off the coast of County Dublin.

Figure 7: Estimated coefficients for land covering for each species.

Figure 8 shows the estimated values for the coefficients for the grouped land
coverings that inform our detection probabilities. The impact of a lower
covariate coefficient is shown when the percentage of the associated land
grouping is high, by lowering the associated probability of detection. The
four groups are classified as Bad, Good, Agriculture, and Coniferous and Peat
(as shown in Table 1), while also including an intercept term. The Good
coefficient for red deer is quite negative, suggesting that the higher the
percentage of good land in a square, the lower the probability of detection of
a red deer is for a given site visit, provided they are present. This is
presumably a reflection of the small numbers of red deer with habitats of this
type. Conversely, the Coniferous and Peat coefficient for red deer is
positive, indicating that the higher the percentage of coniferous forests and
peat bogs in a square, the higher the probability of detecting a red deer is
on a given site visit. This is presumably a reflection of their grazing habits
- conversely the negative coefficient for sika deer is an artefact of changing
grazing patterns given hunting and poaching pressures. There may be an element
of slight statistical confounding between the estimated abundance and presence
parameters though substantially more data would be required to establish this.

Figure 8: Coefficients for detection probabilities. The lower the value of
$\delta$, the smaller the probability of detection.

###  3.4 Sensitivity of population estimates to cull parameters

In this section we explore the impact of the selected cull percentage
intervals on our results. Section 3.4.1 compares these results to analysis
when only the mean cull percentage is used, similar to [30], and Section 3.4.2
shows the impact of changing the assumed cull percentage in one county on the
sika population results.

####  3.4.1 Sensitivity of results to ignoring cull uncertainty at the county
level

Table 4 shows the difference in estimated wild deer populations when we
compare the sampled cull percentage results to those when uncertainty in cull
parameters is ignored. While the median estimated populations are relatively
similar, the range of values covered by the 95% credible intervals is
considerably larger when accounting for more uncertainty in cull rates. Across
all three species, the upper end of the 95% credible interval increases
significantly. This additional range accounts for the uncertainty we have in
our cull percentages.

Table 4: Comparison of wild deer population estimates and associated 95% credible interval for 200 sampled cull rates and mean cull rate. Cull percentages |  Fallow |  Red |  Sika
---|---|---|---
200 Sampled | 155,964 [133,669 - 213,054] | 39,211 [34,033 - 52,264] | 185,056 [135,793 - 298,827]
Mean | 150,384 [147,956 - 153,505] | 38,444 [37,339 - 39,555] | 180,706 [178,010 - 183,883]

From the formulation of the MICAR model for the correlated spatial
distributions, we can calculate the residual correlation between the species
that is included in the spatial surface. These correlation estimates account
for the unexplained or left over correlation, after we take all other
variables in the model, such as land type, into account. Figure 9 shows the
estimates and associated 95% credible intervals for the three between-species
correlations for both the 200 sampled cull percentages and when we ignore cull
uncertainty. The variability differences in the two plots is clear, with the
credible intervals being narrower when cull uncertainty is ignored. Expert
ecological opinion is reflected in the positive between-species correlation
for red and sika deer. Red and sika deer can exist on overlapping habitats, as
observed in the hybridisation of these species in County Wicklow. However,
there is a great deal of uncertainty in these between-species correlation
estimates due to the low quality of data available for model fitting.

Figure 9: Estimated between-species correlations, comparing between 200
sampled cull percentages and using the mean cull percentages. Sampling 200
different cull percentages increases the variability and reduces the median
estimate towards zero.

####  3.4.2 Sensitivity to choice of cull ranges

Sika deer were originally introduced to the Republic of Ireland in 1860 to the
Powerscourt Estate in County Wicklow [51]. From the NPWS hunting returns for
2017/18 deer season, a total of 13,502 sika were shot in County Wicklow. For
the 2017/2018 hunting season, there was 652 licences granted in County
Wicklow, indicates that just under 21 sika were shot per deer hunter.
Naturally, the proportion of the total culled number varies between each
hunter. Most of the 652 hunters shooting on average between 0 and 10 deer in
the season, mainly for home consumption of venison; there are a few deer
hunters associated with estates that would shoot over 80 deer in a season. It
is difficult to even determine if the population of sika in the county is
decreasing, increasing or is being maintained. We have included the cull
percentage range for sika deer in Wicklow to be in the Low grouping, such that
the herd population is expanding. As it is possible the population of sika
deer in Wicklow is being maintained, we assess the impact of changing the mean
cull percentage to that of the Mid interval. That is, we investigate changing
the mean cull percentage for sika in Wicklow is changed from 10% to 15%.

Table 5 shows population estimates when the mean cull percentage for sika in
Wicklow is changed to 15% from 10%. With this change, we see a dramatic
decrease in the estimated population of sika deer, with more minor changes for
both fallow and red estimates. Increasing the mean cull percentage from 10% to
15% may account for part of the natural mortality numbers, but not the high
levels of reproduction in yearling and adult female sika.

Table 5: Comparison of wild deer population estimates and associated 95% credible interval for mean cull rate and a change in one cull rate for sika deer (changed from 15% to 10%). Cull percentages |  Fallow |  Red |  Sika
---|---|---|---
Mean | 150,384 [147,956 - 153,505] | 38,444 [37,339 - 39,555] | 180,706 [178,010 - 183,883]
Altered | 150,441 [148,110 - 153,575] | 38,454 [37,414 - 39,534] | 135,683 [133,756 - 138,419]

County Wicklow is atypical relative to other counties in Ireland due to the
landscape features and the large tracts and density of forestry plantations.
Wicklow is the highest forested county amongst all the counties [16], with
37,470 ha of forested areas in lowland and upland areas, with agricultural
pasture and cereal crops in between. However, there are other forms of
evidence that can be taken into account regarding estimates of population
numbers in this county and if occurring at high densities: (i) the amount of
reported and documented adverse impacts and damage to agricultural crops and
pasture (sika feed on spring grass before livestock is turned out after
spending winter indoors), (ii) very little or non-existent natural
regeneration of flora and trees in semi/native woodlands and broadleaf
forests, (iii) as well as levels of observed damage to commercial forestry
timber yields, would lead us to ascertain that the sika population in County
Wicklow is not being maintained or decreasing but is increasing year-on-year
[52, 42]. The level of hunting, recreationally or professionally, is
inadequate. Indeed, the model median estimates either at 10% or at a 15% cull
rate, still estimates an overabundance of sika in County Wicklow, whereby deer
densities for favourable timber yields and reduction of damage to a tolerable
level is 2.5 deer per km2 [14]; County Wicklow covers an area of 2,027km2,
thus 5,068 deer for the entire county is the ideal number for tolerable levels
of damage and for natural regeneration of flora to occur.

###  3.5 Ecological perspective and population estimate validation

Based on research from the UK, the wild deer population in the UK is estimated
to be between 650,000 and 2 million [15]. This is based on a max culled value
of 350,000, representing a cull rate of between approximately 60% and 15% of
the total population per annum. Scaling this estimate to the Republic of
Ireland, where approximately 42,000 deer were culled in the hunting season
2017/2018, we would expect the range of total deer populations to be between
78,000 and 240,000 based on these UK cull percentages. However, our estimates
from Section 3.3 (between 323,866 and 510,049 for the total population of all
deer species) do not align with these values. This could be due to the
granularity of the source data and deer survey range with the count estimates
adjusted to higher than real estimates due to the difficulty in counting and
recording deer in dense habitats such as thicket forestry. Furthermore, expert
ecological opinion is that the cull numbers in the Republic of Ireland are
substantially underestimated due to poaching and illegal hunting, and thus not
reflected in official statistics.

We derive further confidence in our estimates from historical literature. In
2000, [43] estimated the population of sika deer in Wicklow to be 19,215. With
a conservative 16% increase of population numbers (sika increase their
population numbers between 16-21% annually [29]) year-on-year to 2018, this
yields an estimated population of 277,116 sika for 2017/18. However, this
estimate does not take into account the annual culling rate which varies but
increases from 2000 to 2018. Our model may even be too conservative, but until
such time that there are more accurate count numbers and associated data from
County Wicklow for sika, we cannot elucidate accurate numbers, though we based
our model assumptions on expert opinion and known data available.

The sika deer estimate does skew theestimates of the total deer population
numbers in our model but the sika in Wicklow have a documented high fertility
rate, 87% of yearling females are pregnant and 85% of adult females are
pregnant and may reproduce one calf a year up to ages 12 years and older. Very
high survival rates of calves (85%) also compound higher estimated numbers of
the sika population in Wicklow [43]. Similar high rates of fertility,
pregnancy and calf survival was also recorded in Killarney National Park in
the past [44]. With an annual increase of population numbers of between
16-21%, it is feasible that sika numbers are very high in 2017/18 and within
our model’s ranges, based on initial population estimates of just over 19,000
sika in Wicklow in 2007.

##  4 Discussion

In this article we have introduced a novel N-mixture modelling framework for
joint modelling of multiple species populations for ecological data presented
at multiple differing spatial scales. The sparse data provided from spatially
replicated hunter visits on a 10km $\times$ 10km grid level are used to
effectively estimate population sizes at that spatial resolution through the
targeted downscaling of the aggregated deer cull data available at the county
level, while also accounting for the detection and culling process. Our model
has the attractive aspect of allowing for estimation of inter-species
abundance correlations, ensuring as much information as possible is extracted
from the sparse datasets typically gathered by and available to ecological
researchers. All sources of uncertainty on model parameters are coherently
modelled using a Bayesian framework, and incorporating best expert opinion
from deer ecologists. In terms of contributions to the ecological literature,
this article presents the most complete proof of concept approach to the
estimation of the populations of the three main species of wild deer in the
Republic of Ireland given currently available data, and represents the first
time their populations have been estimated with quantifiable uncertainty.
These population estimates are of interest to many stakeholders including
conservationists, forest managers and policy makers, and our population
estimates align with the expectations, or are deemed as plausible, by
ecological experts within the Republic of Ireland. This research also provides
for a framework for assessing the impact of differing control strategies to
limit deer numbers where overpopulation exists, as is suspected to be the case
for sika deer numbers in the Wicklow county region. Due to the density of sika
deer observed in County Wicklow, an increase in the culling percentage may be
required to curtail damages and bring the population back to a manageable
number.

Finally, we remark on some general research directions. It is straightforward
to extend the presented model to account for additional sources of
information. Data used by [41], including two surveys conducted by the Coillte
forestry agency, were used in their presence-only analysis: (a) a survey
completed by property managers on the presence of deer, and the species in the
forests they maintain and (b) deer density surveys based on faecal pellet
sampling in a subset of their properties between the years 2007 and 2020.
While the data collected by Coillte appears to be point referenced data, this
can be aggregated to the spatial levels used within this study for seamless
integration. The data presented in this paper also typify those collected in
many ecological studies and provide an avenue for the borrowing of information
across multiple species datasets where inter-species correlation on abundances
provide an approach to reduce uncertainties in population estimates. It would
be of interest to extend the approach to other non-ruminant animals, say to
avian settings or similar. Our model can also be extended to include temporal
data, say through an autoregressive process, though this avenue was not
applicable with the motivating data as the 10km $\times$ 10km observational
data is only available for the 2017/2018 hunting season. An additional future
development of the model could be explore the impact of possible preferential
sampling on species estimates [17], as the observational data on deer is
collected in areas where deer are likely to be observed. The results of this
research will be used to further develop strategies in future surveys, while
also allowing informed discussions on wider species management and
conservation programmes to manage the wild deer populations within the
Republic of Ireland.

## Acknowledgements

This publication has emanated from research conducted with the financial
support of Science Foundation Ireland under Grant number 18/CRT/6049. For the
purpose of Open Access, the author has applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising from this
submission.

## ORCID

Aoife K. Hurley | 0000-0002-0407-8724
---|---
Ruth F. Carden | 0000-0002-2829-4667
Sally Cook | 0000-0002-8287-4523
Pieter A. J. Brama | 0000-0003-3405-0551
James Sweeney | 0000-0002-5649-3233

## References

  * Adams et al. [2020] Adams, H. L., Kissell, R. E., Ratajczak, D., Warr, E. L., Applegate, R. D., Barrett, L., Lavacot, T. and Graves, D. (2020) Relationships among white-tailed deer density, harvest, and landscape metrics in TN, USA.  European Journal of Wildlife Research, 66, 19.
  * Aitchison [1982] Aitchison, J. (1982) The statistical analysis of compositional data.  Journal of the Royal Statistical Society: Series B (Methodological), 44, 139–160.
  * Arambepola et al. [2022] Arambepola, R., Lucas, T. C. D., Nandi, A. K., Gething, P. W. and Cameron, E. (2022) A simulation study of disaggregation regression for spatial disease mapping.  Statistics in Medicine, 41, 1–16.
  * Balázs et al. [2021] Balázs, B., Mooney, P., Nováková, E., Bastin, L. and Jokar Arsanjani, J. (2021) Data Quality in Citizen Science, 139–157.  Cham: Springer International Publishing.
  * Barker and MacIsaac [2022] Barker, J. R. and MacIsaac, H. J. (2022) Species distribution models applied to mosquitoes: Use, quality assessment, and recommendations for best practice.  Ecological Modelling, 472, 110073.
  * Beglane et al. [2018] Beglane, F., Baker, K., Carden, R. F., Hoelzel, A. R., Lamb, A. L., Fhionnghaile, R. M., Miller, H. and Sykes, N. (2018) Ireland’s fallow deer: their historical, archaeological and biomolecular records.  Proceedings of the Royal Irish Academy: Archaeology, Culture, History, Literature, 118C, 141–165.  URL: https://www.jstor.org/stable/10.3318/priac.2018.118.01.
  * Benestad et al. [2008] Benestad, R. E., Chen, D. and Hanssen-Bauer, I. (2008) Empirical-statistical downscaling.  World Scientific Publishing Company.
  * Bengsen and Sparkes [2016] Bengsen, A. J. and Sparkes, J. (2016) Can recreational hunting contribute to pest mammal control on public land in Australia?  Mammal Review, 46, 297–310.
  * Böhm et al. [2007] Böhm, M., White, P. C. L., Chambers, J., Smith, L. and Hutchings, M. R. (2007) Wild deer as a source of infection for livestock and humans in the uk.  The Veterinary Journal, 174, 260–276.
  * Bowen and Lidgard [2013] Bowen, W. D. and Lidgard, D. (2013) Marine mammal culling programs: review of effects on predator and prey populations.  Mammal Review, 43, 207–220.
  * Carden et al. [2011] Carden, R. F., Carlin, C. M., Marnell, F., Mcelholm, D., Hetherington, J. and Gammell, M. P. (2011) Distribution and range expansion of deer in Ireland.  Mammal Review, 41, 313–325.
  * Carden et al. [2012] Carden, R. F., McDevitt, A. D., Zachos, F. E., Woodman, P. C., O’Toole, P., Rose, H., Monaghan, N. T., Campana, M. G., Bradley, D. G. and Edwards, C. J. (2012) Phylogeographic, ancient DNA, fossil and morphometric analyses reveal ancient and modern introductions of a large mammal: the complex case of red deer (Cervus elaphus) in Ireland.  Quaternary Science Reviews, 42, 74–84.
  * Cripps et al. [2019] Cripps, J. K., Pacioni, C., Scroggie, M. P., Woolnough, A. P. and Ramsey, D. S. L. (2019) Introduced deer and their potential role in disease transmission to livestock in Australia.  Mammal Review, 49, 60–77.
  * De Nahlik [1992] De Nahlik, A. J. (1992) Management of deer and their habitat: principles and methods. Wilson Hunt.
  * Department for Environment , Food & Rural Affairs [2022] Department for Environment , Food & Rural Affairs (2022) Consultation on the proposed deer management strategy.  Tech. rep., Department for Environment , Food & Rural Affairs.  URL: https://consult.defra.gov.uk/team-trees/consultation-on-proposals-for-the-deer-management/supporting_documents/Deer%20management%20strategy%20consultation%20.pdf.
  * Department of Agriculture, Food & the Marine [2022] Department of Agriculture, Food & the Marine (2022) Forest Statistics Ireland 2022.  Tech. rep., Department of Agriculture, Food & the Marine.
  * Diggle et al. [2010] Diggle, P. J., Menezes, R. and Su, T.-l. (2010) Geostatistical inference under preferential sampling.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 59, 191–232.
  * Doser et al. [2022] Doser, J. W., Finley, A. O. and Banerjee, S. (2022) Joint species distribution models with imperfect detection for high-dimensional spatial data.  arXiv preprint arXiv:2204.02707.
  * European Union [2018] European Union (2018) Copernicus Land Monitoring Service 2018, European Environment Agency (EEA).  URL: https://land.copernicus.eu/pan-european/corine-land-cover/clc2018?tab=download%20(CLC).
  * Gibbs et al. [1975] Gibbs, E., Herniman, K., Lawman, M. and Sellers, R. (1975) Foot-and-mouth disease in British deer: transmission of virus to cattle, sheep and deer.  The Veterinary Record, 96, 558–563.
  * Goldstein and de Valpine [2022] Goldstein, B. R. and de Valpine, P. (2022) Comparing N-mixture models and GLMMs for relative abundance estimation in a citizen science dataset.  Scientific Reports, 12, 12276.
  * Gomez et al. [2018] Gomez, J. P., Robinson, S. K., Blackburn, J. K. and Ponciano, J. (2018) An efficient extension of N-mixture models for multi-species abundance estimation.  Methods in Ecology and Evolution, 9, 340–353.
  * Guilfoyle et al. [2023] Guilfoyle, C., Wilson-Parr, R. and O’Brien, J. (2023) Assessing the ecological suitability of the Irish landscape for the Eurasian lynx (Lynx lynx).  Mammal Research, 68, 151–166.  URL: https://doi.org/10.1007/s13364-022-00670-2.
  * Haines [2016] Haines, L. M. (2016) Maximum likelihood estimation for n-mixture models.  Biometrics, 72, 1235–1245.
  * Hepler et al. [2023] Hepler, S. A., Kline, D. M., Bonny, A., McKnight, E. and Waller, L. A. (2023) An integrated abundance model for estimating county-level prevalence of opioid misuse in Ohio.  Journal of the Royal Statistical Society Series A: Statistics in Society, 186, 43–60.
  * Hickey [2011] Hickey, K. (2011) Wolves in Ireland: A Natural and Cultural History.  Four Courts Press.
  * Jin et al. [2007] Jin, X., Banerjee, S. and Carlin, B. P. (2007) Order-free co-regionalized areal data models with application to multiple-disease mapping.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 69, 817–838.
  * Joseph et al. [2009] Joseph, L. N., Elkin, C., Martin, T. G. and Possingham, H. P. (2009) Modeling abundance using n-mixture models: the importance of considering ecological mechanisms.  Ecological Applications, 19, 631–642.
  * Kaji et al. [2004] Kaji, K., Okada, H., Yamanaka, M., Matsuda, H. and Yabe, T. (2004) Irruption of a colonizing sika deer population.  The Journal of Wildlife Management, 68, 889–899.
  * Kelly et al. [2021] Kelly, D. J., Mullen, E. and Good, M. (2021) Bovine Tuberculosis: The Emergence of a New Wildlife Maintenance Host in Ireland.  Frontiers in Veterinary Science, 8.
  * Lawson [2020] Lawson, A. B. (2020) NIMBLE for Bayesian Disease Mapping.  Spatial and Spatio-temporal Epidemiology, 33, 100323\.
  * Liu et al. [2018] Liu, Y., Nieuwenhuis, M. and McCullagh, A. (2018) The effect of roadside land-use on the occurrence of deer vehicle collisions.  Irish Forestry, 75, 8–25.
  * Lucey [2023] Lucey, A. (2023) Record number of wild deer shot in 2021/2022 season.  The Irish Times.  URL: www.irishtimes.com/ireland/2023/01/29/record-number-of-wild-deer-shot-in-20212022-season/.
  * Madsen and Royle [2023] Madsen, L. and Royle, J. A. (2023) A review of N-mixture models.  WIREs Computational Statistics, e1625.
  * Manjerovic et al. [2014] Manjerovic, M. B., Green, M. L., Mateus-Pinilla, N. and Novakofski, J. (2014) The importance of localized culling in stabilizing chronic wasting disease prevalence in white-tailed deer populations.  Preventive Veterinary Medicine, 113, 139–145.
  * Manning et al. [2009] Manning, A. D., Gordon, I. J. and Ripple, W. J. (2009) Restoring landscapes of fear with wolves in the scottish highlands.  Biological Conservation, 142, 2314–2321.
  * Massei et al. [2015] Massei, G., Kindberg, J., Licoppe, A., Gačić, D., Šprem, N., Kamler, J., Baubet, E., Hohmann, U., Monaco, A., Ozoliņš, J., Cellina, S., Podgórski, T., Fonseca, C., Markov, N., Pokorny, B., Rosell, C. and Náhlik, A. (2015) Wild boar populations up, numbers of hunters down? a review of trends and implications for europe.  Pest Management Science, 71, 492–500.
  * McDevitt et al. [2009] McDevitt, A. D., Edwards, C. J., O’Toole, P., O’Sullivan, P., O’Reilly, C. and Carden, R. F. (2009) Genetic structure of, and hybridisation between, red (Cervus elaphus) and sika (Cervus nippon) deer in Ireland.  Mammalian Biology, 74, 263–273.  URL: https://doi.org/10.1016/j.mambio.2009.03.015.
  * McGee [2023] McGee, H. (2023) Calls for larger deer culls after record 55,000 shot dead last year.  The Irish Times.  URL: www.irishtimes.com/ireland/2023/05/08/calls-for-larger-deer-culls-after-record-55000-shot-dead-last-year/.
  * Mimnagh et al. [2022] Mimnagh, N., Parnell, A., Prado, E. and Moral, R. d. A. (2022) Bayesian multi-species N-mixture models for unmarked animal communities.  Environmental and Ecological Statistics, 29, 755–778.
  * Morera-Pujol et al. [2023] Morera-Pujol, V., Mostert, P. S., Murphy, K. J., Burkitt, T., Coad, B., McMahon, B. J., Nieuwenhuis, M., Morelle, K., Ward, A. I. and Ciuti, S. (2023) Bayesian species distribution models integrate presence-only and presence–absence data to predict deer distribution and relative abundance.  Ecography, 2023, e06451.
  * Murphy et al. [2013] Murphy, V., Carden, R. F., Harrison, S., O’Halloran, J., Irwin, S. and Butler, F. (2013) Deer in Irish commercial forests.  Irish Forestry, 91–103.
  * O’Brien et al. [2007] O’Brien, D. J., Rooney, S. M. and Hayden, T. J. (2007) Reproduction and potential rate of increase of the sika deer herd in Co. Wicklow.  Irish Forestry.
  * O’Donoghue [1991] O’Donoghue, Y. (1991) Growth, reproduction and survival in a feral population of Japanese sika deer (Cervus nippon nippon Temminck, 1836).  Ph.D. thesis, National University of Ireland, Dublin.
  * O’Halloran [2015] O’Halloran, M. (2015) Deer potentially causing TB in cattle, says Coveney.  The Irish Times.  URL: https://www.irishtimes.com/news/ireland/irish-news/deer-potentially-causing-tb-in-cattle-says-coveney-1.2215847.
  * O’Mahony et al. [2023] O’Mahony, J., Vanmechelen, A. and Holloway, P. (2023) Quantifying the distribution and potential biotic interactions between deer and flora using species distribution modelling.  Annals of GIS, 1–16.
  * Pacifici et al. [2019] Pacifici, K., Reich, B. J., Miller, D. A. W. and Pease, B. S. (2019) Resolving misaligned spatial data with integrated species distribution models.  Ecology, 100, e02709.
  * Park and Casella [2008] Park, T. and Casella, G. (2008) The bayesian lasso.  Journal of the American Statistical Association, 103, 681–686.
  * Plummer [2003] Plummer, M. (2003) JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.  In Proceedings of the 3rd International Workshop on Distributed Statistical Computing.
  * Ponisio et al. [2020] Ponisio, L. C., de Valpine, P., Michaud, N. and Turek, D. (2020) One size does not fit all: Customizing mcmc methods for hierarchical models using nimble.  Ecology and Evolution, 10, 2385–2416.
  * Powerscourt [1884] Powerscourt, V. (1884) On the acclimatization of the Japanese deer at Powerscourt.  Proceedings of the Zoological Society of London, 52, 207–209.
  * Purser et al. [2009] Purser, P., Wilson, F. and Carden, R. F. (2009) Deer and forestry in Ireland: A review of current status and management requirements.  Tech. rep. URL: https://www.woodlandsofireland.com/wp-content/uploads/DeerStrategy-1.pdf.
  * R Core Team [2021] R Core Team (2021) R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.  URL: https://www.R-project.org/.
  * Royle [2004] Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts.  Biometrics, 60, 108–115.
  * Rue et al. [2009] Rue, H., Martino, S. and Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71, 319–392.
  * Smith et al. [2014] Smith, S. L., Carden, R. F., Coad, B., Birkitt, T. and Pemberton, J. M. (2014) A survey of the hybridisation status of Cervus deer species on the island of Ireland.  Conservation Genetics, 15, 823–835.  URL: https://doi.org/10.1007/s10592-014-0582-3.
  * Spiegelhalter et al. [2007] Spiegelhalter, D., Thomas, A., Best, N. and Lunn, D. (2007) Openbugs user manual.  Version, 3, 2007.
  * Stan Development Team [2023] Stan Development Team (2023) RStan: the R interface to Stan.  URL: https://mc-stan.org/.  R package version 2.21.8.
  * Sun et al. [2021] Sun, C. C., Hurst, J. E. and Fuller, A. K. (2021) Citizen Science Data Collection for Integrated Wildlife Population Analyses.  Frontiers in Ecology and Evolution, 9.
  * Tildesley et al. [2009] Tildesley, M. J., Bessell, P. R., Keeling, M. J. and Woolhouse, M. E. J. (2009) The role of pre-emptive culling in the control of foot-and-mouth disease.  Proceedings of the Royal Society B: Biological Sciences, 276, 3239–3248.
  * de Valpine et al. [2017] de Valpine, P., Turek, D., Paciorek, C. J., Anderson-Bergman, C., Lang, D. T. and Bodik, R. (2017) Programming with models: writing statistical algorithms for general model structures with nimble.  Journal of Computational and Graphical Statistics, 26, 403–413.