# Spatio-temporal Occupancy Models with INLA

Jafet Belmonta1∗, Sara Martinob∗, Janine Illiana, Håvard Ruec
aSchool of Mathematics and Statistics, University of Glasgow, Glasgow, U.K.
bDepartment of Mathematical Scienqces, NTNU, Trondheim, Norway.
cKing Abdullah University of Science and Technology, Thuwal, SA

1Corresponding author: Jafet Belmont<EMAIL_ADDRESS>
∗Equally contributing authors

###### Abstract

Modern methods for quantifying, predicting, and mapping species distribution
have played a crucial part in biodiversity conservation. Occupancy models have
become a popular choice for analyzing species occurrence data due to their
ability to separate out the observational error induced by imperfect detection
of a species, and the sources of bias affecting the occupancy process.
However, the spatial and temporal variation in occupancy that is not accounted
for by environmental covariates is often ignored or modelled through simple
spatial structures as the computational costs of fitting explicit spatio-
temporal models is too high. In this work, we demonstrate how INLA may be used
to fit complex spatio-temporal occupancy models and how the R-INLA package can
provide a user-friendly interface to make such complex models available to
users.

We show how occupancy models, provided some simplification on the detection
process is assumed, can be framed as latent Gaussian models and, as such,
benefit from the powerful INLA machinery. A large selection of complex
modelling features, and random effect models, including spatio-temporal
models, have already been implemented in R-INLA. These also become available
for occupancy models, providing the user with an efficient, reliable and
flexible toolbox.

We illustrate how INLA provides a flexible and computationally efficient
framework for developing and fitting complex occupancy models using two
complex case studies. Through these, we show how different spatio-temporal
models that include spatial-varying trends, smooth terms, and spatio-temporal
random effects can be fitted to aggregated detection/non-detection data. At
the cost of limiting the complexity of the detection model structure, INLA can
incorporate a range of rather complex structures in the ecological process of
interest and hence, extend the functionality of occupancy models.

The limitations of occupancy models in terms of scalability for large spatio-
temporal data sets remains a challenge and an active area of research. INLA-
based occupancy models provide an alternative inferential and computational
framework to fit complex spatiotemporal occupancy models. The need for new and
more flexible computationally efficient approaches to fit such models makes
INLA an attractive option for addressing complex ecological problems, and a
promising area of research.

##  1 Introduction

The growing concern of the impact of anthropogenic actions on the environment
such as habitat loss, climate change, and land use modification, and of how
the distribution of species might be altered as a result of these has spurred
the development of a wide range of methods for modelling and predicting
species presence and abundance in space (and time) Elith and Leathwick (2009).
The resulting species distribution models may be used to estimate population
sizes for conservation purposes or to understand habitat preferences and
predict species presence or abundance in areas that have not been surveyed or
into the future.

Species distribution models are based on data collected on the presence (or
abundance) of species in space and time. Depending on the species and system
of interest, available resources and other practical issues such as
accessibility, data collection methodology varies among surveys. Different
data collection methods result in different data structures and in turn,
require different modelling approaches. However, one thing that most of these
surveys have in common is the problem of imperfect detection, resulting from
the fact that some individuals may be missed as they are elusive or otherwise
hard to detect, or only temporarily absent. To avoid bias, imperfect detection
has to be accounted for, and detection probabilities have to be estimated as
part of the modelling process (Cressie et al., 2009).

Occupancy data are data that result from a sampling approach that is
relatively simple to implement (as opposed to, e.g., mark capture-recapture)
but still provides a means of accounting for imperfect detection. Here,
specific sites are visited by observers repeatedly over time and at each
instance a species’ presence and absence is recorded. Occupancy models are
then used to model jointly both the specie presence in space-time and relative
to relevant covariates and the detection probabilities, which also might
depend on covariates MacKenzie et al. (2002); Tyre et al. (2003); Cressie et
al. (2009); Dorazio and Royle (2005); Dorazio et al. (2006); Kéry et al.
(2010); Kéry and Royle (2020).

A wide range of occupancy models have been developed and applied, from simple
single species models (MacKenzie et al., 2002; Tyre et al., 2003) to complex
spatio-temporal models (Rushing et al., 2019; Diana et al., 2022), and have
proven to be a flexible tool that ecologists and conservationists use to study
species distributions. Increasingly complex modelling approaches and
associated model fitting software have been used to fit these to data.
Initially, classical maximum likelihood approaches were a popular choice, in
particular during the early 2000’s (Kellner and Swihart, 2014), for example
using packages such as PRESENCE (Hines, 2006) or unmarked (Fiske and Chandler,
2011; Guillera-Arroita, 2017)). These libraries account for imperfect
detection and allow the use of random effects, but lack functionality to
account for spatial and temporal autocorrelation beyond the dependence on
covariates. With increasing model complexity, Bayesian modelling approaches
have become more common, facilitating useful model expansions, such as
incorporating spatial or temporal dependence in the description of the
processes governing occupancy and detection (Devarajan et al., 2020).
Accounting for spatial dependence, in particular, is of crucial importance
when modelling the spatial distribution of species is of interest since
ignoring that measurements taken at nearby locations are dependent can, in
fact, produce biased occupancy estimates (Wright et al., 2019). Since spatial
dependence may be the result of complex unobserved ecological processes
driving species distribution that covariates cannot explain it is necessary to
include model structures that reflect spatio-temporal dynamics to accounts for
autocorrelation not explained by covariates, and that reflect ecological
processes such as dispersal limitation or site fidelity, ultimately leading to
increased model complexity.

Fitting spatially explicit occupancy models requires powerful computational
machinery. Spatially explicit occupancy models can be implemented in widely-
used software packages such as NIMBLE (de Valpine et al., 2017) or Stan
(Carpenter et al., 2017), allowing spatial random effects to be modelled via
conditional autoregressive models (CAR models), where the choice of
neighborhood structure can have a large impact on the results. Despite these
packages have proven to be a flexible tool for Bayesian analyses, they were
not specifically designed for occupancy models, hence users typically
experience long running-times when fitting such models (Clark and Altwegg,
2023). Moreover, while these programming languages are very flexible, they are
often not optimized for efficient computation involving dense covariance
matrices and also may require relatively advanced programming skills. On the
other hand, bespoke packages such as stocc Johnson (2021), ubms (Kellner,
2021), hSDM (Vieilledent, 2019) and Rcppocc (Clark and Altwegg, 2019) provide
a user-friendly interface for fitting spatial occupancy models and can
implement areal spatial models using (intrinsic) CAR models. However, these
packages cannot fit multivariate models and may not be adequate for large
datasets (e.g. data collected at several thousand locations). Recent attempts
to make spatial occupancy models more efficient have focused on using Gaussian
Process approximations, for example, the package spOccupancy (Doser et al.,
2022) uses Nearest Neighbor Gaussian Processes (NNGP) (Diana et al., 2022) to
allow the user to fit models to large spatial datasets.

These bespoke implementations use MCMC algorithms to generate posterior
samples, but are typically limited to the supported model types each package
provides. As a result, users cannot customize their model structure unless
they code their own MCMC algorithm (Ponisio et al., 2020). In addition, MCMC
algorithms can be particularly inefficient for spatial and spatio-temporal
(Blangiardo and Cameletti, 2015; Jona Lasinio et al., 2013). Moreover, large
binary occurrence data in occupancy models can make convergence and chain
mixing difficult, even for non-explicit spatial models (Northrup and Gerber,
2018).

