Preprint DOI to follow Preprint article

Koh and Opitz

[$\ast$]Corresponding author<EMAIL_ADDRESS>
0Year 0Year 0Year

# Extreme-value modelling of migratory bird arrival dates: Insights from
citizen science data

Jonathan Koh0000-0002-7287-0588  Thomas Opitz0000-0002-5863-5020  Institute of
Mathematical Statistics and Actuarial Science, Oeschger Centre for Climate
Change Research, University of Bern, Switzerland  Biostatistics and Spatial
Processes (UR546), INRAE, Avignon, France

(2023; 2019; Date; Date; Date)

###### Abstract

Citizen science mobilises many observers and gathers huge datasets but often
without strict sampling protocols, which results in observation biases due to
heterogeneity in sampling effort that can lead to biased statistical
inferences. We develop a spatiotemporal Bayesian hierarchical model for bias-
corrected estimation of arrival dates of the first migratory bird individuals
at a breeding site. Higher sampling effort could be correlated with earlier
observed dates. We implement data fusion of two citizen-science datasets with
sensibly different protocols (BBS, eBird) and map posterior distributions of
the latent process, which contains four spatial components with Gaussian
process priors: species niche; sampling effort; position and scale parameters
of annual first date of arrival. The data layer includes four response
variables: counts of observed eBird locations (Poisson); presence-absence at
observed eBird locations (Binomial); BBS occurrence counts (Poisson); first
arrival dates (Generalized Extreme-Value). We devise a Markov Chain Monte
Carlo scheme and check by simulation that the latent process components are
identifiable. We apply our model to several migratory bird species in the
northeastern US for 2001–2021. The sampling effort is shown to significantly
modulate the observed first arrival date. We exploit this relationship to
effectively debias predictions of the true first arrival dates.

###### keywords:

Bayesian hierarchical model; Bias correction; Bird phenology; Opportunistic
data; Sampling effort; Species distribution

##  1 Introduction

###  1.1 The rise of citizen science

Defining citizen science and its boundaries is difficult (Haklay et al.,
2021). Broadly speaking, the field involves a collaborative approach to
scientific inquiry that engages volunteers and non-professionals in the
collection, analysis, and interpretation of data. This participatory model is
not new. Friedrich Wilhelm Hersche, a musician by training, discovered the
planet Uranus in 1781 using his own telescopes. Charles Darwin conducted
crowd-sourcing projects in the $19^{\text{th}}$ century by recruiting
acquaintances and travelers to write to him with their own observations Browne
(1996). The longest running citizen science initiative is the annual Christmas
Bird Count (Bock and Root, 1981), a census administered by the National
Audubon Society of wintering birds in the Western Hemisphere by volunteer
birdwatchers since 1900.

During the past two decades, citizen science has gained significant traction,
spurred on in part by growing world literacy levels and the advent of devices
and infrastructures in information technology for gathering, reporting,
sharing and storing data (Newman et al., 2012; Wynn, 2017). These datasets are
often collected at relatively low cost and sometimes with separate and
alternative funding streams (Silvertown, 2009).
Zooinverse111https://www.zooniverse.org/, the world’s largest platform for
volunteer-based research, has seen a surge in registered volunteers to over
1.9 Million since its inception in 2009. Climateprediction.net runs climate
modelling experiments using the home computers of thousands of volunteers.
There is now a peer-reviewed journal (Bonney et al., 2016) dedicated to
disseminating research on citizen science, and the 2019 Citizen Science
Association bienniaol conference attracted 818 registered delegates from 28
countries.

Citizen science is expected to continue to grow in importance. Further
technological progress will improve how machines and citizen science can work
together; for example, volunteer-classified training sets have already been
used to improve the performance of machine learning approaches in astronomy
(Marshall et al., 2015). Fraisl et al. (2020) note that traditional data
sources are not sufficient for measuring the United Nations Sustainable
Development Goals (UN, 2015), but sources from citizen science are required to
better inform policies and actions. Citizen science initiatives can also have
community-level impacts on the participants (Jordan et al., 2012), such as by
empowering them and giving then a voice in local environmental decision-making
Bonney et al. (2016) while increasing public awareness of the scientific
process. Though citizen science’s global economic value has been estimated to
be over 2.5 billion USD annually, the majority of data collected through
citizen-science initiatives do not reach analysis in peer-reviewed literature
(Theobald et al., 2015). There is a need for the development and dissemination
of statistical methods that could facilitate wider scientific use of these
datasets.

###  1.2 Analysis of data for species monitoring

Much of the recent growth in citizen science has taken place in the ecological
and environmental sciences (McKinley et al., 2015; Pocock et al., 2018; Fraisl
et al., 2022). A prominent example in ecology is the eBird project (Cornell
Lab of Ornithology, 2022), launched in 2002, which has led to a database
volume of more than 200 Gigabytes providing information about over 600 million
bird observations. Another important bird monitoring program within the scope
of citizen science is the Breeding Bird Survey (BBS). It was launched in 1966
for official monitoring purposes and engages large numbers of trained
volunteers in collecting standardised data on bird populations according to a
relatively strict sampling protocol along a number of selected routes and
during each year. We will exploit these two data source in this work.
Regarding non avian species, another recent example is the deployment of
camera-trapping projects for mammal monitoring (Hsing et al., 2022).

Many observation data in ecology are _opportunistic_ : they were collected
incidentally or without a predefined research question in mind, such that
information about the criteria applied by the observer for sampling and for
reporting observations is only limited. However, opportunistic data often
allow for a spatiotemporal coverage of species monitoring that would not be
attainable with only protocol-based data collected by professionals at
relatively sparse and deterministically chosen sampling locations. Citizen-
science data is particularly valuable for statistical inference on low-
probablility events (e.g., occurrences of rare species, or unusual or extreme
phenological events) and on occurrence times of events, such as the arrival of
migratory birds at their breeding location in spring. The general data quality
from citizen-science initiatives has been shown to be high (Kosmala et al.,
2016), though the need for statistical methods to account for different data
biases prevails. Isaac et al. (2014) discuss bias-correction approaches for
ecological trend estimates, and conclude that opportunistic data would be
further enhanced if information about the sampling effort at the points of
data collection can be captured. They further identify four key aspects of
biases induced by volunteer sampling, due to uneven sampling in (1) time, (2)
space, (3) effort per visit, and (4) uneven detectability of species
individuals, which can vary by observer, species and land cover. For example,
detectability could be higher in open landscape by contrast with dense forest.
Moreover, the expertise of individual observers for detecting and identifying
species may increase over time due to increasing experience and knowledge-
sharing with other observers, thus reducing observation biases (Kelling et
al., 2015; Johnston et al., 2018).

Using eBird, statistical modelling of sampling effort expressed through
various criteria was discussed by Tang et al. (2021). More generally, a large
body of literature has emerged on how to characterize heterogeneous sampling
effort and quantify it from available data (e.g., Gelfand and Shirota, 2019;
Fink et al., 2020; Johnston et al., 2021, 2023).

Another issue arises when only the presence of species is reported but not
their absence (_presence-only data_), such that areas and times without any
reported occurrences could correspond either to true absence of the species or
to presence of the species but absence of observers. The BBS and eBird
datasets do not explicitly report absences, but both rely on observation
protocols requiring that all species detected at the sampled location and
identified by the observer must be reported. In eBird, data entries satisfying
this protocol are known as _checklist data_ , with a unique space-time
coordinate associated with the checklist, and they account for the vast
majority of eBird entries. Then, a common approach in data analysis that we
also adopt here is to consider species not reported within a checklist as
being absent from the space-time location of the checklist. By systematically
adding such (pseudo-)absence information to the dataset, we obtain a binary
observation for each combination of checklist and species, with $1$
corresponding to presence and $0$ corresponding to (pseudo-)absence of the
species. This allows working with a binomial likelihood to estimate the
probability of presence of the species in terms of available predictors.

###  1.3 Extreme-value analysis for migratory bird arrivals

Extreme-value theory (EVT, Coles, 2001) is a branch of probability theory and
statistics that provides a theoretically justified framework to model extreme
events, i.e., the tails of data distributions. It has been extensively used in
fields such as finance and environmental sciences, especially for climate
data. Applications in ecology have so far mostly focused on extremes of
abiotic environmental processes (Gaines and Denny, 1993; Katz et al., 2005),
which can have strong influence on biotic processes, for example when very low
winter temperatures contribute to limit the extent of outbreaks of forest
pests (Thibaud et al., 2016). However, direct applications of EVT to biotic
variables are scarce, which could be due to the strongly discretised nature of
species observation data, often collected in the form of presence-absence
information or relatively small occurrence counts. By contrast, standard EVT
is more focused on variables measured on a continuous scale. Recent examples
of studies leveraging EVT aim to model species accumulation curves (Borda-de
Água et al., 2021), extremes of species movements Wijeyakulasuriya et al.
(2019) or the first arrival dates of migratory birds (Wijeyakulasuriya et al.,
2023), similar to our goal here. Moreover, the statistical approaches
implemented in these studies do not explicitly take into consideration the
specific observation biases inherent to the ecological datasets that were
used, which is the focus of our work here.

First arrival dates of migratory birds at their destination are an example of
phenological events (i.e., of recurring life history events of species, such
as migration or breeding). Data arising on continuous scales, such as seasonal
timing of an event, are more common for phenological events, such that the
analysis of extreme events using EVT is appropriate.

Bird migration is a complex phenomenon taking place at large geographic scales
(Somveille et al., 2015), and with possible temporal trends in migration
patterns resulting from changes in climate, land cover, as well as other
ecological and evolutionary factors (Cotton, 2003; Conklin et al., 2021). The
study of migration patterns, and of arrival times at the breeding location in
particular, has attracted strong interest (Linden, 2011; Youngflesh et al.,
2021). Approximately 20% of all bird species are concerned by migration, which
allows birds to adapt to seasonal cycles in climatic stresses and in the
availability of resources such as food, especially during their breeding
period.

