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

Jonathan Koh, Thomas Opitz

► To cite this version:

Jonathan Koh, Thomas Opitz. Extreme-value modelling of migratory bird arrival dates: insights from citizen science data. Journal of the Royal Statistical Society: Series A Statistics in Society, 2025, 188 (3), pp.674-699. 10.1093/jrsssa/qnae108. hal-05198967

HAL Id: hal-05198967

https://hal.inrae.fr/hal-05198967v1

Submitted on 20 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient statistical inference methods for assessing changes in species' populations using citizen science data

Emily B. Dennis<fnref n="1" />,<fnref n="2" /> ID, Alex Diana<fnref n="3" />, Eleni Matechou<fnref n="2" /> and Byron J.T. Morgan<fnref n="2" />

*Address for correspondence:* Byron J.T. Morgan, School of Mathematics, Statistics and Actuarial Science, University of Kent, Canterbury, Kent CT2 7NF, UK. Email: b.j.t.morgan@kent.ac.uk

Read before The Royal Statistical Society at the Discussion Meeting on the 'Analysis of citizen science data' held at the Society's 2024 annual conference in Brighton on Tuesday, 3 September 2024, the President, Dr Andrew Garrett, in the Chair.

## Abstract

The global decline of biodiversity, driven by habitat degradation and climate breakdown, is a significant concern. Accurate measures of change are crucial to provide reliable evidence of species' population changes. Meanwhile citizen science data have witnessed a remarkable expansion in both quantity and sources and serve as the foundation for assessing species' status. The growing data reservoir presents opportunities for novel and improved inference but often comes with computational costs: computational efficiency is paramount, especially as regular analysis updates are necessary. Building upon recent research, we present illustrations of computationally efficient methods for fitting new models, applied to three major citizen science data sets for butterflies. We extend a method for modelling abundance changes of seasonal organisms, firstly to accommodate multiple years of count data efficiently, and secondly for application to counts from a snapshot mass-participation survey. We also present a variational inference approach for fitting occupancy models efficiently to opportunistic citizen science data. The continuous growth of citizen science data offers unprecedented opportunities to enhance our understanding of how species respond to anthropogenic pressures. Efficient techniques in fitting new models are vital for accurately assessing species' status, supporting policy-making, setting measurable targets, and enabling effective conservation efforts.

**Keywords:** biodiversity change, citizen science, concentrated likelihood, generalized abundance index, occupancy models, variational Bayes

## 1 Introduction

An existential crisis of our time is the alarming decrease of biodiversity, due to anthropogenic factors such as climate breakdown and loss of habitat. Producing robust measures of change is vital for evaluating species' status, understanding rates of change, and monitoring responses to pressures, as well as progress of conservation actions, such as towards biodiversity targets (Butchart et al., 2010). Both data and appropriate statistical models are critical for measuring biodiversity change. Citizen (or community) science (CS) data, where information is gathered by voluntary participants, are increasingly used for this purpose (Chandler et al., 2017; Silvertown, 2009).

CS data from systematic, designed surveys and monitoring schemes, collected by skilled and committed volunteers, have long been used to produce estimates of species' status, whereas

<footnote n="1">¹Butterfly Conservation, Manor Yard, East Lulworth, Dorset BH20 5QP, UK</footnote>

<footnote n="2">²School of Mathematics, Statistics and Actuarial Science, University of Kent, Canterbury, Kent CT2 7NF, UK</footnote>

<footnote n="3">³School of Mathematics, Statistics and Actuarial Science, University of Essex, Wivenhoe Park, Colchester, Essex CO4 3SQ, UK</footnote>
observations of species from less structured, opportunistic or mass-participation sampling, attracting contributions from wider society, are increasingly gathered (Pocock et al., 2017), broadening the scope to measure changes in populations from extensive geographic areas and for a variety of taxa. The many sources of CS data provide vast opportunities for biodiversity monitoring, but require suitable analytical approaches, for example to deal with sources of bias (Isaac & Pocock, 2015). Johnston et al. (2023) outlined a diversity of challenges for biodiversity monitoring using CS data, relating to dealing with observer behaviour, data structures, statistical models, and communication. This paper focuses upon one such key challenge, which is in the computational cost of analysing CS data, particularly with increasing data and complex models, leading to the 'necessity to identify and develop suitable modifications to improve computational efficiency and scalability, adapting traditional (and developing new) methods to big data' (McCrea et al., 2023).

In this paper, we present illustrations of computationally efficient methods for analysing CS data. Our work is motivated by applications to data for British butterflies, to model changes in both abundance and distribution, but the methods and overall need for efficiency apply to CS data for a range of taxa and locations. The global decline of insects has been widely reported, particularly in western Europe and North America (Wagner et al., 2021), yet there is an ongoing need for robust data and rigorous analysis methods (Didham et al., 2020; Thomas et al., 2019). Many taxa and geographic regions are lacking in sufficient data to appropriately assess trends, but there is a wealth of data gathered on butterflies in the UK. Butterflies are the most comprehensively monitored invertebrate taxon and their population status provides a valuable indicator for changes in biodiversity as they respond rapidly to environmental change.

We demonstrate efficient analysis methods for three sources of CS data for UK butterflies: (i) The UK Butterfly Monitoring Scheme (UKBMS) began in 1976 and is one of the longest-running insect monitoring schemes in the world. Counts of butterflies are made each year by recorders walking transects for six months from 1 April, according to a strict protocol (Pollard & Yates, 1993). (ii) The annual Big Butterfly Count (BBC) launched in 2010 and is a UK-wide mass-participation CS project in which members of the public record how many individual butterflies are seen for 15 min periods during 23-24 days in late July and early August. It is the largest CS project of its kind in the world, with almost 95,000 participants in 2023. (iii) The Butterflies for the New Millenium (BNM) data base collates opportunistic records of where butterflies are seen, and consists of over 14 million records (Fox et al., 2023). These data are unstandardized, and more than 80 similar recording schemes exist for various taxa in the UK (BRC, 2022; Pocock et al., 2015) and are typically used for describing changes in species' distributions.

Analyses of UKBMS and BNM data feed in to regular reporting such as the 'State of the UK's Butterflies' (Fox et al., 2023), Red List assessments of extinction risk (Fox et al., 2022), government biodiversity indicators (JNCC, 2022), and multitaxa outputs such as the State of Nature (Burns et al., 2023), all of which contribute to providing robust evidence for conservation, policy development, and the wider state of biodiversity.

During a time of biodiversity change, frequent analysis updates are essential for monitoring species' populations, and for understanding responses to both pressures and conservation action. Efficient computational methods are therefore vital, especially in CS surveys, which typically involve data sets that are large and continuously growing in size.

We are motivated by the need to measure population changes, by modelling abundance and distribution using the data sets mentioned above. In this paper, we summarize our recent research, describe new methods and applications, and suggest avenues for future development. We start with models for abundance, where we demonstrate how in certain circumstances it is possible to greatly reduce the dimensionality of the effective model parameter space, considerably reducing computation time. We then consider opportunistic data, where we show how variational inference (VI) can provide an efficient alternative to Markov chain Monte Carlo (MCMC) for Bayesian inference when fitting occupancy models.

## 2 Generalized abundance index

Models for abundance data, such as from the UKBMS, need to account for two key features: seasonality and missing data. Seasonality results in counts that vary throughout the season according to the emergence of one or more generations of adult insects over time in any year. Missing data
arise in two ways. Firstly, some visits are missed by volunteers. UKBMS transects are sampled by committed and skilled volunteers, walking transects weekly across six months of the year, but inevitably some weeks are missed: Dennis et al. (2013) estimated that this occurs for roughly 8 of 26 weeks of the transect season. Secondly, there is turnover in transects sampled each year; for example transects may cease to be monitored, and new transects continue to be introduced to the scheme.

A variety of methods have produced analyses of these data. Dennis et al. (2013) presented a method based on using generalized additive models (GAMs) to model seasonality, followed by a generalized linear model (GLM) analysis to account for the annual variation in the transects sampled. This approach produced annual indices of abundance for each species, which can be used to estimate time trends in abundance. Alternatively, Matechou et al. (2014) proposed a stopover model approach, which enables estimation of within-season survival probabilities. Dennis et al. (2016) then proposed the generalized abundance index (GAI) approach, which provides a framework for both of these types of models, as well as a third alternative, which describes seasonality parametrically using an appropriate mixture of distributions, such as Normal (see Dennis et al., 2022, for examples).

In outline, for a given species, in any year we suppose that counts are obtained at S sites, each visited on at most V occasions. Each count, $y_{s,v}$, for site s and visit v is regarded as the realization of a random variable, such as Poisson, with expectation $\lambda_{s,v} = N_s a_{s,v}$, where the likelihood then takes the form

$$L(N, \theta; y) \propto \prod_{s=1}^{S} \prod_{v=1}^{V} \exp(-N_s a_{s,v}) (N_s a_{s,v})^{y_{s,v}}. \quad (1)$$

Here the $N_s$ are parameters describing site abundance, and $a_{s,v}$ denotes a function determined by parameters, $\theta$, describing seasonal variation. Here and later the product over visits only includes terms corresponding to when visits are made. Background is provided by Dennis et al. (2016), who derive a concentrated-likelihood approach that substantially reduces the number of parameters to estimate via maximum likelihood. Briefly, using maximum likelihood, the site parameters can be estimated by

$$N_s = \frac{y_{s,v}}{a_{s,v}}, \quad (2)$$

where we use the dot notation to indicate summation corresponding to visited sites. The total observed count for each site is therefore rescaled to account for incomplete sampling within the season. Substitution of the estimates of equation (2) into equation (1) results in a Poisson likelihood which can be maximized with respect to only the parameters $\theta$. Counts can be expected to be overspread relative to the Poisson and/or contain additional zeros, for example due to small counts at the ends of the season. Alternative discrete distributions such as negative binomial and zero-inflated Poisson, respectively, may then be appropriate, as described in Dennis et al. (2016).

# 3 Two adaptations of the GAI

We now provide two new adaptations of the GAI. Firstly, the model is extended to include appropriate site and year effects using an annual model integrated within the GAI, combined with the use of concentrated likelihood. Secondly, we adapt the GAI to improve the analysis of BBC data, a ‘snapshot’, mass-participation CS scheme.

## 3.1 An extended GAI with site and time effects

The basic GAI, described above, is a static model, where data for each year are analysed separately. However, to deal with the turnover of sites surveyed in the UKBMS each year, the Poisson GLM stage of Dennis et al. (2013), based on the annual model of ter Braak et al. (1994), is typically used. Here the abundance estimates ($\hat{N}_{s,r}$), for sites s and now year r, are the dependent values for an additive model with year and site effects and the estimated year effects are then used to form abundance indices. Use of log-linear regression in this way is widely used in ecology, for example using TRIM (Trends and Indices for Monitoring Data, Bogaart et al., 2020; van Strien et al., 2004). This
two-stage GAI approach is now a valuable tool for analysing seasonal count data in the UK (Fox et al., 2021, 2023; JNCC, 2022; UKBMS, 2023) and beyond (Schmucki et al., 2016; Van Swaay et al., 2020). In practice, the proportion of the species' flight period surveyed for a given site and year is typically included as a weighting in the GLM stage, such that better sampled sites have a higher contribution to the estimated species' abundance indices (Brereton et al., 2018).

A drawback of the two-stage GAI approach is the need to bootstrap to account for variance propagation between the two model stages. Bootstrapping is time consuming for large data sets, and can also be problematic when resampling small data sets, for example for rare species. We now describe an alternative approach which effectively incorporates the annual model within the GAI, extending the model to consider all years at once, and thus solving the issue of variance propagation.

In brief, we expand the Poisson likelihood of equation (1) to counts $y_{s,v,r}$ where $r$ denotes one of $Y$ successive years, with a corresponding expansion to $a_{s,v,r}$. We then incorporate the expression for the annual model $N_{s,r} = e^{a_s + \beta_r}$ - see ter Braak et al. (1994) - which results in the following expression for the log-likelihood, ignoring an additive constant.

$$ \ell(\alpha, \beta, \theta; y) = \sum_{r} \sum_{s} \sum_{v} \{-e^{(a_s + \beta_r)} a_{s,v,r} + y_{s,v,r}(a_s + \beta_r) + y_{s,v,r}\log(a_{s,v,r})\} \quad (3) $$

Here {$a_s$} and {$\beta_r$} are respectively site and year effects to be estimated. As above, we use a concentrated likelihood approach to form maximum-likelihood parameter estimates efficiently, by concentrating out the parameters $\alpha$, resulting in the log-likelihood

$$ \ell(\beta, \theta; y) = \sum_{r} \sum_{s} \sum_{v} \left[ -\frac{y_{s,...} e^{\beta_r} a_{s,v,r}}{\sum_{j} e^{\beta_j} a_{s,j}} + y_{s,v,r} \{\beta_r - \log(\sum_{j} e^{\beta_j} a_{s,j})\} \right] \quad (4) $$

which we can maximize efficiently with respect to the parameters $\theta$ and $\beta$. Full details are given in Section S1 of the online supplementary material.

The aim here is to produce a new model which can be fitted efficiently and produce estimates of uncertainty. From the above modelling, estimates of error result from inverting the estimated Hessian at the maximum-likelihood estimates, which is far more efficient than bootstrapping. Profile confidence intervals are also easily obtained.

We illustrate this extended GAI approach with application to UKBMS data spanning 1976–2022 for two species: Chalk Hill Blue *Polyommatus coridon* and Gatekeeper *Pyronia tithonus*. Chalk Hill Blue is a species confined to calcareous grassland in southern England, with UKBMS counts of around 900,000 individuals from ~ 460 sites, whereas Gatekeeper is a widespread butterfly species across southern Britain, with UKBMS counts from more than 4,300 sites, counting over 3.3 million individuals.

Options for describing seasonal variation (parameters associated with $\theta$) are flexible, as in the original GAI, but here we use a Normal distribution, since both species are univoltine (one generation per year), and thus parameters {$\mu_r$} and {$\sigma_r$} are estimated for each year, which describe the mean flight dates and length of the flight period, respectively (see Figure S1, online supplementary material where we see the tendency for earlier emergence over time). Flight periods are thus assumed to be fixed over sites within each year.

Figure 1 shows the estimates of {$\beta_r$} from the extended GAI, with 95% confidence intervals (CI) produced from the estimated Hessian at the maximum-likelihood estimates. Error estimates were scaled to account for overdispersion in the Poisson GAI formulation.

In order to estimate a linear trend over time, the parameters $\beta$ may also be expressed by a linear form, where we set $\beta_r = \gamma r, \forall r$, in equation (4) and maximize the resulting log likelihood with respect to the parameters $\theta$ and $\gamma$. The intercept parameter for a linear regression on year cannot be estimated due to confounding. This is also the case for the annual model (McCullagh & Nelder, 1989, p. 63, van Strien et al., 2004), where we therefore set $\beta_1 = 0$. Trend lines based on the extended GAI with the linear form are shown in Figure 1(i), and due to the lack of intercept, for ease of comparison we equate the ordinate at the middle year to that from a posthoc linear regression (also shown).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Results from applying the extended (generalized abundance index) GAI to two species: Chalk Hill Blue a) and Gatekeeper b). Plots (i) show the estimates of {$\beta_r$} (black) with 95% confidence intervals (CI). Linear trend estimates for {$\beta_r$} from the extended GAI are shown in green (with 95% CI). The trend line from a posthoc linear regression through the estimates of {$\beta_r$} is also shown (purple dashed line). Plots (ii) compare abundance indices from the extended GAI (black) and two-stage GAI (blue), with 95% CI produced from the estimated Hessian at the maximum-likelihood estimates and from nonparametric bootstrapping respectively. Abundance indices correspond to {$\beta_r$} converted to the log$_{10}$ scale with a mean of 2, as is standard practice for UK Butterfly Monitoring Scheme indices.</figcaption></figure>

Abundance indices produced from the extended GAI are virtually indistinguishable from indices produced from the two-stage GAI (Figure 1(ii)). The two-stage GAI was based on the implementation used for annual analyses of UKBMS data, where a GAM is used for the first stage. Despite the similarity in the indices, the 95% CI produced from bootstrapping the two-stage GAM are much wider than those estimated from the Hessian for the extended GAI. This difference is due to differences between the models, rather than due to differences in the methods of error estimation. For example the extended GAI does not make a distributional assumption (Poisson) about the site abundance estimates, N, as is the case in the GLM of two-stage GAI. The standard errors from the extended GAI may also be underestimated if the data are overdispersed relative to the Poisson distribution assumed, which could be explored in future work using alternative discrete distributions.

The implementation of the two-stage GAI includes a weighting in the GLM, accounting for the proportion of the flight period sampled, whereas the extended GAI does not. Interestingly, the close resemblance of the two indices in Figure 1(ii) implies that the new extended GAI effectively has this weighting built-in, by modelling the counts directly so that better-sampled sites have a greater contribution to the likelihood.

The extended model took 42 min to run for the Chalk Hill Blue, whereas the nonparametric bootstrap (with 1,000 replicates) for the two-stage GAI took 98 min. For the more widespread Gatekeeper, the extended GAI took 1.3 hr, whereas the two-stage GAI bootstrap took approximately 7 hr, but based on only 200 replicates. Analyses were parallelized across four CPUs on a computer equipped with Intel Core i7-8,700@3.2 GHz with 32 GB of RAM. Clearly, there are
substantial time savings to be had for the 59 UK butterfly species, and potentially far more for more numerous taxa. We have seen that the advantages of the extended GAI are not confined to efficiency, but may also involve greater parameter precision provided the model and distribu- tional assumptions are appropriate. In the Discussion, we outline the potential contributions of future work on the extended GAI.

## 3.2 Analysis of big butterfly count data

The Big Butterfly Count (BBC; https://bigbutterflycount.org) is an annual survey of common butterfly species which encourages wide participation particularly from members of the general public. The BBC has a high media profile and attracts a large number of participants, many of whom may have limited or no prior experience of biodiversity monitoring. The sampling protocol is minimal: participants simply count numbers of individuals seen of widespread butterfly (and day-flying moth) species for 15 min from any location.

The BBC generates a large amount of data for the UK's widespread butterfly species - more than 11 million butterflies have been counted since 2011- but to date few analyses have been undertaken and annual reporting of BBC results is based on simple comparisons with respect to the previous year only. Using data for 2011–2014, Dennis, Morgan, Brereton, et al. (2017) demonstrated that estimates of change in abundance from BBC were comparable to those estimated from standardized monitoring (UKBMS data), but that the short snapshot sampling period of three weeks results in bias caused by the inter-annual variation in species' flight periods.

This is demonstrated in Figure 2a which shows the estimated flight period for Marbled White Melanargia galathea for four example years. This species is a univoltine (single-generation), summer-flying species, which is therefore likely to be particularly susceptible to phenological bias (in the estimation of timing) with respect to the BBC sampling period. Flight periods were estimated from UKBMS data using the GAI with a spline formulation for {$a_v$}, which is fixed across sites within a given year. The timing of peak emergence varies year-to-year, and thus the proportion of the flight period sampled by the BBC varies annually (see Figure 2b). For example in 2011 the BBC only captures the tail end of the flight period, whereas in 2012 more than half of the flight period is captured.

Here we describe a modification of the GAI for producing abundance indices and trends from BBC data. The GAI was developed for standardized monitoring data, typically collected along transects, where sites are clearly defined and revisited many times within and across years. Mass-participation CS data such as from the BBC does not have such a structure; it consists of many locations, which often only have one count undertaken. Hence we define a BBC ‘site’ to be a 1km × 1km square and pool BBC data to this spatial scale. BBC sites may then have counts across multiple days within the BBC sampling period, and potentially multiple entries (15 min counts) on a given day.

To accommodate these multiple entries, the GAI likelihood for a given year now takes the form

$$L \propto \prod_{m=1}^{M} \prod_{v=1}^{V} \exp(-N_m a_{m,v} \kappa_{m,v}) (N_m a_{m,v})^{y_{m,v}}, \quad (5)$$

where $\kappa_{m,v}$ is the number of entries for site $m$ and visit $v$, and $y_{m,v} = \sum_{\kappa} y_{m,v,\kappa}$ is the sum of the counts over those entries.

The number of sites in a mass-participation CS dataset such as BBC is high, and we adopt a concentrated-likelihood approach to reduce the functional size of the parameter space, as in the previous section, where equation (2) becomes

$$N_m = \frac{y_{m,..}}{\sum_{v=1}^{V} a_{m,v} \kappa_{m,v}} \quad (6)$$

The BBC sampling period represents just a snapshot of most species' flight periods, therefore accurate flight period estimation from applying the GAI to BBC data is not possible for all species. Instead, for each species and year, we produce estimates of {$N_m$} for each BBC site
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Demonstration of the phenology adjustment approach for the Marbled White butterfly: a) flight period curves estimated from UK Butterfly Monitoring Scheme (UKBMS) data for 4 yr. The blue shaded areas represent the Big Butterfly Count (BBC) sampling period each year. b) the proportion of the Marbled White flight period covered by the BBC sampling period each year. c) relative abundance indices produced from the generalized abundance index applied to UKBMS data (black), from BBC data without phenology adjustment (i, blue squares), and from BBC data with phenology adjustment (ii, green triangles). Indices are on the log₁₀ scale with a mean value of 2 (indicated by the horizontal dashed lines).</figcaption></figure>