An alternative approach to MCMC-based algorithms is Integrated Nested Laplace
Approximation (INLA) introduced in Rue et al. (2009). Unlike MCMC, which can
virtually be applied to any Bayesian model, INLA is a more specialized tool
and only applies to the, albeit very wide, class of Latent Gaussian models
(LGMs). A LGM consists of three elements: a likelihood model, a latent
Gaussian field and a vector of non-Gaussian hyperparameters. The data are
considered conditionally independent given the latent field, and hence a
univariate likelihood model describes the marginal distribution of the
observations. As in the generalized linear model framework, the mean (or
another measure of central tendency) of the observations is linked to a
Gaussian linear predictor through a known link function. The core of LGMs lies
in its second element: spatio-temporal structures (both discrete and
continuous), covariate effects (either linear and non-linear), and other
random effects can all be elements of the latent Gaussian field. Generalized
linear mixed models (GLMM), generalized additive mixed models (GAMM), spatial
and spatio-temporal models are all examples of LGM Rue et al. (2009).

INLA exploits the Gaussianity of the latent field and uses sparse methods for
Gaussian Markov random fields (GRMF) (Rue and Held, 2005) that allow for fast
computations even for large models. Together with the so called SPDE approach
(Lindgren et al., 2011), which represents spatially continuous Mátern Gaussian
field as a GMRF, INLA has become a very popular tool that has facilitated the
increased use of Bayesian modelling within the community of applied users (for
a recent review see Rue et al. (2017); Bakka et al. (2018))

The popularity of the INLA/SPDE approach stems from three sources: the high
computational efficiency of the algorithm, the flexibility that the LGM
framework offers, and the fact that the INLA algorithm is implemented in the
R-library R-INLA. Thanks to these, INLA has made complex Bayesian modeling
more easily available to practitioners. In recent years, many ecological
studies have used the INLA machinery, and its flexibility and ease of
implementation has facilitated the fitting of complex spatio-temporal model
(Paradinas et al., 2015; Arce Guillen et al., 2023), point processes
(Williamson et al., 2022; Laxton et al., 2023), as well as joint models where
different data sources are merged together Martino et al. (2021); this
includes integrated species distribution models (IDMs) (Simmonds et al., 2020;
Paradinas et al., 2023; Altwegg and Nichols, 2019).

Occupancy models have, until now, been excluded from the class of models
amenable to INLA because of the specific form of the likelihood function,
expressed as a function of both occupancy and detection probabilities.

In this paper we show that, if we are willing to assume a simplified form for
one of the two unknown probabilities, the occupancy model can be written in
the form of an LGM and, therefore, benefit from the entire existing INLA
machinery. The advantage is that a large class of random effects, already
implemented in R-INLA, also become available for occupancy model. The
practitioner no longer needs to either implement bespoke MCMC algorithms for
each new study or to be limited to the model choices made by developers of
specific occupancy modelling software.

An added advantage of fitting model within the INLA framework is that measures
of model fit such as DIC, WAIC and marginal likelihood (Spiegelhalter et al.,
2002; Watanabe, 2013) are readily available. In addition, measures based on
leave one out cross-validation (LOOCV) and the recently introduced measures
based on leave group out cross-validation (LGOCV) can also be computed without
the need to refit the model. The LGOCV is a novel idea introduced in Liu and
Rue (2022) and is a better option than LOOCV for evaluating the predictive
performance of models incorporating structured random effects (Adin et al.,
2023).

In this work, we illustrate, through case studies, how INLA can accommodate a
wide variety of already existing occupancy models and extend some of these to
more complex scenarios. Using this approach can be useful both as a
theoretical framework that can cover many of the specific models and as a
practical tool for practitioners that can relate to a single software for
their analysis.

##  2 Statistical models

In ecological studies, information on the distribution of a species in space
may be gleaned by the presence or absence of individuals of the species across
sites. Whilst presence can be confirmed by the detection of individuals, non-
detections may be the result of either the species being undetected during the
survey or being truly absent at a site. Thus, occupancy models account for
non-detection of species even though it is present at a site, by jointly
modelling the species’ occurrence and detection processes through occupancy
and detection probabilities, respectively (MacKenzie et al., 2002). In order
for these probabilities to be estimated, sites are surveyed at multiple
occasions, across which the true occupancy state is assumed to remain
constant. For a set of $M$ discrete sites, occupancy data are structured in
terms of the detection history of whether a species is detected or not at a
given site across multiple visits, i.e. $\bm{d}_{i}=d_{i1},\ldots,d_{iK}$ for
the $i$-th site with $K$ visits. Thus, the joint likelihood of occupancy
states and detection through time is:

$\displaystyle
L(\bm{\psi},\bm{p}|\bm{d}_{1},\ldots,\bm{d}_{M})=\prod_{i=1}^{m}\left[\psi_{i}\prod_{j}^{K_{i}}p_{ij}^{d_{ij}}(1-d_{ij})^{1-d_{ij}}\right]\prod_{i=m+1}^{M}\left[\psi_{i}\prod_{j}^{K_{i}}(1-p_{ij})+(1-\psi_{i})\right],$
(1)

where $\bm{\psi}$ and $\bm{p}$ are the occupancy and detection probabilities,
respectively. The first part of Equation (1) refers to those $m$ sites where
the species is detected at least once throughout the $K$ visits. The second
part of Equation (1) corresponds to the remaining sites, where the species is
not detected at all (i.e. $\sum_{j}y_{ij}=0$). If detection probabilities are
assumed constant over the $K$ visits (i.e. $p_{ij}=p_{i}\ \forall\ j$), the
observational model can then be seen as an aggregated process across the
visits. Under this assumption, let $y_{i}=\sum_{j}d_{ij}$, i.e. the total
number of times the species was detected; now the likelihood becomes:

$\displaystyle L(\bm{\psi},\bm{p}|\bm{y})$
$\displaystyle=\left[\prod_{i=1}^{m}\psi_{i}p_{i}^{y_{i}}(1-p_{i})^{K_{i}-y_{i}}\right]\left[\prod_{i=m+1}^{M}\psi_{i}(1-p_{i})^{K_{i}}+(1-\psi_{i})\right]$
$\displaystyle=\prod_{i=1}^{M}\left[\psi_{i}p_{i}^{y_{i}}(1-p_{i})^{K_{i}-y_{i}}+I_{[y_{i}=0]}(1-\psi_{i})\right].$
(2)

The likelihood in Equation (2) is now factorized into univariate models for
the marginal distribution of the observations. This likelihood is a mixture of
a binomial density and point mass at $y_{i}=0$, i.e. there is a probability
$\psi_{i}$ of drawing from a Binomial $(K_{i},p_{i})$ density, and a
probability (1-$\psi_{i}$) of drawing a zero (i.e. a species absence). If we
link either the detection $p_{i}$ or the occupancy $\psi_{i}$ probability to a
linear predictor $\eta_{i}$, using some known link function $g(\cdot)$, while
letting the other probability be part of the hyperparameter vector, then the
likelihood in Equation (2) has a form that is amenable to the INLA framework.
If, in addition, we assume a Gaussian distribution for the linear predictor
$\eta_{i}$ then the model can be seen as an LGM. Note that the link function
$g(\cdot)$ can be chosen to be a logit, probit or cloglog link.