Previous studies often use rather complex approaches to define and estimate a
date that can be viewed as representative for the arrival of birds at their
breeding location. For example, for a given area and year, the representative
date could be based on both the first arrival date of an individual bird and
the first local maximum of detection probabilities of the species estimated
from eBird data (Youngflesh et al., 2021). In our approach, we focus only on
the first arrival date for each pixel of a mesh covering the study area, i.e.,
we model the minimum of all dates of occurrence of the species during the
year. We use the generalized extreme-value distribution suggested by EVT to
model this extreme value of the sample of all observations of the year,
similar to Wijeyakulasuriya et al. (2023).

In the northeastern US, the study region we have chosen for this work, birds
arrive during spring from more southern regions where they have spent
overwintering. Our goal is to analyse how the probability distribution of the
date of the first arrivals of birds varies across species and space, and in
response to factors related to climate and land cover. The mapping of the
arrival at breeding locations should be restricted to locations where birds
actually are present during the year, i.e., locations that provide a suitable
environment for birds to breed. Therefore, the considered locations must be
part of the spatial area forming the ecological niche occupied by the bird
species for breeding. Therefore, we will design our statistical model to allow
us to predict both the niche and the times of first arrival. At locations
where the probability of presence of a species (i.e., of _occupancy_) is very
low according to the niche model, we will not map arrival times.

To our knowledge, the modelling approach we will develop here is the first
that aims to appropriately capture the interplay of the niche, the sampling
effort and the observed first arrival dates with the goal to provide debiased
predictions for the true first arrival dates. In extension to the work of
Wijeyakulasuriya et al. (2023), we take into account the role of the niche and
the sampling effort when predicting first arrival dates. The approach of
Youngflesh et al. (2021) is based on a different metric for the arrival dates,
with a less direct and intuitive interpretation, and it does not account for
sampling effort except for choosing a study area with an relatively high
overall sampling effort in eBird.

###  1.4 Bayesian modelling of complex ecological data

Observational data in ecology often do not directly measure the latent
processes of interest, such as the ecological niche, the sampling effort and
the timing of a phenological event in our case. This is due to observation
biases and complex interactions among such processes. Moreover, different data
sources (e.g., eBird and BBS) can contribute complementary information about
the same latent process (e.g., the ecological niche). The data we use here
provide concrete illustration. The eBird checklists are available for a large
number of spatiotemporal locations chosen by the observers and provide a
generally good spatiotemporal coverage of the study area, especially for the
most recent years in the study period. By contrast, the BBS observations are
only available along predefined routes of length around 40 km, with up to 50
stops separated by around $800$ m at which observers can report occurrence
numbers of detected bird species. Time intervals for observation are also
imposed by the study protocol. Therefore, observation always takes place at
prespecified spatiotemporal locations in BBS, with the land-cover type at
those locations marked by the presence of a usually relatively large road. As
a consequence, BBS data may be less representative of all possible land-cover
types in comparison to eBird, and also do not provide direct information about
events taking place outside the prescribed observation interval, such as first
arrivals. On the other hand, within BBS many routes have been sampled during
each year since the 1970s, such that the temporal coverage of BBS over the
full study period is more homogeneous and complete as compared with eBird. In
summary, combining information from both datasets offers the possibility of
improved inferences for properties of the niche and of phenological events of
bird species.

With Bayesian hierarchical models (BHMs, Banerjee et al., 2003), data are
assumed to have been generated conditional on latent processes, which in turn
are conditioned on hyperparameters (e.g., variance or spatialcorrelation
range). This makes it possible to account for spatial patterns, complex
relationships and uncertainties. Using Bayes’ Theorem, the combination of
prior information based on expert knowledge and the data likelihood results in
posterior distributions for latent processes and hyperparameters that reflect
updated beliefs about the processes of interest given the observed data van de
Schoot et al. (2021).

An interesting application of BHMs in ecology is the fusion of presence-only
data, available through large datasets but with strongly heterogeneous and
unknown sampling effort, with presence-absence data, available through smaller
datasets but with known sampling effort, to infer species distribution maps
(e.g., Gelfand and Shirota, 2019). An important element of such approaches is
the inclusion of a latent Gaussian field that represents spatial heterogeneity
in preferential sampling, similar to the foundational work of Diggle et al.
(2010).

To infer our complex BHM for mapping the niche, observation effort and first
dates of arrival, we devise a Markov Chain Monte Carlo scheme with Vecchia
likelihood approximations and Metropolis-adjusted Langevin Algorithm updates
for the latent Gaussian fields.

###  1.5 Outline of the paper

In the remainder of the paper, Section 2 presents the datasets we use and
extensive preprocessing steps. In Section 3, we recall basics of EVT and
introduce the Bayesian hierarchical model we use to identify the spatial
variability with respect to three aspects (sampling effort; niche of a given
species; date of arrival of the first individuals during spring migration for
a given species). Details of latent Gaussian process models and Markov-Chain
Monte-Carlo inference are presented in Section 4, together with a simulation
study that confirms statistical identifiability of the model components we
seek to infer from data. Results for estimated arrival dates and other model
components are presented for a selection of species in Section 5. A discussion
of results and an outlook towards relevant research topics related to our
approach and to citizen-science data in general concludes the paper in Section
6, where we also highlight the potential utility of EVT in ecology.

##  2 Datasets and preprocessing steps

###  2.1 Bird observation data

We extracted bird observation data from eBird and BBS databases for the study
period 2001–2021. Both databases provide separate data files for each state in
the selected study area in the Northeastern US composed of the following
states: Connecticut, Maine, Massachusetts, New Hampshire, New Jersey, New
York, Pennsylvania, Rhode Island, Vermont. We merged extracted data into a
single dataset during preprocessing. We extracted the full observation record
of the eBird Basic Dataset (EBD), where each entry reports the observation of
one or several individuals of a species, and the following attributes were
retained: name of species, count of observed individuals, longitude, latitude,
date, duration of observation (i.e., temporal sampling effort), and a flag
indicating if the observation is part of a so-called checklist. Recall that a
checklist is an observation event where observers report all species they
observe, provided they succeed in identifying the species. For BBS, we
extracted information about observed species at the level of the survey route
and year.

First arrivals of migratory birds in the study area are known to take place at
the earliest towards end of March. We here identify migratory species as
having no observations during the winter months December and January in eBird.
For the migratory species characterized by this property, we remove a few
February occurrences that are understood as carryover from the preceding
autumn season. Only migratory species that were also observed in BBS are kept.
This selection procedure led us to consider around 50 species as migratory,
and we apply our modelling approach to the ten species with the largest number
of reported presences in eBird: Red-eyed Vireo, Eastern Wood-Pewee, Chimney
Swift, Great Crested Flycatcher, Warbling Vireo, Veery, Chestnut-sided
Warbler, Magnolia Warbler, Purple Martin, Blackburnian Warbler, in decreasing
order of occurrence entries.

The following two data tables were generated from eBird:

  * •

_eBird checklists:_ each row corresponds to one checklist in eBird, with
attributes longitude, latitude, year and duration of observation.

  * •

_eBird species occurrence:_ for each migratory species, the table contains the
same number of rows as there are checklists, and each row contains the
attributes of the checklist, the name of the species, and the binary presence-
absence flag to indicate if the species was observed. Species observations
flagged as “flyover” in eBird (i.e., where birds clearly did not develop any
nesting or breeding activity near the observation location) were declared as
absences.

We define four response variables for the components of our regression model,
based on BBS species counts, and on eBird data appropriately aggregated to a
regular pixel grid, here with pixel width fixed to 20km. This pixel size is
similar to the spatial mapping unit used for most official eBird
communications. Response variables are as follows:

  * •

_BBS species counts:_ the response variable is the number of observed birds,
available for each combination of migratory species, route and year, where
only routes surveyed during the specific year are considered; we use the
Poisson distribution for this variable.

  * •

_eBird checklist count:_ the response variable is the number of available
checklists, available for each configuration of year and pixel; we use the
Poisson distribution for this variable;

  * •

_eBird occurrence count:_ the response variable is the number of checklists
(among all available checklists) for which the species was reported present,
available for each configuration of species, year and pixel; we use the
binomial distribution for this variable;

  * •

_eBird first arrival dates:_ the minimum date of observation, expressed as the
number of days since January 1st, available for each configuration of species,
year and pixel where at least one observation of the species occurred; we use
the Generalized Extreme-Value distribution (GEV) for the negated variable,
which therefore represents a maximum.

Illustrations of some of the above datasets are provided in the schematic in
Figure 3.

###  2.2 Effect of observation effort on first arrival dates

Figure 1 shows observed dates of first arrival plotted against positive
checklist counts, for each pixel-year combination and several species.
Observed dates tend to occur clearly later during the year when less
checklists are available, that is, when sampling effort was lower in space.
This provides us with strong motivation for modelling the influence of
sampling effort on observed dates, and then use the estimated relationship to
provide bias-corrected predictions for true dates by setting the sampling
effort to a very high value during the prediction step.

Figure 2 further explores whether the duration of observation during a
checklist helps explain the residuals obtained after fitting a local
regression curve in Figure 1. We detect a slight effect of duration for the
two example species highlighted in Figure 2. This suggests that both the
number of checklists and the duration of observation per checklist are
relevant components of the overall sampling effort that could contribute to
modulate the distribution of the date of first arrival of birds at their
breeding location, with higher effort associated with earlier observation of
the first arrival.

To explain the remaining variability in the residuals in Figures 1 and 1, we
conjecture that other factors related to spatial variation in bird phenology
and to spatial variation in occupancy probabilities (i.e., the spatial “range”
of the species) are relevant. We will construct our model to allow us to
estimate spatial effects for these properties and to use them to correct
predicted dates of first arrival.

Figure 1:  Scatterplots of the logarithm of positive checklist counts (x-axes)
and dates of first arrival (y-axes) observed at each pixel-year combination
for the following species (top left to bottom right): Red-eyed Vireo, Chimney
Swift, Warbling Vireo and Chestnut-sided Warbler. The larger dots show the 20
binned estimates, the black line a smooth fitted curve, and the shaded region
the $95\%$ pointwise confidence intervals of the curve.

Figure 2:  Scatterplots of the inverse duration of observation during
checklists (x-axis) and the residuals obtained from fitting a smooth local
regression curve to the dates of first arrival as a function of the logarithm
of positive checklist counts (y-axis), for two species: Willow Flycatcher
(left); Orchard Oriole (right).