(1km × 1km square) from equation (6) using daily flight period estimates of {$a_v$} from the GAI applied to UKBMS data, where flight periods are assumed to be the same across sites. The average of {$N_m$} provides an overall measure of BBC abundance per year and species (Dennis et al., 2016).

Figure 2c demonstrates the benefit of the phenology adjustment approach for the Marbled White butterfly. Fluctuations in the BBC abundance index produced without phenology adjustment (Figure 2c(i)) largely reflect year-to-year variation in the proportion of the flight period captured by the BBC each year, whereas adjusting for phenology produces an index that shows a pattern more similar to the UKBMS abundance index (Figure 2c(ii))). Applying the phenology adjustment approach is less influential for a multivoltine species such as the Comma *Polygonia c-album* (Figure S2, online supplementary material), for which the BBC sampling period covers a smaller, and less variable, proportion of the overall flight period each year.

Dennis et al. (2024) describe the new phenology adjustment approach for snapshot citizen science data in full, including the use of bootstrapping to estimate uncertainty. The approach is applied to BBC data for 17 species and explored further via simulation. The method enables data from snapshot CS schemes such as BBC to contribute to monitoring biodiversity. BBC results
receive high-profile media coverage and are already beginning to reflect the advantages of the new analyses outlined here.

## 4 Efficient occupancy model fitting

Occupancy models (Altwegg & Nichols, 2019; MacKenzie et al., 2018) are widely used for inferring species distributions from presence/absence data at multiple sites, across a single or multiple seasons. They have been employed to extract meaningful species distribution trends from opportunistic citizen science data, addressing challenges associated with nonsystematic sampling and variable observation effort (Isaac et al., 2014; Kéry et al., 2010). However, as the size of the corresponding presence/absence data increases, occupancy models can be computationally demanding, especially in a Bayesian framework.

We start by defining the standard occupancy model that we fit in this section. We assume that there are n sampling units, where each sampling unit corresponds to a site in a specific year. In each sampling unit, we have a set of observations $y_i$, equal to 1/0 if the species was detected/not detected. The sampling unit to which observation i belongs is denoted by $k_i$. The hierarchical model representation is