Placing the occupancy model into the LGM framework gives access to the entire
existing INLA machinery that is already implemented in the R-INLA software.
The model for the parameter that is linked to the linear predictor $\eta_{i}$
can include spatial effects (both continuous and areal) as well as smooth
effects of covariates and other available random effects (see
names(inla.models()$latent) for all available random effects models). Usually,
in the INLA framework, hyperparameters are constant, i.e. they do not depend
on covariates. That would imply a constant detection probability over the
whole study domain (both in the time and space component). As this seem a too
restrictive assumption, when implementing the likelihood in Equation (2) into
R-INLA we have added some extra flexibility allowing for a transformation of
the probability $\theta_{i}=h(p_{i})$ (here $h(\cdot)$ can be the logit,
probit or cloglog link) to depend linearly on site covariates. No random
effect can be included in in this part of the model. In the context of
occupancy models, the main interest typically focuses on the occupancy
probability as it associated with the true latent occupancy state. Hence, a
most natural choice is to link the linear predictor to the probability of
occurrence. The reverse modelling choice is still available and some
practitioners might find this option relevant. Thus, while maintaining a
relatively simple structure within one model sub-component, the other sub-
component can flexibly accommodate a wide range of complex terms.

In the next section, we use two case studies to illustrate some of the spatio-
temporal models that can be fitted with this approach. First, we consider two
spatial models with independent temporal and spatial (structured and
unstructured) random effects. Such models can also be fitted with the
spOccupancy library Doser et al. (2022). However, we showcase how to include
smooth terms in the linear predictor of an occupancy model in addition to the
spatial and temporal random effects. A feature that, to our knowledge, is only
available in INLA. Then, we provide a spatio-temporal extension of both types
of models where a separable-space time model is fitted, i.e. a dynamic spatio-
temporal model, where the structured spatial component evolves through time
with autoregressive dynamics (Cameletti et al., 2013). The second case study
illustrates how spatially-varying coefficients in occupancy models can be
estimated using INLA. Spatially-varying coefficient models (SVCMs) are a
flexible class of models used when the relationship between the response and
some observed covariates is not uniform across space (Gelfand et al., 2003;
Finley, 2011). These models have been by used by Meehan et al. (2019) to
quantify the spatial variation in birds populations trends across discretized
regions. In this work, we analyse data from the North American Breeding Bird
Survey (BBS) to demonstrate how spatially-varying coefficients can be
estimated in an occupancy model for a set established birds conservation
regions. Furthermore, we formulate this model in a continuous-space framework
and use the proposed INLA-based occupancy model to identify significant
changes in the continuously and spatially varying trends.

While the structures of all of the aforementioned models vary widely in
complexity, the syntax for fitting such models in INLA varies only slightly.
Users, who are already familiar with INLA can easily fit these models or
tailor our code to their specific requirement, while new users can benefit
from the accessible R code and simulated examples available in the
supplementary online material. This online resource offers a comprehensive
explanation of the models under consideration and illustrates how various
datasets can be simulated for their subsequent analysis.

##  3 Case studies

###  3.1 Modelling the spatio-temporal distribution of the Red-eyed Vireo in
the Hubbard Brook Experimental Forest, New Hampshire, USA.

The Hubbard Brook Experimental Forest (HBEF), located in central New Hampshire
within the White Mountain National Forest, is a 3600 ha hydrological research
area that has been used extensively to monitor long-term changes in avian
populations (Holmes, 2011). In this study, we modelled the spatio-temporal
distribution of Vireo olivaceus (a.k.a Red-eyed Vireo) a migrant species and
established breeder in the HBEF. We used V. olivaceus occupancy data obtained
from the Valleywide Bird Survey (VBS) (Rodenhouse and Sillett, 2019). The VBS
is a monitoring scheme, that has been running in HBEF since 1999, where
transect surveys have been conducted each year to collect land bird point-
count data during the breeding seasons (May-July). A total of 373 locations
are monitored, for each location $i=1,\dots,373$ and time $t=2010,\dots,2018$,
the number $K_{it}$ of visits varies between 1 and 3 (see Van Tatenhove et al.
(2019) for more details).

We model _V. olivaceus_ occurrences as:

$\displaystyle z_{it}$ $\displaystyle\sim\mathrm{Bernoulli}(\psi_{it})$
$\displaystyle\mathrm{logit}(\psi_{it})=\eta_{it}$ $\displaystyle
y_{it}|z_{it}$ $\displaystyle\sim\mathrm{Binomial}(K_{it},z_{it}\times
p_{it})$ $\displaystyle\mathrm{logit}(p_{it})$
$\displaystyle=\alpha_{0}+\alpha_{1}\mbox{ mean date of
survey}_{it}+\alpha_{2}\mbox{mean date of survey}^{2}_{it},$ (3)

where $z_{it}$ indicates the occupancy status at site $i$ at time $t$, while
$y_{it}$ indicates how many of the $K_{it}$ visits resulted in a detection.
The occupancy state $z_{it}$ for the $i$-th site in year $t$ is determined by
the occupancy probability $\psi_{it}$, defined on the logit scale through the
linear predictor $\eta_{it}$. We fitted three separate space-time models for
$\eta_{it}$ all described in Table 1.

Table 1: Spatio-temporal structures and covariate effects of the linear predictor fitted in a space-time occupancy model for the Red-eyed Vireo VBS data. Model | Linear predictor $\eta_{it}$
---|---
Model 1: IID spatial effect | $\beta_{0}+f_{1}(\mbox{elevation}_{i})+f_{2}(t)+f_{3}(i)$
Model 2: Continuous spatial model | $\beta_{0}+f_{1}(\mbox{elevation}_{i})+f_{2}(t)+\omega(i)$
Model 3: Continuous space-time model | $\beta_{0}+f_{1}(\mbox{elevation}_{i})+\omega(i,t)$

All three models have an intercept $\beta_{0}$ (i.e. the logit-scale baseline
occupancy probability), and a smooth effect, $f_{1}(\mbox{elevation})$, of the
standardized elevation at each site. They differ in how the space-time
component is treated: Model (1) includes an autoregressive AR1 effect of time
$f_{2}(t)$ and an unstructured spatial effect, $f_{3}(i)$, that simply
consists of independently and identically distributed (IID) random effects of
the locations. Model (2) has the same temporal effect ($f_{2}(t)$) but
explicitly accounts for spatial correlation through the incorporation of a
Matérn Gaussian spatial field, $\omega(i)$. Finally, Model (3) includes a full
space-time model $\omega(i,t)$, where changes in time happen according to a
first order autoregressive process with spatially correlated innovations (see
for example Cameletti et al. (2013); Krainski et al. (2005) for details about
this model). In all three cases, the detection probability $p_{it}$ is assumed
to vary according to the date of each survey averaged over each year. All
three models can easily be implemented in the R-INLA library without the need
for ad hoc functions. Further modifications of the model, for example the use
of a linear instead of a smooth effect of altitude, can easily be implemented.

Running times where approximately 2 minutes for the more complex spatio-
temporal Model 3, and less than 10 seconds for the simple iid spatial model
(Table 2). Figure 1 shows the estimated non-linear relationship between
(logit) occupancy probability and elevation together with 95% credibility
bands. Table 2 reports the posterior means and 95% credibility intervals for
the intercept. All three parameters give similar estimates for this parameter.
The parameters $\alpha_{0},\alpha_{1},\alpha_{2}$ of the detection probability
were treated as hyperparameters in the implementation. Posterior means and 95%
posterior credible intervals are reported in Table 2. Both covariates
considered in this part of the model do not appear to be statistically
significant.