###  2.3 Climate and land cover data

The North Atlantic Oscillation (NAO) Index describes changes in the strength
of two recurring pressure patterns in the atmosphere over the North Atlantic:
a low near Iceland, and a high near the Azores Islands. Positive NAO index
values indicate these features are strong, creating a big pressure difference
between them. Strongly positive values are linked to warm conditions across
the Eastern US and Northern Europe, and cold conditions across Southern
Europe. Negative NAO index values indicate these features are relatively weak,
and the pressure difference between them is smaller. Strongly negative values
are linked to cold conditions in the Eastern US and Northern Europe, and warm
conditions in Southern Europe. As a covariate at annual temporal resolution,
we here use the average of NAO for the month of March, which was identified as
a relevant climate predictor for bird migration timing by Wijeyakulasuriya et
al. (2023).

We use the 2021 National Land Cover Database of the US to extract land cover
information. To facilitate model construction and identifiability, we will not
use land cover variables as covariates within our BHM. Instead, we will
compare posterior maps of niche, observation effort and first date of arrival
with land cover maps to detect and interpret significant correlations between
model components and land cover. Land cover proportions for the pixel grid are
shown in the Appendix in Figure 10 for four land-cover types, defined by
appropriately merging original NLCD categories: Developed areas (including all
areas with buildings and infrastructure such as roads); Forest; Vegetation
(excluding forest but including planted and cultivated land); Water (including
wetlands).

Figure 3:  Schematic summarising the data sources, the data preprocessing
steps and the four components in the data layer of the BHM (panels labelled
with numbers 1 to 4). The dashed lines illustrate our modelling choices and
how model components interact, but they do not necessarily show the chronology
of how these components are estimated, since the whole model is estimated
jointly using a fully Bayesian approach. The plots in the panels numbered ‘1’,
‘2’ and ‘4’ show the number of available eBird checklists, the median duration
spent per checklist in each pixel, and the number of stops along a route from
the BBS data, respectively, in 2001 (left display in each panel) and 2021
(right display).

Figure 3 shows a schematic summarising the different data sources and
preprocessing steps of our approach, and the general workflow of the model.
The plot in Panel 1 illustrates a strong increase in eBird checklists and in
spatial coverage between the beginning and the end of the study period; the
plot in Panel 2 highlights that the median duration per checklist varies
substantially in both space and time. In fact, in 2001 we see longer durations
in specific spatial regions, while for 2021 the durations are more homogeneous
across larger areas, and there is a general overall increase in activity. The
display in Panel 3 shows the occurrence dates of observations of a species
pooled together for a pixel-year, and we extract the minimum day of the year,
which will serve as an observation for the GEV response variable in our BHM.
Finally, the figure in Panel 4 shows the number of BBS stops in 2001 and 2021.

##  3 Bayesian hierchical model

###  3.1 General setting and notations

Figure 3 also highlights the four components in the data layer of the BHM we
propose, and we will further detail its notations and structure in the
following.

We write $\mathcal{S}\times\mathcal{T}$ for the space-time study domain, and
denote by $A_{i}$, $i=1,\ldots,D\times T$, the set of non-intersecting space-
time cells based on a division of the study area $\mathcal{S}$ into $D$
pixels, replicated over $T=21$ years. We identify each cell $A_{i}$ with a
representative location $(\bm{s}_{i},t_{i})$, such as its barycenter, where
$\bm{s}_{i}\in\mathcal{S}$ and $t_{i}=1,\dots,T$.

Let $N^{\text{ckl}}_{i}$ be the number of eBird checklists in $A_{i}$. If
$N^{\text{ckl}}_{i}>0$, let $N^{\text{spc}}_{i}$ denote the sum of the binary
indicators for the detection of a given species from each of the
$N^{\text{ckl}}_{i}$ checklists. Binary indicators take value $1$ if the
species was detected and $0$ otherwise. Given $N^{\text{spc}}_{i}>0$, we write
$\bm{Y}_{i}=(Y_{i,1},\ldots,Y_{i,N^{\text{ckl}}_{i}})\in[1,365]^{N^{\text{spc}}_{i}}$
to denote the corresponding arrival dates in $\bm{s}_{i}$ within the year
$t_{i}$ from the $N^{\text{spc}}_{i}$ observed species occurrences. Based on
these arrival dates, we calculate
$Z^{\prime}_{i}=\min_{j}^{N_{i}}Y_{i,j}=-\max_{j}^{N_{i}}-Y_{i,j}$, the date
of first arrival in year $t_{i}$ within $\bm{s}_{i}$. We transform
$Z_{i}=\log\\{(Z_{i}^{\prime}+365)/-Z_{i}^{\prime}\\}$ to reflect that
$Z_{i}^{\prime}$ has a natural lower bound at $-365$, and we model
$Z_{i}\in\mathbb{R}$ instead.

To simplify notation in this section, we sometimes concatenate our data
responses into vector form, so for example we write
$\bm{N}^{\text{ckl}}=(N^{\text{ckl}}_{1},\dots,N^{\text{ckl}}_{D\times T})$.

###  3.2 Extreme-value distribution for first arrival dates

Provided that the partitioning of $\mathcal{S}$ provides enough ‘independent’
arrival times within each pixel-year $A_{i}$, Extreme-Value Theory (EVT,
Coles, 2001) motivates using a Generalised Extreme-Value distribution (GEV) to
model the transformed first date of arrival in $A_{i}$. Invoking EVT, we
assume convergence in distribution of each $Z_{i}$ to a limiting GEV
distribution as $N_{i}^{\text{spc}}\rightarrow\infty$. This suggests using the
GEV distribution as the model for the finite-sample maxima $Z_{i}$.

To state the relevant theory more precisely, let $M_{n}$ denote the maxima of
independent and identically distributed random variables $X_{1},\dots,X_{n}$.
If $M_{n}$ converges in distribution to a non-degenerate distribution after
linear renormalisation, then this distribution must be a GEV. In practice, for
fixed $n$ large enough (where $n=N_{i}^{\text{spc}}$ in our setting), one uses
the approximation

$\text{Pr}(M_{n}\leq
z)\approx\begin{cases}\exp\left\\{-\exp\left(-\frac{z-\mu}{\sigma}\right)\right\\},&\xi=0,\\\\[5.0pt]
\exp\left\\{-\left[1+\xi\left(\frac{z-\mu}{\sigma}\right)\right]^{-1/\xi}\right\\},&\xi\neq
0,\end{cases}$ (1)

where the right-hand side of (1) is the distribution function of the GEV,
defined on $\\{z:1+\xi(z-\mu)/\sigma>0\\}$ with location $\mu\in\mathbb{R}$,
scale $\sigma>0$ and shape $\xi\in\mathbb{R}$. The shape parameter determines
the tail behaviour of the GEV: for $\xi>0$, the GEV has a finite lower
endpoint of its support and an upper heavy tail of power-law form; for
$\xi<0$, the GEV has a finite upper endpoint but is unbounded below; for
$\xi=0$, the support of the GEV has no finite bounds, and the upper tail
exhibits an exponential decay rate. In our application, some sample sizes
$N_{i}^{\text{spc}}$ can be very small and dates $Y_{i,j}$,
$j=1,\ldots,N_{i}^{\text{spc}}$ can be dependent, but the three-parameter GEV
still provides a flexible model to accommodate relevant distributional
properties of the maxima $Z_{i}$, even if the theoretical asymptotics are not
reached. We incorporate linear predictors into the GEV regression equation
used to model $\bm{Z}$ in our BHM: one for the location $\mu$ and another for
the scale $\sigma$.

###  3.3 Data layer of the Bayesian hierarchical model

We construct a system of four regression equations with latent Gaussian
processes used in the linear predictors. Recall that we denote the response
variables, as described in Section 2.1, as follows: $N_{i}^{\text{ckl}}$ –
eBird checklist count per pixel-year i; $N_{j}^{\text{BBS}}$ – BBS counts of
the species of interest per year-route $j$; $N_{i}^{\text{spc}}$ – eBird
occurrence count for the species of interest per pixel-year i; $Z_{i}$ –
transformed negated first arrival date for the species of interest in eBird
per pixel-year $i$.

For the response variable $N_{i}^{\text{ckl}}$ reporting the checklist counts,
we weight the Poisson intensities according to the surface areas
$\mathcal{A_{i}}$ of the pixels (i.e., we use an offset of
$\log\mathcal{A_{i}}$ in the linear predictor), since some pixels at the
boundary of the study region have smaller area. For the response variable
$N_{j}^{\text{BBS}}$ reporting the BBS occurrence counts, observations are
made along a subset of 50 equidistant stops along each of the routes, and we
weight the Poisson intensities accordingly with the number of stops (between
$1$ and $50$) that were visited. To alleviate notations in the predictor
formulas below, we do not explicitly write down these weights related to pixel
surface areas and road lengths. Moreover, routes can span across several
pixels; therefore, we construct the Poisson intensity parameters as a weighted
mean of the parameters for the involved pixels with weights $w_{k}$ defined as
the percentage length inside each of the pixels $A_{k}$.

The hierarchical system of the four regression equations can be written as
follows:

$\displaystyle
N^{\text{BBS}}_{j}\mid\lambda^{\text{BBS}},\bm{\theta}_{\text{bbs}}$
$\displaystyle\sim\mathrm{Pois}\left\\{n_{j}^{\text{stop}}\times\sum_{k\in\text{route}_{j}}\omega_{k}\lambda^{\text{BBS}}(\bm{s}_{k};\bm{\theta}_{\text{bbs}})\right\\},$
$\displaystyle
N^{\text{ckl}}_{i}\mid\lambda^{\text{ckl}},\bm{\theta}_{\text{ckl}}$
$\displaystyle\sim\mathrm{Pois}\left\\{\lambda^{\text{ckl}}(\bm{s}_{i},t_{i};\bm{\theta}_{\text{ckl}})\right\\},$
$\displaystyle N^{\text{spc}}_{i}\mid p^{\text{spc}},\bm{\theta}_{\text{spc}}$
$\displaystyle\sim\mathrm{Bin}\\{N^{\text{ckl}}_{i},p^{\text{spc}}(\bm{s}_{i},t_{i};\bm{\theta}_{\text{spc}})\\},$
$\displaystyle Z_{i}\mid\mu,\bm{\theta}^{\mu},\sigma,\bm{\theta}_{\sigma}$
$\displaystyle\sim\mathrm{GEV}\\{\mu(\bm{s}_{i},t_{i};\bm{\theta}_{\mu}),\sigma(\bm{s}_{i};\bm{\theta}_{\sigma}),\xi\\},$