$$ 
\left\{
\begin{array}{l}
\operatorname{logit}(\psi_j) = X_j^{\psi} \beta_{\psi} \\
z_j \sim \operatorname{Be}(\psi_j) \\
\operatorname{logit}(p_i) = X_i^{p_i} \beta_p \\
y_i \sim \begin{cases} \operatorname{Be}(p_i) & \text{if } z_{k_i} = 1 \\ 0 & \text{if } z_{k_i} = 0 \end{cases}
\end{array}
\right.
\qquad (7)
$$

where $z_j$ is the occupancy state of sampling unit j; $\beta_{\psi}$ are the covariate coefficients of the occupancy probability $\psi$; $\beta_p$ are the covariate coefficients of the detection probability $p$; $X_j^{\psi}$ and $X_i^{p_i}$ are the available covariates for sampling unit j and observation i, respectively.

Similar versions of this model have been considered in Diana et al. (2023) and in Doser et al. (2023), who employed Gaussian processes to account for spatio-temporal autocorrelation in the occupancy probability. Bayesian inference for these occupancy models can be easily performed using Markov chain Monte Carlo (MCMC) (see Diana et al., 2023; Doser et al., 2023, who employed a Pólya-Gamma sampling scheme for logistic regression models Polson et al., 2013). In this MCMC framework, the z terms from equation (7), indicating species presence/absence at each site, are treated as latent variables and hence inferred and updated, typically at each MCMC iteration. Inferring the latent variables z allows us to easily write the complete data likelihood—see for example King (2014) and Newman et al. (2023)—for the observations in y, as shown in equation (7). However, when the number of sites is large, this leads to a computationally intensive MCMC, even when efficient model-fitting approaches, such as the Pólya-Gamma scheme, are employed for updating the model parameters. In addition when the complete data likelihood is used there may be correlations between estimators which can also slow down MCMC—see Newman et al. (2023) and Borowska and King (2022), where a subset of latent states are not treated as auxiliary variables but are integrated out numerically to reduce the correlation between latent variables. Therefore, MCMC-based inference can be prohibitively slow, which limits its application to large data sets such as CS data. For example, for the data sets considered in Diana et al. (2023), obtaining acceptable effective sample sizes from the posterior distributions of all parameters requires running times of around 19 hr (on an Intel Core i7-10610U@1.8 GHz). An obvious alternative is to use classical inference to fit occupancy models since it can be much faster (see for example the approach of Dennis, Morgan, Freeman, et al., 2017). In this case, the likelihood function is written by marginalizing over the z variables and hence the observed data likelihood is used for inference. Expressions for the complete and observed data likelihood for occupancy models are given in Section S2 of the online supplementary material. However, quantifying uncertainty around functions of parameters can sometimes be computationally intensive, relying on bootstrap methods when closed form expressions of variances are not available. Bayesian inference also offers the potential to more readily account for spatio-temporal autocorrelation in the occupancy probabilities.
Variational inference (VI) has been proposed as an alternative tool to overcome the computational issues of MCMC (Jordan et al., 1999). VI is traditionally faster than MCMC-based approaches because it transforms the problem from sampling (from a posterior distribution) to optimization. Therefore, VI combines the speed of classical inference, with the interpretability of Bayesian inference. However, while MCMC-based inference always recovers the true posterior (given enough MCMC iterations), in VI the true posterior distribution $p(\theta|y)$ is approximated using an appropriate flexible family of distributions, which is called the variational family. We denote the variational distributions by $q_\lambda(\theta)$, where $\lambda$ is the set of variational parameters, the observed data likelihood by $p(y|\theta)$ and the prior distribution by $p(\theta)$. The parameters $\lambda$ corresponding to the optimal variational distribution can be found by minimizing the Kullback-Leibler (KL) divergence between the true posterior distribution $p(\theta|y)$ and the variational distribution $q_\lambda(\theta)$. It can be proved that this is equivalent to finding the $q_\lambda(\theta)$ that minimizes the quantity

$$E_{\theta \sim q_{\lambda}(\theta)}[\log p(y, \theta) - \log q_{\lambda}(\theta)] \quad (8)$$

which is known as the Evidence Lower Bound (ELBO), and forms the basis of VI inference. We note that $\log p(y, \theta) = \log \{p(y|\theta)p(\theta)\}$ and more details on VI can be found in Blei et al. (2017).

In VI, it is common to assume that parameters are a-posteriori independent, which is equivalent to assuming a variational family of the form $q(\theta) = \prod_i q_i(\theta_i)$ (the mean-field assumption) since this assumption considerably simplifies the inference. However, if the assumption is not valid, then posterior variance is underestimated (Wang & Titterington, 2005). In the case of occupancy models, occupancy and detection probability are independent conditionally on $z$, but not independent of $z$. Hence, assuming (according to the mean field assumption) that they are independent of $z$ implies that $\psi$ and $p$ are independent of each other a-posteriori, which clearly does not hold. If the dependence structure of the model is ignored, then it leads to underestimation of the posterior variance of the occupancy and detection probability parameters, $(\beta^{\psi}, \beta^p)$ (Clark et al., 2016).

However, if the observed data likelihood is used for inference, instead of the complete data likelihood, with the latter being common practice in a Bayesian framework, then we do not need to assume that $(\beta^{\psi}, \beta^p)$ are independent of $z$, or of each other. In the observed data-likelihood case, we assume, as is standard in VI (Titsias & Lázaro-Gredilla, 2014), that $q_\lambda(\theta)$ is a multivariate normal distribution, that is, if $\theta = (\beta^{\psi}, \beta^p)$ are the model parameters, we assume $\theta \sim q_\lambda(\theta) = N(\mu, \Sigma)$, where the variational parameters are $\lambda = (\mu, C)$, with $\mu$ the mean of the variational distribution and $C$ the Cholesky factor of the covariance matrix $\Sigma$.

We perform inference using stochastic gradient descent. Computing the gradient of equation (8) with respect to $\lambda$ is complicated by the fact that $\lambda$ itself appears in the expectation. To overcome this issue, we use the reparameterization trick (Rezende et al., 2014), which is a cornerstone of variational inference as it allows us to easily obtain these types of gradients. More details on the inference are presented in Section S2 of the online supplementary material.

To investigate the efficacy of our novel VI approach, we have performed a small preliminary simulation study to assess the coverage of Bayesian credible intervals generated using our procedure. Across the simulations, we varied the number of sites and the average occupancy and detection probability. We have chosen $n$, the number of sites, to be 500, 1,000, and 2,000 and $\psi$ and $p$ to be 0.25, 0.5, and 0.75. For each simulation, we ran 500 replications. The coverage was computed across 4 covariate coefficients for occupancy and 4 covariate coefficients for detection, with the coefficients randomly set to be either -1 or 1. Results are reported in Table 1 where it is seen that the 95% posterior credible intervals have the nominal coverage.

We also analysed the dataset of Ringlet butterflies collated through the Butterflies for the New Millennium (BNM) recording scheme run by Butterfly Conservation, using records collected between 1970 and 2014, which is also used in Diana et al. (2023). The data set consists of > 2 million records from ~ 140,000 unique 1 km² (defined as sites), of which > 218,000 detections of Ringlet have been made, and nondetections were produced using observations of other butterfly species (Kéry et al., 2010). In this case, we do not account for spatial autocorrelation, but we model year as a factor variable and assume independence between sites, a point which we discuss in the next section. We also use relative list length, obtained by dividing the list length, which is the number of species recorded for a given site/date, divided by the maximum recorded list length
<table><caption>**Table 1.** Coverage of the 95% Bayesian credible intervals generated using variational Bayes, for varying *n*, the number of sites, and values of *ψ* and *p*</caption><thead><tr><th>n / ψ = p</th><th>0.25</th><th>0.5</th><th>0.75</th></tr></thead><tbody><tr><td>500</td><td>0.950</td><td>0.949</td><td>0.961</td></tr><tr><td>1,000</td><td>0.961</td><td>0.953</td><td>0.952</td></tr><tr><td>2,000</td><td>0.951</td><td>0.948</td><td>0.951</td></tr></tbody></table>

in a neighbouring square area of 50x50 km, as a covariate for detection probability and model the effect of the day of the year using a second-order polynomial of Julian date.

In this case, convergence was determined by assessing when parameter updates were smaller than a prespecified tolerance, which required approximately 1 hr of computation time (results were again obtained on an Intel Core i7-10610U@1.80 GHz). In Figure 3, we present the posterior distribution of the occupancy index, which is the year-specific occupancy probability, here assumed the same for all sites. The increase in precision over time reflects the growth of the underlying opportunistic CS data set, which has shown ongoing expansion, particularly since the mid-1990s (Fox et al., 2023).

The occupancy probabilities show similarities with results presented in Diana et al. (2023), in particular an increasing trend in recent years reflecting an expansion in the range of Ringlet in the UK and similar drops in the species' prevalence in mid-1970s and early and late-1990s are identified. However, in this case, the pattern is less smooth, as expected, since the effect of year is not constrained in any way. Furthermore, a direct comparison cannot be made, since Diana et al. (2023) fit a more complex model accounting for spatio-temporal correlation, and plot an occupancy index rather than the annual occupancy probabilities.

In Figure 4, we present the posterior distribution of the detection probability throughout the year which can be interpreted as an estimate of species' flying time. Unlike in Diana et al. (2023), in this example detection probability does not vary with year, and thus represents an average for 1970–2014. We also note here that we do not account for spatial effects on the probability of detection. However, the species' flight patterns are likely to vary in space, and hence the curve in Figure 4 represents a mixture of spatially varying curves. Additionally, we note that the peak of detection probability in this case is lower than that obtained by Diana et al. (2023), which is due to the different model structure for detection and occupancy probabilities, as discussed above. The 95% posterior credible interval of the coefficient of relative list length is (0.881, 0.883).

We have shown that using a VI framework to fit occupancy models efficiently to CS data shows promising results, with various avenues for future development—some of which we discuss in the next section.

# 5 Discussion

CS data have an important rôle in the future of biodiversity monitoring, but the ongoing growth of such data requires novel statistical models along with efficient inference methods and available computer code. In this paper, we have built upon recent work and demonstrated that concentrated likelihoods can greatly reduce the dimensionality of the effective model parameter space and that variational inference can substantially reduce computational time for large CS data sets. The approaches proposed in this paper address computational problems, but also lead to model developments for CS data.

The GAI is efficient due to the use of concentrated likelihood (Dennis et al., 2016). There is a broad analogy here with the efficiency that can result from adopting a hidden Markov model, when appropriate, and forming a forward algorithm for likelihood construction; see for example Cowen et al. (2017). We have shown concentrated likelihood to be useful in the new extension of the GAI to incorporate the annual model which conveniently allows for formal variance propagation, negating the need for time-consuming bootstrapping. A similar approach is given in Bravington et al. (2021), where the two stages for density surface models, one involving a detection probability and the other a GAM, are combined into one.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** 95% posterior credible intervals of the occupancy probabilities of each year for Ringlet. The dots represent the posterior medians.</figcaption></figure>

The extended GAI has been introduced in this paper, with application to two univoltine species, but future work will test wider application to more species, in particular with varying quantities of data and life histories, including species with two or more generations per year, requiring more complex functions to describe seasonal variation (Dennis et al., 2016). As mentioned, alternative distributions to the Poisson may also be explored, in particular to assess for effects of overdispersion and the associated impact on estimates of uncertainty.

In both adaptations of the GAI presented, species' flight periods were assumed to be the same across sites within each year, as is typical in the production of abundance trends for UK butterflies (UKBMS, 2023). However, greater flexibility in the seasonal variation function a can be readily accounted for through the inclusion of appropriate spatial covariates. For example this is done by Schmucki et al. (2016) and Dennis et al. (2022). This is similarly a future direction to explore for occupancy models, accounting for spatial variation in detection probabilities, for example due to spatial variation in species' phenologies, via spatial covariates or spatial effects. See also for example Clark and Altwegg (2019) and Dennis et al. (2019) for models with spatial variation in occupancy.

The extended GAI also provides a basis for further extensions within this computationally efficient approach, such as formal data integration by maximizing a joint likelihood; see Besbeas et al. (2002) and Schaub and Kéry (2021). For example, UKBMS (count data from standardized monitoring) and BBC (mass-participation CS data) could be used to produce new urban butterfly indicators. Integrated modelling approaches can optimize the use of available CS data sets, but present outstanding questions and challenges (Isaac et al., 2020; Johnston et al., 2023; Zipkin et al., 2021), including the need for computationally efficient approaches for data integration.

We have introduced the basic occupancy model within a VI framework and have discussed the use of the observed, instead of the complete, data likelihood for avoiding the standard issue of underestimated posterior variances when parameters are assumed to be independent. In this case, the dependence structure in the model parameters and latent variables was such that the use of the observed data likelihood allowed us to use VI without having to assume that parameters are a-posteriori independent. Intuitively, the same approach could be employed in other ecological models within a VI framework, although at the moment this is only an intuition and future work would need to explore the quality of inference for different types of data and corresponding models.

We believe that VI provides a powerful and versatile framework for efficiently fitting a wide range of ecological models. The advantage of VI is that it combines the speed of classical inference
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Posterior median and 95% posterior credible intervals of the detection probability in each week for Ringlet.</figcaption></figure>

with the interpretability of Bayesian inference. VI relies on the ability to obtain the gradient of the likelihood function, which might seem like an obstacle. However, it is possible to take advantage of recent developments of deep learning methodologies such as the use of automatic differentiation, which can automatically compute gradients of this type (as long as the likelihood function is tractable), such as for example using the package TMB (Kristensen et al., 2016). We envisage that analyses currently based on hidden Markov models for likelihood computation—see for example Cowen et al. (2017), Besbeas and Morgan (2019, 2020) and McClintock et al. (2020)—can be a fruitful research avenue, since automatic differentiation engines would allow us to compute gradients by differentiating through the forward recursion used to compute the likelihood in this case. The use of automatic differentiation also enables us to consider highly nonlinear extensions, such as neural networks, to be introduced in the model. However, as discussed, particular attention needs to be paid to the chosen variational family, since that will determine the accuracy of the approximation.

We have considered a simple occupancy model, in comparison to models that have been employed previously for data of this type, as a means of discussing VI and the use of observed data likelihood within the VI framework. The flexible spatio-temporal models, for example based on Gaussian processes (Diana et al., 2023; Doser et al., 2022), which have been developed within an MCMC framework could be considered within this VI framework in future work. However, since this leads to a regression with as many covariates as the number of support of points, assuming a variational approximation with a full covariance matrix is computationally prohibitive, since the number of parameters of the Cholesky factor C scales quadratically with the number of covariates. One option to overcome this problem is to induce sparsity on the inverse of the covariance matrix by zeroing elements of the Cholesky factor C (Tan & Nott, 2018). For example, it is possible to assume a variational approximation where the covariate coefficients for the spatial approximation are independent in the posterior. Although this step can potentially reintroduce bias in the model, since it assumes a-posteriori independence of parameters, it leads to a substantial reduction in the number of parameters, making the model feasible to estimate. Finally, in cases where the observed data likelihood cannot be obtained, for example in complex data-generating processes or models with individual random effects, then the complete data likelihood and efficient MCMC inference (see for example King et al., 2023, who devise an importance sampling approach for ecological models with random effects) may provide the only viable alternative, at least at the moment.
To some extent formal design considerations do not arise with the data that we have considered, however there is an issue of nonrandom sampling (Boyd, Powney, et al., 2023; Johnston et al., 2023) and the need to account for issues such as preferential sampling and spatial and temporal biases (Altwegg & Nichols, 2019; Boersch-Supan et al., 2019; Conn et al., 2017; Pati et al., 2011). See for example papers by King et al. (2023) and by Lahoz-Monfort et al. (2014), respectively on sampling the data, and on how to design studies when resources are limited.

We have demonstrated efficient analysis methods for sources of CS data for UK butterflies, but efficient statistical inference methods are needed for understanding population changes from CS data for a wide range of taxa and locations. For example the GAI approach, or related models, has been applied to moths, bees and beetles (Dennis et al., 2021; Fox et al., 2021; Matechou et al., 2018). Development of efficient inference for occupancy models is also vital given their application to various taxa (for example Boyd, August, et al., 2023; Burns et al., 2023; Outhwaite et al., 2019). The need for efficiency will continue to increase with the growth of data sets such as the Global Biodiversity Information Facility (GBIF), which has amassed more than 2.5 billions occurrences of more than one million species (GBIF.org, 2023).

CS data are increasingly ‘big’, not only in terms of volume, but also involving characteristics such as variety. Farley et al. (2018) and McCrea et al. (2023) discuss the ‘Four Vs Framework’ in which data may be characterized as ‘big’. Analysing CS data for biodiversity monitoring presents various challenges (Johnston et al., 2023), but methods also need to be suitably scalable for these increasingly large data sets.

The examples in this paper are based upon analyses of data featuring observations of species submitted by citizen scientists, but computational challenges also arise from other data types, for example with the growth of data from technological advances such as automated interpretation of images submitted by citizen scientists (Terry et al., 2020; van Klink et al., 2022).

It may be argued that high performance computing (HPC) and cloud computing can be used to address the challenge of fitting computationally demanding models to CS data (Farley et al., 2018). For example, a supercomputer has enabled the production of occupancy trend estimates for thousands of UK species using Bayesian occupancy models (Boyd, August, et al., 2023; Outhwaite et al., 2019). However ultimately, as data sets continue to grow, and models become increasingly complex, we argue that a trade-off is needed, and that using more and more computing resources is not a simple solution. Fitting computationally demanding methods can make appropriate model validation and inference difficult: for example variable selection may become impractical (Johnston et al., 2023), as well as suitable goodness-of-fit assessment. Achieving model convergence for all parameters in MCMC may also become difficult (Boyd, August, et al., 2023; Outhwaite et al., 2019).

Furthermore, to maximize the use of CS data in biodiversity monitoring, there is a need for statistical approaches that are appropriately disseminated and accessible for use in practice. Johnston et al. (2023) suggest that ‘accessible communication of novel methods could democratize analysis of these data and thus enable CS data to reach their broadest potential’. Statistical methods that depend upon HPC may be a barrier for analysis to those without easy or affordable access to such resources, thus hindering the potential of biodiversity monitoring with CS data globally (Pocock et al., 2018).

Efficient methods are also crucial for producing frequent analysis updates for reporting on the status of species, particularly as time lags in data availability continue to reduce. The need for accurate reporting is ever necessary in monitoring species’ status, measuring against biodiversity targets, supporting policy-making and guiding effective conservation effort. This paper has presented examples for fitting computationally efficient models to CS data, but, as also suggested by Johnston et al. (2023), with the growth of such data and its importance for biodiversity monitoring (Callaghan et al., 2021; Pocock et al., 2018), there is an ongoing need to develop efficient statistical inference methods, with the potential to learn from developments in mainstream statistics.

## Acknowledgments

We thank the Associate Editor and two referees for their positive and helpful reviews. The UKBMS is organized and funded by Butterfly Conservation, the British Trust for Ornithology (BTO), and the Joint Nature Conservation Committee (JNCC). The UKBMS is indebted to all volunteers who
contribute data to the scheme. We are also very grateful to all of the volunteers who have contrib-
uted to the Butterflies for the New Millennium project, which is run by Butterfly Conservation
with support from Natural England, and also the citizen scientists who submitted Big Butterfly
Count data.

*Conflicts of interest:* None declared.

Funding

B.M. was supported by a Leverhulme Emeritus Fellowship, EM-2017-061. This work was partly funded by NERC grant NE/T010045/1.

Data availability

Data and code associated with the analyses in this paper are available as follows: (i) example of the extended GAI approach applied to UKBMS data https://github.com/EBDennis/Extended_GAI_UKBMS_example (ii) example of the GAI applied to BBC data https://github.com/EBDennis/GAI_BBC_example (iii) example of VI applied to BNM data for Ringlet https://github.com/AlexDiana/VBOccupancy.

Supplementary material

Supplementary material is available online at *Journal of the Royal Statistical Society: Series A*.

Author contributions statement

All authors devised the paper. E.D. curated the data sets. E.D. did the computing for the new GAI modelling and A.D. did the computing for the VI analyses. All authors collaborated in the writing of the paper, and E.D. and A.D. produced the online supplementary material.

References

Altwegg R., & Nichols J. D. (2019). Occupancy models for citizen-science data. *Methods in Ecology and Evolution, 10*(1), 8-21. https://doi.org/10.1111/mee3.2019.10 issue-1

Besbeas P., Freeman S. N., Morgan B. J. T., & Catchpole E. A. (2002). Integrating mark–recapture–recovery and census data to estimate animal abundance and demographic parameters. *Biometrics*, *58*(3), 540–547. https://doi.org/10.1111/j.0006-341X.2002.00540.x

Besbeas P., & Morgan B. J. T. (2019). Exact inference for integrated population modelling. *Biometrics*, 75(2), 475-484. https://doi.org/10.1111/biom.13045

Besbeas P., & Morgan B. J. T. (2020). A general framework for modelling population abundance data. *Biometrics*, 76(1), 281-292. https://doi.org/10.1111/biom.v76.1

Blei D. M., Kucukelbir A., & McAuliffe J. D. (2017). Variational inference: A review for statisticians. *Journal of the American Statistical Association*, 112(518), 859-877. https://doi.org/10.1080/01621459.2017.1285773

Boersch-Supan P. H., Trask A. E., & Baillie S. R. (2019). Robustness of simple avian population trend models for semi-structured citizen science data is species-dependent. *Biological Conservation*, 240, 108286. https://doi.org/10.1016/j.biocon.2019.108286

Bogaart P., van der Loo M., & Pannekoek J. (2020). *rtrim: Trends and Indices for Monitoring Data*. https://CRAN.R-project.org/package=rtrim. R package version 2.1.1.

Borowska A., & King R. (2022). Semi-complete data augmentation for efficient state space model fitting. *Journal of Computational and Graphical Statistics*, 32(1), 19–35. https://doi.org/10.1080/10618600.2022.2077350

Boyd R. J., August T. A., Cooke R., Logie M., Mancini F., Powney G. D., Roy D. B., Survey K., & Isaac N. J. (2023). An operational workflow for producing periodic estimates of species occupancy at national scales. *Biological Reviews*, 98(5), 1492–1508. https://doi.org/10.1111/brv.98.5

Boyd R. J., Powney G. D., & Pescott O. L. (2023). We need to talk about nonprobability samples. *Trends in Ecology & Evolution*, 38(6), 521–531. https://doi.org/10.1016/j.tree.2023.01.001

Bravington M. V., Miller D. L., & Hedley S. L. (2021). Variance propagation for density surface models. *Journal of Agricultural, Biological, and Environmental Statistics*, 26(2), 306–323. https://doi.org/10.1007/s13253-021-00438-2

BRC (2022). Biological Records Centre Home Page. www.brc.ac.uk Date accessed November 13, 2023.
Brereton T. M., Botham M. S., Middlebrook I., Randle Z., Noble D., & Harris S. E. A. (2018). *United Kingdom butterfly monitoring scheme report for 2017*. (Technical Report). Centre for Ecology & Hydrology, Butterfly Conservation, British Trust for Ornithology and Joint Nature Conservation Committee.

Burns F., Mordue S., al Fulaij N., Boersch-Supan P. H., Boswell J., Boyd R. J., Bradfer-Lawrence T., de Ornellas P., de Palma A., de Zylva P., Dennis E. B., Foster S., Gilbert G., Halliwell L., Hawkins K., Haysom K. A., Holland M. M., Hughes J., Jackson A. C., ...Gregory R. D. (2023). State of Nature 2023, the State of Nature partnership, www.stateofnature.org.uk.

Butchart S. H. M., Walpole M., Collen B., van Strien A., Scharlemann J. P. W., Almond R. E. A., Baillie J. E. M., Bomhard B., Brown C., Bruno J., Carpenter K. E., Carr G. M., Chanson J., Chenery A. M., Csirke J., Davidson N. C., Dentener F., Foster M., Galli A., ...Watson R. (2010). Global biodiversity: Indicators of recent declines. *Science*, 328(5982), 1164–1168. https://doi.org/10.1126/science.1187512

Callaghan C. T., Poore A. G., Mesaglio T., Moles A. T., Nakagawa S., Roberts C., Rowley J. J., Vergés A., Wilshire J. H., & Cornwell W. K. (2021). Three frontiers for the future of biodiversity research using citizen science data. *BioScience*, 71, 55–63. https://doi.org/10.1093/biosci/biaa131

Chandler M., See L., Copas K., Bonde A. M. Z., López B. C., Danielsen F., Legind J. K., Masinde S., Miller-Rushing A. J., Newman G., Rosemartin A., & Turak E. (2017). Contribution of citizen science towards international biodiversity monitoring. *Biological Conservation*, 213, 280–294. https://doi.org/10.1016/j.biocon.2016.09.004

Clark A. E., & Altwegg R. (2019). Efficient Bayesian analysis of occupancy models with logit link functions. *Ecology and Evolution*, 9(2), 756–768. https://doi.org/10.1002/ece3.2019.9 issue-2

Clark A. E., Altwegg R., & Ormerod J. T. (2016). A variational Bayes approach to the analysis of occupancy models. *PloS One*, 11(2), e0148966. https://doi.org/10.1371/journal.pone.0148966

Conn P. B., Thorson J. T., & Johnson D. S. (2017). Confronting preferential sampling in wildlife surveys: Diagnosis and model-based triage. *Methods in Ecology and Evolution*, 8(11), 1535–1546. https://doi.org/10.1111/mee3.2017.8 issue-11

Cowan L., Besbeas P. T., Morgan B. J. T., & Schwarz C. (2017). Hidden Markov models for extended batch data. *Biometrics*, 73(4), 1321–1331. https://doi.org/10.1111/biom.12701

Dennis E. B., Brereton T. M., Morgan B. J. T., Fox R., Shortall C. R., Prescott T., & Foster S. (2019). Trends and indicators for quantifying moth abundance and occupancy in Scotland. *Journal of Insect Conservation*, 23(2), 369–380. https://doi.org/10.1007/s10841-019-00135-z

Dennis E. B., Fagard-Jenkin C., & Morgan B. J. T. (2022). rGAI: An R package for fitting the generalized abundance index to seasonal count data. *Ecology and Evolution*, 12(8), e9200. https://doi.org/10.1002/ece3.9200

Dennis E. B., Freeman S. N., Brereton T., & Roy D. B. (2013). Indexing butterfly abundance whilst accounting for missing counts and variability in seasonal pattern. *Methods in Ecology and Evolution*, 4(7), 637–645. https://doi.org/10.1111/mee3.2013.4 issue-7

Dennis E. B., Kéry M., Morgan B. J. T., Coray A., Schaub M., & Baur B. (2021). Integrated modelling of insect population dynamics at two temporal scales. *Ecological Modelling*, 441, 109408. https://doi.org/10.1016/j.ecolmodel.2020.109408

Dennis E. B., Morgan B. J. T., Brereton T. M., Roy D. B., & Fox R. (2017). Using citizen science butterfly counts to predict species population trends. *Conservation Biology*, 31(6), 1350–1361. https://doi.org/10.1111/cobi.2017.31 issue-6

Dennis E. B., Morgan B. J. T., Freeman S. N., Brereton T. M., & Roy D. B. (2016). A generalized abundance index for seasonal invertebrates. *Biometrics*, 72(4), 1305–1314. https://doi.org/10.1111/biom.12506

Dennis E. B., Morgan B. J. T., Freeman S. N., Ridout M. S., Brereton T. M., Fox R., Powney G. D., & Roy D. B. (2017). Efficient occupancy model-fitting for extensive citizen-science data. *PLoS One*, 12(3), e0174433. https://doi.org/10.1371/journal.pone.0174433

Dennis E. B., Morgan B. J. T., Harrower C. A., Bourn N. A. D., & Fox R. (2024). Incorporating phenology to estimate species' population trends from snapshot citizen-science data. *Journal of Agricultural, Biological, and Environmental Statistics, To appear*.

Diana A., Dennis E. B., Matechou E., & Morgan B. J. T. (2023). Fast Bayesian inference for large occupancy data sets, using the Pólya-Gamma scheme. *Biometrics*, 79, 2503–2515. https://doi.org/10.1111/biom.13816

Didham R. K., Basset Y., Collins C. M., Leather S. R., Littlewood N. A., Menz M. H. M., Müller J., Packer L., Saunders M. E., Schönrogge K., Stewart A. J. A., Yanoviak S. P., & Hassall C. (2020). Interpreting insect declines: Seven challenges and a way forward. *Insect Conservation and Diversity*, 13(2), 103–114. https://doi.org/10.1111/icad.v13.2

Doser J. W., Finley A. O., & Banerjee S. (2023). Joint species distribution models with imperfect detection for high-dimensional spatial data. *Ecology*, 104, e4137. https://doi.org/10.1002/ecy.4137

Doser J. W., Finley A. O., Kéry M., & Zipkin E. F. (2022). spoccupancy: An R package for single-species, multi-species, and integrated spatial occupancy models. *Methods in Ecology and Evolution*, 13, 1670–1678. https://doi.org/10.1111/mee3.v13.8
Farley S. S., Dawson A., Goring S. J., & Williams J. W. (2018). Situating ecology as a big-data science: Current advances, challenges, and solutions. *BioScience*, **68**, 563-576. https://doi.org/10.1093/biosci/biy068

Fox R., Dennis E. B., Brown A. F., & Curson J. (2022). A revised red list of British butterflies. *Insect Conservation and Diversity*, **15**, 485-495. https://doi.org/10.1111/icad.v15.5

Fox R., Dennis E. B., Harrower C. A., Blumgart D., Bell J. R., Cook P., Davis A. M., Evans-Hill L. J., Haynes F., Hill D., Isaac N. J. B., Parsons M. S., Pocock M. J. O., Prescott T., Randle Z., Shortall C. R., Tordoff G. M., Tucson D., & Bourn N. A. D. (2021). The State of Britain's Larger Moths 2021. Butterfly Conservation, Rothamsted Research and UK Centre for Ecology & Hydrology, Wareham, Dorset, UK.

Fox R., Dennis E. B., Purdy K. M., Middlebrook I., Roy D. B., Noble D. G., Botham M. S., & Bourn N. A. D. (2023). The State of the UK's Butterflies 2022. Butterfly Conservation, Wareham, UK.

GBIF.org (2023). GBIF Home Page. https://www.gbif.org Date accessed November 13, 2023.

Isaac N. J. B., Jarzyna M. A., Keil P., Dambly L. I., Boersch-Supan P. H., Browning E., Freeman S. N., Golding N., Guillera-Arroita G., Henrys P. A., Jarvis S., Lahoz-Monfort J., Pagel J., Pescott O. L., Schmucki R., Simmonds E. G., & O'Hara R. B. (2020). Data integration for large-scale models of species distributions. *Trends in Ecology & Evolution*, **35**(1), 56-67. https://doi.org/10.1016/j.tree.2019.08.006

Isaac N. J. B., & Pocock M. J. O. (2015). Bias and information in biological records. *Biological Journal of the Linnean Society*, **115**(3), 522-531. https://doi.org/10.1111/bij.12532

Isaac N. J. B., van Strien A. J., August T. A., de Zeeuw M. P., & Roy D. B. (2014). Statistics for citizen science: Extracting signals of change from noisy ecological data. *Methods in Ecology and Evolution*, **5**, 1052-1060. https://doi.org/10.1111/mee3.2014.5 issue-10

JNCC (2022). UK Biodiversity Indicators 2022 - C6 Insects of the wider countryside (butterflies). https://jncc.gov.uk/our-work/ukbi-c6-insects-of-the-wider-countryside/ Date accessed November 8, 2023.

Johnston A., Matechou E., & Dennis E. B. (2023). Outstanding challenges and future directions for biodiversity monitoring using citizen science data. *Methods in Ecology and Evolution*, **14**, 103-116. https://doi.org/10.1111/mee3.v14.1

Jordan M. I., Ghahramani Z., Jaakkola T. S., & Saul L. K. (1999). An introduction to variational methods for graphical models. *Machine Learning*, **37**, 183-233. https://doi.org/10.1023/A:1007665907178

Kéry M., Gardner B., & Monnerat C. (2010). Predicting species distributions from checklist data using site-occupancy models. *Journal of Biogeography*, **37**, 1851-1862. https://doi.org/10.1111/j.1365-2699.2010.02345.x

King R. (2014). Statistical ecology. *Annual Review of Statistics and Its Application*, **1**, 401-426. https://doi.org/10.1146/annurev-statistics-022513-115633

King R., Sarzo B., & Elvira V. (2023). When ecological individual heterogeneity models and large data collide: An importance sampling approach. *The Annals of Applied Statistics*, **17**, 3112-3132. https://doi.org/10.1214/23-AOAS1753

Kristensen K., Nielsen A., Berg C. W., Skaug H., & Bell B. M. (2016). TMB: Automatic differentiation and Laplace approximation. *Journal of Statistical Software*, **70**, 1-21. https://www.jstatsoft.org/index.php/jss/article/view/v070i05. https://doi.org/10.18637/jss.v070.i05

Lahoz-Monfort J. J., Harris M. P., Morgan B. J. T., Freeman S. N., & Wanless S. (2014). Exploring the consequences of reducing survey effort for detecting individual and temporal variability in survival. *Journal of Applied Ecology*, **51**, 534-543. https://doi.org/10.1111/1365-2664.12214

MacKenzie D. I., Nichols J. D., Royle J. A., Pollock K. H., Bailey L. L., & Hines J. E. (2018). *Occupancy estimation and modeling: Inferring patterns and dynamics of species occurrence* (2nd ed.). Academic Press.

Matechou E., Dennis E. B., Freeman S. N., & Brereton T. (2014). Monitoring abundance and phenology in (multivoltine) butterfly species: A novel mixture model. *Journal of Applied Ecology*, **51**, 766-775. https://doi.org/10.1111/jpe.2014.51 issue-3

Matechou E., Freeman S. N., & Comont R. (2018). Caste-specific demography and phenology in bumblebees: Modelling Beewalk data. *Journal of Agricultural, Biological and Environmental Statistics*, **23**, 427-445. https://doi.org/10.1007/s13253-018-0332-y

McClintock B. T., Langrock R., Gimenez O., Cam E., Borchers D. L., Glennie R., & Patterson T. A. (2020). Uncovering ecological state dynamics with hidden Markov models. *Ecology Letters*, **23**, 1878-1903. https://doi.org/10.1111/ele.v23.12

McCrea R., King R., Graham L., & Börger L. (2023). Realising the promise of large data and complex models. *Methods in Ecology and Evolution*, **14**, 4-11. https://doi.org/10.1111/mee3.v14.1

McCullagh P., & Nelder J. A. (1989). *Generalized linear models* (2nd ed.). Chapman & Hall.

Newman K., King R., Elvira V., de Valpine P., McCrea R. S., & Morgan B. J. T. (2023). State-space models for ecological time series: Practical model-fitting. *Methods in Ecology and Evolution*, **14**, 26-42. https://doi.org/10.1111/mee3.v14.1

Outhwaite C. L., Powney G. D., August T. A., Chandler R. E., Rorke S., Pescott O. L., Harvey M., Roy H. E., Fox R., Roy D. B., Alexander K., Ball S., Bantock T., Barber T., Beckmann B. C., Cook T., Flanagan J., Fowles A., Hammond P., ...Isaac N. J. B. (2019). Annual estimates of occupancy for bryophytes, lichens and invertebrates in the UK, 1970-2015. *Scientific Data*, **6**, 1-12. https://doi.org/10.1038/s41597-019-0269-1
Pati D., Reich B. J., & Dunson D. B. (2011). Bayesian geostatistical modelling with informative sampling locations. *Biometrika*, 98, 35–48. https://doi.org/10.1093/biomet/asq067

Pocock M. J. O., Chandler M., Bonney R., Thornhill I., Albin A., August T., Bachman S., Brown P. M. J., Cunha D. G. F., Grez A., Jackson C., Peters M., Rabarijaon N. R., Roy H. E., & Danielsen F., 2018). A vision for global biodiversity monitoring with citizen science. In *Advances in Ecological Research* (Vol. 59, pp. 169–223). Elsevier.