Parameter | IID Spatial effect | Continuous spatial model | Continuous space-time model
---|---|---|---
Occupancy
$\quad\beta_{0}$ | 0.31 ( -0.53, 1.14) | 0.17 (-1.08, 1.36) | 0.46 (-0.21,1.13)
Detection
$\quad\alpha_{0}$ |  0.19 (0.13, 0.25) |  0.18 (0.11, 0.24) |  0.15 (0.09, 0.22)
$\quad\alpha_{1}$ | 0.26 (-0.68,1.07) | 0.26 (-1.25, 1.49) | -0.09 (-0.48, 0.28)
$\quad\alpha_{2}$ | -0.53 (-1.34, 0.41) | -0.52 ( -1.74, 0.99) | -0.16 (-0.53, 0.24)
WAIC | 6694.11 | 6679.09 | 6865.82
DIC | 6714.50 | 6688.25 | 6865.85
Model likelihood | -3590.46 | -3560.36 | -3615.87
LGOCV log-score | -1.12 | -1.07 | -1.09
Run time (sec) | 5.91 | 77.80 | 174.92
Table 2: Posterior mean with 95% Credible Intervals (CrI) for the fixed
effects (bold numbers are point-wise estimates for whose CrI do not contain
zero) along with the WAIC, DIC, likelihood, log-score and run time for (i) a
simple spatial occupancy model with iid site-random effects (Model 1); (ii) a
spatially explicit occupancy model with continuous spatial processes (Model
2); (iii) spatio-temporal occupancy model with spatially continuous process
and AR1 component (Model 3). Figure 1: Smooth elevation effect on the logit-
scaled occupancy probabilities estimated from the space-time occupancy models
for the Red-eyed Vireo VBS data.

Predictions are also easily available in R-INLA as independent samples from
the fitted posterior model can be easily generated (Cristian Chiuchiolo and
Rue, 2023). Figure 2 shows the predicted surface of occupancy probabilities
for 2010 and 2018, for each of the three models in Table 1. The posterior
samples can be used to compute different metrics; for example, Figure 3
illustrates the uncertainty associated with model predictions by showing the
differences between the 0.975 and 0.025 quantiles of the predicted occupancy
probabilities for each of the three models. These maps are useful to visualize
how prediction uncertainty varies in space: for example it is easy to see that
Model (1) is, in general, associated with higher prediction uncertainty than
the other two models.

Figure 2: Predicted occupancy probability mean for Red-eyed Vireo across the
HBEF in 2010 and 2018 estimated from three different space-time occupancy
models. Figure 3: Difference between the 0.975 and 0.025 quantiles of the
estimated occupancy probabilities of the Red-eyed Vireo across the HBEF in
2010 and 2018 estimated from three different space-time occupancy models.

As mentioned, one convenient side effect of the use of INLA is that several
model comparison criteria (e.g., WAIC and DIC as shown in Table 2) are readily
available from the model output. In addition, the recently developed leave-
group-out cross-validation (LGOCV) framework (Liu and Rue, 2022), provides
cross-validation scores without the need to rerun the model. We computed the
log-score function based on the LGOCV and defined as:

$\mathrm{LOGS}_{LGOCV}=\frac{1}{n}\sum_{i=1}^{n}\pi(Y_{i}=y_{i}|\mathbf{y}_{-\mathbb{I}_{i}}),$
(4)

where the posterior predictive density
$\pi(Y_{i}=y_{i}|\mathbf{y}_{-\mathbb{I}_{i}})$ is computed for every
observation (testing point) $i$ and an appropriate group $\mathbb{I}_{i}$ that
defines the training set $\mathcal{D}_{-\mathbb{I}_{i}}$. In this work, we
used INLA’s correlation-based automatic construction of the leave out group
$\mathbb{I}_{i}$ (we use the inla.group.cv function and set num.level.sets = 3
which controls the degree of independence ), but see the supplementary online
material for an example, where the leave out groups are manually constructed.
All these measures, for the three considered models, are reported in Table 2.
All of them disfavour Model 1 while the ranking between Model 2 and 3 is not
unique and depends on the chosen criteria.

###  3.2 Occupancy model with spatially-varying coefficients for the Gray
Catbird

The North American Breeding Bird Survey (BBS) is a comprehensive monitoring
program devised to track the status and trends of bird populations across
North America (Ziolkowski et al., 2023). In this work, we use BBS data on
distribution of the species _Dumetella carolinensis_ (a.k.a. Gray Catbird) to
illustrate how INLA can be used to quantify spatially-explicit trends while
correcting for imperfect detection. The data comprise of presence/absence
records of the Gray Catbird at $i=1,\dots,1846$ sites sampled at least once
between $t=2000,\dots,2019$. The model structure is similar to Equation (3),
and again, we showcase three alternative models for the linear predictor
$\eta_{it}$, linked to the probability of occurrence $\psi_{it}$.

Table 3: Space-varying time trends in occupancy models for the Gray Catbird BBS Data. Model | linear predictor
---|---
Model 1: Fixed time trend | $\eta_{it}=\beta_{0}+\omega(i)+\epsilon_{t}+\beta_{t}t$
Model 2: Continuously varying time trend | $\eta_{it}=\beta_{0}+\omega(i)+\epsilon_{t}+\beta_{t}(i)t$
Model 3: Regional varying time trend | $\eta_{it}=\beta_{0}++\omega(i)+\epsilon_{t}+\beta_{t}(\mbox{bcr}_{i})t$

The three models are summarised in Table 3. Again, all models contain an
intercept $\beta_{0}$, a spatially continuous Gaussian random field with
Matérn covariance structure, $\omega(i)$, and an iid effect of the year
$\epsilon_{t}$. They differ in the way a linear time trend is treated. In
Model 1, the time trend $\beta\ t$ is constant over the whole space, while in
the other two models the time trend varies in space. In Model 2 the time trend
varies continuously with time while in Model 3 the time trend varies
discretely, by bird conservation region (brc) as defined by the North American
Bird Conservation Initiative (Bird Studies Canada and NABCI, 2014)). We use a
Matérn continuous Gaussian model for $\beta_{t}(i)$ in Model 2 and the Besag-
York-Mollié model (BYM) (Besag, 1991) (an extension to the intrinsic CAR model
that contains an i.i.d. model component for the non-spatial heterogeneity) for
$\beta_{t}(\text{bcr}_{i})$ in Model 3.

Finally, we define the logit of the detection probabilities as a linear and
quadratic function of the Julian date of the visit, i.e.
$\mathrm{logit}(p_{i})=\alpha_{0}+\alpha_{1}\mbox{day}_{i}+\alpha_{2}\mbox{day}^{2}_{i}$.

Table 4 summarizes the posterior means and 95% credibility intervals for the
baseline occupancy probability ($\beta_{0}$) and the parameters of the
detection model ($\alpha_{0},\alpha_{1},\alpha_{2}$) for all three models.
These estimates are rather similar for the three models. The same Table also
reports the DIC, marginal likelihood and log-LGOCV score: all three scores
marginally prefer the continuously varying coefficients.