where

$\bm{\theta}_{\text{bbs}},\bm{\theta}_{\text{ckl}},\bm{\theta}_{\text{spc}},\bm{\theta}_{\mu},\bm{\theta}_{\sigma}\sim\text{Hyperpriors}$

are hyperparameters for the random predictors $\lambda^{\text{BBS}}$,
$\lambda^{\text{ckl}}$, $p^{\text{spc}}$, $\mu$ and $\sigma$ that govern the
different model components; we discuss their specifics next. Checklist
occurrences form a point pattern, and due to the structure of the ckl-
regression equation, we model them as a spatiotemporal log-Gaussian Cox point
process (Møller et al., 1998), discretised according to the pixels $A_{i}$.
Specifics of $\bm{\theta}$ are discussed in the next subsection.

###  3.4 Structure of latent Gaussian processes

We write $X(\bigcdot)$ to denote any latent Gaussian effect, here either
indexed by pixels $\bm{s}\in\mathcal{S}$ or years $t\in\mathcal{T}$. We use
Gaussian process priors for all latent effects. For fixed effects, such as
intercepts, we assign a centered normal prior $\mathcal{N}(0,10)$ with high
variance. We further use the exponential covariance function to parametrise
the spatial effects, given as follows for two locations $\bm{s}_{1}$ and
$\bm{s}_{2}$:

$\mathrm{Cov}\\{X(\bm{s}_{1}),X(\bm{s}_{2})\\}=\sigma^{2}\exp(-||\bm{s}_{1}-\bm{s}_{2}||/\kappa),\quad\sigma,\kappa>0,$

with Euclidean distance $||\cdot||$, standard deviation $\sigma>0$ and range
parameter$\kappa>0$.

The four spatial Gaussian process priors included in our model are

$\displaystyle
X^{\text{effort}}(\bigcdot)\sim\mathcal{GP}(\bm{\omega}_{1}),\quad
X^{\text{niche}}(\bigcdot)\sim\mathcal{GP}(\bm{\omega}_{2}),$ $\displaystyle
X^{\text{GEV-}\mu}(\bigcdot)\sim\mathcal{GP}(\bm{\omega}_{3}),\quad
X^{\text{GEV-}\sigma}(\bigcdot)\sim\mathcal{GP}(\bm{\omega}_{4}),$ (2)

where $\bm{\omega}_{1}$, $\bm{\omega}_{2}$, $\bm{\omega}_{3}$ and
$\bm{\omega}_{4}$ contain individual range and standard deviation parameters
that control the spatial dependence of each process. These hyperparameters are
assigned identical joint PC priors (Fuglstad et al., 2018). We also impose
sum-to-zero constraints on all spatial effects to aid with identifiability.
Spatial patterns of the niche of the species are captured by
$X^{\mathrm{niche}}$, whereas spatial patterns of the sampling effort are
modeled through $X^{\mathrm{effort}}$.

The structure of link functions and linear predictors, explained in detail
subsequently, is:

$\displaystyle\log\lambda^{\text{BBS}}(\bm{s}_{i})$
$\displaystyle=\beta_{0}^{\text{BBS}}+X^{\text{niche}}(\bm{s}_{i}),$
$\displaystyle\log\lambda^{\text{ckl}}(\bm{s}_{i},t_{i})$
$\displaystyle=\beta_{0}^{\text{ckl}}+X^{\text{year}}(t_{i})+X^{\text{pref}}(\bm{s}_{i}),$
$\displaystyle\mathrm{cloglog}\\{p_{i}^{\text{spc}}(\bm{s}_{i},t_{i})\\}$
$\displaystyle=\beta_{0}^{\text{spc}}+X^{\text{niche}}(\bm{s}_{i})+\dfrac{1}{\beta^{\text{act}}d_{\bm{s}_{i},t_{i}}},$
$\displaystyle\mu(\bm{s}_{i},t_{i})$
$\displaystyle=g\\{\beta_{0}^{\text{GEV-}\mu}+X^{\text{GEV-}\mu}(\bm{s}_{i})+\beta_{1}^{\text{GEV-}\mu}\mathrm{NAO_{t_{i}}}$
$\displaystyle\quad\quad+\theta^{\text{niche-
GEV}}X^{\text{niche}}(\bm{s}_{i}),x_{\text{effort}}(\bm{s}_{i},t_{i})\\},$
$\displaystyle\log\sigma(\bm{s}_{i})$
$\displaystyle=\beta_{0}^{\text{GEV-}\sigma}+X^{\text{GEV-}\sigma}(\bm{s}_{i}),$

where the temporal latent effect $X^{\text{year}}$ captures the variability in
checklist counts across years and has prior
$X^{\text{year}}(\cdot)\sim\mathcal{GP}(\bm{\omega}_{5})$. The climate
covariable, i.e. the NAO value for year $t_{i}$, is denoted as
$\mathrm{NAO_{t_{i}}}$, and we include it inside the link function $g$
(detailed below) as a fixed effect modulating the upper bound of the GEV
location parameter. The scaling parameters $\theta^{\text{niche-
spc}}\in\mathbb{R}$ and $\theta^{\text{niche-GEV}}$ control the scale of
sharing of the niche spatial effect. The complementary log-log link function
$\text{cloglog}(p)=\log(-\log(1-p))$ for binomial probabilities allows
interpreting them as the probabilities of observing a count of at least $1$ in
a Poisson model (which is equal to $1-\exp(-\lambda))$ with Poisson intensity
$\lambda=\exp(-x))$). This ensures consistency with the Poisson model for the
BBS counts.

We differentiate between two dimensions of sampling effort that we combine
into $x_{\text{effort}}$. The first one is modelled through the space-time
intensity $\lambda^{\text{ckl}}$ of checklists, which we call the preference.
We hypothesise that it influences the first arrival dates $Z_{i}$, where a
higher preference likely leads to more and thus earlier sightings in the year.
The second dimension of sampling effort that we consider here is the median
duration of observation per checklist in a pixel-year $d_{\bm{s}_{i},t_{i}}$,
which we call the activity. This effect could influence both the binomial
probability of observing the species of interest in a checklist, and again the
first arrival dates. A higher activity could increase the probability of
observing the species in a checklist, which could also induce an earlier
observed arrival date. We combine the two effects preference and activity into

$x_{\text{effort}}(\bm{s}_{i},t_{i})=\theta^{\text{pref}}\log\lambda^{\text{ckl}}(\bm{s}_{i},t_{i})+\dfrac{\theta^{\text{act}}}{d_{\bm{s}_{i},t_{i}}}$
(3)

with scaling hyperparameters $\theta^{\text{pref}}$ and $\theta^{\text{act}}$.
With this structure, setting an infinite duration $d_{\bm{s}_{i},t_{i}}$
causes the activity term to vanish.

The next subsection explains the choice we make for the nonlinear link
function $g$.

###  3.5 Saturated sharing towards the GEV predictor

Sharing latent processes across carefully chosen regression equations, as
implemented in Section 3.4, enables us to leverage data from several response
variables for appropriately identifying these four latent processes.
Specifically, we want to identify the effect of spatiotemporally heterogeneous
sampling effort and correct it in predictions, such as those of the true first
arrival dates.

The most common approach of sharing latent effects between regression
equations is via a linear structure where the linear predictor of one
regression equation has an additive component $X(\bigcdot)$ and the linear
predictor of another regression equation has additive component $\beta\times
X(\bigcdot)$, with a scaling factor $\beta\in\mathbb{R}$. This yields models
akin to the preferential-sampling framework of Diggle et al. (2010); Pati et
al. (2011), and to Koh et al. (2023); Yadav et al. (2023) in the context of
joint modelling of occurrences and sizes of wildfires or landslides. However,
such linear sharing is not flexible enough for sharing the sampling effort
towards the linear predictor of the $\mu$ parameter of the GEV distribution
for first arrival dates. Indeed, with linear sharing, letting the effort tend
to infinity would lead to first arrival dates tending to January 1st, the
smallest possible value, which is wrong and would prevent us from setting the
observation effort to a very high level for predicting the true first arrival
dates. As illustrated by the exploratory data analysis in Figure 1, we reach a
saturation of the effect of observation effort on first arrival dates at very
high levels. For this reason, we share $x_{\text{effort}}$ with the $\mu$
component of the GEV distribution through the following upper-bounded function
that always attains a finite saturation threshold:

$g(x_{\text{bound}},\,x_{\text{effort}})=\frac{x_{\text{bound}}}{1+\exp(-x_{\text{effort}})},$
(4)

where the first argument modulates the upper bound in the GEV location
parameter, which can be attained with unlimited (i.e., infinite) effort. The
second argument represents the sampling effort. If $x_{\text{effort}}$ is
close to $0$, then the $\mu$ component of the GEV is relatively insensitive to
its variations. Figure 4 illustrates that the class of functions (4) is
bounded and monotonically increasing, and it therefore allows the effect of
the sampling effort on the first arrival dates to saturate with high effort.

The structure we put into (4) is also inspired by species accumulation curves
in ecology (Colwell et al., 1994), but here we replace the species abundance
with the dates of first observed arrivals. As far as we are aware, such an
approach has not yet been considered in the literature.

Figure 4:  Plots of the nonlinear sharing function $g$ in (4) for varying
$x_{\text{effort}}$ (x-axis).

##  4 Simulation-based Bayesian inference

Fast and relatively accurate off-the-shelf implementations of Bayesian
hierarchical inference exist, such as the INLA framework (Rue et al., 2017).
However, due to certain specifics of our model, such as nonlinear sharing of
latent spatial effects and the inclusion of linear predictors in several
parameters of the GEV distribution, it is not feasible to use them here.
Instead, we develop a Markov-Chain Monte-Carlo sampler for simulation-based
inference. For the purpose of keeping the numerical linear algebra tractable
when handling the multivariate Gaussian vectors of dimension $D=1268$ for the
pixel-based latent Gaussian fields, we implement the Vecchia approximation for
the exponential covariance model.