Pocock M. J. O., Roy H. E., Preston C. D., & Roy D. B. (2015). The biological records centre: A pioneer of citizen science. *Biological Journal of the Linnean Society*, 115, 475–493. https://doi.org/10.1111/bij.12548

Pocock M. J. O., Tweddle J. C., Savage J., Robinson L. D., & Roy H. E. (2017). The diversity and evolution of ecological and environmental citizen science. *PloS One*, 12, e0172579. https://doi.org/10.1371/journal.pone.0172579

Pollard E., & Yates T. J. (1993). *Monitoring butterflies for ecology and conservation: The British butterfly monitoring scheme*. Chapman & Hall.

Polson N. G., Scott J. G., & Windle J. (2013). Bayesian inference for logistic models using Pólya-Gamma latent variables. *Journal of the American Statistical Association*, 108, 1339–1349. https://doi.org/10.1080/01621459.2013.829001

Rezende D. J., Mohamed S., & Wierstra D. (2014). Stochastic backpropagation and approximate inference in deep generative models. In *International conference on machine learning* (pp. 1278–1286). PMLR.

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

Schmucki R., Pe'er G., Roy D. B., Stefanescu C., Van Swaay C. A. M., Oliver T. H., Kuussaari M., Van Strien A. J., Ries L., Settele J., Musche M., Carnicer J., Schweiger O., Brereton T. M., Harpke A., Heliölä J., Kühn E., & Julliard R. (2016). A regionally informed abundance index for supporting integrative analyses across butterfly monitoring schemes. *Journal of Applied Ecology*, 53, 501–510. https://doi.org/10.1111/jpe.2016.53 issue-2

Silvertown J. (2009). A new dawn for citizen science. *Trends in Ecology & Evolution*, 24, 467–471. https://doi.org/10.1016/j.tree.2009.03.017

Tan L. S. L., & Nott D. J. (2018). Gaussian variational approximation with sparse precision matrices. *Statistics and Computing*, 28(2), 259–275. https://doi.org/10.1007/s11222-017-9729-7

ter Braak C. J. F., van Strien A. J., Meijer R., & Verstrael T. J. (1994). Analysis of monitoring data with many missing values: which method? In E. J. M. Hagemeijer, & T. J. Verstrael (Eds.), *Bird Numbers 1992. Distribution, Monitoring and Ecological Aspects: Proceedings of the 12th International Conference of IBCC and EOAC, Noordwijkerhout, The Netherlands* (pp. 663–673). Statistics Netherlands, Voorburg/Heerlen & SOVON, Beek-Ubbergen.

Terry J. C. D., Roy H. E., & August T. A. (2020). Thinking like a naturalist: Enhancing computer vision of citizen science images by harnessing contextual data. *Methods in Ecology and Evolution*, 11(2), 303–315. https://doi.org/10.1111/mee3.112

Thomas C., Jones T. H., & Hartley S. E. (2019). ‘Insectageddon’: A call for more robust data and rigorous analyses. *Global Change Biology*, 25(6), 1891–1892. https://doi.org/10.1111/gcb.2019.25 issue-6

Titsias M., & Lázaro-Gredilla M. (2014). Doubly stochastic variational Bayes for non-conjugate inference. In *International conference on machine learning* (pp. 1971–1979). PMLR.

UKBMS (2023). UKBMS official statistics. https://ukbms.org/official-statistics Date accessed November 8, 2023.

van Klink R., August T., Bas Y., Bodesheim P., Bonn A., Fossøy F., Hoye T. T., Jongejans E., Menz M. H. M., Miraldo A., Roslin T., Roy H. E., Ruczyński I., Schigel D., Schäffler L., Sheard J. K., Svenningsen C., Tschan G. F., Wäldchen J., ...Bowler D. E. (2022). Emerging technologies revolutionise insect ecology and monitoring. *Trends in Ecology & Evolution*, 37, 872–885. https://doi.org/10.1016/j.tree.2022.06.001

van Strien A., Pannekoek J., Hagemeijer W., & Verstrael T. (2004). A loglinear Poisson regression method to analyse bird monitoring data. In A. Anselin (Ed.), *Bird Numbers 1995, Proceedings of the International Conference and 13th Meeting of the European Bird Census Council, Pärnu, Estonia*. Bird Census News 13 (2000) (pp. 33–39).

Van Swaay C. A. M., Dennis E. B., Schmucki R., Sevillega C. G., Aghababyan K., Åström S., Balalaikins M., Bonelli S., Botham M., Bourn N., Brereton T., Cancela J. P., Carlisle B., Collins S., Dopagne C., Dziekanska I., Escobés R., Faltynek Fric Z., Feldmann R., ...Roy D. B. (2020). Assessing Butterflies in Europe - Butterfly Indicators 1990-2018 Technical report. Butterfly Conservation Europe & ABLE/eBMS (www.butterfly-monitoring.net).

Wagner D. L., Grames E. M., Forister M. L., Berenbaum M. R., & Stopak D. (2021). Insect decline in the anthropocene: Death by a thousand cuts. *Proceedings of the National Academy of Sciences*, 118(2), e2023989118. https://doi.org/10.1073/pnas.2023989118

Wang B., & Titterington D. M. (2005). Inadequacy of interval estimates corresponding to variational Bayesian approximations. In *International Workshop on Artificial Intelligence and Statistics* (pp. 373–380). PMLR.

Zipkin E. F., Zylstra E. R., Wright A. D., Saunders S. P., Finley A. O., Dietze M. C., Itter M. S., & Tingley M. W. (2021). Addressing data integration challenges to link ecological processes across scales. *Frontiers in Ecology and the Environment*, 19(1), 30–38. https://doi.org/10.1002/fee.v19.1
# Frequentist prediction sets for species abundance using indirect information

Elizabeth Bersson ID and Peter D. Hoff

Department of Statistical Science, Duke University, Durham, NC 27701, USA

*Address for correspondence:* Elizabeth Bersson, Department of Statistical Science, Duke University, Durham, NC 27701, USA. Email: elizabeth.bersson@duke.edu