Parameter | Constant spatial trend | Continuously spatial-varying trend | Regionally spatial-varying trend
---|---|---|---
Occupancy
$\quad\beta_{0}$ | 0.064 (-1.38, 1.51) | 0.30 (-0.84, 1.44) | 0.11 (-1.28,1.50)
$\quad\beta_{t}$ |  0.17 ( -0.13,-0.21) | — | —
Detection
$\quad\alpha_{0}$ |  0.46 (0.45,0.48) |  0.46 (0.45,0.48) |  0.45 (0.46,0.47)
$\quad\alpha_{1}$ |  3.24 (2.95,23.5) |  3.73 (3.44,3.99) |  3.26 (2.95,3.58)
$\quad\alpha_{2}$ |  -3.12 (-3.37,-3.12) |  -3.60 (-3.84,-3.31) |  -3.13 (-3.45,-2.83)
WAIC | 84062.2 | 83848.44 | 84018.57
DIC | 84160.72 | 83978.49 | 84121.27
Marginal likelihood | -42464.14 | -42394.32 | -42451.11
LGOCV log-score | -1.591 | -1.589 | -1.591
Run time (sec) | 28.53 | 134.11 | 43.77
Table 4: Posterior mean with 95% Credible Intervals for the fixed effects
(bold numbers are point-wise estimates for whose CrI do not contain zero)
along with the WAIC, DIC, likelihood and run times of the three SVC occupancy
models fitted to the Gray Catbird data.

The posterior mean of the estimated temporal trend in space is shown in Figure
4(upper panel), while areas in space where the credibility interval for the
estimated trend does not include 0 are shown in Figure 4 (lower panel). Model
1 estimates a constant and significant negative trend across the whole domain
(the estimated constant linear trend for Model 1 is also reported in Table 4).
When we allow the trend to vary in space we get a more nuanced picture: Model
2 identifies ares in the central part of the domain with significantly
negative trend but also small areas with significantly positive trend, while
Model 3 identifies only areas with a significantly negative trend. Model 2 and
Model 3 agree on identifying such areas. Visualization of the SVC models
predictions over space can be used to inform the monitoring of long-term
changes in occupancy trends, which is a major goal in many ecological studies
as this allow hot spots relevant for conservation and managements to be
identified (Altwegg and Nichols, 2019; Outhwaite et al., 2020; Coomber et al.,
2021).

Figure 4: Estimated spatiotemporal trends for the Gray Catbird from 2000-2019.
The first row shows the estimated mean trend for SVCM 1 (fixed constant
trend), SVCM 2 (continuously spatial-varying trend), and SVCM 3 (regional
spatial-varying trend with BYM structured spatial random effects ). The second
row shows significant changes in the trend, i.e. a negative trend if the 95 %
credible interval ¡ zero, a positive trend if the credible interval was ¿
zero, and no significant change otherwise.

Maps of posterior means of the predicted occupancy probabilities for Gray
Catbird in 2020 for all three models are reproduced in Figure 5. The maps
appear to be quite similar, however, spatial differences in logit-scale
occupancy probabilities can be identified by comparing the predictions from
the fixed trend model (Model 1) against the continuously and regionally
spatial-varying trend models (Model 2 and Model 3, respectively; Figure 6 ).

Figure 5: Gray Catbird estimated occupancy probabilities from three SVCMs in
2020. Figure 6: Differences between the estimated logit-scaled mean occupancy
probabilities across the SVC occupancy models for the Gray Catbird. On the
left, the mean differences between the fixed trend model and the continuously
spatially varying trend model. On the right, the mean differences between the
model with a fixed trend and a regionally varying trend.

##  4 Discussion

In recent years, the rise of new technologies and data submission platforms
has provided access to many and often multiple data sources (e.g. derived from
citizen science programs, atlases, museums, and planned surveys) for species
distributions modelling (August et al., 2015). As the volume of recorded
ecological data continues to grow, so does the need for computationally
efficient methods for fitting complex spatio-temporal ecological models.
Recent developments in spatio-temporal occupancy models have focused on using
Gaussian processes (GP) or nearest‐neighbor Gaussian processes (NNGP)
approximations while taking the computational advantage of using the Pólya-
Gamma data augmentation scheme (Polson et al., 2013) to improve computational
running times. The Pólya-Gamma data augmentation approach has proven to be a
useful method for constructing computationally efficient Gibbs samplers for
estimating occupancy models parameters (Clark and Altwegg, 2019), and has
recently been used for introducing spatio-temporal random effects in occupancy
models fitted to large data sets (Diana et al., 2022; Doser et al., 2022).
However, users interested in employing this approach must either code their
own sampling algorithm or use the predefined set of models provided by a
specific package. Additionally, the Pólya-Gamma scheme relies on the logit-
link function for defining the linear predictor, hence imposing limitations on
certain applications (e.g., when aiming to employ the cloglog function to link
the linear predictor to the log of the intensity of the underlying Poisson
process (Simmonds et al., 2020)).

INLA is an alternative to MCMC that has already proven to open avenues for
groundbreaking progress in many ecological studies (Laxton et al., 2023;
Paradinas et al., 2023; Williamson et al., 2022, among others). In this study,
we show that, under some conditions, occupancy models can also be fitted
within the INLA framework and therefore benefit from the wide range of models
already implemented in the R-INLA software. The main limitation of the INLA
approach is that one cannot use survey-specific covariates in the detection
model and that random effects can only be included either in the detection or
in the occupancy sub-component. In cases where random effects are required in
both model components, or if survey specific covariates need to be accounted
for, an MCMC approach might be more appropriate. Nonetheless, Meehan et al.
(2020) showed that averaging survey level covariates across sites, yields very
similar results for N-Mixture models when comparing INLA with classical or
MCMC-based approaches. Our approach is more flexible than Meehan et al. (2020)
as it allows: (1) detection probabilities to vary across sites and (2) random
effects to be included in the linear predictor associated with the ecological
process of interest.

Further, it is also important to notice that the log-likelihood for the
occupancy models is not log-concave, therefore, special care has to be taken
when formulating the models and choosing the priors. We believe that, for many
applied situations, the flexibility and the ease of implementation that the
INLA framework offers to the user, will compensate for the slightly more rigid
model that one has to adopt for the detection probability.

In recent years, INLA has become a major computational inferential framework
used by practitioners and statisticians alike (Van Niekerk et al., 2023; Van
Niekerk et al., 2021). Its popularity is owed not only to the efficient
computational machinery, but also to its ongoing development and derivative
projects that enhance and expand upon INLA’s capabilities (Van Niekerk et al.,
2021). Such developments can greatly enhance current occupancy models theory
and applications. For example, prior choice is an important step in Bayesian
modelling, penalized complexity (PC) priors (Simpson et al., 2017), as
implemented in INLA has proven to be a useful tool in ecological modelling as
a selection tool (Sadykova et al., 2017) or as a way of defining the
parameters of the covariance structure of the Matérn field (Laxton et al.,
2023). Although PC priors can also be used in contexts outside the INLA
framework, the fact that they are implemented in the R-INLA packages, makes
them especially convenient for INLA users. Furthermore, INLA provides a range
of useful tools for model checking and model selection with model comparison
criteria as DIC and WAIC. It also provides the recently introduced LGOCV
strategy, which enables both manual and automatic construction of validation
groups while avoiding the need for refitting the model multiple times for each
leave-out group to assess the predictive performance of different models in a
computationally efficient manner. In this work, we have presented how INLA can
be used to fit a single species occupancy model. However, since the model has
been made amenable to the INLA framework, it is easy to build more complex
models for multiple or for combining different data sources, where different
species or data sets can be represented by different likelihoods that share
some of the components in the linear predictor.