###  4.1 Vecchia approximation for Gaussian components

Due to the large number of space-time cells in our case study, keeping spatial
latent Gaussian effects and their conditional distributions numerically
tractable is challenging. Various solutions have been proposed to tackle this
fundamental problem in spatial statistics with large-dimensional data,
including Gauss–Markov random fields used to approximate certain types of
covariance functions based on Stochastic Partial Differential Equations
(Lindgren et al., 2011, SPDE approach,), or Nearest-Neighbor Gaussian
processes (NNGPs, Datta et al., 2016). We here use another flexible and
popular option to construct approximate Gauss–Markov representations of any
Gaussian covariance function through the Vecchia approximation (Vecchia,
1988), which sparsifies the Cholesky factor of the precision matrix. The
Vecchia approximation is known to be quite accurate in practice and still
represents a valid Gaussian process. Writing $x_{i}$ to denote a Gaussian
field $X$ evaluated at $s_{i}\in\mathcal{S}$, the exact joint density of the
set of observations $\\{x_{1},\dots,x_{D}\\}$ can be written as

$f(x_{1},\dots,x_{D})=f(x_{1})\prod^{D}_{i=2}f(x_{1})f(x_{i}\mid
x_{i-1},\dots,x_{1}).$

The Vecchia approximation assumes that

$\displaystyle f(x_{1},\dots,x_{D})$
$\displaystyle\approx\hat{f}(x_{1},\dots,x_{D})$
$\displaystyle=f(x_{1})\prod^{D}_{i=2}f(x_{1})f(x_{i}\mid
x_{i-1},\dots,x_{\max\\{i-k,1\\}}),$ (5)

i.e., it deliberately reduces the conditioning history in the conditional
densities, which decreases the computational complexity of evaluating the
density from $\mathcal{O}(D^{3})$ to $\mathcal{O}(Dk^{3})$. This reduction is
considerable if $k\ll D$, as in our setting with number of spatial pixels
$D>1000$. An active research area (Katzfuss and Guinness, 2021) is to
understand how one should order the conditioning variables in (4.1) to obtain
approximations that satisfy certain criteria of optimality. In his seminal
work, Vecchia suggested to order the conditioning variables lexicographically
based on their spatial coordinates, but this approach has been shown to be
inefficient (Guinness, 2018). Instead, we follow the ordering proposed by the
latter study, which ensures that points are ordered in a quasi-random fashion.
We select $k=5$ conditioning locations in our case study to make inference
feasible. The simulation study detailed below confirmed that this is
reasonable.

###  4.2 Proposal scheme for Markov-Chain Monte-Carlo

We draw posterior inference about the parameters using a general Markov-Chain
Monte-Carlo (MCMC) sampling scheme with customised Metropolis-Hastings (MH)
updates. The hyperparameters are updated using Gibbs sampling. The
computationally most intensive part of our scheme involves sampling from the
latent Gaussian components, which are updated jointly through a Metropolis
Adjusted Langevin Algorithm (MALA; Roberts and Rosenthal, 2001) designed to
speed up mixing of the chains. We write $\bm{x}$ for the latent Gaussian
component of interest, $\bm{x}^{\star}$ for a proposal for that component, and
$\bm{y}$ for the data vector. A discrete approximation of the Langevin
diffusion implies that MALA is a MH algorithm with proposal distribution

$(\bm{x}|\bm{x}^{\star})\sim
N_{p}\left(\bm{x}^{\star}+\dfrac{\delta^{2}}{2}\nabla_{\bm{\theta}}\log\pi(\bm{x}^{\star}|\bm{y}),\delta^{2}\bf{W}_{p}\right),$
(6)

where tuning parameters are the pre-whitening matrix $\bf{W}_{p}$, used to
account for the correlation between parameters within the component, and the
stepsize $\delta$. For example, if we set
$\bm{x}=\\{X^{\text{niche}}(\bm{s}_{1}),\dots,X^{\text{niche}}(\bm{s}_{D})\\}$,
the gradient in (6), when updating this component, is

$\displaystyle\nabla_{\bm{x}}\log\pi(\bm{x}|\bm{y})\approx$
$\displaystyle\nabla_{\bm{x}}\log\pi_{\text{BBS}}(\bm{N}^{\text{BBS}}\mid\bm{x})+\nabla_{\bm{x}}\log\pi_{\text{spc}}(\bm{N}^{\text{spc}}\mid\bm{x})$
$\displaystyle+\nabla_{\bm{x}}\log\pi_{\text{GEV}}(\bm{Z}\mid\bm{x})+\nabla_{\bm{x}}\log\pi_{\text{Vecchia}}(\bm{x}),$

where the first three terms represent the gradients of log-likelihoods for the
model components that $X^{\text{niche}}$ is included in, and

$\nabla_{\bm{x}}\log\pi_{\text{Vecchia}}(\bm{x})=-\tilde{\Sigma}_{\bm{\omega}_{2}}^{-1}\bm{x},$

where the $D\times D$ precision matrix $\tilde{\Sigma}_{\bm{\omega}_{2}}^{-1}$
is the Vecchia-approximated precision of the Gaussian field prior assigned to
$X^{\text{niche}}$, as detailed in §4.1 and implied by (4.1). The precision
matrix is parametrised by the vector $\bm{\omega}_{2}=(\sigma,\kappa)$, see
§3.4, to which we assign a hyperprior.

For our data application, we parallelise and run separate chains for each
species. We generate $120^{\prime}000$ posterior samples and discard the first
$80^{\prime}000$ as the burn-in period. We then perform a standard thinning
operation of the Markov chains and keep one of four consecutive samples, thus
obtaining $10^{\prime}000$ samples overall to perform posterior inference. We
divide the number of posterior samples and the burn-in period by two for our
simulation study for computational reasons. We monitor convergence and mixing
of the chains through trace plots and by assessing the effective sample size.
Code implemented for the MCMC procedure is available at
https://github.com/kohrrelation/mcmc_birds. To obtain a shorter burn-in period
and to avoid having initial values too far from the region where posterior
mass concentrates, we ran test chains that model each model component
separately first (where we can estimate all parameters except the sharing
parameters). Then, we used the last values of those chains as starting values
when fitting the full model. Calculations were performed on
UBELIXhttps://www.id.unibe.ch/hpc, the HPC cluster at the University of Bern,
and the computational time of our algorithm on there is roughly 6 seconds per
iteration.

###  4.3 Simulation study for sanity checking

We check by simulation that the model parameters (latent Gaussian fields and
hyperparameters) are identifiable and can be estimated appropriately through
posterior mean point estimates obtained from the above Bayesian inference
procedure. We especially wish to check whether the spatiotemporal field of
first arrival dates and the niche of the species are identifiable. This
verification seems important for the proposed models owing to their high
structural complexity with a relatively large number of parameters. For the
simulation study, we mimic our data application and choose the same spatial
domain and temporal replicates, i.e., we simulate datasets at the $D=1268$
pixels over the $T=21$ years from 2001 to 2021 with the same covariates and
data structures, including the climate covariate (NAO), the number of BBS
stops and the median duration in each pixel-year. The values we set for latent
Gaussian processes and hyperparameters are similar to those obtained when
having fitted this model to a real species. We estimate all the parameters for
$100$ independent replicates of such data.

The different model components are generally well identified and reproduced by
our estimation procedure. The results in the Appendix (Figure 11) suggest that
the posterior means of the median first arrival dates estimated with this
procedure are indeed essentially unbiased, and the truth is mostly well
recovered. Importantly, the temporal evolution and spatial patterns are also
captured satisfactorily. For example, in our simulation setting we have
spatially more homogeneous and later first arrivals for the earlier years such
as 2003, and the estimated posterior means reproduce this property. The
boxplots in Figure 11 suggest that first arrivals for certain pixels are
harder to estimate than others, but the truth falls within the predicted
interquartile ranges in most cases.

##  5 Results

We now apply the model constructed and checked in Sections 3 and 4 to the data
described in Section 2, where we run separate models for the $10$ species with
largest occurrence numbers in eBird. Due to space constraints, we report a
selection of results for four species: Chimney Swift, Great Crested
Flycatcher, Magnolia Warbler and Purple Martin.

The MCMC trace plots shown in the Appendix (Figure 12) for one of the
aforementioned species indicate that the chain is mixing relatively well,
though the shape parameter $\xi$ still displays autocorrelation after the
burn-in period and thinning.

Figure 5 illustrates that there are more checklists in the more recent years,
and the model appropriately captures these temporal differences and also the
spatial pattern in observed checklist numbers.