[Read before The Royal Statistical Society at the Discussion Meeting on the 'Analysis of citizen science data' held at the Society's 2024 annual conference in Brighton on Tuesday, 3 September 2024, the President, Dr Andrew Garrett, in the Chair.]

## Abstract

Citizen science databases that consist of volunteer-led sampling efforts of species communities are relied on as essential sources of data in ecology. Summarizing such data across counties with frequentist-valid prediction sets for each county provides an interpretable comparison across counties of varying size or composition. As citizen science data often feature unequal sampling efforts across a spatial domain, prediction sets constructed with indirect methods that share information across counties may be used to improve precision. In this article, we present a nonparametric framework to obtain precise prediction sets for a multinomial random sample based on indirect information that maintain frequentist coverage for each county. We detail a simple algorithm to obtain prediction sets for each county using indirect information where the computation time does not depend on the sample size and scales nicely with the number of species considered. The indirect information may be estimated by a proposed empirical Bayes procedure based on information from auxiliary data. Our approach makes inference for under-sampled counties more precise, while maintaining area-specific frequentist validity for each county. Our method is used to provide a useful description of avian species abundance in North Carolina, USA based on citizen science data from the eBird database.

**Keywords:** categorical data, conformal prediction, empirical Bayes, exchangeability, frequentist coverage, nonparametric

## 1 Introduction

Understanding species abundance across heterogeneous spatial areas is an important task in ecology. Citizen science databases that consist of observations of species counts gathered by volunteers are increasingly regarded as one of the richest sources of data for such a task. One of the largest such data sources is the eBird database in which citizen scientists throughout the world input counts of bird sightings (B. Sullivan et al., 2009). In addition to its use for describing avian species abundance, eBird is a principal resource for understanding global biodiversity and is widely used in constructing and implementing conservation action plans (B. L. Sullivan et al., 2017).

More generally, analyses from such databases may be used for informing policy, conservation efforts, habitat preservation, and more, for which understanding species prevalence for non-overlapping geographic areas, such as counties across a state or country, is important. In practice, species abundance from citizen science data are commonly summarized within areas such as counties by empirical proportions from a sample, as in, e.g. Arnold et al. (2021) and Camerini and Groppali (2014). Such proportions can be used to construct a prediction set for each county that provides a description of species prevalence for that county with guaranteed frequentist coverage.
Given the impact on policy design, corresponding uncertainty quantification is of particular import (Lele, 2020), and so it is desirable that precise prediction sets maintain a target coverage rate regardless of the county's size or composition. This is challenging as a common feature of citizen science data is unequal sampling efforts that results in some counties with large amounts of data information and others with very little. Using direct procedures that only make use of within-county information, a prediction set may be imprecise in these counties with low sampling efforts. This suggests using indirect information such as data from neighbouring counties to improve prediction set precision for a given county.

In this article, we describe species abundance across sampling areas such as counties with frequentist-valid prediction sets that are constructed to contain an unobserved bird with $1 - \alpha$ probability. That is, a valid prediction set for a given county is a set of avian species such that an unobserved bird will belong to one of those species with $1 - \alpha$ probability in a frequentist sense. We develop a valid nonparametric prediction method that allows for information to be shared across counties. Specifically, our approach results in prediction sets with guaranteed frequentist coverage for each county that are constructed with the incorporation of indirect or prior information. We detail and provide code for an empirical Bayes procedure to estimate such prior information from auxiliary data such as neighbouring counties. If this indirect information used to construct the prediction sets is accurate, the prediction sets will be smaller than direct sets that only make use of within-county information.

Identifying avian species that are present across a spatial domain is of general interest in ecology (Lebrun et al., 2012; Shanahan & Possingham, 2009; Twedt et al., 2010). To this aim, in Section 4, we detail the usefulness of the proposed approach in summarizing the eBird citizen science data. For this task, frequentist-valid prediction sets provide useful summaries of the data that may be used to compare information across subregions and better inform policy. In particular, a frequentist-valid prediction set consists of species that are likely to be present in a given area and can be used to draw statistically robust conclusions regarding future observations. Such a prediction set is constructed based on observed species counts, and, by maintaining a specified frequentist coverage rate, reflects the uncertainty in the data. Commonly reported summaries of observed species abundance data, such as a list of observed proportions or counts do not reflect the uncertainty in the data. Moreover, the approach we propose constructs frequentist-valid prediction sets with the incorporation of auxiliary information. This is desirable as sharing information across counties generally results in smaller prediction sets as compared to direct prediction approaches, particularly so in counties with low sampling efforts. In contrast, indirect approaches constructed with auxiliary information such as a Bayesian prediction set may reflect uncertainty in the data, but are not calibrated.

# 2 Methodology

## 2.1 Background and notation

For county $j \in \{1, \dots, J\}$, let $X_j$ be a vector of length $K$ where $X_{j,i} = x_{j,i}$ is the observed count of species $i$ over some set sampling period that may vary across counties. We model $X_j$ with a $K$-dimensional multinomial distribution with $N_j = \sum_{i=1}^K x_{j,i}$ trials and population proportions vector $\theta_j$,

$$X_j \sim MN_K(\theta_j, N_j). \qquad (1)$$

We construct a prediction set $A_\alpha(X_j)$ for an observation of a new bird arising from the same distribution, $Y_j \sim MN_K(\theta_j, 1)$, where $Y_j \in \mathcal{Y}$ for $\mathcal{Y} = \{(y_1, \dots, y_K) : \sum_{i=1}^K y_i = 1, y_i \in \{0, 1\} \ (i = 1, \dots, K)\}$. Let $y_j^{(k)} \in \mathcal{Y}$ denote a prediction of category $k$, that is, let $y_j^{(k)}$ be a vector of length $K$ with a one at index $k$ and zeros elsewhere. In particular, we are interested in a prediction set for $Y_j$ that maintains frequentist validity for some error rate $\alpha$. Formally, we refer to this as an $\alpha$-valid prediction set:
**Definition** (α-Valid Prediction Set). An α-valid prediction set for a predictand $Y_j \in \mathcal{Y}$ is any subset $A_\alpha$ of the sample space $\mathcal{Y}$ that contains $Y_j$ with probability greater than or equal to $1 - \alpha$,

$$P_{\theta}(Y_j \in A_{\alpha}(X_j)) \geq 1 - \alpha, \quad \forall \theta, \qquad (2)$$

where the probability is taken with respect to $Y_j$ and $X_j$.

Additionally, small or precise α-valid prediction sets are of particular interest, where prediction set size is measured by expected cardinality, that is, expected number of the K categories in the sample space included in the prediction set.

## 2.2 Order-based prediction for a single area

A standard approach to construct α-valid prediction sets for each county or area is with a direct method that only makes use of within-area information. As such, we first consider construction of a prediction set for a single area j, using only data from area j. For ease of notation, we drop the area-identifying subscript in this subsection.

For multinomial data in general, if the event probability vector $\theta$ is known, an α-valid prediction set is any combination of categories such that their event probabilities cumulatively sum to be greater than or equal to $1 - \alpha$. Equivalently stated, an α-valid prediction set may be constructed by excluding categories such that the cumulative sum of the excluded categories' event probabilities is less than $\alpha$. Such a prediction set may be constructed by admitting categories in some pre-specified order into the prediction set until the cumulative sum of their event probabilities is at least $1 - \alpha$. The resulting prediction set will have $1 - \alpha$ coverage regardless of the ordering used to admit categories. In fact, the class of all α-valid prediction sets may be constructed by following this procedure for non-strict total orderings of categories.

Perhaps intuitively, constructing such a prediction set by including categories with the largest event probabilities will result in the smallest α-valid prediction set. In the terminology of ordering, this corresponds to constructing a prediction set based on an ordering of categories that matches the ordering of the elements in $\theta$. We refer to this optimal ordering as the oracle ordering:

**Theorem 1** (Oracle order-based prediction). Let $Y \sim MN_K(\theta, 1)$ for $\theta$ known. Then,

1. the class of all α-valid prediction sets for a given $\theta$ consists of prediction sets of the form,

$$A_{\alpha}^{\theta, o} = \{ \mathbf{y}^{(k)} \in \mathcal{Y} : \left[ \sum_{l=1}^{K} \mathbb{1}(o_k \ge o_l) \theta_l \right] > \alpha \}, \qquad (3)$$

for some vector $o \in \mathbb{R}^K$, and

2. the **oracle ordering** is that which corresponds to the increasing order statistics of $\theta$,

$$o^{\theta} = \{o : \theta_m < \theta_n \Rightarrow o_m < o_n \quad \forall m, n \in \{1, \dots, K\}, m \neq n\},$$

and $A_{\alpha}^{\theta, o^{\theta}}$ has the smallest cardinality among all orderings.

In practice, $\theta$ is unknown, but a prediction set may be constructed based on an observed sample $X = x$. It turns out, in fact, that any conditional α-valid prediction set can be written similarly to the previous construction (equation (3)), where the cumulative sum is computed with respect to the empirical proportions given by $x$ and $\mathcal{Y}$. This is a generalization of the conformal prediction framework, a popular machine learning approach to construct prediction regions based on measuring conformity (or non-conformity) of a predictand to an observed sample (Vovk et al., 2005).
**Theorem 2** (α-valid order-based prediction). Let $X \sim MN_K(\theta, N)$, $Y \sim MN_K(\theta, 1)$. Then, every conformal α-valid prediction set based on observed data $\boldsymbol{x}$ can be written

$$A_{\alpha}(\boldsymbol{x}) = \left\{ \boldsymbol{y}^{(k)} \in \mathcal{Y} : \left[ \sum_{l=1}^{K} \mathbb{1}(o_k \ge o_l) \frac{x_l + \boldsymbol{y}_l^{(k)}}{N+1} \right] > \alpha \right\}, \quad (4)$$

for some vector $\boldsymbol{o} \in \mathbb{R}^K$.

Note that the prediction set depends on the vector $\boldsymbol{o}$ only through the order of its elements. Such a vector may be defined based on the data, or otherwise. We elaborate in the remainder of this section.

For any ordering of the K categories, constructing a prediction set following Theorem 2 results in a prediction set with guaranteed finite-sample $1 - \alpha$ frequentist coverage. The choice of ordering, however, will impact prediction set precision, that is, the set's cardinality. For inference for a single area, a natural approach is to order the categories with respect to their empirical proportions. The empirical proportions are unbiased for population proportions, so, if the area has a large sample size, an ordering based on the empirical proportions will approximate the oracle ordering well. It turns out this approach is well-motivated by classical prediction approaches. Specifically, a standard direct prediction method constructs a prediction set separately for an area based on an area-specific conditional pivotal quantity (Faulkenberry, 1973; Tian et al., 2022). For a multinomial population, $Y | X + Y$ is such a quantity that follows a multivariate hypergeometric distribution which does not depend on the event probability vector. See Thatcher (1964) for work on prediction sets of this type for binomial data. A prediction set constructed to contain species belonging to a highest mass region of this pivotal distribution is obtained by including species with the largest empirical counts until their cumulative proportion sum exceeds $1 - \alpha$,

$$A_{\alpha}^{D}(\boldsymbol{x}) = \left\{ \boldsymbol{y}^{(k)} \in \mathcal{Y} : \left[ \sum_{l=1}^{K} \mathbb{1}\left( (x_k + \boldsymbol{y}_k^{(k)}) \ge (x_l + \boldsymbol{y}_l^{(k)}) \right) \frac{x_l + \boldsymbol{y}_l^{(k)}}{N+1} \right] > \alpha \right\}. \quad (5)$$

This direct prediction set based on an ordering of the empirical proportions is appealing as it is easy to interpret and has finite-sample guaranteed $1 - \alpha$ frequentist coverage. For an area with low sampling effort, though, the empirical proportions will not precisely estimate the true proportions. As a result, a prediction set may have prohibitively large cardinality such that it is not practically useful. For such an area, incorporating indirect information from neighbouring counties can improve the estimates of the county proportions and thereby increase the precision of a prediction set.

## 2.3 Order-based prediction for multiple areas

In general, in analysing small area data, that is, areal data featuring small within-area sample sizes in some areas, it is common to utilize indirect methods that share information across areas (Rao & Molina, 2015). The eBird database is a rich data source, and inference in any given county may be improved upon by taking advantage of auxiliary data using an indirect method. In this subsection, we detail how information from neighbouring counties may be used in estimating an ordering of categories to improve prediction set precision.

As opposed to a direct prediction set based on an ordering corresponding to within-county empirical proportions, an indirect prediction set can be constructed similarly whereby species are admitted into the prediction set based on an ordering corresponding to empirical posterior proportions estimated from a hierarchical model. Such an estimate may be obtained based on a conjugate Dirichlet prior distribution parameterized with a common concentration hyperparameter for the J areas,

$$\theta_1, \dots, \theta_J \sim \text{Dirichlet}_K(\boldsymbol{\gamma}). \quad (6)$$
Given a hyperparameter $\gamma \in \mathbb{R}^K$, the posterior expectation of the proportions $\theta_j$ in county $j$ is $\tilde{x}_j/(N_j + \sum_{i=1}^K \gamma_i)$, where $\tilde{x}_j = x_j + \gamma$. In this way, $\tilde{x}_j$ may be interpreted as a posterior vector of counts for county $j$. Then, an $\alpha$-valid prediction set based on $\tilde{x}_j$ is,

$$A_{\alpha}^{I}(\mathbf{x}_j) = \left\{ \mathcal{Y}^{(k)} \in \mathcal{Y} : \left[ \sum_{l=1}^{K} \mathbb{1}\left( (\tilde{x}_{j,k} + \mathbf{y}_{k}^{(k)}) \geq (\tilde{x}_{j,l} + \mathbf{y}_{l}^{(k)}) \right) \frac{\mathbf{x}_{j,l} + \mathbf{y}_{l}^{(k)}}{N_j + 1} \right] > \alpha \right\}. \quad (7)$$

By Theorem 2, $A_{\alpha}^{I}(\mathbf{x}_j)$ is an $\alpha$-valid procedure, and it is constructed based on prior information. Specifically, it differs from the direct set given in equation (5) in that categories are admitted into the prediction set based on an ordering determined by posterior counts that incorporate indirect information $\gamma$, as opposed to an ordering based on the observed sample. Moreover, it has been shown that if the indirect information used is accurate, $A_{\alpha}^{I}(\mathbf{x}_j)$ may be more precise than a direct prediction set with the same coverage rate (Bersson & Hoff, 2024; Hoff, 2023).

In total, $A_{\alpha}^{D}(\mathbf{x}_j)$ and $A_{\alpha}^{I}(\mathbf{x}_j)$ are both $\alpha$-valid prediction procedures. They differ in the order in which species are admitted into the prediction sets, as species are admitted into the direct set in terms of decreasing empirical proportions and into the indirect set in terms of decreasing posterior counts. As a result, for an area with a small sample size, incorporating accurate prior information can result in an ordering used to construct a prediction set that more accurately approximates the oracle ordering as the empirical proportions might be too unstable. Of note, these two approaches are equivalent for a uniform prior $\gamma = c1$, for any constant $c$. This includes, for example, a standard noninformative prior $c = 1$, a standard objective Bayes Jeffrey's prior $c = 1/2$, and an improper prior $c = 0$.

## 2.4 Empirical Bayes estimation of indirect information

To obtain an $\alpha$-valid indirect prediction set for county $j \in \{1, \dots, J\}$, all that is required is an estimate of the prior concentration parameter $\gamma$. We propose an empirical Bayesian approach whereby values of $\gamma$ to be used for county $j$ are estimated from data collected in neighbouring counties. Specifically, we use the maximum likelihood estimate of the marginal likelihood based on the conjugate hierarchical model given by equations (1) and (6),

$$\begin{aligned} \gamma_j &= \arg \max_{\gamma} \log p\left(\bigcup_{l \in L_j} \mathbf{X}_l | \gamma\right) \\ &= \arg \max_{\gamma} \log \prod_{l \in L_j} \left[ \frac{\Gamma\left(\sum_{i=1}^K \gamma_i\right)}{\Gamma\left(\sum_{i=1}^K \mathbf{x}_{l,i} + \gamma_i\right)} \times \prod_{i=1}^K \frac{\Gamma(\mathbf{x}_{l,i} + \gamma_i)}{\Gamma(\gamma_i)} \right], \end{aligned} \quad (8)$$

where $L_j \subseteq \{1, \dots, J\} \setminus \{j\}$ is a non-empty set containing the indices of counties neighbouring county $j$. Information is shared across neighbouring counties to inform an estimate of the prior for county $j$, and, when estimated in this way, the prior concentration represents an across-county pooled prior concentration. This optimization problem can be solved numerically with a Newton-Raphson algorithm. See Appendix B for details and derivation of such an algorithm. Code to implement this procedure in the R Statistical Programming language is available online, see Section 5.

When $\gamma_j$ is estimated using data independent of area $j$ and used to construct $A_{\alpha}^{I}(\mathbf{x}_j)$, the finite-sample coverage guarantee of $A_{\alpha}^{I}(\mathbf{x}_j)$ holds regardless of the accuracy of the estimated prior hyper-parameter. If the estimated vector $\gamma_j$ is accurate, then $A_{\alpha}^{I}(\mathbf{x}_j)$ may also be more precise than direct prediction approaches.

## 3 Simulation study

To illustrate how the incorporation of indirect information can affect precision of prediction sets, we compare expected set cardinality obtained from the indirect and direct prediction methods for a single simulated area. In contrast to the eBird data, for example, the analysis of this section corresponds to that of one county. Because citizen science data such as these often feature unequal
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Monte Carlo approximations (±1 standard deviation) of the expected cardinality ratio of indirect method based on an accurate prior with moderate precision to direct method, and expected cardinality ratio of α-valid prediction set given the oracle ordering to direct method.</figcaption></figure>

sampling efforts across counties, we are particularly interested in demonstrating the difference in cardinality between these two approaches for a range of sample sizes *N* = 10, 100, 1,000. Moreover, we compare results for varying number of categories *K*. Throughout, we consider a low entropy regime in which ⌈*K*/4⌋ categories unequally split nearly all of the probability mass, and the rest of the categories have nearly probability 0. While we do not necessarily expect real populations in practice to have such a distribution, it is chosen to clearly demonstrate the benefit of including indirect information in the construction of prediction sets that maintain frequentist coverage.

In one construction of indirect prediction sets, we consider a prior based on full information with moderate prior precision $\gamma = \theta \times 10$. We compare with direct prediction sets given by equation (5), or, equivalently, indirect prediction sets constructed with a uniform prior $\gamma = 1$. Finally, we compare the approaches to α-valid order-based prediction sets obtained based on an oracle ordering. Results comparing Monte Carlo approximations of the expected prediction set cardinality ratios between the various approaches obtained from 25,000 replications are displayed in Figure 1.

As all methods considered are α-valid procedures, the crucial difference between them is the incorporation of indirect information. Utilizing accurate prior information in the construction of prediction sets generally results in prediction sets distinctly smaller than direct sets, particularly so if there are a large number of categories relative to the sample size. This is evidenced by the red dashes in Figure 1 showing the expected cardinality ratios of the indirect to direct prediction sets are always at or below a value of 1. An accurate prior may be one that approximates the true probability mass vector well with large precision relative to sample size, as seen in the left plot of Figure 1 for sample size *N* = 10. More generally, though, all that is needed is a prior that results in posterior counts that accurately approximate the oracle ordering of categories. We discuss the three sample size regimes in detail below.

For a small sample size of *N* = 10, the prior $\gamma$ used to construct the indirect prediction sets is an informative prior with strong precision in that the scale used is equal to the sample size in this case. As a result, the posterior distributions contain notably more information than what is in each simulated dataset. As a result, the ordering of categories induced by the posterior counts, used to construct the indirect prediction sets, are accurately approximating the oracle ordering of categories. This is evidenced by the nearly identical behaviour of the two cardinality ratios explored. In conjunction with the instability of the direct method in the presence of such a small sample size, this results in notably smaller cardinality of the indirect set as compared to the direct set, even for relatively small total numbers of categories. At its best, the indirect prediction set is about 80% smaller than the direct set.

For a moderate sample size of *N* = 100, the prior precision used to construct the indirect prediction sets is not overwhelming as compared to the sample size, and hence the posterior counts do not approximate the oracle ordering as well as in the regime with a smaller sample size. This is evidenced by the divergence of the red and blue dashes in the middle plot of Figure 1. Still, particularly as the number of categories increases for fixed *N*, the benefit of utilizing prior information of
this type is highlighted by the decline of the cardinality ratio of the indirect to direct prediction sets (red lines). For example, in the case of $N = 100$ and $K = 150$, the indirect prediction set constructed with $\gamma$ is about 15% smaller than the direct prediction set.

A similar but less pronounced pattern is seen in the presence of a larger sample size of $N = 1,000$. For this sample size with $K \le 150$, all methods considered perform relatively similarly. However, as the number of categories increases, there is a distinct gain in prediction set precision given the input of indirect information in prediction set construction.

## 3.1 Comparison with Bayesian prediction sets

For multi-group count data such as we are considering, a standard Bayesian prediction set may be constructed from the posterior predictive distribution (see, for example, Gelman et al., 2014 Section 1.3) based on the hierarchical model given by equations (1) and (6). Specifically, the posterior predictive distribution for group $j$ is a multinomial distribution for a single trial with proportions vector $\tilde{x}_j/(N_j + \sum_{i=1}^K \gamma_i)$. Then, a Bayesian prediction set $A_\alpha^B$ may be taken to be the categories corresponding to the highest mass region of this predictive distribution such that the categories' cumulative predictive probability sum exceeds $1 - \alpha$,

$$A_{\alpha}^{B}(\mathbf{x}_{j}) = \left\{ \mathcal{Y}^{(k)} \in \mathcal{Y} : \left[ \sum_{l=1}^{K} 1(\tilde{x}_{j,l} \ge \tilde{x}_{j,k}) \frac{\tilde{x}_{j,l}}{N_j + \sum_{i=1}^{K} \gamma_i} \right] > \alpha \right\} \quad (9)$$

Such a prediction set may have smaller cardinality than the indirect or direct approach, but it will not be $\alpha$-valid. In fact, the frequentist coverage may arbitrarily be nearly zero under a misspecified prior parameter with a large scale.

To illustrate the relationship between cardinality and coverage, we compare the Bayesian prediction method to the proposed indirect approach (7) in simulation. We follow the same simulation set-up as before and, for both methods, consider a prior based on full information with moderate prior precision $\gamma = \theta \times 10$. Results for a moderate sample size of $N = 100$ are plotted in Figure 2. By construction, the indirect approach is $\alpha$-valid and maintains frequentist coverage rates at or above the nominal level. In contrast, the Bayesian prediction sets may be narrower than the indirect sets, but, correspondingly, the frequentist coverage rate of the Bayesian sets may fall below the nominal rate.

<figure><img src="image_7.png" /><figcaption>**Figure 2.** (Left) Monte Carlo approximations (±1 standard deviation) of the expected cardinality of the indirect $\alpha$-valid prediction set (indirect) and Bayesian posterior predictive set (Bayes), both constructed with prior concentration $\gamma = \theta \times 10$. (Right) Monte Carlo approximations (with 95% Clopper-Pearson intervals) of the frequentist coverage rate of each method, both constructed for an error rate of $\alpha = 0.15$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Within-county log sample size.</figcaption></figure>

# 4 Summarizing eBird species abundance data

In this section, we describe avian species abundance in North Carolina, USA from eBird data obtained from citizen-uploaded complete checklists of species observations in the first week of May 2023. Across the 99 counties, 393 unique species were identified. Some species such as the Northern Cardinal, Carolina Wren, and American Robin were identified frequently. Many others like the Northern Saw-whet Owl and the Solitary Sandpiper were rarely seen; in fact, 50% of species were seen fewer than 100 times each across the entire state. Moreover, within-county sample sizes vary drastically (Figure 3) from approximately 50,000 individual birds identified in Wake County, one of the most populous counties in NC that contains the state's capital, to only 8 in Pasquotank County, a small coastal county consisting of about 1/30th of the human population of Wake County.

As motivated in Section 1, describing such data with $\alpha$-valid prediction sets for each county provides a useful summary with unambiguous statistical interpretation. That is, with at least probability $1 - \alpha$, an unobserved bird in a given county will belong to a species contained in the specified prediction set, where the probability is taken with respect to the random sample and the prediction. Here, we demonstrate the usefulness of this approach in gaining better understanding of species abundance. Moreover, we elaborate on the benefit of utilizing indirect information in the construction of practically useful sets that are precise, particularly for counties with small within-county sample sizes.

For each county in NC, we construct an indirect prediction set based on a prior hyperparameter estimated from data in the five nearest neighbouring counties, following the procedure described in Section 2.4. The eBird data consist of independent samples collected across the state, so samples are independent across counties. As a result of this independence, finite-sample coverage of the indirect prediction approach is guaranteed. We compare the cardinality of these indirect prediction sets to that of direct prediction sets, both of which maintain at least 95% coverage for each county. The cardinality ratios of the indirect to direct prediction sets across the counties in NC are plotted in Figure 4. To highlight the impact of within-county sample size, the lower quantile sample sizes are overlaid on their respective county.

In general, the incorporation of indirect information in the construction of prediction sets results in notably smaller cardinality of the indirect prediction sets as compared with that of the direct prediction sets. Of the 99 counties in NC, indirect sets have smaller cardinality in 65, and the two approaches result in the same cardinality in 20 counties. The improvement in cardinality is particularly conspicuous in counties with small to moderate sample sizes, as evidenced by the sample sizes of counties corresponding to the smallest cardinality ratios, seen in Figure 4. Moreover, 10 counties have trivial direct sets consisting of all K species, while only two counties with the smallest within-county sample sizes, 8 and 14, have trivial indirect prediction sets. For the county with the third smallest sample size (24), the indirect prediction set only includes 80 species, or about 20% of all possible species, while the direct prediction set is the trivial set.

Overall, even in counties with larger sample sizes, it is most common for the indirect and direct prediction sets to contain a different set of species. In fact, the indirect and direct prediction sets
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Cardinality ratio of indirect to direct prediction sets. Prior hyperparameters estimated with an empirical Bayesian procedure based on five nearest neighbours for each county. The lowest quantile sample sizes are overlaid on their respective counties.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Empirical probability masses of species included in only the indirect (Prothonotary Warbler, Eastern Towhee, Summer Tanager, Red-eyed Vireo, Pine Warbler, and Chipping Sparrow), only the direct (Double-crested Cormorant, Eastern Meadowlark, Downy Woodpecker (Eastern), Eastern Kingbird, and Red-headed Woodpecker), or both (all remaining species) prediction sets, sorted by the posterior proportion (Post.Pred). The empirical proportion (MLE) and the posterior proportion based on a prior estimated from data in nearest five neighbouring counties are plotted.</figcaption></figure>

disagree for nearly every county in NC. They are equivalent for only six counties where they are not both trivial sets. Commonly, this discrepancy corresponds with smaller indirect sets, and hence highlights the benefit of inclusion of indirect information in the construction of prediction sets.

## 4.1 Order-based prediction in Robeson County

To further compare the two approaches and elucidate the role of the ordering of the species, we elaborate on the construction of indirect and direct prediction sets for Robeson County. Robeson is located near the southeastern border of NC and features a moderately small within-county sample size of 247 birds observed, with species-specific observation counts ranging from zero to ten. The two prediction sets have nearly the same cardinality but contain differences in species inclusion. Specifically, the indirect prediction set contains 33 species, and the direct set contains 32, with an overlap of 27 species.

To illustrate the role of the ordering used in the construction of α-valid prediction sets, the empirical proportions based on the observed sample [maximum likelihood estimate (MLE)] and posterior proportions (Post.Pred) are plotted in Figure 5 for the union of included species in the two sets. In the figure, the species are sorted by increasing posterior proportions. The indirect and direct sets include species based on the posterior and empirical distributions, respectively. Discrepancies between the indirect and direct sets occur when these two distributions disagree.
<table><caption>**Table 1.** Percentage of all birds observed within each respective county, for select species included in either the indirect set (Pine Warbler and C. Sparrow) or the direct set (D. Cormorant and E. Kingbird)</caption><thead><tr><th></th><th>D. Cormorant</th><th>E. Kingbird</th><th>Pine Warbler</th><th>C. Sparrow</th></tr></thead><tbody><tr><td>Robeson</td><td>0.81%</td><td>0.4%</td><td>0.00%</td><td>0.00%</td></tr><tr><td>NC-017</td><td>0.00%</td><td>2.85%</td><td>1.97%</td><td>2.63%</td></tr><tr><td>NC-047</td><td>0.00%</td><td>0.00%</td><td>1.29%</td><td>0.64%</td></tr><tr><td>NC-051</td><td>0.00%</td><td>0.68%</td><td>2.62%</td><td>5.59%</td></tr><tr><td>NC-093</td><td>0.00%</td><td>0.00%</td><td>1.09%</td><td>0.55%</td></tr><tr><td>NC-165</td><td>0.00%</td><td>0.86%</td><td>2.58%</td><td>5.16%</td></tr><tr><td>γ</td><td>0.00</td><td>1.33</td><td>3.42</td><td>4.69</td></tr></tbody></table>

Note. Estimated prior hyperparameter $\gamma$ for Robeson County recorded in the last row.

From Figure 5, it is easy to see the indirect prediction set consists of the species with the 33 largest posterior predictive proportions. In contrast, the direct set consists of species with the largest sample probability mass. Naturally, the ordering of these two estimates agree for species common to the region, and, as such, there is a fair amount of overlap of species inclusion.

As a result of our estimation procedure for the prior hyperparameter $\gamma$ for Robeson County, the disparity between inclusion or exclusion of a species among the two prediction set methods is further elucidated by examining species presence in neighbouring counties. In short, species with more frequent occurrence in neighbouring counties will have a larger estimated prior count than those seen rarely in neighbouring counties. Species occurrences in neighbouring counties are displayed in Table 1 for a select few species along with the estimated $\gamma$ for Robeson County, obtained by solving equation (8) using data in these neighbouring counties.

Intuitively, species that are seen in neighbouring counties with some relative frequency, such as the Chipping Sparrow or Pine Warbler, are probably also present in Robeson County, and hence should be included in a prediction set. In practice, these species have a comparatively high estimated prior of about 5 and 4, respectively, and hence are included in the indirect prediction set even though they were not recorded as being observed in Robeson County in the dataset. Alternatively, consideration of indirect information yields the conclusion that species like the Eastern Kingbird and Cormorant may be rare in the area in general, as reflected by small $\gamma$ values, and thus these species are not included in the indirect prediction set.

## 4.2 Inference among species with tied observed counts in Haywood County

In species abundance data, particularly for areas or counties with small sample sizes, it is common for multiple species to have the same observed count. A feature of the construction of the direct order-based prediction approach as presented is that species with the same observed counts will either be jointly included or excluded from the prediction set. As a result, a direct prediction set constructed from a sample with tied species counts may have increased cardinality over an indirect prediction set that does not necessarily jointly admit all species with tied observed counts. If the direct set has increased cardinality for this reason, the direct set will also have increased coverage over the indirect set.

When constructing a prediction set based on the empirical proportions without consideration of indirect information, as in the construction of the direct set, this may commonly occur, and there is no clear approach to choose among the species with tied counts without further information than what is provided in the sample in that county. One could randomly choose to include one of the species from the set of species with tied counts, for example, but a more principled manner is to utilize indirect information to determine which species should be included. This is the mechanism used by the indirect prediction approach when the prior hyperparameter is a real-valued vector estimated from indirect information. As such, a more nuanced benefit of utilizing indirect information in the construction of a prediction set is the capacity to include a select few categories with tied empirical proportions.
<table><caption>**Table 2.** Percentage of all birds observed within each respective county for species included in either both prediction sets (Bobolink) or only the direct set (L. Flycatcher, R. Hawk, C. Yellowthroat, and E. Kingbird)</caption><thead><tr><th></th><th>L. Flycatcher</th><th>R. Hawk</th><th>C. Yellowthroat</th><th>E. Kingbird</th><th>Bobolink</th></tr></thead><tbody><tr><td>Haywood</td><td>0.21%</td><td>0.24%</td><td>0.21%</td><td>0.24%</td><td>0.24%</td></tr><tr><td>NC-021</td><td>0.06%</td><td>0.29%</td><td>0.18%</td><td>0.43%</td><td>0.14%</td></tr><tr><td>NC-099</td><td>0.00%</td><td>0.08%</td><td>0.16%</td><td>0.00%</td><td>0.00%</td></tr><tr><td>NC-115</td><td>0.07%</td><td>0.14%</td><td>0.21%</td><td>0.28%</td><td>0.00%</td></tr><tr><td>NC-173</td><td>0.18%</td><td>0.18%</td><td>0.66%</td><td>0.09%</td><td>1.41%</td></tr><tr><td>NC-175</td><td>0.00%</td><td>0.30%</td><td>1.00%</td><td>0.54%</td><td>0.84%</td></tr><tr><td>γ</td><td>0.4</td><td>1.29</td><td>2.43</td><td>1.49</td><td>2.15</td></tr></tbody></table>

Note. Estimated prior hyperparameter $\gamma$ for Haywood County recorded in the last row. Species are sorted by posterior proportion.

To demonstrate, we elaborate on species inclusion in the indirect and direct prediction sets in Haywood County. Haywood is popular destination in the Blue Ridge Mountains, located near the western border of North Carolina. It features a moderately large within-county sample size of roughly 4,000 birds observed. In Haywood County, the indirect prediction set contains 70 species, and the larger direct set contains 74. In the construction of these prediction sets, the ordering of species with regards to the posterior proportions and the empirical proportions agree for most species. As a result, all 70 species included in the indirect set are also included in the direct set. The disparity in species inclusion occurs primarily as a result of tied counts of species occurrence in the sample.

Empirical proportions in Haywood and neighbouring counties are reported in Table 2 for the five species included in Haywood County's prediction sets with the smallest posterior proportions. The species with the four smallest posterior proportions are included only in the direct set, and the other species, the Bobolink, is included in both the indirect and direct sets. The Bobolink was observed nine times in the sample from Haywood County, or about 0.24% of the Haywood sample. For an ordering determined by either the empirical counts or the posterior counts, this species is required to be included in the order-based prediction set to guarantee $1 - \alpha$ coverage. Two of the other species, the Red-shouldered Hawk and Eastern Kingbird, were each also observed nine times in the sample from Haywood, and, by construction of the order-based prediction approach, must also be included in the direct set. When admitting the species into a prediction set by posterior counts based on the real-valued prior hyperparameter $\gamma$ estimated from data in neighbouring counties, as in the indirect approach considered, the 'tie' among these three species is broken, and only one, the Bobolink, is included in the indirect prediction set.

## 4.3 Comparison with Bayesian prediction sets

The $\alpha$-valid prediction approaches may also be compared to the Bayesian prediction method detailed in equation (9). As the Bayesian method and the indirect method are constructed based on the same hierarchical working model, the Dirichlet prior parameter in (9) may be estimated following the same empirical Bayes approach used for the indirect prediction method. While the Bayesian sets may be smaller than the direct or indirect sets, the Bayesian sets are not $\alpha$-valid, and thus the frequentist coverage of a set for any given county may fall below the nominal level.

As conveyed in the cardinality comparison between the indirect and direct prediction sets, incorporating auxiliary information in the construction of the prediction sets results in improved cardinality overall. Specifically, the Bayesian prediction sets are smaller than the direct sets in 59 of the 99 counties in North Carolina. In comparing the two methods that utilize auxiliary information, the Bayesian sets are smaller than indirect sets in 33 counties, and they have the same cardinality in 15 counties. Overall, the indirect method outperforms both the Bayesian method and the direct method in terms of minimizing prediction set cardinality. Moreover, the Bayesian
approach is not α-valid, and thus the frequentist coverage rates of the Bayesian prediction sets may fall below the nominal level.

# 5 Discussion

Species abundance data collected across heterogeneous areas is increasingly important in understanding biodiversity. Some of the largest sources of such data are citizen science databases for which volunteers spearhead the data collection. As a result of the civilian-led scientific effort, such data often feature unequal sampling across a spatial domain where some areas have large within-area sample sizes and others have much smaller within-area sample sizes.

In this article, we propose summarizing species abundance data of this type with valid prediction sets that are constructed by sharing information across areas. Utilizing indirect information may result in smaller prediction sets than otherwise achievable with direct methods. Meanwhile, maintaining validity of the prediction sets for each area allows for an accessible interpretation that enables a straightforward comparison across areas. In particular, maintaining interpretable statistical guarantees on a descriptor of such data is important as analyses from such data often have far reaching policy implications. Smaller prediction sets may be attainable based on Bayesian inference of a spatial hierarchical model such as that presented in Tang et al. (2023), for example, but these approaches introduce bias and a resulting prediction set would not retain the nominal frequentist coverage rate guarantee for each county.

The usefulness of our approach for summarizing citizen science data is motivated in part to combat the common problem of varying sampling efforts across areas. We detail how α-valid prediction sets can be constructed with the incorporation of indirect information to improve within-county prediction set precision and propose an empirical Bayes procedure to do so. Incorporation of accurate indirect information results in a narrower prediction set for a given county than a direct prediction set by exploiting data in nearest neighbouring counties. The proposed empirical Bayes procedure is based on a standard hierarchical model that is straightforward to understand, and the authors provide code for implementation.

There may, however, be a benefit to utilizing a more structured prior that incorporates indirect information in a more complex manner such as a prior that weights data from different parts of the state differently. For example, a model based on a learned intrinsic distance between counties was shown in Christensen and Hoff (2022) to fit a subset of the eBird data better than standard methods based on geographic adjacency structure. In the sample analysed in Section 4, we found an indirect prediction set constructed with a hyperparameter estimated from five nearest neighbours results in overall narrower prediction sets than a direct approach, but it would be valuable to explore if this can be further improved upon with a more detailed prior. More broadly, different applications may warrant an alternative information sharing prior if, for example, there is no notion of spatial distance across the different areas. For example, it may be of interest to compare species abundance variation across different time frames for a given county.

All replication codes for this article, including functions to implement the empirical Bayes estimation procedure for the prior hyperparameter, are available at https://github.com/betsybersson/FreqPredSets_Indirect.

*Conflicts of interest:* No competing interest is declared.

# Funding

None to declare.

# Data availability

The eBird project is managed by the Cornell Lab of Ornithology, and the raw data are publicly available online at https://science.ebird.org/en/use-ebird-data. The sample analysed in Section 4 is available at https://github.com/betsybersson/FreqPredSets_Indirect/tree/main/data.
## Author contributions statement

E.B. and P.H. conceived the experiments, E.B. conducted the experiments, E.B. and P.H. analysed the results, E.B. and P.H. wrote and reviewed the manuscript.

## Appendix A: Proofs

**Remark 1** (Concerning Theorem 1). We first elaborate on the construction of an order-based prediction set following equation (3). To test if an element $y^{(k)}$ in the sample space $\mathcal{Y}$ is included in a prediction set for a given vector $o$ and known event probability vector $\theta$, the cumulative sum of event probabilities of the categories corresponding to the minimum element of $o$ up to element $k$, following the ordering of $o$, is computed. If this cumulative sum is greater than the error rate $\alpha$, then element $k$ is included in the prediction set. As a result, all elements with such cumulative sums greater than $\alpha$ are included in the prediction set. The elements with such cumulative sums less than or equal to $\alpha$ are not included. Therefore, by construction, $P(Y \notin A_{\alpha}^{\theta,o} | \theta) \le \alpha$. Consequently,

$$P(Y \in A_{\alpha}^{\theta,o} | \theta) = 1 - P(Y \notin A_{\alpha}^{\theta,o} | \theta) \ge 1 - \alpha,$$

so $A_{\alpha}^{\theta,o}$ is $\alpha$-valid.

**Proof of Theorem 1.** (1): Note first that $o$ enters equation (3) only through the ordering of its elements. As such, without loss of generality, consider vectors of the form $o \in \{0, 1\}^K$. When constructing a prediction set following equation (3) based on such a vector $o$, the category space is effectively divided into two disjoint subsets,

$$\mathcal{Y}_0 = \{y^{(k)} \in \mathcal{Y} : o_k = 0\}$$

$$\mathcal{Y}_1 = \{y^{(k)} \in \mathcal{Y} : o_k = 1\},$$

such that, by construction,

$$A_{\alpha}^{\theta,o} = \begin{cases} \mathcal{Y} & \text{if } \sum_{j: y^{(j)} \in \mathcal{Y}_0} \theta_j > \alpha \\ \mathcal{Y}_1 & \text{else} \end{cases}.$$

For a given error rate $\alpha$, clearly any $\alpha$-valid prediction set may be constructed under considerations of permutations of the vector $o$ of the form $o \in \{0, 1\}^K$.

(2): We wish to show no other ordering results in an $\alpha$-valid prediction set with strictly smaller cardinality than $o^{\theta}$. In words, consider switching the ordering of one element at a time. This will always result in a prediction set with the same cardinality or greater cardinality than that under $o^{\theta}$. More formally, let $\tilde{o} = \{o_1^{\theta}, o_2^{\theta}, \dots, o_{k^*}^{\theta}, o_{k^*-1}^{\theta}, o_{k^*+1}^{\theta}, \dots, o_K^{\theta}\}$, that is, equivalent to $o^{\theta}$ with the $(k^*)$th and $(k^*-1)$th ordering flipped. But, by construction, $\theta_{k^*} \ge \theta_{k^*-1}$, so $|A_{\alpha}^{\theta,o^{\theta}}| \le |A_{\alpha}^{\theta,\tilde{o}}|$. $\square$

**Proof of Theorem 2.** Let $Y_1, \dots, Y_{N+1} \sim \text{i.i.d. MN}_K(\theta, 1)$. Then, we wish to construct a conformal prediction set for $Y_{N+1}$ based on an observation of $X = \sum_{i=1}^N Y_i$, and some conformity measure C.

To determine if a candidate category $k \in \{1, \dots, K\}$ is included in a $1-\alpha$ conformal prediction set, the conformal algorithm proceeds
as follows, see Section 2.1 of Bersson and Hoff (2024) for more details:

1. Set $y_{N+1} = y^{(k)}$ where $y^{(k)}$ is a vector of length K with a 1 in the kth index and 0s elsewhere.

2. For $j = 1, ..., N+1$, compute conformity scores $c_j = C(x - y_j + y^{(k)}, y_j)$.

3. Set

$$p_k = \frac{\{\#j \in \{1, ..., N+1\} : c_{N+1} \ge c_j\}}{N+1}$$

More compactly, and by symmetry in the problem, this conformal p-value may be equivalently written as:

$$p_k = \sum_{l=1}^{K} \mathbb{1}(o_k^* \ge o_l^*) \frac{x_l + y_l^{(k)}}{N+1},$$

where $o_k^* = c_{N+1}$ and $o_l^* = c_j$ for a $j \in \{1, ..., N\}$ such that $y_j = y^{(l)}$. Then, for prediction mis-coverage rate $\alpha$, the category k is included in the prediction set if $p_k > \alpha$. A prediction set constructed from this procedure may be concisely written as follows:

$$A_{\alpha}(\mathbf{x}) = \left\{ y^{(k)} \in \mathcal{Y} : \left[ \sum_{l=1}^{K} \mathbb{1}(o_k \ge o_l) \frac{x_l + y_l^{(k)}}{N+1} \right] > \alpha \right\},$$

for some $o \in \mathbb{R}^K$.

□

## Appendix B: Maximization of the Marginal Multinomial-Dirichlet Likelihood

In this section, we detail a Newton-Raphson algorithm to maximize the log marginal likelihood of a conjugate multinomial-Dirichlet model, sometimes referred to as the Dirichlet-multinomial compound distribution:

$$X_j \sim MN_K(\theta_j, N_j), \text{ independently for } j = 1, ..., J$$
$$\theta_1, ..., \theta_J \sim \text{Dirichlet}_K(\gamma).$$

The log likelihood of the marginal likelihood is as follows:

$$\mathcal{L}(\gamma) \propto \sum_{j=1}^{J} \left[ \log \Gamma \left( \sum_{i=1}^{K} \gamma_i \right) - \log \Gamma \left( N_j + \sum_{i=1}^{K} \gamma_i \right) + \sum_{i=1}^{K} \log \Gamma(x_{i,i} + \gamma_i) - \sum_{i=1}^{K} \log \Gamma(\gamma_i) \right].$$
Define

$$ \Psi(s) = \frac{d}{ds} \log \Gamma(s) = -\xi + \sum_{n=0}^{\infty} \left[ \frac{1}{n+1} - \frac{1}{n+s} \right], $$

where $\xi$ is the Euler-Mascheroni constant. Then, it is straightforward to obtain the first and second derivatives of the marginal log likelihood,

$$ \begin{aligned} \frac{d}{d\gamma_k'} &= \sum_{j=1}^{J} \left[ \Psi\left(\sum_{i=1}^{K} \gamma_i\right) - \Psi\left(N_j + \sum_{i=1}^{K} \gamma_i\right) + \Psi(x_{j,k} + \gamma_k) - \Psi(\gamma_k) \right] \\ \frac{d}{d\gamma_k^2} &= \sum_{j=1}^{J} \left[ \Psi'\left(\sum_{i=1}^{K} \gamma_i\right) - \Psi'\left(N_j + \sum_{i=1}^{K} \gamma_i\right) + \Psi'(x_{j,k} + \gamma_k) - \Psi'(\gamma_k) \right] \\ \frac{d}{d\gamma_k d\gamma_k'} &= \sum_{j=1}^{J} \left[ \Psi'\left(\sum_{i=1}^{K} \gamma_i\right) - \Psi'(N_j + \sum_{i=1}^{K} \gamma_i) \right], \end{aligned} $$

where $\Psi'$ is the trigamma function. Let $g$ be the gradient vector of length $K$ and $H$ the Hessian matrix. Finally, Newton's method updates $\gamma$ as follows:

$$ \gamma^{(t+1)} = \gamma^{(t)} - H^{-1}(\gamma^{(t)})g(\gamma^{(t)}), $$

where the algorithm is iterated until convergence.

Solving for the MLE of the unknown parameter in the Dirichlet-multinomial compound distribution is covered in literature dating back to at least (Mosimann, 1962), and it is stated in Wallach (2008) that the likelihood is concave. For a detailed overview of the Dirichlet-multinomial compound distribution, see Ng et al. (2011). For derivations and comparisons of various algorithms to solve the maximization problem, including computational expense and initialization sensitivity, we refer the reader to Minka (2000) and Wallach (2008).

## References

* Arnold Z. J., Wenger S. J., & Hall R. J. (2021). Not just trash birds: Quantifying avian diversity at landfills using community science data. *PLoS One*, 16(9), 1–14. https://doi.org/10.1371/journal.pone.0255391

* Bersson E., & Hoff P. D. (2024). Optimal conformal prediction for small areas. *Journal of Survey Statistics and Methodology*, *n/a*(2325-0992), smae010. https://doi.org/10.1093/jssam/smae010

* Camerini G., & Groppali R. (2014). Landfill restoration and biodiversity: A case of study in Northern Italy. *Waste Management and Research*, 32(8), 782–790. https://doi.org/10.1177/0734242X14545372

* Christensen M. F., & Hoff P. D. (2022). A flexible and interpretable spatial covariance model for data on graphs. *Environmetrics*, *n/a*(n/a), e2879. https://doi.org/10.1002/env.2879

* Faulkenberry G. D. (1973). A method of obtaining prediction intervals. *Journal of the American Statistical Association*, 68(342), 433–435. https://doi.org/10.1080/01621459.1973.10482450

* Gelman A., Carlin J. B., Stern H. S., Dunson D. B., Vehtari A., & Rubin D. B. (2014). *Bayesian data analysis* (3rd ed.). CRC Press, Taylor & Francis Group.

* Hoff P. (2023, May). Bayes-optimal prediction with frequentist coverage control. *Bernoulli*, 29(2), 901–928. https://doi.org/10.3150/22-BEJ1484

* Lebrun J. J., Thogmartin W. E., & Miller J. R. (2012, August). Evaluating the ability of regional models to predict local avian abundance. *Journal of Wildlife Management*, 76(6), 1177–1187. https://doi.org/10.1002/jwmg.374

* Lele S. R. (2020). How should we quantify uncertainty in statistical inference? *Frontiers in Ecology and Evolution*, 8, 1–18. https://doi.org/10.3389/fevo.2020.00035

* Minka T. P. (2000). *Estimating a Dirichlet distribution* (Technical report).

* Mosimann J. E. (1962). On the compound multinomial distribution, the multivariate beta-distribution, and correlations among proportions. *Biometrika*, 49(2), 65. https://doi.org/10.2307/2333468
Ng K. W., Tian G.-L., & Tang M.-L. (2011, May). Dirichlet-multinomial distribution. In *Dirichlet and related distributions: Theory, methods and applications* (1st ed., Vol. 895, chapter 6, pp. 199-225). John Wiley & Sons, Incorporated.

Rao J. N. K., & Molina I. (2015). *Small area estimation* (2nd ed.). John Wiley and Sons, Inc.

Shanahan D. F., & Possingham H. P. (2009, October). Predicting avian patch occupancy in a fragmented landscape: Do we know more than we think? *Journal of Applied Ecology*, *46*(5), 1026-1035. https://doi.org/10.1111/j.1365-2664.2009.01694.x

Sullivan B., Wood C., Iliff M., Bonney R., FInk D., & Kelling S. (2009). eBird: A citizen-based bird observation network in the biological sciences. *Biological Conservation*, *142*(10), 2282-2292. https://doi.org/10.1016/j.biocon.2009.05.006

Sullivan B. L., Phillips T., Dayer A. A., Wood C. L., Farnsworth A., Iliff M. J., Davies I. J., Wiggins A., Fink D., Hochachka W. M., Rodewald A. D., Rosenberg K. V., Bonney R., & Kelling S. (2017, April). Using open access observational data for conservation action: A case study for birds. *Biological Conservation*, *208*, 5-14. https://doi.org/10.1016/j.biocon.2016.04.031

Tang B., Clark J. S., Marra P. P., & Gelfand A. E. (2023, March). Modeling community dynamics through environmental effects, species interactions and movement. *Journal of Agricultural, Biological, and Environmental Statistics*, *28*(1), 178-195. https://doi.org/10.1007/s13253-022-00520-3

Thatcher A. R. (1964). Relationships between Bayesian and confidence limits for predictions. *Journal of the Royal Statistical Society. Series B (Methodological)*, *26*(2), 176-210. https://doi.org/10.1111/j.2517-6161.1964.tb00551.x

Tian Q., Nordman D. J., & Meeker W. Q. (2022, November). Methods to compute prediction intervals: A review and new results. *Statistical Science*, *37*(4), 580-597. https://doi.org/10.1214/21-STS842

Twedt D. J., Tirpak J. M., Jones-Farrand D. T., Thompson F. R., Uihlein W. B., & Fitzgerald J. A. (2010, August). Change in avian abundance predicted from regional forest inventory data. *Forest Ecology and Management*, *260*(7), 1241-1250. https://doi.org/10.1016/j.foreco.2010.07.027

Vovk V., Gammerman A., & Shafer G. (2005). *Algorithmic learning in a random world*. Springer.

Wallach H. M. (2008). *Structured topic models for language*. University of Cambridge.
# Extreme-value modelling of migratory bird arrival dates: insights from citizen science data

Jonathan Koh<fnref n="1" /> ID and Thomas Opitz<fnref n="2" /> ID

*Address for correspondence:* Jonathan Koh, Institute of Mathematical Statistics and Actuarial Science, Oeschger Centre for Climate Change Research, University of Bern, Alpeneggstrasse 22, 3012 Bern, Switzerland. Email: jonathan.koh@stat.math.ethz.ch

[To be read before The Royal Statistical Society at the Discussion Meeting on the Analysis of citizen science data' held at the Societys 2024 annual conference in Brighton on Tuesday, 3 September 2024, the President, Dr Andrew Garrett, in the Chair]

## Abstract

Citizen science mobilizes many observers and gathers huge datasets but often without strict sampling protocols, resulting in observation biases due to heterogeneous sampling effort, which can lead to biased predictions. We develop a spatio-temporal Bayesian hierarchical model for bias-corrected estimation of arrival dates of the first migratory bird individuals at their breeding sites. Higher sampling effort could be correlated with earlier observed dates. We implement data fusion of two citizen-science datasets with fundamentally different protocols (Breeding Bird Survey, eBird) and obtain posterior distributions of the latent process, which contains four spatial components endowed with Gaussian process priors: species niche; sampling effort; position and scale parameters of annual first arrival date. The data layer consists of 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 and find that the sampling effort significantly modulates the observed first arrival dates. We exploit this relationship to effectively bias-correct predictions of the true first arrivals.

**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, the field involves a collaborative approach to scientific inquiry that engages volunteers and nonprofessionals in the collection, analysis, and interpretation of data. This participatory model is not new. Friedrich Wilhelm Herschel, a musician by training, discovered the planet Uranus in 1781 using his telescopes. Charles Darwin conducted crowd-sourcing projects in the 19th century by recruiting acquaintances and travellers to write to him about their observations (Browne, 1996). The longest-running citizen science initiative is the annual Christmas Bird Count (Bock & Root, 1981), a census administered by the National Audubon Society of wintering birds in the Western Hemisphere by volunteer birdwatchers since 1900.

<footnote n="1">¹Institute of Mathematical Statistics and Actuarial Science, Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland</footnote>

<footnote n="2">²Biostatistics and Spatial Processes (UR546), INRAE, Avignon, France</footnote>
Citizen science has gained significant traction over the past two decades, spurred on by growing world literacy levels and the advent of technological devices and infrastructures 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 unconventional funding sources (Silvertown, 2009). Zooniverse,<fnref n="1" /> the world's largest platform for volunteer-based research, has seen a surge to over 1.9 Million registered volunteers 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, Cooper, et al., 2016) dedicated to disseminating research on citizen science, and the 2019 Citizen Science Association biannual conference attracted 818 registered delegates from 28 countries.

Citizen science is expected to grow in importance. Technological progress continually improves how machines and citizen science 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), and 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, giving them a voice in local environmental decision-making (Bonney, Phillips, et al., 2016), and increasing public awareness of the scientific process. Though citizen science's global economic value has been estimated to exceed 2.5 billion USD annually, the majority of data collected through citizen-science initiatives has not yet reached analysis in peer-reviewed literature (Theobald et al., 2015). There is a need to develop and disseminate statistical methods that facilitate wider scientific use of these datasets.

## 1.2 Analysis of data for species monitoring

Much of the recent growth in citizen science has occurred in the ecological and environmental sciences (Fraisl et al., 2022; McKinley et al., 2015; Pocock et al., 2018). A prominent example in ecology is the eBird project (Cornell Lab of Ornithology, 2022), launched in 2002, which has led to a database of over 200 Gigabytes providing information on over 600 million bird observations. Another important bird monitoring program within the scope of citizen science is the North American Breeding Bird Survey (BBS). It was launched in 1966 for official monitoring purposes and engages large numbers of trained volunteers in collecting standardized data on bird populations according to a rather strict sampling protocol along several selected routes and during each year. BBS data can be considered of high quality but sampling is relatively sparse in space and occurs only during fixed periods in the year, whereas the sampling of eBird data is much more heterogeneous but provides a generally denser spatio-temporal coverage, as we will detail in Section 2. We exploit these two data sources in this work. Regarding nonavian species, another recent example, among many others, is the deployment of camera traps for mammal monitoring (Hsing et al., 2022).

Many observation data in ecology are opportunistic, i.e. they were collected incidentally or without a predefined research question in mind, so information on the criteria applied by the observer for sampling and reporting observations is limited. However, opportunistic data often allow for a spatio-temporal coverage of species monitoring that would not be attainable with protocol-based data collected by professionals at relatively sparse and preselected sampling locations. Citizen-science data are particularly valuable for statistical inference on low-probability events (e.g. the presence of rare species individuals, or unusual or extreme phenological events) and on occurrence times of events, such as the arrival of migratory birds at their breeding site in spring. The general data quality from citizen-science initiatives is high (Kosmala et al., 2016), although 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 on the sampling effort at the data collection points could be captured. They further identify four key dimensions of biases induced by volunteer sampling: uneven sampling in time and space, uneven effort per visit, and uneven detectability of species individuals, which can vary by observer, species, and land cover. For example, detectability could be higher in open landscapes than in dense forests. Moreover, the expertise of individual observers for detecting and identifying species may increase over time due to

<footnote n="1">¹ https://www.zooniverse.org/</footnote>
increasing training, experience, and knowledge-sharing with other observers, thus reducing observation biases (Johnston et al., 2018; Kelling et al., 2015).

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 to characterize heterogeneous sampling effort and quantify it from available data (e.g. Fink et al., 2020; Gelfand & Shirota, 2019; Johnston et al., 2021, 2023).

Another problem arises when only the presence of species is reported in a nonexhaustive way but not their absence (*presence-only data*). Then, areas and times without any reported occurrences could correspond either to the true absence of the species or to the presence of the species but the absence of any observer. 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 be reported, leading to exhaustive sampling for detected and identified bird individuals. In eBird, data entries satisfying this protocol are known as *checklist data*, with a unique space-time coordinate associated with each checklist; these account for the vast majority of eBird entries. Checklist data in eBird are quality-controlled by experts, and nonexpert observers can train in birding with resources offered by eBird, such as tutorials, online courses and bird identification apps. Observer skill for detecting and correctly identifying various bird species still varies across observers, even for experts. Due to the large number of observers contributing to the eBird dataset, our model will capture an observation effort that represents the average observation skill.

We adopt the common approach of considering 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 the presence and 0 corresponding to (pseudo-) absence of the species. We can then estimate the species presence probability using available predictors with a binomial likelihood.

## 1.3 Extreme-value analysis for migratory bird arrivals

Extreme-value theory (EVT, Coles, 2001) is a branch of probability theory and statistics that originated in the first half of the 20th century. It 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 & Denny, 1993; Katz et al., 2005) which can have a strong influence on biotic processes, for example when very low winter temperatures contribute to limiting the extent of outbreaks of forest pests (Thibaud et al., 2016). However, direct applications of EVT to biotic variables in ecology are scarce, due perhaps to the strongly discretized 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. Regarding human biology, some extreme-value studies exist and have considered extreme sports performances (e.g. for professional swimming competitions, Spearing et al., 2021) and longevity (Belzile et al., 2022). Recent ecological 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., 2024). However, the aforementioned approaches do not explicitly consider the specific observation biases inherent to ecological datasets, which is the focus of our work.

The first arrivals of migratory birds at their destination for a given year are an example of phenological events, i.e. of recurring life history events of species, such as migration or breeding. Bird migration is a complex phenomenon taking place over large geographic scales (Somveille et al., 2015), with possible temporal trends in migration patterns resulting from changes in climate, land cover and other ecological and evolutionary factors (Conklin et al., 2021; Cotton, 2003). The study of migration patterns, and particularly of arrival times at the breeding site, has attracted strong interest (e.g. 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 used rather complex approaches to define and estimate a date that can be viewed as representative of the arrival of birds at their breeding site. Throughout this paper, we
<figure><img src="image_11.png" /><figcaption>**Figure 1.** Scatterplots of the logarithm of positive checklist counts (x-axes) and first arrivals (y-axes, with later dates corresponding to lower values) observed at each pixel-year combination for *Chimney Swift* (left) and *Chestnut-sided Warbler* (right). The larger dots show 20 binned estimates, the black line a smooth fitted curve, and the shaded region the 95% pointwise confidence intervals of the curve.</figcaption></figure>

will simply write 'first arrival' to refer to the calendar date of first arrival. Youngflesh et al. (2021) based the representative date for a given area and year on both the first arrival of an individual bird and the first local maximum of the species' detection probabilities estimated from eBird data. In our approach, we focus only on the first arrival for each pixel of a spatial 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 (GEV) distribution motivated by EVT to model this sample extreme of all observed occurrence times in the year, similar to Wijeyakulasuriya et al. (2024). We treat dates as a continuous variable, which is sensible since the distribution of observed first arrivals spans over a long enough period so that the discretization effects inherent to dates have a weak influence; Figure 1 suggests that observed dates span over several months. Although distributions for discrete extremes have been proposed in the literature (e.g. Hitz et al., 2024; Prieto et al., 2014; Ranjbar et al., 2022), they usually come with overhead for numerical computation and modelling. Moreover, asymptotic theory is limited to heavy-tailed variables, by contrast with our setting where tails have a natural finite bound and are therefore light-tailed.

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

To our knowledge, the modelling approach developed here is the first that aims to appropriately capture the interplay of the niche, the sampling effort and the observed first arrivals to provide bias-corrected predictions of the true first arrivals. Wijeyakulasuriya et al. (2024) used only eBird data and focused on modelling the spatial dependence among the first arrivals for each year using the class of max-infinitely divisible models, so that first arrivals at locations without ob- servations during certain years can be predicted. We extend their approach and consider the niche and the sampling effort when predicting first arrivals, and we focus on revealing spatial patterns that remain stable across the whole study period. 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 a relatively high overall sampling effort in eBird.

## 1.4 Bayesian hierarchical modelling of complex ecological data

Observation data in ecology often do not directly measure the latent (i.e. not directly observed) processes of interest, such as the ecological niche, the sampling effort and the timing of a pheno- logical 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
a concrete illustration. The eBird checklists are available for spatio-temporal locations chosen by the observers and provide generally good spatio-temporal 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 around 40 km long, 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 pre-specified spatio-temporal 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 provide no direct information about events taking place outside the prescribed observation time interval, such as first arrivals. On the other hand, within BBS many routes have been sampled during each year since the 1970s, so the temporal coverage of BBS over the full study period is more homogeneous and complete than eBird. 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 be generated conditional on latent processes, which in turn are conditioned on hyperparameters (e.g. the variance or the spatial correlation 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 & Shirota, 2019) using so-called integrated species distribution models. An often 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 the complex BHM, we design for mapping the niche, observation effort, and first arrivals, we devise a Markov chain Monte Carlo (MCMC) scheme with Vecchia likelihood approximations (Katzfuss & Guinness, 2021; Vecchia, 1988) and Metropolis-adjusted Langevin algorithm (MALA, Roberts & Rosenthal, 2001) updates for the latent Gaussian fields. Our approach marks the first use of the Vecchia approximation in Bayesian hierarchical modelling for spatial extremes, extending other recent approaches in EVT (Huser et al., 2023; Majumder et al., 2024).

## 1.5 Outline of the paper

Section 2 presents the datasets and extensive preprocessing steps. In Section 3, we recall the basics of EVT and introduce the BHM, we use to identify the spatial variability of 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 MCMC inference are presented in Section 4, together with a simulation study that confirms the statistical identifiability of the model components. Results for estimated first arrivals and other model components are presented for several species in Section 5. A discussion and an outlook towards future research related to our approach and to citizen-science data more generally concludes the paper in Section 6.

# 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 nine following states: Connecticut, Maine, Massachusetts, New Hampshire, New Jersey, New York, Pennsylvania, Rhode Island, and Vermont. We merged the extracted data into a single dataset. We used the full observation record of the eBird Basic Dataset, where each entry reports the observation of one or several individuals of a species, and retained the following attributes: name of species, count of observed individuals, longitude, latitude, date, duration of observation (i.e. what we consider as 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.

The first arrivals of migratory birds in the study area are known to take place at the earliest towards the end of March. We identify migratory species as having no observations during the winter months of December and January in eBird. For the migratory species characterized by this property, we remove a few February occurrences that are understood as carryovers from the preceding autumn season. Only migratory species that had at least 200 occurrences in eBird for 2021 and 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 10 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 did not nest or breed 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 eBird data aggregated to a regular pixel grid of 20 km width. This pixel size is similar to the spatial mapping unit used for most official eBird communications and provides a good compromise between relatively noisy data behaviour at very small scales and the intended modelling of spatial variability at relatively small scales. Our response variables are:

* *BBS species counts:* 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 number of available checklists, calculated for each configuration of year and pixel; we use the Poisson distribution for this variable;

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

* *eBird first arrivals:* the minimum date of observation, expressed as the number of days since 1st January, calculated for each configuration of species, year, and pixel where at least one observation of the species occurred; we use the GEV distribution for the negated variable, which therefore represents a maximum.

Some of the above datasets and their components are illustrated in Figure 3.

## 2.2 Effect of observation effort on first arrivals

Figure 1 shows observed first arrivals plotted against positive checklist counts, for each pixel-year combination and two species. Observed dates tend to occur later during the year when fewer checklists are available, i.e. when sampling effort was lower in space. This gives strong motivation for modelling the influence of sampling effort on observed dates, and then using the estimated relationship to bias-correct predictions of true dates by setting the sampling effort to a very high, saturating 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
<figure><img src="image_12.png" /><figcaption>**Figure 2.** Scatterplots of the residuals (in days) of the smooth local regression fitted to first arrivals in Figure 1 (with later dates corresponding to lower values) as a function of the logarithm (with base 10) of duration (in minutes) of observation (x-axis, for *Chimney Swift* (left) and *Chestnut-sided Warbler* (right)). To improve readability, the range of values shown for the residuals does not include all points.</figcaption></figure>

duration for the two example species in Figure 2, with relatively earlier arrivals arising for longer durations. This suggests that both the number of checklists and the duration of observation per checklist are relevant components of the overall sampling effort modulating the distribution of the observed first arrival of birds at their breeding sites, with higher effort associated with earlier observation of the first arrival.

To explain the remaining variability in the residuals in Figures 1 and 2, we conjecture that other factors related to spatial variation in bird phenology and in species presence are relevant. We construct our model to estimate spatial effects for these properties and then use them to correct predicted first arrivals.

## 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 indicates that 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 indicates that these features are relatively weak, and the pressure difference between them is smaller. Strongly negative NAO is 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 use the average NAO for March, which was already identified as a relevant climate predictor for bird migration timing by Wijeyakulasuriya et al. (2024).

We use the 2021 National Land Cover Database (NLCD) 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 the niche, observation effort and components of the first arrival distribution with land cover maps to detect and interpret significant correlations between these model components and land cover. Land cover proportions for the pixel grid are reported in the Appendix (Figure A1) for four land-cover types, defined by 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 shows a schematic summarizing the different data sources and preprocessing steps of our approach, and the overall 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 2001, we see relatively longer durations in specific spatial regions, while in 2021 the durations are more homogeneous across larger areas and also generally longer. 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.
<figure><img src="image_13.png" /><figcaption>**Figure 3.** Schematic summarizing the data sources, the data preprocessing steps and the four components in the data layer of the BHM (panels labelled with numbers 1–4). The dashed lines illustrate our modelling choices and how model components interact, but they need not 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).</figcaption></figure>

# 3 Bayesian hierarchical model

## 3.1 General setting and notations

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

We write $S \times T$ for the space-time study domain, and denote by $A_i, i = 1, \dots, D \times T$, the set of nonintersecting space-time cells based on a division of the study area $S$ into $D = 1,268$ pixels, replicated over $T = 21$ years. We identify each cell $A_i$ with a representative location $(s_i, t_i)$, such as its barycenter, where $s_i \in S$ and $t_i = 1, \dots, T$.

Let $N_i^{\text{ckl}}$ be the number of eBird checklists in $A_i$. If $N_i^{\text{ckl}} > 0$, then $N_i^{\text{spc}} \in \{0, 1, \dots, N_i^{\text{ckl}}\}$ denotes the sum of the binary indicators of detecting a given species from each of the $N_i^{\text{ckl}}$ checklists. Binary indicators equal 1 if the species was detected and 0 otherwise. Given $N_i^{\text{spc}} > 0$, we write $Y_i = (Y_{i,1}, \dots, Y_{i,N_i^{\text{spc}}}) \in [1, 365]^{N_i^{\text{spc}}}$ to denote the corresponding first arrivals in $s_i$ within the year $t_i$ from the $N_i^{\text{spc}}$ observed species occurrences. Based on these arrival dates, we calculate $Z'_i = -\min_j N_i^{\text{spc}} Y_{i,j} = \max_j N_i^{\text{spc}} (-Y_{i,j})$, the first arrival in year $t_i$ within $s_i$. We set $Z_i = -\log(-Z'_i/366)$ to reflect that $Z'_i$ has a natural lower bound at $-366$, and model $Z_i > 0$ instead.