Occupancy models are better suited for protocolized surveys where the sampling
design can effectively correct for imperfect detection while minimizing other
sources of sampling biases (although occupancy models have also proven to be
useful in certain opportunistic surveys (van Strien et al., 2013)). However,
unstructured ecological data without a standardized collection scheme have
become more common over recent years, offering a much greater spatio-temporal
and taxonomic coverage than protocolized data sets. Unfortunately, such data
sources often consist solely of presence-only records, lacking information on
species absence, and are often spatially and temporally biased towards the
choice of the sampling locations (van Strien et al., 2013). Thus, integrated
distribution models (IDMs) have emerged as a way to combine the information
from different data sources, such as those arising from structured and
unstructured surveys, to enhance the accuracy of the estimates of species
distribution models (Koshkina et al., 2017; Simmonds et al., 2020). INLA has
proven to be powerful platform to fit complex integrated models to multiple
data sources with different biases and observational processes (Martino et
al., 2021; Villejo et al., 2023). Thus, future work may focus on the
integration of the occupancy model framework to develop accurate and efficient
IDMs. This can also help in cases where site-survey covariates are relevant
for the analysis, but not controlled due to the lack of a survey design. For
example, data from (i) planned surveys where site-level covariates can account
for observational process biases can be analysed with INLA’s occupancy model
formulation, and (ii) data from unstructured surveys could be analyzed as a
thinned point process as in Yuan et al. (2017); Koshkina et al. (2017), while
jointly estimating the parameters related to the true species distributions.

The challenge of scalability for large and complex spatiotemporal datasets is
a well-known limitation of occupancy models. Thus, INLA offers an alternative
inferential and computationally efficient framework for handling intricate
spatiotemporal occupancy models within an established modelling framework. The
appeal INLA has as a modelling framework for occupancy models, lies in its
potential to address complex ecological problems by using the machinery
already available within the R-INLA software.

## conflict of interest

We warrant that this manuscript is the original work of the authors listed and
has not previously been published. All the authors listed made substantial
contributions to the manuscript and qualify for authorship, and no authors
have been omitted. We warrant that none of the authors has any conflict of
interest in regard to this manuscript

## Supporting Information

All data sources used for this work are publicly available and have been
properly cited in the manuscript. Red-eyed Vireo occurrences data are
available from the Valleywide Bird Survey knb-lter-hbr.178. Gray Catbrid data
from the North American BBS can be obtained from
https://doi.org/10.5066/P9GS9K64. R code and simulated examples are available
in the online supplementary material.