The posterior means of the sharing parameters are all positive for
$\theta^{\text{pref}}$ and all negative for $\theta^{\text{act}}$ (see Table
1, which corroborates our exploratory findings from Section 2 that higher
sampling effort leads to earlier observed dates of first arrival.

Both $\theta^{\text{niche-GEV}}$ and $\beta_{1}^{\text{GEV-}\mu}$ behave
similarly (not shown here). Their posterior means are small, e.g., the latter
effect has posterior means hovering around $[0.02,0.04]$, which shifts the
location parameter of the GEV by one to two days for a year with high NAO of
$+1^{\circ}$C. Higher values for NAO and temperature anomaly correspond to
warmer weather, which could explain earlier arrivals. However, for the
majority of the species we investigate these slightly positive NAO
coefficients appear not to be significant, as their credible intervals include
zero. We conclude that the slight NAO effect on first arrival dates could be
present for some species but there is no strong signal of a significant NAO
effect across all species; Wijeyakulasuriya et al. (2023) obtained a similar
result when applying their analyses to the Magnolia Warbler.

The GEV shape parameters are negative for all species and have posterior means
around $-0.6$ to $-0.9$, which means that we estimate a fixed lower bound for
the possible first arrival dates.

Figure 5: Observed numbers of checklists, $N^{\text{ckl}}_{i}$ (left) against
model-based posterior mean estimates of the Poisson mean number of checklists,
$\lambda^{\text{ckl}}(\bm{s}_{i},t_{i})$ (right), using the model fitted to
Chimney Swift in 2001 (top) and 2021 (bottom).

###  5.1 Excursion sets of latent spatial fields

We study more precisely the spatial patterns in the posterior estimation of
the four latent spatial fields in (2) for the sampling effort, the species
niche, and the $\mu$ and $\sigma$ parameters of the distribution of the first
arrival dates. We highlight areas with very low or high values through
excursion set analysis.

To identify areas with very high or low values of the spatial random effects,
we utilise credible sets for excursion regions (Bolin and Lindgren, 2015). We
evaluate where the fields exceed or fall below certain thresholds. To
visualise excursion sets simultaneously for all values of the probability
level $\alpha$, Bolin and Lindgren (2015) introduced the positive and negative
excursion functions $F^{+}_{u}(s)=1-\inf\\{\alpha\mid
s\in\mathrm{E}^{+}_{u,\alpha}\\}\in[0,1]$ and
$F^{-}_{u}(s)=1-\inf\\{\alpha\mid s\in\mathrm{E}^{-}_{u,\alpha}\\}\in[0,1]$.
Figure 6 highlights these excursion functions for the species _Purple Martin_.
Sampling effort is high in the southeastern ares near the ocean and drops in
the further northwestern areas. The niche of the species tends to concentrates
in southwestern areas. The fields related to first arrival dates show more
fragmented spatial patterns for this species.

Figure 6: Excursion set analysis for the four latent spatial fields
$X^{\text{effort}}$ (top left), $X^{\text{niche}}$ (top right),
$X^{\text{GEV-}}\mu$ (bottom left) and $X^{\text{GEV-}}\sigma$ (bottom right)
estimated by the model fitted to the Purple Martin. Thresholds are fixed at
$2$ ,$1$, $0.1$ and $0.1$ for the four fields, using the same order.

###  5.2 Correlation of latent spatial fields with land cover

To further explain the patterns in the four latent spatial fields, we
calculate the correlation of the their posterior means with the land-cover
maps described in Section 2 and shown in the Appendix (Figure 10. This can
help interpretation by providing better understanding about the land-cover
types that drive the four latent processes for a given species.

Figure 7 shows these correlations for four species. The correlations between
the estimated latent field for sampling effort and land cover proportions is
very similar across species, which is expected since the checklist count data
for the corresponding regression equation are the same in the four models.
Sampling effort has strong positive correlation with Developed (built) areas,
weaker correlation with Water-dominated or Vegetation-dominated areas, and
negative correlation with Forest. Regarding the latent field for niche, it
shows different patterns for different species, for example with positive
correlations with Forest for species that are known to breed in forests (e.g.,
_Magnolia Warbler_), and with Developed land cover for species known to breed
in Developed areas (e.g., _Chimney Swift_). The latent field describing the
spatial variability in first arrival dates show generally weaker correlation
with specific land-cover types. For example, _Chimney Swift_ tends to settle
earlier in Developed and Vegetation areas in comparison to Water- and Forest-
dominated areas.

Figure 7: Correlations of land-cover proportions (x-axis) with posterior means
of the latent fields $X^{\text{effort}}$, $X^{\text{niche}}$,
$X^{\text{GEV-}}\mu$ and $X^{\text{GEV-}}\sigma$ (y-axis), for four species:
Chimney Swift (top left), Great Crested Flycatcher (top right), Magnolia
Warbler (bottom left) and Purple Martin (bottom right).

###  5.3 Predicting true first arrivals

To further illustrate the utility of our model, we predict the true first
arrival times in all pixels.

First, Figure 8 shows that the model predicts well the spatial variability in
the niche, i.e. the probability of occupancy. We use this component of the
model to avoid predicting first arrival dates at locations where the presence
of the species is uncertain in terms of the estimated niche. Therefore, we do
not map predictions of first arrival dates at pixel-years where the estimated
occupancy probability $p^{\text{spc}}\\{(\bm{s}_{i},t_{i})\\}$ is smaller than
$1\%$. This feature of our model avoids making predictions in areas where data
do not give a clear signal of breeding activity of birds.

For the prediction of the true first arrival dates, we set the values of
$\lambda^{\text{ckl}}(\bm{s}_{i},t_{i})$ and of $d_{\bm{s}_{i},t_{i}}$ to
values close to infinity in order to mimic infinite sampling effort, such that
we expect to retrieve the true and unbiased dates of first arrival. This
approach is illustrated in Figure 9 for the Great Crested Flycatcher, where we
also show that the model appropriately captures the spatial pattern in the
observed first arrival dates.

To give some more detail about how the prediction works, we arrow our focus to
two pixels, marked with triangles in Figure 8. The first pixel is chosen near
the southwestern boundary of our study region and has $83\%$ forest cover but
only $6\%$ developed land cover. It is located in a rather secluded region
with generally low sampling effort. In 2022, there were no reported sightings
of the four species of interest in this pixel. Table 1 (middle three columns)
highlights that the model can extrapolate to this pixel. It would also be
possible to draw posterior predictive samples for that pixel-year. The second
pixel is chosen in the middle of the study region. For example, a first
arrival was observed in 2022 for the Purple Martin on the $26^{\text{th}}$ of
June, but Table 1 (last three columns) shows that our model flags this pixel
as being uncertain to be located in the species’ niche, and that is therefore
possible that in general no breeding activity takes place in this pixel. This
illustrates the capability of our model to prevent inaccurate extrapolation to
spatial regions with low occupancy probability.

Figure 8: Empirical mean probability (left displays) of observing the species
Great Crested Flycatcher (top) and Purple Martin (bottom) in checklists
aggregated over all years; posterior mean probability (right displays) of
observing the species from the occupancy model by using the mean checklist
duration aggregated over all years as covariate. The two triangles are the two
pixels used in the predictions of first arrival dates detailed in Section 5.3.
Figure 9:  Spatial plots for the species Great Crested Flycatcher: Observed first arrivals of Great Crested Flycatcher in 2022 (top left), posterior mean of the GEV median first arrivals in 2022 (top right), bias-corrected mean of the GEV median first arrivals in 2022 (bottom left), and the same plot as before but with indications (in white rectangles over the corresponding pixels) for where the model estimates an occupancy probability of less than $1\%$ for the species (bottom right). The grey pixels in the top left display indicate where no species occurrences were recorded in 2022, or where no checklists and therefore no median durations were available in that pixel. Species | $\hat{\theta}^{\text{pref}}$ | $\hat{\theta}^{\text{act}}$ | Observed | Predicted | Debiased | Observed | Predicted | Debiased
---|---|---|---|---|---|---|---|---
Chimney Swift | 0.247 (0.236, 0.258) | -0.295 (-0.361,-0.204) | NA | 23/05 | 23/04 | 01/05 | 08/05 | 23/04
Great Crested Flycatcher | 0.126 (0.121, 0.131) | -0.735 (-0.820 ,-0.652) | NA | 20/05 | 25/04 | 04/05 | 11/05 | 26/04
Magnolia Warbler | 0.167 (0.163, 0.170) | -0.620 (-0.706, -0.543) | NA | 27/05 | 06/05 | 06/05 | 16/05 | 02/05
Purple Martin | 0.273 (0.297, 0.240) | -0.367 (-0.625, 0.033) | NA | 24/05 | 18/04 | 29/06 | 07/05 | 24/04
Table 1: Examples of predictions related to first arrival dates for two pixel-
years in year 2022. The first three columns report the species and the
posterior mean estimates of the sharing parameters, with $95\%$ credible
intervals in brackets. The next columns give the observed, predicted and
debiased first arrivals (in day/month format) at the two pixels indicated by
the two triangles in Figure 8); the first three columns correspond to the
pixel in the southeastern area, the last three columns to the pixel located in
the middle of the study region. An “NA” value in the Observed column indicates
that there were no observations within that pixel in 2022. The strikethroughs
over the dates indicate that our model estimates an occupancy probability of
the species of less than $1\%$ for that pixel-year.

##  6 Discussion and outlook

###  6.1 Statistical learning from participatory data

Participatory data collection and citizen science programs can provide huge
added value for advancing and sharing scientific knowledge, thanks in
particular to the large volume of data and their diversity in terms of
spatiotemporal coverage. However, sampling protocols from these sources are
usually not very strict and there may be different levels of observer
expertise, such that statistical techniques are required to correct for
heterogeneous sampling effort and observation biases.

In this work, we have developed a relatively complex Bayesian hierarchical
model where sampling effort and its influence on the observation bias of a
phenological event is quantified, and then corrected for during prediction.
Our results show that the necessary corrections are often substantial.
Moreover, by singling out the sampling effort, we were also able to estimate
spatial occupancy probabilities, which allowed us to precisely filter out
areas of space where it does not make sense to provide model-based predictions
of the phenological event since it might never occur.

###  6.2 Bayesian statistics and machine learning

We advocate the use of carefully constructed statistical models allowing for
inference and interpretation of latent components. Bayesian hierarchical
models (BHMs) have many advantages: they permit to identify various latent
processes (such as sampling effort) and their interactions, to keep track of
uncertainties and assess them through posterior sampling, to incorporate
expert knowledge, and to properly separate the process of how data are
observed from the latent ecological processes we seek to infer. The Bayesian
framework further offers the advantage of being relatively robust with noisy
data and overdispersion in count data through the inclusion of appropriately
designed random effects.

Regarding alternative approaches, various machine learning and deep learning
tools have become widely used in ecology in recent years and are extremely
powerful for predicting the data as they have been observed. Further
developments may be useful to allow such algorithms to predict the processes
of interest if these are not directly observed, to extrapolate beyond the
range of observed data (as we have done by setting the observation effort to a
very high level when performing prediction), and to avoid propagating or even
reinforcing biases from the input data when generating predictive outputs
(Dunson, 2018). An interesting goal of research in the field of widely used
machine learning and deep learning approaches would be to devise mechanisms
similar to the sharing of random effects in BHMs that can be implemented
within the architectures of such alternative models and algorithms.

###  6.3 Modelling sampling effort

The sampling effort as assessed in eBird has various dimensions, and we have
here taken into account the annual number of visits to an area (i.e., using
checklists) but also the median duration of those visits. Understanding how
the combination of the two aspects influences the number of detected species
individuals is a complex task (see also the work of Tang et al., 2021), and
would merit further exploration going beyond the additive structure of
$x_{\text{effort}}$ that we posed in (3). Other relevant attributes that we
did not consider here are partially reported for eBird checklists, such as the
surface of the area where observations took place, or the distance over which
the observers moved during the checklist event. Such information could be
further studied for a more complete characterisation of the sampling effort
and its impact on observation biases with respect to species distributions.

###  6.4 Detectability of species

Observation biases related to the detectability of species by the observer
remain a major challenge since they are difficult to identify with available
data, i.e., it is possible that no species occurrences were reported although
individuals of the species were present near the sampling location. This risk
of non detection can vary with land cover, species behaviour, observer
experience and other factors. In our model, detectability can be viewed as
being included in the niche process (which could more precisely be called the
_detectable niche_), but it seems hard to come up with a statistical approach
to further disentangle the niche and the detectability.

###  6.5 Validation of results

Validation of modelling results is inherently challenging with data obtained
through heterogeneous sampling effort and ridden with biases not well known
beforehand. Certainly, we can validate whether the model appropriately
reproduces data as they have been observed. However, the latent processes of
interest, as well as sampling effort and observation biases, can be entangled
in complex ways. For example, we performed such checks for observed first
arrival dates in Section 5. On the other hand, direct validation of inferences
on the latent processes of interest is often difficult, such that careful
construction of models and inference algorithms is paramount.

Sometimes, external data sources not suffering from these drawbacks can be
used, but often the available data have much smaller coverage of space, time
and different species. We point out that the fusion of data sources, such as
done in our study with BBS and eBird data to infer the niche of species, can
help improve the estimation of certain model components, but could also be
used for validation purposes at space-time locations covered by both datasets.

Nowadays, more and more data collection projects are initiated with the goal
to use advanced sensing technologies that do not require the presence of a
human observer and therefore provide exciting opportunities to collect large
amounts of data obeying strict sampling protocols. For instance, one could
track bird arrivals with radar data (e.g., Nussbaumer et al. (2021), or the
vogelwarte.ch project). Camera traps combined with artificial intelligence for
species identification can provide near-continuous temporal coverage, although
with rather limited spatial coverage.

###  6.6 Possible extensions of our approach

Here, we have focused on identifying spatiotemporal evolutions in ecological
processes for the past and present. In the context of species distribution
modelling, important research questions concern the potential impact of
climate change and its imbrication with other global change processes (e.g.,
land use). Our modelling approach could be further refined to establish the
relationships between ecological dynamics and climate variables in more
detail, and could then be used for projecting dates of first arrival under
future climate simulations, similar to Wijeyakulasuriya et al. (2023).

We have fixed the spatial pixel mesh using a mesh size similar to that
routinely used in modelling results published by the eBird project. Regarding
our work here, it would be interesting to further explore the influence of the
mesh size on the properties of the observed dates of first arrivals. For
example, a larger mesh size would lead to larger numbers of observations
within each pixel and could thus potentially reduce the observation bias or
facilitate its quantification. On the other hand, it would also decrease the
spatial resolution of predictive maps and the hamper the identification of
relationships of ecological processes with land cover. Therefore, further
research efforts could be oriented to investigate the role of the choice of
the mesh, and how to optimize this choice, or how to optimally combine results
from different mesh sizes.

Another interesting extension would consist of studying interactions among
different species in terms of their niche and phenological events, for example
to infer which groups of species show similar spring migration patterns.

###  6.7 Ecological extreme-value analysis

Finally, we believe that Extreme-Value Theory could have a much more important
role to play in ecological science. Extreme climatic and environmental events
are known to have major impacts on biodiversity. Therefore, a first important
domain of application is the analysis of the impacts of such extreme events on
ecological processes. A second domain of application, of which the present
work is part of, concerns the analysis of extreme phenological events using
Extreme-Value Theory. Applications of this type are still in their early
infancy in the existing literature, but we hope that our work shows the
utility of a more widespread application of extreme-value techniques.

## References

  * Banerjee et al. (2003) Banerjee, S., B. P. Carlin, and A. E. Gelfand (2003).  Hierarchical modeling and analysis for spatial data.  Chapman and Hall/CRC.
  * Bock and Root (1981) Bock, C. E. and T. L. Root (1981).  The christmas bird count and avian ecology.  Studies in Avian Biology 6, 17–23.
  * Bolin and Lindgren (2015) Bolin, D. and F. Lindgren (2015).  Excursion and contour uncertainty regions for latent Gaussian models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77(1), 85–106.
  * Bonney et al. (2016) Bonney, R., C. Cooper, and H. Ballard (2016, May).  The theory and practice of citizen science: Launching a new journal.  Citizen Science: Theory and Practice.
  * Bonney et al. (2016) Bonney, R., T. B. Phillips, H. L. Ballard, and J. W. Enck (2016).  Can citizen science enhance public understanding of science?  Public Understanding of Science 25, 2–16.
  * Borda-de Água et al. (2021) Borda-de Água, L., S. Alirezazadeh, M. Neves, S. P. Hubbell, P. A. V. Borges, P. Cardoso, F. Dionísio, and H. M. Pereira (2021).  Species Accumulation Curves and Extreme Value Theory, pp. 211–226.  Ecology, Biodiversity and Conservation. Cambridge University Press.
  * Browne (1996) Browne, J. (1996).  Charles darwin: A biography, vol. i, voyaging.  Journal of the History of Biology 29(2), 314–316.
  * Coles (2001) Coles, S. (2001).  An introduction to statistical modeling of extreme values.  Springer.
  * Colwell et al. (1994) Colwell, R. K., J. A. Coddington, and D. L. Hawksworth (1994).  Estimating terrestrial biodiversity through extrapolation.  Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 345(1311), 101–118.
  * Conklin et al. (2021) Conklin, J. R., S. Lisovski, and P. F. Battley (2021).  Advancement in long-distance bird migration through individual plasticity in departure.  Nature Communications 12(1), 4780.
  * Cornell Lab of Ornithology (2022) Cornell Lab of Ornithology (2022).  eBird Basic Dataset, Version: EBD relNov-2022.
  * Cotton (2003) Cotton, P. A. (2003).  Avian migration phenology and global climate change.  Proceedings of the National Academy of Sciences 100(21), 12219–12222.
  * Datta et al. (2016) Datta, A., S. Banerjee, A. O. Finley, and A. E. Gelfand (2016).  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association 111(514), 800–812.
  * Diggle et al. (2010) Diggle, P. J., R. Menezes, and T.-l. Su (2010).  Geostatistical inference under preferential sampling.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 59(2), 191–232.
  * Dunson (2018) Dunson, D. B. (2018).  Statistics in the big data era: Failures of the machine.  Statistics & Probability Letters 136, 4–9.
  * Fink et al. (2020) Fink, D., T. Auer, A. Johnston, V. Ruiz-Gutierrez, W. M. Hochachka, and S. Kelling (2020).  Modeling avian full annual cycle distribution and population trends with citizen science data.  Ecological Applications 30(3), e02056.
  * Fraisl et al. (2020) Fraisl, D., J. Campbell, L. See, U. Wehn, J. Wardlaw, M. Gold, I. Moorthy, R. Arias, J. Piera, J. L. Oliver, J. Masó, M. Penker, and S. Fritz (2020).  Mapping citizen science contributions to the un sustainable development goals.  Sustainability Science 15(6), 1735–1751.
  * Fraisl et al. (2022) Fraisl, D., G. Hager, B. Bedessem, M. Gold, P.-Y. Hsing, F. Danielsen, C. B. Hitchcock, J. M. Hulbert, J. Piera, H. Spiers, M. Thiel, and M. Haklay (2022).  Citizen science in environmental and ecological sciences.  Nature Reviews Methods Primers 2(1), 64.
  * Fuglstad et al. (2018) Fuglstad, G.-A., D. Simpson, F. Lindgren, and H. Rue (2018).  Constructing priors that penalize the complexity of Gaussian random fields.  Journal of the American Statistical Association 114(525), 445–452.
  * Gaines and Denny (1993) Gaines, S. D. and M. W. Denny (1993).  The largest, smallest, highest, lowest, longest, and shortest: extremes in ecology.  Ecology 74(6), 1677–1692.
  * Gelfand and Shirota (2019) Gelfand, A. E. and S. Shirota (2019).  Preferential sampling for presence/absence data and for fusion of presence/absence data with presence-only data.  Ecological Monographs 89(3), e01372.
  * Guinness (2018) Guinness, J. (2018).  Permutation and grouping methods for sharpening gaussian process approximations.  Technometrics 60(4), 415–429.  PMID: 31447491.
  * Haklay et al. (2021) Haklay, M., D. Fraisl, B. Greshake Tzovaras, S. Hecker, M. Gold, G. Hager, L. Ceccaroni, B. Kieslinger, U. Wehn, S. Woods, C. Nold, B. Balázs, M. Mazzonetto, S. Ruefenacht, L. A. Shanley, K. Wagenknecht, A. Motion, A. Sforzi, D. Riemenschneider, D. Dorler, F. Heigl, T. Schaefer, A. Lindner, M. Weißpflug, M. Mačiulienė, and K. Vohland (2021).  Contours of citizen science: a vignette study.  Royal Society Open Science 8(8), 202108.
  * Hsing et al. (2022) Hsing, P.-Y., R. A. Hill, G. C. Smith, S. Bradley, S. E. Green, V. T. Kent, S. S. Mason, J. Rees, M. J. Whittingham, J. Cokill, M. C. Scientists, and P. A. Stephens (2022).  Large-scale mammal monitoring: The potential of a citizen science camera-trapping project in the united kingdom.  Ecological Solutions and Evidence 3(4), e12180.
  * Isaac et al. (2014) Isaac, N. J. B., A. J. van Strien, T. A. August, M. P. de Zeeuw, and D. B. Roy (2014).  Statistics for citizen science: extracting signals of change from noisy ecological data.  Methods in Ecology and Evolution 5(10), 1052–1060.
  * Johnston et al. (2018) Johnston, A., D. Fink, W. M. Hochachka, and S. Kelling (2018).  Estimates of observer expertise improve species distributions from citizen science data.  Methods in Ecology and Evolution 9(1), 88–97.
  * Johnston et al. (2021) Johnston, A., W. M. Hochachka, M. E. Strimas-Mackey, V. Ruiz Gutierrez, O. J. Robinson, E. T. Miller, T. Auer, S. T. Kelling, and D. Fink (2021).  Analytical guidelines to increase the value of community science data: An example using ebird data to estimate species distributions.  Diversity and Distributions 27(7), 1265–1277.
  * Johnston et al. (2023) Johnston, A., E. Matechou, and E. B. Dennis (2023).  Outstanding challenges and future directions for biodiversity monitoring using citizen science data.  Methods in Ecology and Evolution 14(1), 103–116.
  * Jordan et al. (2012) Jordan, R. C., H. L. Ballard, and T. B. Phillips (2012).  Key issues and new approaches for evaluating citizen-science learning outcomes.  Frontiers in Ecology and the Environment 10(6), 307–309.
  * Katz et al. (2005) Katz, R. W., G. S. Brush, and M. B. Parlange (2005).  Statistics of extremes: modeling ecological disturbances.  Ecology 86(5), 1124–1134.
  * Katzfuss and Guinness (2021) Katzfuss, M. and J. Guinness (2021).  A General Framework for Vecchia Approximations of Gaussian Processes.  Statistical Science 36(1), 124 – 141.
  * Kelling et al. (2015) Kelling, S., A. Johnston, W. M. Hochachka, M. Iliff, D. Fink, J. Gerbracht, C. Lagoze, F. A. La Sorte, T. Moore, A. Wiggins, et al. (2015).  Can observation skills of citizen scientists be estimated using species accumulation curves?  PloS one 10(10), e0139600.
  * Koh et al. (2023) Koh, J., F. Pimont, J.-L. Dupuy, and T. Opitz (2023).  Spatiotemporal wildfire modelling through point processes with moderate and extreme marks.  The Annals of Applied Statistics 17(1), 560–582.
  * Kosmala et al. (2016) Kosmala, M., A. Wiggins, A. Swanson, and B. Simmons (2016).  Assessing data quality in citizen science.  Frontiers in Ecology and the Environment 14(10), 551–560.
  * Linden (2011) Linden, A. (2011).  Using first arrival dates to infer bird migration phenology.  Boreal Environment Research 16, 49–60.
  * Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Lindström (2011).  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.
  * Marshall et al. (2015) Marshall, P. J., C. J. Lintott, and L. N. Fletcher (2015).  Ideas for citizen science in astronomy.  Annual Review of Astronomy and Astrophysics 53(1), 247–278.
  * McKinley et al. (2015) McKinley, D. C., A. J. Miller-Rushing, H. L. Ballard, R. Bonney, H. Brown, D. M. Evans, R. A. French, J. K. Parrish, T. B. Phillips, S. F. Ryan, L. A. Shanley, J. L. Shirk, K. F. Stepenuck, J. F. Weltzin, A. Wiggins, O. D. Boyle, R. D. Briggs, S. F. Chapin III, D. A. Hewitt, P. W. Preuss, and M. A. Soukup (2015).  Investing in citizen science can improve natural resource management and environmental protection.  19.
  * Møller et al. (1998) Møller, J., A. R. Syversveen, and R. P. Waagepetersen (1998).  Log gaussian cox processes.  Scandinavian journal of statistics 25(3), 451–482.
  * Newman et al. (2012) Newman, G., A. Wiggins, A. Crall, E. Graham, S. Newman, and K. Crowston (2012).  The future of citizen science: emerging technologies and shifting paradigms.  Frontiers in Ecology and the Environment 10(6), 298–304.
  * Nussbaumer et al. (2021) Nussbaumer, R., S. Bauer, L. Benoit, G. Mariethoz, F. Liechti, and B. Schmid (2021).  Quantifying year-round nocturnal bird migration with a fluid dynamics model.  Journal of The Royal Society Interface 18(179), 20210194.
  * Pati et al. (2011) Pati, D., B. J. Reich, and D. B. Dunson (2011, 03).  Bayesian geostatistical modelling with informative sampling locations.  Biometrika 98(1), 35–48.
  * Pocock et al. (2018) Pocock, M. J., M. Chandler, R. Bonney, I. Thornhill, A. Albin, T. August, S. Bachman, P. M. Brown, D. G. F. Cunha, A. Grez, C. Jackson, M. Peters, N. R. Rabarijaon, H. E. Roy, T. Zaviezo, and F. Danielsen (2018).  Chapter six - a vision for global biodiversity monitoring with citizen science.  In D. A. Bohan, A. J. Dumbrell, G. Woodward, and M. Jackson (Eds.), Next Generation Biomonitoring: Part 2, Volume 59 of Advances in Ecological Research, pp. 169–223. Academic Press.
  * Roberts and Rosenthal (2001) Roberts, G. O. and J. S. Rosenthal (2001).  Optimal scaling for various Metropolis-Hastings algorithms.  Statistical Science 16(4), 351 – 367.
  * Rue et al. (2017) Rue, H., A. Riebler, S. H. Sørbye, J. B. Illian, D. P. Simpson, and F. K. Lindgren (2017).  Bayesian computing with INLA: a review.  Annual Review of Statistics and Its Application 4, 395–421.
  * Silvertown (2009) Silvertown, J. (2009).  A new dawn for citizen science.  Trends in Ecology & Evolution 24(9), 467–471.
  * Somveille et al. (2015) Somveille, M., A. S. Rodrigues, and A. Manica (2015).  Why do birds migrate? A macroecological perspective.  Global Ecology and Biogeography 24(6), 664–674.
  * Tang et al. (2021) Tang, B., J. S. Clark, and A. E. Gelfand (2021).  Modeling spatially biased citizen science effort through the ebird database.  Environmental and Ecological Statistics 28(3), 609–630.
  * Theobald et al. (2015) Theobald, E., A. Ettinger, H. Burgess, L. DeBey, N. Schmidt, H. Froehlich, C. Wagner, J. HilleRisLambers, J. Tewksbury, M. Harsch, and J. Parrish (2015).  Global change and local solutions: Tapping the unrealized potential of citizen science for biodiversity research.  Biological Conservation 181, 236–244.
  * Thibaud et al. (2016) Thibaud, E., J. Aalto, D. S. Cooley, A. C. Davison, and J. Heikkinen (2016).  Bayesian Inference For The Brown-Resnick Process, With An Application To Extreme Low Temperatures.  Annals of Applied Statistics 10(4), 2303–2324.
  * UN (2015) UN (2015).  A/res/70/1 un general assembly transforming our world: the 2030 agenda for sustainable development.  Seventieth session of the General Assembly on 25 Sept 2015.
  * van de Schoot et al. (2021) van de Schoot, R., S. Depaoli, R. King, B. Kramer, K. Märtens, M. G. Tadesse, M. Vannucci, A. Gelman, D. Veen, J. Willemsen, et al. (2021).  Bayesian statistics and modelling.  Nature Reviews Methods Primers 1(1), 1.
  * Vecchia (1988) Vecchia, A. V. (1988).  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society. Series B (Methodological) 50(2), 297–312.
  * Wijeyakulasuriya et al. (2023) Wijeyakulasuriya, D. A., E. M. Hanks, and B. A. Shaby (2023).  Modeling first arrival of migratory birds using a hierarchical max-infinitely divisible process.
  * Wijeyakulasuriya et al. (2019) Wijeyakulasuriya, D. A., E. M. Hanks, B. A. Shaby, and P. C. Cross (2019).  Extreme value-based methods for modeling elk yearly movements.  Journal of Agricultural, Biological and Environmental Statistics 24(1), 73–91.
  * Wynn (2017) Wynn, J. (2017).  Citizen science in the digital age : rhetoric, science, and public engagement.  Rhetoric, culture, and social critique. Tuscaloosa: The University of Alabama Press Tuscaloosa.
  * Yadav et al. (2023) Yadav, R., R. Huser, T. Opitz, and L. Lombardo (2023, 09).  Joint modelling of landslide counts and sizes using spatial marked point processes with sub-asymptotic mark distributions.  Journal of the Royal Statistical Society Series C: Applied Statistics, qlad077.
  * Youngflesh et al. (2021) Youngflesh, C., J. Socolar, B. R. Amaral, A. Arab, R. P. Guralnick, A. H. Hurlbert, R. LaFrance, S. J. Mayor, D. A. Miller, and M. W. Tingley (2021).  Migratory strategy drives species-level variation in bird sensitivity to vegetation green-up.  Nature Ecology & Evolution 5(7), 987–994.

blabla

##  7 Appendix

###  7.1 Land-cover plots

Figure 10 shows plots of the four land-cover covariates defined as proportions
of certain land-cover categories in the pixels of the spatial region.

|
---|---
|
Figure 10: Land-cover proportions in the study area (NLCD, 2021) according to
four types described in Section 2.

###  7.2 Plots from the simulation study

Figure 11 shows plots assessing the identifiability of the observed first
arrivals in our simulation study.

Figure 11:  Simulation study: Truth vs. posterior means for the median first
arrival dates in 2003 (left) and 2020 (right), shown with spatial maps (for 3
different seeds in the simulation study), and with boxplots of the differences
in days between true and predicted arrival dates, based on 100 simulations.
The boxplot shows results for 150 randomly chosen pixels, marked with black
boxes in the spatial plots.

###  7.3 Further plots of MCMC results

Figure 12 shows trace plots from the MCMC run of the model on the species
Magnolia Warbler.

Figure 13 shows spatial plots of the four latent spatial fields from the model
estimated for the species Chimney Swift.

Figure 12:  MCMC trace plots for some fixed effect and sharing parameters (top
left to bottom right): $\beta_{0}^{\text{ckl}}$, $\beta_{0}^{\text{spc}}$,
$\beta_{0}^{\text{ckl}}$, $\beta_{0}^{\text{spc}}$ from our model for the
species Magnolia Warbler. The red dashed line indicates the burn-in period.
Figure 13: Posterior means plots of the four spatial random effects in (2) for
the model fitted to the species Chimney Swift.