## 3.2 Extreme-value distribution for first arrivals

Provided that the partitioning of $S$ provides enough ‘effectively independent’ arrival times within each pixel-year $A_i$, EVT (Coles, 2001) motivates using a GEV to model the transformed first
arrivals in $A_i$. Invoking EVT, we assume convergence in distribution of each $Z_i$ to a limiting GEV as $N_i^{\text{spc}} \to \infty$.

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 nondegenerate distribution after linear renormalization, then this distribution must be a GEV; this theoretical result was further extended and shown to hold for maxima over dependent data, although under technical assumptions that are difficult to formally validate in applications. In practice, for fixed $n$ large enough (where $n = N_i^{\text{spc}}$ in our setting), one uses the approximation

$$ \Pr(M_n \le z) \approx \begin{cases} \exp\{-\exp\left(-\frac{z-\mu}{\sigma}\right)\}, & \zeta = 0, \\ \exp\left\{-\left[1 + \zeta\left(\frac{z-\mu}{\sigma}\right)\right]^{-1/\zeta}\right\}, & \zeta \ne 0, \end{cases} \quad (1) $$

where the right-hand side of (1) is the distribution function of the GEV, defined on $\{z: 1 + \zeta(z - \mu)/\sigma > 0\}$ with location $\mu \in \mathbb{R}$, scale $\sigma > 0$ and shape $\zeta \in \mathbb{R}$. The shape parameter determines the tail behaviour of the GEV: for $\zeta > 0$, the GEV has support with a finite lower endpoint and an upper heavy tail of power-law form; for $\zeta < 0$, the GEV has a finite upper endpoint but is unbounded below; for $\zeta = 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, \dots, N_i^{\text{spc}}$ can be dependent, but the three-parameter GEV remains a flexible model to accommodate relevant distributional properties of the maxima $Z_i$, even in cases where the theoretical asymptotics are not reached. The transformation from $Z'_i$ to $Z_i$ is very close to linear in a large neighbourhood of $Z'_i = 100$ (corresponding to approximately mid-April), so it does not substantially modify the shape $\zeta$ but the parameters $\mu$ and $\sigma$. We incorporate linear predictors into the GEV regression equation used to model $Z_i$ in our BHM: one for the location $\mu$ and another for the scale $\sigma$.