## References

  * Adin et al. [2023] A Adin, E Krainski, A Lenzi, Z Liu, J Martínez-Minaya, and H Rue.  Automatic cross-validation in structured models: Is it time to leave out leave-one-out?  November 2023.
  * Altwegg and Nichols [2019] Res Altwegg and James D Nichols.  Occupancy models for citizen-science data.  _Methods in Ecology and Evolution_ , 10(1):8–21, 2019.
  * Arce Guillen et al. [2023] Rafael Arce Guillen, Finn Lindgren, Stefanie Muff, Thomas W. Glass, Greg A. Breed, and Ulrike E. Schlägel.  Accounting for unobserved spatial variation in step selection analyses of animal movement via spatial random effects.  _Methods in Ecology and Evolution_ , 14(10):2639–2653, 2023.  doi: https://doi.org/10.1111/2041-210X.14208.  URL https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.14208.
  * August et al. [2015] Tom August, Martin Harvey, Paula Lightfoot, David Kilbey, Timos Papadopoulos, and Paul Jepson.  Emerging technologies for biological recording.  _Biological Journal of the Linnean Society_ , 115(3):731–749, 2015.
  * Bakka et al. [2018] Haakon Bakka, Håvard Rue, Geir-Arne Fuglstad, Andrea Riebler, David Bolin, Janine Illian, Elias Krainski, Daniel Simpson, and Finn Lindgren.  Spatial modeling with r-inla: A review.  _WIREs Computational Statistics_ , 10(6):e1443, 2018.  doi: https://doi.org/10.1002/wics.1443.  URL https://wires.onlinelibrary.wiley.com/doi/abs/10.1002/wics.1443.
  * Besag [1991] J Besag.  Bayesian image restoration with two applications in spatial statistics (with discussion).  _Ann I Stat Math 43_ , 1:59, 1991.
  * Bird Studies Canada and NABCI [2014] Bird Studies Canada and NABCI.  Bird Conservation Regions.  Published by Bird Studies Canada on behalf of the North American Bird Conservation Initiative, 2014.  https://birdscanada.org/bird-science/nabci-bird-conservation-regions Accessed: 16.08.2023.
  * Blangiardo and Cameletti [2015] Marta Blangiardo and Michela Cameletti.  _Spatial and spatio-temporal Bayesian models with R-INLA_.  John Wiley & Sons, 2015.
  * Cameletti et al. [2013] Michela Cameletti, Finn Lindgren, Daniel Simpson, and Håvard Rue.  Spatio-temporal modeling of particulate matter concentration through the spde approach.  _AStA Advances in Statistical Analysis_ , 97:109–131, 2013\.
  * Carpenter et al. [2017] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus A Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: A probabilistic programming language.  _Journal of statistical software_ , 76, 2017.
  * Clark and Altwegg [2023] A. Clark and R. Altwegg.  Gibbs sampler for multi-species occupancy models.  _Environmental and Ecological Statistics_ , 2023.
  * Clark and Altwegg [2019] Allan E. Clark and Res Altwegg.  Efficient bayesian analysis of occupancy models with logit link functions.  _Ecology and Evolution_ , 9(2):756–768, 2019\.  doi: https://doi.org/10.1002/ece3.4850.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/ece3.4850.
  * Coomber et al. [2021] Frazer G Coomber, Bethany R Smith, Tom A August, Colin A Harrower, Gary D Powney, and Fiona Mathews.  Using biological records to infer long-term occupancy trends of mammals in the uk.  _Biological Conservation_ , 264:109362, 2021.
  * Cressie et al. [2009] Noel Cressie, Catherine A Calder, James S Clark, Jay M Ver Hoef, and Christopher K Wikle.  Accounting for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical modeling.  _Ecological Applications_ , 19(3):553–570, 2009\.
  * Cristian Chiuchiolo and Rue [2023] Janet van Niekerk Cristian Chiuchiolo and Håvard Rue.  Joint posterior inference for latent gaussian models with r-inla.  _Journal of Statistical Computation and Simulation_ , 93(5):723–752, 2023.  doi: 10.1080/00949655.2022.2117813.
  * de Valpine et al. [2017] Perry de Valpine, Daniel Turek, Christopher J Paciorek, Clifford Anderson-Bergman, Duncan Temple Lang, and Rastislav Bodik.  Programming with models: writing statistical algorithms for general model structures with nimble.  _Journal of Computational and Graphical Statistics_ , 26(2):403–413, 2017.
  * Devarajan et al. [2020] Kadambari Devarajan, Toni Lyn Morelli, and Simone Tenan.  Multi-species occupancy models: review, roadmap, and recommendations.  _Ecography_ , 43:1612–1624, 2020.
  * Diana et al. [2022] Alex Diana, Emily Beth Dennis, Eleni Matechou, and Byron John Treharne Morgan.  Fast bayesian inference for large occupancy datasets.  _Biometrics_ , 2022.
  * Dorazio and Royle [2005] Robert M Dorazio and J Andrew Royle.  Estimating size and composition of biological communities by modeling the occurrence of species.  _Journal of the American Statistical Association_ , 100(470):389–398, 2005.
  * Dorazio et al. [2006] Robert M Dorazio, J Andrew Royle, Bo Söderström, and Anders Glimskär.  Estimating species richness and accumulation by modeling species occurrence and detectability.  _Ecology_ , 87(4):842–854, 2006.
  * Doser et al. [2022] Jeffrey W Doser, Andrew O Finley, Marc Kéry, and Elise F Zipkin.  spoccupancy: An r package for single-species, multi-species, and integrated spatial occupancy models.  _Methods in Ecology and Evolution_ , 13(8):1670–1678, 2022.
  * Elith and Leathwick [2009] Jane Elith and John R Leathwick.  Species distribution models: ecological explanation and prediction across space and time.  _Annual Review of Ecology, Evolution, and Systematics_ , 40:677–697, 2009.
  * Finley [2011] Andrew O. Finley.  Comparing spatially-varying coefficients models for analysis of ecological data with non-stationary and anisotropic residual dependence.  _Methods in Ecology and Evolution_ , 2(2):143–154, 2011.  doi: https://doi.org/10.1111/j.2041-210X.2010.00060.x.  URL https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/j.2041-210X.2010.00060.x.
  * Fiske and Chandler [2011] Ian Fiske and Richard Chandler.  unmarked: An R package for fitting hierarchical models of wildlife occurrence and abundance.  _Journal of Statistical Software_ , 43(10):1–23, 2011.  URL http://www.jstatsoft.org/v43/i10/.
  * Gelfand et al. [2003] Alan E. Gelfand, Hyon-Jung Kim, C. F. Sirmans, and Sudipto Banerjee.  Spatial modeling with spatially varying coefficient processes.  _Journal of the American Statistical Association_ , 98(462):387–396, 2003.  ISSN 01621459.  URL http://www.jstor.org/stable/30045248.
  * Guillera-Arroita [2017] Gurutzeta Guillera-Arroita.  Modelling of species distributions, range dynamics and communities under imperfect detection: advances, challenges and opportunities.  _Ecography_ , 40(2):281–295, 2017.
  * Hines [2006] JE Hines.  Presence: Software to estimate patch occupancy and related parameters, version 11.5, 2006.
  * Holmes [2011] Richard T Holmes.  Avian population and community processes in forest ecosystems: Long-term research in the hubbard brook experimental forest.  _Forest Ecology and Management_ , 262(1):20–32, 2011.
  * Johnson [2021] Devin S. Johnson.  _stocc: Fit a Spatial Occupancy Model via Gibbs Sampling_ , 2021.  URL https://CRAN.R-project.org/package=stocc.  R package version 1.31.
  * Jona Lasinio et al. [2013] Giovanna Jona Lasinio, Gianluca Mastrantonio, and Alessio Pollice.  Discussing the “big n problem”.  _Statistical Methods & Applications_, 22:97–112, 2013\.
  * Kellner [2021] Ken Kellner.  _ubms: Bayesian Models for Data from Unmarked Animals using ’Stan’_ , 2021.  URL https://CRAN.R-project.org/package=ubms.  R package version 1.0.2.
  * Kellner and Swihart [2014] Kenneth F Kellner and Robert K Swihart.  Accounting for imperfect detection in ecology: a quantitative review.  _PLoS One_ , 9:e111436, 2014.
  * Kéry and Royle [2020] Marc Kéry and J Andrew Royle.  _Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 2: Dynamic and Advanced Models_.  Academic Press, 2020.
  * Kéry et al. [2010] Marc Kéry, J Andrew Royle, Hans Schmid, Michael Schaub, Bernard Volet, Guido Haefliger, and Niklaus Zbinden.  Site-occupancy distribution modeling to correct population-trend estimates derived from opportunistic observations.  _Conservation Biology_ , 24(5):1388–1397, 2010\.
  * Koshkina et al. [2017] Vira Koshkina, Yan Wang, Ascelin Gordon, Robert M Dorazio, Matt White, and Lewi Stone.  Integrated species distribution models: combining presence-background data and site-occupancy data with imperfect detection.  _Methods in Ecology and Evolution_ , 8(4):420–430, 2017.
  * Krainski et al. [2005] Elias Krainski, Virgilio Gómez-Rubio, Haakon Bakka, Amanda Virgilio, Lenzi, Daniela Castro-Camilo, Daniel Simpson, Finn Lindgren, and Håvard Rue.  _Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA_.  Chapman & Hall/CRC Press, 2005.
  * Laxton et al. [2023] Megan R Laxton, Óscar Rodríguez De Rivera, Andrea Soriano-Redondo, and Janine B Illian.  Balancing structural complexity with ecological insight in spatio-temporal species distribution models.  _Methods in Ecology and Evolution_ , 14(1):162–172, 2023.
  * Lindgren et al. [2011] Finn Lindgren, Håvard Rue, and Johan Lindström.  An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 73(4):423–498, 2011.
  * Liu and Rue [2022] Zhedong Liu and Haavard Rue.  Leave-group-out cross-validation for latent gaussian models.  _arXiv preprint arXiv:2210.04482_ , 2022.
  * MacKenzie et al. [2002] Darryl I MacKenzie, James D Nichols, Gideon B Lachman, Sam Droege, J Andrew Royle, and Catherine A Langtimm.  Estimating site occupancy rates when detection probabilities are less than one.  _Ecology_ , 83(8):2248–2255, 2002.
  * Martino et al. [2021] Sara Martino, Daniela Silvia Pace, Stefano Moro, Edoardo Casoli, Daniele Ventura, Alessandro Frachea, Margherita Silvestri, Antonella Arcangeli, Giancarlo Giacomini, Giandomenico Ardizzone, et al.  Integration of presence-only data from several sources: a case study on dolphins’ spatial distribution.  _Ecography_ , 44(10):1533–1543, 2021.
  * Meehan et al. [2020] T. Meehan, N. Michel, and H. Rue.  Estimating animal abundance with n-mixture models using the r-inla package for r.  _Journal of Statistical Software_ , 95(2):1–26, 2020.  doi: 10.18637/jss.v095.i02.  URL https://www.jstatsoft.org/index.php/jss/article/view/v095i02.
  * Meehan et al. [2019] Timothy D. Meehan, Nicole L. Michel, and Håvard Rue.  Spatial modeling of audubon christmas bird counts reveals fine-scale patterns and drivers of relative abundance trends.  _Ecosphere_ , 10(4):e02707, 2019.  doi: https://doi.org/10.1002/ecs2.2707.  URL https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1002/ecs2.2707.
  * Northrup and Gerber [2018] Joseph M Northrup and Brian D Gerber.  A comment on priors for bayesian occupancy models.  _PloS one_ , 13:e0192819, 2018.
  * Outhwaite et al. [2020] Charlotte L Outhwaite, Richard D Gregory, Richard E Chandler, Ben Collen, and Nick JB Isaac.  Complex long-term biodiversity change among invertebrates, bryophytes and lichens.  _Nature ecology & evolution_, 4(3):384–392, 2020.
  * Paradinas et al. [2015] Iosu Paradinas, David Conesa, M Grazia Pennino, Facundo Muñoz, Angel M Fernández, Antonio López-Quílez, and José María Bellido.  Bayesian spatio-temporal approach to identifying fish nurseries by validating persistence areas.  _Marine Ecology Progress Series_ , 528:245–255, 2015.
  * Paradinas et al. [2023] Iosu Paradinas, Janine B Illian, Alexandre Alonso-Fernändez, Maria Grazia Pennino, and Sophie Smout.  Combining fishery data through integrated species distribution models.  _ICES Journal of Marine Science_ , 80(10):2579–2590, 2023.
  * Polson et al. [2013] Nicholas G Polson, James G Scott, and Jesse Windle.  Bayesian inference for logistic models using pólya–gamma latent variables.  _Journal of the American statistical Association_ , 108(504):1339–1349, 2013.
  * Ponisio et al. [2020] Lauren C Ponisio, Perry de Valpine, Nicholas Michaud, and Daniel Turek.  One size does not fit all: Customizing mcmc methods for hierarchical models using nimble.  _Ecology and evolution_ , 10(5):2385–2416, 2020\.
  * Rodenhouse and Sillett [2019] Nick Rodenhouse and Scott Sillett.  Valleywide bird survey, hubbard brook experimental forest, 1999-2016 (ongoing).  2019\.
  * Rue and Held [2005] Havard Rue and Leonhard Held.  _Gaussian Markov Random Fields: theory and applications_.  CRC press, 2005.
  * Rue et al. [2009] Håvard Rue, Sara Martino, and Nicolas Chopin.  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 71(2):319–392, 2009.
  * Rue et al. [2017] Håvard Rue, Andrea Riebler, Sigrunn H. Sørbye, Janine B. Illian, Daniel P. Simpson, and Finn K. Lindgren.  Bayesian computing with inla: A review.  _Annual Review of Statistics and Its Application_ , 4(1):395–421, 2017.  doi: 10.1146/annurev-statistics-060116-054045.  URL https://doi.org/10.1146/annurev-statistics-060116-054045.
  * Rushing et al. [2019] Clark S Rushing, J Andrew Royle, David J Ziolkowski, and Keith L Pardieck.  Modeling spatially and temporally complex range dynamics when detection is imperfect.  _Scientific Reports_ , 9(1):1–9, 2019.
  * Sadykova et al. [2017] Dinara Sadykova, Beth E Scott, Michela De Dominicis, Sarah L Wakelin, Alexander Sadykov, and Judith Wolf.  Bayesian joint models with inla exploring marine mobile predator–prey and competitor species habitat overlap.  _Ecology and Evolution_ , 7(14):5212–5226, 2017\.
  * Simmonds et al. [2020] Emily G Simmonds, Susan G Jarvis, Peter A Henrys, Nick JB Isaac, and Robert B O’Hara.  Is more data always better? a simulation study of benefits and limitations of integrated distribution models.  _Ecography_ , 43(10):1413–1422, 2020.
  * Simpson et al. [2017] Daniel Simpson, Håvard Rue, Andrea Riebler, Thiago G. Martins, and Sigrunn H. Sørbye.  Penalising model component complexity: A principled, practical approach to constructing priors.  _Statistical Science_ , 32(1):1–28, 2017.  ISSN 08834237, 21688745.  URL http://www.jstor.org/stable/26408114.
  * Spiegelhalter et al. [2002] David J. Spiegelhalter, Nicola G. Best, Bradley P. Carlin, and Angelika Van Der Linde.  Bayesian measures of model complexity and fit.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 64(4):583–639, 2002.  doi: https://doi.org/10.1111/1467-9868.00353.  URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/1467-9868.00353.
  * Tyre et al. [2003] Andrew J Tyre, Brigitte Tenhumberg, Scott A Field, Darren Niejalke, Kirsten Parris, and Hugh P Possingham.  Improving precision and reducing bias in biological surveys: estimating false-negative error rates.  _Ecological Applications_ , 13(6):1790–1801, 2003\.
  * Van Niekerk et al. [2021] Janet Van Niekerk, Haakon Bakka, Håvard Rue, and Olaf Schenk.  New frontiers in bayesian modeling using the inla package in r.  _Journal of Statistical Software_ , 100(2):1–28, 2021.  doi: 10.18637/jss.v100.i02.  URL https://www.jstatsoft.org/index.php/jss/article/view/v100i02.
  * Van Niekerk et al. [2023] Janet Van Niekerk, Elias Krainski, Denis Rustand, and Håvard Rue.  A new avenue for bayesian inference with inla.  _Computational Statistics & Data Analysis_, 181:107692, 2023.  ISSN 0167-9473.  doi: https://doi.org/10.1016/j.csda.2023.107692.  URL https://www.sciencedirect.com/science/article/pii/S0167947323000038.
  * van Strien et al. [2013] Arco J van Strien, Chris AM van Swaay, and Tim Termaat.  Opportunistic citizen science data of animal species produce reliable estimates of distribution trends if analysed with occupancy models.  _Journal of Applied Ecology_ , 50(6):1450–1458, 2013.
  * Van Tatenhove et al. [2019] Aimee Van Tatenhove, Emily Filiberti, T Scott Sillett, Nicholas Rodenhouse, and Michael Hallworth.  Climate-related distribution shifts of migratory songbirds and sciurids in the white mountain national forest.  _Forests_ , 10(2):84, 2019.
  * Vieilledent [2019] G. Vieilledent.  _hSDM: Hierarchical Bayesian species distribution models._ , 2019\.  URL https://CRAN.R-project.org/package=hSDM.  R package version 1.4.4.
  * Villejo et al. [2023] Stephen Jun Villejo, Janine B Illian, and Ben Swallow.  Data fusion in a two-stage spatio-temporal model using the inla-spde approach.  _Spatial Statistics_ , 54:100744, 2023.  ISSN 2211-6753.  doi: https://doi.org/10.1016/j.spasta.2023.100744.  URL https://www.sciencedirect.com/science/article/pii/S2211675323000192.
  * Watanabe [2013] Sumio Watanabe.  A widely applicable bayesian information criterion.  _J. Mach. Learn. Res._ , 14(1):867–897, mar 2013\.  ISSN 1532-4435.
  * Williamson et al. [2022] Laura D. Williamson, Beth E. Scott, Megan Laxton, Janine B. Illian, Victoria L.G. Todd, Peter I. Miller, and Kate L. Brookes.  Comparing distribution of harbour porpoise using generalized additive models and hierarchical bayesian models with integrated nested laplace approximation.  _Ecological Modelling_ , 470:110011, 2022.  ISSN 0304-3800.  doi: https://doi.org/10.1016/j.ecolmodel.2022.110011.  URL https://www.sciencedirect.com/science/article/pii/S0304380022001223.
  * Wright et al. [2019] Wilson J Wright, Kathryn M Irvine, and Megan D Higgs.  Identifying occupancy model inadequacies: can residuals separately assess detection and presence?  _Ecology_ , 100(6):e02703, 2019.
  * Yuan et al. [2017] Yuan Yuan, Fabian E Bachl, Finn Lindgren, David L Borchers, Janine B Illian, Stephen T Buckland, Håvard Rue, and Tim Gerrodette.  Point process models for spatio-temporal distance sampling data from a large-scale survey of blue whales.  2017\.
  * Ziolkowski et al. [2023] D.J. Ziolkowski, M. Lutmerding, V.I. English, W.B.and Aponte, and M-A.R. Hudson.  North American Breeding Bird Survey Dataset 1966 - 2022: U.S. Geological Survey data release, 2023.