## 3.3 Data layer of the BHM

We construct a system of four regression equations with latent Gaussian processes used in the linear predictors. 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_i^{\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 negative first arrivals 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 $A_i$ of the pixels (i.e. we use an offset of $\log A_i$ in the linear predictor), since some pixels at the boundary of the study region have smaller area. For the response variable $N_i^{\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, so 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 cells $A_k$.

The hierarchical system of the four regression equations is

$$ \begin{aligned} N_i^{\text{BBS}} | \lambda^{\text{BBS}}, \theta_{\text{bbs}} &\sim \text{Pois}\left\{\sum_{k \in \text{route}_i} \omega_k \lambda^{\text{BBS}}(s_k; \theta_{\text{bbs}})\right\}, \\ N_i^{\text{ckl}} | \lambda^{\text{ckl}}, \theta_{\text{ckl}} &\sim \text{Pois}\left\{\lambda^{\text{ckl}}(s_i, t_i; \theta_{\text{ckl}})\right\}, \\ N_i^{\text{spc}} | N_i^{\text{ckl}}, p^{\text{spc}}, \theta_{\text{spc}} &\sim \text{Bin}\{N_i^{\text{ckl}}, p^{\text{spc}}(s_i, t_i; \theta_{\text{spc}})\}, \\ Z_i | \mu, \theta_{\mu}, \sigma, \theta_{\sigma} &\sim \text{GEV}\{\mu(s_i, t_i; \theta_{\mu}), \sigma(s_i; \theta_{\sigma}), \zeta\}, \end{aligned} $$
where

$\theta_{bbs}, \theta_{ckl}, \theta_{spc}, \theta_{\mu}, \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 $N_i^{\text{ckl}}$-related model component, we model them as a spatio-temporal log-Gaussian Cox point process (Møller et al., 1998), discretized according to the cells $A_i$.

## 3.4 Structure of latent Gaussian processes

We write $X(\bullet)$ to denote any latent Gaussian effect, here either indexed by pixels $s \in S$ or years $t \in T$. We use Gaussian process priors for all latent effects. We further use the exponential covariance function to parametrize the spatial effects, given as follows for two locations $s_1$ and $s_2$:

$$ \text{Cov}\{X(s_1), X(s_2)\} = \sigma^2 \exp(-\|s_1 - 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

$$ \begin{aligned} X^{\text{pref}}(\bullet) &\sim \mathcal{GP}(\omega_1), & X^{\text{niche}}(\bullet) &\sim \mathcal{GP}(\omega_2), \\ X^{\text{GEV-}\mu}(\bullet) &\sim \mathcal{GP}(\omega_3), & X^{\text{GEV-}\sigma}(\bullet) &\sim \mathcal{GP}(\omega_4), \end{aligned} \quad (2) $$

where $\omega_1, \omega_2, \omega_3$, and $\omega_4$ contain individual range and standard deviation parameters that control the spatial dependence of each process. These hyperparameters are assigned identical joint penalized complexity 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^{\text{niche}}$, whereas spatial patterns of the sampling effort are modelled through $X^{\text{pref}}$. Moreover, we include a temporal latent effect $X^{\text{year}}$ to capture the variability in checklist counts across years, and set the prior $X^{\text{year}}(\bullet) \sim \mathcal{GP}(\omega_5)$, where $\omega_5$ contains the range and standard deviation parameters of a temporal exponential covariance function.

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

$$ \begin{aligned} \log \lambda^{\text{BBS}}(s_i) &= \beta_0^{\text{BBS}} + X^{\text{niche}}(s_i), \\ \log \lambda^{\text{ckl}}(s_i, t_i) &= \beta_0^{\text{ckl}} + X^{\text{year}}(t_i) + X^{\text{pref}}(s_i), \\ \text{cloglog} \{ p^{\text{spc}}(s_i, t_i) \} &= \beta_0^{\text{spc}} + X^{\text{niche}}(s_i) + \frac{\beta^{\text{act}}}{d_{s_i, t_i}}, \\ \mu(s_i, t_i) &= g\{\beta_0^{\text{GEV-}\mu} + X^{\text{GEV-}\mu}(s_i) + \beta_1^{\text{GEV-}\mu} \text{NAO}_{t_i} \\ &\quad + \theta^{\text{niche-GEV}} X^{\text{niche}}(s_i), x_{\text{effort}}(s_i, t_i)\}, \\ \log \sigma(s_i) &= \beta_0^{\text{GEV-}\sigma} + X^{\text{GEV-}\sigma}(s_i). \end{aligned} $$

The climate covariable, i.e. the NAO value for year $t_i$, is denoted as $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 parameter $\theta^{\text{niche-GEV}}$ controls the sharing of the niche spatial effect. All scaling hyperparameters and intercepts in the model, along with $\xi$, are given flat normal prior distributions with mean 0 and variance 100. The complementary log-log link function $\text{cloglog}(p) = \log(-\log(1-p))$ for binomial probabilities allows interpreting them as the probabilities of observing at least one count 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 hypothesize that it influences the first arrivals $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_{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 arrivals. A higher *activity* could increase the probability of observing the species in a checklist, which could also induce an earlier observed first arrival. We combine the two effects *preference* and *activity* into

$$x_{\text{effort}}(s_i, t_i) = \theta^{\text{eff}} + \theta^{\text{pref}} \log \lambda^{\text{ckl}}(s_i, t_i) + \frac{\theta^{\text{act}}}{d_{s_i,t_i}} \quad (3)$$

with real-valued scaling hyperparameters $\theta^{\text{pref}}$ (typically expected to be positive) and $\theta^{\text{act}}$ (typically expected to be negative), and an intercept $\theta^{\text{eff}}$ representing a baseline sampling effort that is modulated according to the checklist abundance and observer activity. With this structure, setting an infinite duration $d_{s_i,t_i}$ causes the activity term to vanish.

The next subsection explains the choice we make for the nonlinear link function g of the GEV location parameter.

## 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 to identify the four latent processes. Specifically, we want to identify the effect of spatio-temporally heterogeneous sampling effort and correct it in predictions, such as those of the true first arrivals.

The most often used 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(\bullet)$, while another has an additive component $\beta \times X(\bullet)$, with a scaling factor $\beta \in \mathbb{R}$. This yields models akin to the preferential-sampling framework of Diggle et al. (2010) and Pati et al. (2011), and to Koh et al. (2023) and 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 for first arrival dates. Letting the effort tend to infinity with linear sharing would lead to first arrivals tending to 1st January, the smallest possible value. As illustrated by the exploratory data analysis in Figure 1, the observation effort effect on first arrivals saturates at very high levels. Therefore, we share $x_{\text{effort}}$ with the $\mu$ component of the GEV through the following upper-bounded function:

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

where the first argument determines the upper bound in the GEV location parameter, which is modulated by the sampling effort in the second argument and can be attained with unlimited (i.e. infinite) effort. The g-function is always positive, which is appropriate here when considering as response the transformed dates $Z_i > 0$. Figure 4 illustrates that the class of functions (4) is bounded and monotonically increasing, therefore allowing the effect of the sampling effort on the first arrivals 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 first observed arrival dates. To our knowledge, such an approach has not yet been considered in the literature. It establishes a representation that is both parsimonious and ‘physically meaningful’, a feature that is desirable in statistical modelling in various contexts such as landslide prediction (Opitz et al., 2022) or temperature extremes (see the 2024 preprint "Integration of physical bound constraints to alleviate shortcomings of statistical models for extreme temperatures" from Noyelle et al.).
<figure><img src="image_14.png" /><figcaption>Figure 4. Plots of the nonlinear sharing function g in (4) for varying $x_{effort}$ (x-axis) and $x_{bound}$ (colours).</figcaption></figure>

# 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 the nonlinear sharing of latent spatial effects and the inclusion of linear predictors in several parameters of the GEV, it is not feasible to use them here. Instead, we develop an MCMC sampler for simulation-based inference.

## 4.1 Vecchia approximation for Gaussian components

Due to the large number of spatial cells in our case study ($D = 1,268$), 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 (SPDE approach, Lindgren et al., 2011), or nearest-neighbour Gaussian processes (NNGPs, Datta et al., 2016). Here, we use another flexible and popular option to construct approximate Gauss-Markov representations of any Gaussian covariance function through the Vecchia (1988) approximation, which leads to a sparse Cholesky factor of the precision matrix. We write {$x_1, ..., x_D$} for the set of observations of a Gaussian field $X$ evaluated at locations $s_1, ..., s_D \in S$, and we consider a permutation $m: \{1, ..., D\} \to \{1, ..., D\}$ defining a reordering of the observations. Based on this reordering, we can write $H(i; m) = [j \in \{1, ..., D\}: m(j) < m(i)]$ for the 'history' of the $i$th index based on the permutation $m$, with $x_{H(i;m)}$ denoting the corresponding subvector of observations. For a chosen permutation $m$, the exact joint density of the observations can be written as the product of conditional densities, i.e.

$$f(x_1, ..., x_D) = f(x_{m(1)}) \prod_{i=2}^{D} f(x_{m(i)} | x_{H(i;m)})$$

The Vecchia approximation still represents a valid Gaussian process, but assumes that

$$f(x_1, ..., x_D) \approx \hat{f}(x_1, ..., x_D) = f(x_{m(1)}) \prod_{i=2}^{D} f(x_{m(i)} | x_{S(i;m)}), \quad (5)$$

where $S(i; m) \subseteq H(i; m)$ and $|S(i; m)| = k$. The approximation reduces the conditioning history in the conditional densities, which decreases the computational complexity of evaluating the density from $O(D^3)$ to $O(Dk^3)$. This reduction is considerable if $k \ll D$, as in our setting with $D > 1,000$. An active research area (Katzfuss & Guinness, 2021) is to understand how one should choose the
permutation *m* and order the conditioning variables in (5) to obtain approximations that satisfy certain criteria of optimality. Vecchia (1988) 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 in Section 4.3 confirms that estimation works well.

## 4.2 Proposal scheme for MCMC

We draw posterior inference on the parameters using a general MCMC sampling scheme with customized Metropolis-Hastings (MH) updates. The hyperparameters are updated using Gibbs sampling. The most computationally intensive part of our scheme involves sampling from the latent Gaussian components, which are updated jointly through a MALA (Roberts & Rosenthal, 2001) designed to speed up mixing of the chains. We write *x* for the current state of the latent Gaussian component of interest, *x*\* for a proposal for that component, and *y* for the data vector. A discrete approximation of the Langevin diffusion implies that MALA is an MH algorithm with proposal distribution

$$ (x^* | x) \sim N_p \left( x + \frac{\delta^2}{2} \nabla_{\theta} \log \pi(x | y), \delta^2 \mathbf{W}_p \right), \quad (6) $$

where tuning parameters are the prewhitening matrix W<sub>p</sub>, used to account for the correlation between parameters within the component, and the stepsize δ. For example, if we set *x* = {X<sup>niche</sup>(s<sub>1</sub>), ..., X<sup>niche</sup>(s<sub>D</sub>)}, the gradient in (6), when updating this component, is

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

where the first three terms represent the gradients of log-likelihoods for the model components that X<sup>niche</sup> is included in, and

$$ \nabla_x \log \pi_{\text{Vecchia}}(x) = -\tilde{\Sigma}_{\omega_2}^{-1} x, $$

where the *D* × *D* precision matrix $\tilde{\Sigma}_{\omega_2}^{-1}$ is the Vecchia-approximated the precision of the Gaussian field prior assigned to X<sup>niche</sup>, as detailed in §4.1 and implied by (5). The precision matrix is parametrized by the vector $\omega_2 = (\sigma, \kappa)$, see §3.4, to which we assign a hyperprior.

For our data application, we parallelize and run separate chains for each species. We generate 80, 000 posterior samples and discard the first 60, 000 as the burn-in period. We then perform a thinning operation of the Markov chains and keep one of four consecutive samples, thus obtaining 5, 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 the 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). When possible, we then used the last values of those chains as starting values when fitting the full model. Calculations were performed on UBELIX,<fnref n="2" /> the HPC cluster at the University of Bern, and the computational time of our algorithm there is roughly 6 s per iteration.

<footnote n="2">² https://www.id.unibe.ch/hpc</footnote>
## 4.3 Simulation study

We check by simulation that the model parameters are identifiable and can be estimated appropriately through posterior mean estimates obtained from the above Bayesian inference procedure. We especially wish to check whether the spatio-temporal field of first arrivals and the niche of the species are identifiable. This verification is 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 = 1,268$ 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 the species *Chimney Swift*. To test model robustness, we allow for potential model misspecification and overdispersion in the checklist counts. We thus add another layer of randomness to the count component, giving

$$N_i^{\text{ckl}} | \Lambda_i^{\text{ckl}}, \theta_{\text{ckl}} \sim \text{Pois}(\Lambda_i^{\text{ckl}}), \\ \Lambda_i^{\text{ckl}} \sim \Gamma\{r, r\lambda_i^{\text{ckl}}(s_i, t_i; \theta_{\text{ckl}})\},$$

where $\Gamma(a, b)$ is the Gamma distribution with shape $a > 0$ and scale $b > 0$. As a result, $N_i^{\text{ckl}}$ has a negative binomial distribution with mean $\lambda_i^{\text{ckl}}(s_i, t_i)$ and an inflated variance $\lambda_i^{\text{ckl}}(s_i, t_i) + \lambda_i^{\text{ckl}}(s_i, t_i)^2/r$. We set $r = 10$ and 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 A2) suggest that the posterior means of the median first arrivals 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 A2 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. We obtained similar simulation results in the setting where there is no model misspecification (not shown), i.e. $r = \infty$.

## 5 Results

We apply the model constructed and checked in Sections 3 and 4 to the bird data described in Section 2, where we run separate models for the 10 species with largest occurrence numbers in eBird. We report a selection of results for four species: *Chimney Swift*, *Great Crested Flycatcher*, *Chestnut-sided Warbler*, and *Purple Martin*. The MCMC trace plots in the Appendix (Figure A3) for one of the aforementioned species indicate that the chains are mixing relatively well, although some parameters like $\theta^{\text{act}}$ and $\beta_0^{\text{GEV-}\mu}$ still display autocorrelation after the burn-in period and thinning.

To compare our approach with a simpler one where only the first arrivals are considered but without modelling the niche or sampling effort, we also report results for a model using only the $Z_i$-variables but not $N_i^{\text{BBS}}$, $N_i^{\text{ckl}}$, and $N_i^{\text{spc}}$. We fit this model separately and call it the GEV-only model, for which $\theta^{\text{niche-GEV}}$ is null and $x_{\text{effort}}$ set to a constant.

Figure 5 illustrates that there have been more checklists in recent years, and the model 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 mostly negative (or positive with 95% posterior credible interval including zero) 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.

Table 1 indicates that the sharing parameter $\theta^{\text{niche-GEV}}$ often has posterior credible intervals that do not contain zero, with most posterior means in the range [0, 0.15]. This means that birds tend to be observed earlier in areas where they are relatively abundant. For $\beta_1^{\text{GEV-}\mu}$, posterior means are generally hovering around [0.01, 0.03] but often not containing 0 in its credible intervals, indicating that first arrivals can take place slightly earlier in years with high NAO. Higher NAO and temperature anomaly correspond to warmer weather, which could explain earlier arrivals. However,
<figure><img src="image_15.png" /><figcaption>**Figure 5.** Observed numbers of checklists, $N_j^{\text{ckl}}$ (left) against model-based posterior mean estimates of the Poisson mean number of checklists, $\lambda^{\text{ckl}}(\mathbf{s}_i, t_i)$ (right), using the model fitted to *Chimney Swift* in 2001 (top) and 2021 (bottom).</figcaption></figure>

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 arrivals could be present for some species but there is no strong signal of a significant NAO effect across all species; Wijeyakulasuriya et al. (2024) obtained a similar result when applying their analyses to the *Magnolia Warbler*.

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

## 5.1 Excursion sets of latent spatial fields

We study 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 arrivals.

To identify areas with very high or low values of the spatial random effects, we utilize credible sets for excursion regions (Bolin & Lindgren, 2015). We evaluate where the fields exceed or fall below certain thresholds. To visualize excursion sets simultaneously for all values of the probability level $\alpha \in (0, 1)$, Bolin and Lindgren (2015) introduced the positive and negative excursion functions $F_u^+(s) = 1 - \inf\{\alpha | s \in E_{u,a}^+\} \in [0, 1]$ and $F_u^-(s) = 1 - \inf\{\alpha | s \in E_{u,a}^-\} \in [0, 1]$, where $E_{u,a}^+$ is the subregion with maximal surface area in which the spatial field fully exceeds the threshold $u$ with probability $\alpha$ (with an analogous definition of $E_{u,a}^-$ for threshold deficits below $u$). Figure 6 highlights these excursion functions for the species *Purple Martin*. Sampling effort is high in the southeastern areas near the ocean and strongly decreases in the further northwestern areas. The niche of the species tends to concentrate in southwestern areas. The fields related to first arrivals show more fragmented spatial patterns for this species.

## 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 their posterior means with the land-cover maps described in Section 2 and shown in the Appendix
<table><caption>**Table 1.** Examples of predictions related to first arrivals for two pixel-years in 2022</caption><thead><tr><th>Species</th><th>Chimney Swift</th><th>Great Crested Flycatcher</th><th>Chestnut-sided Warbler</th><th>Purple Martin</th></tr></thead><tbody><tr><td>β̂<sup>pref</sup></td><td>0.191 (0.184,0.202)</td><td>0.204 (0.199,0.21)</td><td>0.187 (0.183,0.191)</td><td>0.2 (0.178,0.217)</td></tr><tr><td>β̂<sup>act</sup></td><td>-0.15<br>(-0.217,-0.061)</td><td>-0.818<br>(-0.911,-0.696)</td><td>-0.548<br>(-0.619,-0.454)</td><td>-0.03<br>(-0.269,0.236)</td></tr><tr><td>β̂<sup>niche-GEV</sup> (× 10<sup>-2</sup>)</td><td>4.9 (4.664,5.134)</td><td>4 (3.894,4.133)</td><td>0.2 (0.17,0.278)</td><td>6 (5.541,6.443)</td></tr><tr><td>β<sup>GEV-μ</sup><sub>1</sub> (× 10<sup>-1</sup>)</td><td>0.1 (-0.083,0.204)</td><td>0.2 (0.125,0.331)</td><td>0.3 (0.295,0.409)</td><td>0 (-0.384,0.526)</td></tr><tr><td>Observed</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>Predicted (GEV only)</td><td>28/05</td><td>27/05</td><td>01/07</td><td>18/06</td></tr><tr><td>Predicted</td><td>09/05</td><td>03/05</td><td>21/05</td><td>07/06</td></tr><tr><td>Debiased</td><td>03/04</td><td>13/04</td><td>03/05</td><td>28/03</td></tr><tr><td>Observed</td><td>01/05</td><td>04/05</td><td>04/05</td><td>29/06</td></tr><tr><td>Predicted (GEV only)</td><td>25/05</td><td>20/05</td><td>15/05</td><td>19/05</td></tr><tr><td>Predicted</td><td>09/05</td><td>15/05</td><td>12/05</td><td>12/05</td></tr><tr><td>Debiased</td><td>22/04</td><td>05/05</td><td>03/05</td><td>07/04</td></tr></tbody></table>

Note. The first five rows report the species and the posterior mean estimates of selected parameters, with 95% credible intervals in brackets. The next rows give the observed, predicted (with the GEV-only model, and with the full model), and debiased (with the full model) median first arrivals (in day/month format) at the two pixels indicated by the two triangles in Figure 8); the first four rows here correspond to the pixel in the southeastern area, the last four rows to the pixel located in the middle of the study region. An ‘NA’ value in the Observed row indicates that there were no observations within that pixel in 2022. The strikethroughs over the dates indicate that our full model estimates a species presence probability of less than 1% for that pixel-year.

(Figure A1). This can improve our understanding of 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 are 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. Chestnut-sided Warbler), and with Developed land cover for species known to breed in Developed areas (e.g. Chimney Swift). The latent fields describing the spatial variability in first arrivals show generally weaker correlations with specific land-cover types. For example, *Chimney Swift* tends to settle earlier in Developed and Vegetation areas than in Water- and Forest-dominated areas.

## 5.3 Predicting true first arrivals

To illustrate the utility of our model, we predict the true first arrivals in all pixels. Firstly, Figure 8 shows that the model reproduces generally well the spatial variability in the probability of observing species presence in eBird, with some differences in empirical and model-predicted probabilities due to the additional information provided by BBS data. We use this component of the model to avoid predicting first arrivals 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 presence probability $p^{spc}((s_i, t_i))$ is smaller than 1%. This feature of our model avoids making predictions in areas where data give no clear signal of a species’ breeding activity.

For the prediction of the true first arrivals, we set $λ^{ckl}(s_i, t_i)$ and $d_{s_i, t_i}$ to be close to infinity to mimic infinite sampling effort, so we expect to retrieve the true and unbiased dates of the first arrival. This approach is illustrated in Figure 9 for the *Great Crested Flycatcher*, where we also show that the model captures the spatial pattern in the observed first arrivals.
<figure><img src="image_16.png" /><figcaption>**Figure 6.** Excursion set analysis for the four latent spatial fields $X^{\text{effort}}$ (top left), $X^{\text{niche}}$ (top right), and $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.</figcaption></figure>

To expand on our prediction procedure, we focus on 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 highlights that the model can extrapolate and draw posterior predictive samples for that pixel-year. The second pixel is chosen in the central area of the study region. For example, a first arrival was observed in 2022 for the *Purple Martin* on the 29th of June, but Table 1 (last two rows) shows that our model flags this pixel as being uncertain to be located in the species' niche; it is thus 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 very low species presence probability. Table 1 also shows that the GEV-only model generally estimates a later first arrival than that from the full model after correcting biases. The former predictions usually fall outside of what is expected of the four species' phenology in our study region, i.e. late spring/early summer rather than early spring.

# 6 Discussion and outlook

## 6.1 Statistical learning from participatory data

Participatory data and citizen-science programs advance and share scientific knowledge, thanks particularly to the data's large volume, diversity, and spatio-temporal coverage. However, sampling protocols for these data sources are usually lenient, with different levels of observer expertise. Statistical techniques are required to correct for heterogeneous sampling effort and observation biases.

In this work, we developed a BHM in which the sampling effort and its influence on the observation bias of a phenological event can be quantified and then corrected during prediction. Our results show that the required corrections are often substantial. We could also estimate spatial probabilities of species presence, and then filter out implausible regions when predicting the phenological event, since it might never occur.