<figure><img src="image_17.png" /><figcaption>**Figure 7.** Correlations of land-cover proportions (x-axis) with posterior means of the latent fields $X^{\text{effort}}$, $X^{\text{niche}}$, $X^{\text{GEV-}\mu}$, and $X^{\text{GEV-}\sigma}$ (y-axis), for four species: *Chimney Swift* (top left), *Great Crested Flycatcher* (top right), *Chestnut-sided Warbler* (bottom left), and *Purple Martin* (bottom right).</figcaption></figure>

## 6.2 Ecological data fusion using latent processes

We carefully constructed statistical models to infer and interpret latent components based on in- complete and biased observations of the true ecological processes. BHMs have many advantages: they allow us to identify various latent processes and their interactions, to separate the process de- scribing how data are observed from the latent ecological processes we seek to reveal, to keep track of uncertainties and assess them through posterior sampling, and to incorporate expert knowl- edge. The Bayesian framework is also relatively robust to noisy data and overdispersion in count data thanks to the inclusion of appropriately designed random effects. Data fusion in Integrated Species Distribution Models (e.g. Miller et al., 2019, and the 2023 preprint "The Point Process Framework for Integrated Modelling of Biodiversity Data" of Adjei et al.) combining opportun- istic and protocol-based data, is facilitated by shared latent components. Further extensions to bet- ter consider temporal dynamics are very actively explored in the current literature but are challenging since very complex ecological and observational dynamics must be disentangled.

Machine learning and deep learning tools are used in ecology for their skill in predicting ob- served data, though they may be less effective at capturing complex latent processes. Further de- velopments are needed to allow such algorithms to predict the processes of interest when these are not directly observed, to extrapolate beyond the range of observed data (as we have done by set- ting the observation effort to a very high level in predictions), and to avoid propagating or even reinforcing biases from the input data when generating predictive outputs (Dunson, 2018). Future research could devise mechanisms similar to the sharing of random effects in BHMs but implemented within the architectures of general machine learning and deep learning approaches.

## 6.3 Modelling sampling effort

The eBird sampling effort has various dimensions. We have here considered the annual number of visits to an area (using checklists) and the median duration of each of those visits. Understanding how the combination of these two aspects influences the number of detected species individuals is a complex task (e.g. see Tang et al., 2021), and would merit exploration beyond the additive
<figure><img src="image_18.png" /><figcaption>**Figure 8.** Empirical mean probability (left displays) of observing the species *Great Crested Flycatcher* (top) and *Purple Martin* (bottom) in checklists aggregated over all years; posterior mean probability (right displays) of observing the species from the presence model by using the mean checklist duration aggregated over all years as covariate. The two triangles are the two pixels used in the predictions of first arrivals detailed in Section 5.3.</figcaption></figure>

structure of $x_{\text{effort}}$ that we posed in (3). Other relevant attributes that we did not consider here are partially available for eBird checklists, such as the surface of the area where observations took place, or the distance over which the observers moved during the checklist event. Such information could be further studied for a more complete characterization of the sampling effort and its impact on observation biases for species distributions. Identifying sampling effort through BHMs could also be of interest in other domains, for example for reports of hail events and damages, and more generally for citizen-science data in climate science.

In cases where the performance of data fusion of opportunistic data with protocol-cased data could be improved by collecting additional exhaustive field data, it would be instrumental to conduct simulation studies using BHMs to explore how new sampling locations can be placed in an optimal and cost-efficient way.

## 6.4 Detectability of species

Observation biases related to the detectability of species remain a major challenge since they are difficult to identify with available data, i.e. it is possible that individuals of a species were present but not observed. This risk of nondetection can vary with land cover, species behaviour, observer experience, and other factors. In our model, detectability can be viewed as being included in the niche process, which could more precisely be called the *detectable niche*. It is difficult to devise a general statistical approach to further disentangle the niche and the detectability. Expert knowledge about the conditions of detectability could provide further insights and improve models, but including it was beyond the scope of this work. It would also require designing models adapted to species-specific properties, and could hinder simple intercomparison of models for different species.
<figure><img src="image_19.png" /><figcaption>**Figure 9.** Maps for the species *Great Crested Flycatcher*. Observed first arrivals of *Great Crested Flycatcher* in 2022 (top left, posterior mean of the GEV median first arrivals in 2022 (top right), bias-corrected mean of the GEV median first arrivals in 2022 (bottom left), and the same plot as before though with indications (in white rectangles over the corresponding pixels) for where the model estimates a presence probability of less than 1% for the species (bottom right). The grey pixels in the top left display indicate where no species occurrences were recorded in 2022, or where no checklists and therefore no median durations were available in that pixel.</figcaption></figure>

## 6.5 Validation of results

Validation of modelling results is inherently challenging with data that are obtained through heterogeneous sampling efforts and ridden with unknown biases. We can only validate directly whether the model appropriately reproduces data as they have been observed. For example, we performed such checks for checklist counts and observed first arrivals in Section 5. However, the latent processes of interest, as well as sampling effort and observation biases, can be entangled in complex ways. Direct validation of inferences on the latent processes of interest is thus difficult. Careful construction of interpretable models and appropriate inference algorithms is paramount.

External protocol-based data sources not suffering from these drawbacks can be used, but they often have much smaller spatial, temporal, and species coverage. The fusion of data sources, such as in our study with BBS and eBird data, can not only improve the estimation of certain model components, but could also be used for validation purposes at space-time locations covered by both datasets.

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

## 6.6 Possible extensions of our approach

We focused on identifying past and present spatio-temporal trends in ecological processes. In species distribution modelling, important research concerns the potential impact of climate change
and its interplay with other dynamic processes (e.g. land use change). Our modelling approach could be refined to establish more detailed relationships between ecological dynamics and climate variables, to then be used for projecting first arrivals under future climate simulations, similar to Wijeyakulasuriya et al. (2024).

We have fixed the spatial pixel mesh using a size similar to that routinely used in modelling results published by the eBird project. One could explore the influence of the mesh size on the properties of the observed first arrivals. For example, a larger mesh size leads to larger numbers of observations within each pixel. This could reduce the observation bias. On the other hand, it would also decrease the spatial resolution of predictive maps and hamper the identification of relationships between ecological processes and land cover. Further research could investigate the role of the mesh choice, or how to optimally combine results from different mesh sizes.

Regarding temporal trends, the seasonal variation of sampling effort and of the species presence probability, along with how these properties interact with first arrivals, could also be modelled more precisely.

If some spatio-temporal domains are visited at high frequency, the corresponding sampling effort could be regarded as being exhaustive, and the structure of the model, e.g. the g-function in (4), could be modified accordingly to improve the identifiability of model components.

Another extension consists of studying interactions among different species in terms of their niche and phenological events, e.g. to infer which species groups show similar spring migration patterns. One could implement joint species distribution models with each having a set of regression equations and with the possibility to share certain spatial random effects between several species; this could be feasible with our current approach for a small number of species. Alternatively, one could fit models separately for each species as we have done here, and then apply classification approaches on the posterior estimates of latent model components to group species that show similar model behaviour.

## 6.7 Ecological extreme-value analysis

EVT could play a more prominent role in ecological science. An important domain of application is the analysis and probabilistic prediction of extreme climatic and environmental events, which can strongly influence population dynamics and biodiversity. Here, the tools from EVT could pave the way towards ecological analyses focusing not only on long-term climate averages but also on the influence of specific extreme weather events. A second domain of application, of which the present work is part, concerns the analysis of extreme phenological events using EVT. Although analyses of this type are still in their early infancy, we hope that our work will motivate more widespread applications of extreme-value techniques.

## Acknowledgements

The authors thank the Editor and two anonymous referees for their numerous valuable comments and suggestions to improve the quality of the manuscript.

*Conflicts of interest:* None declared.

## Funding

No funding sources are declared.

## Data availability

The data in this article are shared publicly on https://github.com/kohrrelation/mcmc_birds.

## Appendix

### A.1 Land-Cover Plots

Figure A1 shows plots of the four land-cover covariates defined as proportions of certain land-cover categories in the pixels of the study area.
<figure><img src="image_20.png" /><figcaption>**Figure A1.** Land-cover proportions in the study area (NLCD, 2021) according to four types described in Section 2.</figcaption></figure>

## A.2 Plots from the Simulation Study

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

<figure><img src="image_21.png" /><figcaption>**Figure A2.** Simulation study: Truth vs. posterior means for the median first arrivals in 2003 (left) and 2020 (right), shown with spatial maps (for three different seeds in the simulation study), and with boxplots of the differences in days between true and predicted arrival dates, based on 100 simulations. The boxplot shows results for 150 randomly chosen pixels, marked with black boxes in the spatial plots.</figcaption></figure>
## A.3 Further Plots of MCMC Results

Figure A3 shows trace plots from the MCMC run of the model on the species *Great Crested Flycatcher*. Figure A4 shows spatial plots of the four latent spatial fields from the model estimated for the species *Chimney Swift*.

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

**Figure A3.** MCMC trace plots for some fixed effect and sharing parameters (top left to bottom right): $\beta^{\text{act}}$, $\beta_0^{\text{spc}}$, $\beta_0^{\text{BBS}}$, $\theta^{\text{niche-GEV}}$, $\theta^{\text{pref}}$, $\theta^{\text{act}}$, $\beta_0^{\text{GEV-\mu}}$, $\beta_0^{\text{GEV-\sigma}}$, and $\epsilon$ from our model for the species *Great Crested Flycatcher*. The dashed line indicates the burn-in period.

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

**Figure A4.** Posterior means plots of the four spatial random effects in (2) for the model fitted to the species *Chimney Swift*.
## References

Banerjee S., Carlin B. P., & Gelfand A. E. (2003). *Hierarchical modeling and analysis for spatial data*. Chapman and Hall/CRC.

Belzile L. R., Davison A. C., Gampe J., Rootzén H., & Zholud D. (2022). Is there a cap on longevity? A statistical review. *Annual Review of Statistics and Its Application*, 9(1), 21–45. https://doi.org/10.1146/statistics.2022.9 issue-1

Bock C. E., & Root T. L. (1981). The Christmas bird count and avian ecology. *Studies in Avian Biology*, 6, 17–23.
http://hdl.handle.net/1969.3/24710

Bolin D., & Lindgren F. (2015). Excursion and contour uncertainty regions for latent Gaussian models. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 77(1), 85–106. https://doi.org/10.1111/rssb.12055

Bonney R., Cooper C., & Ballard H. (2016). The theory and practice of citizen science: Launching a new journal. *Citizen Science: Theory and Practice*, 1(1). https://doi.org/10.5334/cstp.56

Bonney R., Phillips T. B., Ballard H. L., & Enck J. W. (2016). Can citizen science enhance public understanding of science? *Public Understanding of Science*, 25(1), 2–16. https://doi.org/10.1177/0963662515607406

Borda-de Água L., Alírezazadeh S., Neves M., Hubbell S. P., Borges P. A. V., Cardoso P., Dionísio F., & Pereira H. M. (2021). *Species accumulation curves and extreme value theory*, Ecology, biodiversity and conservation (pp. 211–226). Cambridge University Press.

Browne J. (1996). Charles Darwin: A biography, vol. i, voyaging. *Journal of the History of Biology*, 29(2), 314–316. https://doi.org/10.1126/science.268.5214.1196

Coles S. (2001). *An introduction to statistical modeling of extreme values*. Springer.

Colwell R. K., Coddington J. A., & Hawksworth D. L. (1994). Estimating terrestrial biodiversity through extrapolation. *Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences*, 345(1311), 101–118. https://doi.org/10.1098/rstb.1994.0091

Conklin J. R., Lisovski S., & Battley P. F. (2021). Advancement in long-distance bird migration through individual plasticity in departure. *Nature Communications*, 12(1), 4780. https://doi.org/10.1038/s41467-021-25022-7

Cornell Lab of Ornithology (2022). eBird Basic Dataset, Version: EBD relNov-2022.

Cotton P. A. (2003). Avian migration phenology and global climate change. *Proceedings of the National Academy of Sciences*, 100(21), 12219–12222. https://doi.org/10.1073/pnas.1930548100

Datta A., Banerjee S., Finley A. O., & Gelfand A. E. (2016). Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. *Journal of the American Statistical Association*, 111(514), 800–812. https://doi.org/10.1080/01621459.2015.1044091

Diggle P. J., Menezes R., & Su T.-L. (2010). Geostatistical inference under preferential sampling. *Journal of the Royal Statistical Society: Series C (Applied Statistics)*, 59(2), 191–232. https://doi.org/10.1111/j.1467-9876.2009.00701.x

Dunson D. B. (2018). Statistics in the big data era: Failures of the machine. *Statistics & Probability Letters*, 136, 4–9. https://doi.org/10.1016/j.spl.2018.02.028

Fink D., Auer T., Johnston A., Ruiz-Gutierrez V., Hochachka W. M., & Kelling S. (2020). Modeling avian full annual cycle distribution and population trends with citizen science data. *Ecological Applications*, 30(3), e02056. https://doi.org/10.1002/eap.2056

Fraisl D., Campbell J., See L., Wehn U., Wardlaw J., Gold M., Moorthy I., Arias R., Piera J., Oliver J. L., Masó J., Penker M., & Fritz S. (2020). Mapping citizen science contributions to the un sustainable development goals. *Sustainability Science*, 15(6), 1735–1751. https://doi.org/10.1007/s11625-020-00833-7

Fraisl D., Hager G., Bedessem B., Gold M., Hsing P.-Y., Danielsen F., Hitchcock C. B., Hulbert J. M., Piera J., Spiers H., Thiel M., & Haklay M. (2022). Citizen science in environmental and ecological sciences. *Nature Reviews Methods Primers*, 2(1), 64. https://doi.org/10.1038/s43586-022-00144-4

Fuglstad G.-A., Simpson D., Lindgren F., & Rue H. (2018). Constructing priors that penalize the complexity of Gaussian random fields. *Journal of the American Statistical Association*, 114(525), 445–452. https://doi.org/10.1080/01621459.2017.1415907

Gaines S. D., & Denny M. W. (1993). The largest, smallest, highest, lowest, longest, and shortest: Extremes in ecology. *Ecology*, 74(6), 1677–1692. https://doi.org/10.2307/1939926

Gelfand A. E., & Shirota S. (2019). Preferential sampling for presence/absence data and for fusion of presence/absence data with presence-only data. *Ecological Monographs*, 89(3), e01372. https://doi.org/10.1002/ecm.1372

Guinness J. (2018). Permutation and grouping methods for sharpening Gaussian process approximations. *Technometrics*, 60(4), 415–429. PMID: 31447491. https://doi.org/10.1080/00401706.2018.1437476

Haklay M., Fraisl D., Tzovaras B. G., Hecker S., Gold M., Hager G., Ceccaroni L., Kieslinger B., Wehn U., Woods S., Nold C., Balázs B., Mazzonetto M., Ruefenacht S., Shanley L. A., Wagenknecht K., Motion A., Sforzi A., Riemenschneider D., ... Vohland K. (2021). Contours of citizen science: A vignette study. *Royal Society Open Science*, 8(8), Article 202108. https://doi.org/10.1098/rsos.202108
Hitz A. S., Davis R. A., & Samorodnitsky G. (2024). Discrete Extremes. *Journal of Data Science*, 524–536. https://doi.org/10.6339/24-JDS1120

Hsing P.-Y., Hill R. A., Smith G. C., Bradley S., Green S. E., Kent V. T., Mason S. S., Rees J., Whittingham M. J., Cokill J., Scientists M. C., & Stephens P. A. (2022). Large-scale mammal monitoring: The potential of a citizen science camera-trapping project in the United Kingdom. *Ecological Solutions and Evidence*, 3(4), e12180. https://doi.org/10.1002/2688-8319.12180

Huser R., Stein M. L., & Zhong P. (2023). Vecchia likelihood approximation for accurate and fast inference with intractable spatial max-stable models. *Journal of Computational and Graphical Statistics*, 33(3), 978–990. https://doi.org/10.1080/10618600.2023.2285332

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(10), 1052–1060. https://doi.org/10.1111/mee3.2014.5 issue-10

Johnston A., Fink D., Hochachka W. M., & Kelling S. (2018). Estimates of observer expertise improve species distributions from citizen science data. *Methods in Ecology and Evolution*, 9(1), 88–97. https://doi.org/10.1111/mee3.2018.9 issue-1

Johnston A., Hochachka W. M., Strimas-Mackey M. E., Gutierrez V. R., Robinson O. J., Miller E. T., Auer T., Kelling S. T., & Fink D. (2021). Analytical guidelines to increase the value of community science data: An example using eBird data to estimate species distributions. *Diversity and Distributions*, 27(7), 1265–1277. https://doi.org/10.1111/ddi.v27.7

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(1), 103–116. https://doi.org/10.1111/mee3.v14.1

Jordan R. C., Ballard H. L., & Phillips T. B. (2012). Key issues and new approaches for evaluating citizen-science learning outcomes. *Frontiers in Ecology and the Environment*, 10(6), 307–309. https://doi.org/10.1890/110280

Katz R. W., Brush G. S., & Parlange M. B. (2005). Statistics of extremes: Modeling ecological disturbances. *Ecology*, 86(5), 1124–1134. https://doi.org/10.1890/04-0606

Katzfuss M., & Guinness J. (2021). A general framework for Vecchia approximations of Gaussian processes. *Statistical Science*, 36(1), 124–141. https://doi.org/10.1214/19-STS755

Kelling S., Johnston A., Hochachka W. M., Iliff M., Fink D., Gerbracht J., Lagoze C., La Sorte F. A., Moore T., Wiggins A., & Wong W. K. (2015). Can observation skills of citizen scientists be estimated using species accumulation curves? *PLoS One*, 10(10), e0139600. https://doi.org/10.1371/journal.pone.0139600

Koh J., Pimont F., Dupuy J.-L., & Opitz T. (2023). Spatiotemporal wildfire modelling through point processes with moderate and extreme marks. *The Annals of Applied Statistics*, 17(1), 560–582. https://doi.org/10.1214/22-AOAS1642

Kosmala M., Wiggins A., Swanson A., & Simmons B. (2016). Assessing data quality in citizen science. *Frontiers in Ecology and the Environment*, 14(10), 551–560. https://doi.org/10.1002/fee.2016.14 issue-10

Linden A. (2011). Using first arrival dates to infer bird migration phenology. *Boreal Environment Research*, 16, 49–60. http://hdl.handle.net/10138/232775

Lindgren F., Rue H., & Lindström J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 73(4), 423–498. https://doi.org/10.1111/j.1467-9868.2011.00777.x

Majumder R., Reich B. J., & Shaby B. A. (2024). Modeling extremal streamflow using deep learning approximations and a flexible spatial process. *The Annals of Applied Statistics*, 18(2), 1519–1542. https://doi.org/10.1214/23-AOAS1847

Marshall P. J., Lintott C. J., & Fletcher L. N. (2015). Ideas for citizen science in astronomy. *Annual Review of Astronomy and Astrophysics*, 53(1), 247–278. https://doi.org/10.1146/astro.2015.53 issue-1

McKinley D. C., Miller-Rushing A. J., Ballard H. L., Bonney R., Brown H., Evans D. M., French R. A., Parrish J. K., Phillips T. B., Ryan S. F., Shanley L. A., Shirk J. L., Stepenuck K. F., Weltzin J. F., Wiggins A., Boyle O. D., Briggs R. D., Chapin III S. F., Hewitt D. A., Preuss P. W., & Soukup M. A. (2015). Investing in citizen science can improve natural resource management and environmental protection. 19

Miller D. A., Pacifici K., Sanderlin J. S., & Reich B. J. (2019). The recent past and promising future for data integration methods to estimate species' distributions. *Methods in Ecology and Evolution*, 10(1), 22–37. https://doi.org/10.1111/mee3.2019.10 issue-1

Møller J., Syversveen A. R., & Waagepetersen R. P. (1998). Log Gaussian Cox processes. *Scandinavian Journal of Statistics*, 25(3), 451–482. https://doi.org/10.1111/1467-9469.00115

Newman G., Wiggins A., Crall A., Graham E., Newman S., & Crowston K. (2012). The future of citizen science: Emerging technologies and shifting paradigms. *Frontiers in Ecology and the Environment*, 10(6), 298–304. https://doi.org/10.1890/110294
Nussbaumer R., Bauer S., Benoit L., Mariethoz G., Liechti F., & Schmid B. (2021). Quantifying year-round noc- turnal bird migration with a fluid dynamics model. *Journal of The Royal Society Interface*, 18(179), Article 20210194. https://doi.org/10.1098/rsif.2021.0194

Opitz T., Bakka H., Huser R., & Lombardo L. (2022). High-resolution Bayesian mapping of landslide hazard with unobserved trigger event. *The Annals of Applied Statistics*, 16(3), 1653–1675. https://doi.org/10.1214/21-AOAS1561

Pati D., Reich B. J., & Dunson D. B. (2011, March). Bayesian geostatistical modelling with informative sampling locations. *Biometrika*, 98(1), 35–48. https://doi.org/10.1093/biomet/asq067

Pocock M. J., Chandler M., Bonney R., Thornhill I., Albin A., August T., Bachman S., Brown P. M., Cunha D. G. F., Grez A., Jackson C., Peters M., Rabarijaon N. R., Roy H. E., Zaviezo T., & Danielsen F. (2018). Chapter six - a vision for global biodiversity monitoring with citizen science. In D. A. Bohan, A. J. Dumbrell, G. Woodward, & M. Jackson (Eds.), *Next generation biomonitoring: Part 2*, Volume 59 of *advances in ecological research* (pp. 169–223). Academic Press.

Prieto F., Gómez-Déniz E., & Sarabia J. M. (2014). Modelling road accident blackspots data with the discrete generalized Pareto distribution. *Accident Analysis & Prevention*, 71, 38–49. https://doi.org/10.1016/j.aap.2014.05.005

Ranjbar S., Cantoni E., Chavez-Demoulin V., Marra G., Radice R., & Jaton K. (2022). Modelling the extremes of seasonal viruses and hospital congestion: The example of flu in a Swiss hospital. *Journal of the Royal Statistical Society Series C: Applied Statistics*, 71(4), 884–905. https://doi.org/10.1111/rssc.12559

Roberts G. O., & Rosenthal J. S. (2001). Optimal scaling for various Metropolis-Hastings algorithms. *Statistical Science*, 16(4), 351–367. https://doi.org/10.1214/ss/1015346320

Rue H., Riebler A., Sørbye S. H., Illian J. B., Simpson D. P., & Lindgren F. K. (2017). Bayesian computing with INLA: A review. *Annual Review of Statistics and Its Application*, 4(1), 395–421. https://doi.org/10.1146/ statistics.2017.4 issue-1

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

Somveille M., Rodrigues A. S., & Manica A. (2015). Why do birds migrate? A macroecological perspective. *Global Ecology and Biogeography*, 24(6), 664–674. https://doi.org/10.1111/geb.2015.24 issue-6

Spearing H., Tawn J., Irons D., Paulden T., & Bennett G. (2021). Ranking, and other properties, of elite swimmers using extreme value theory. *Journal of the Royal Statistical Society Series A: Statistics in Society*, 184(1), 368–395. https://doi.org/10.1111/rssa.12628

Tang B., Clark J. S., & Gelfand A. E. (2021). Modeling spatially biased citizen science effort through the eBird database. *Environmental and Ecological Statistics*, 28(3), 609–630. https://doi.org/10.1007/s10651-021-00508-1

Theobald E., Ettinger A., Burgess H., DeBey L., Schmidt N., Froehlich H., Wagner C., HilleRisLambers J., Tewksbury J., Harsch M., & Parrish J. (2015). Global change and local solutions: Tapping the unrealized potential of citizen science for biodiversity research. *Biological Conservation*, 181, 236–244. https://doi.org/10.1016/j.biocon.2014.10.021

Thibaud E., Aalto J., Cooley D. S., Davison A. C., & Heikkinen J. (2016). Bayesian inference for the Brown–Resnick process, with an application to extreme low temperatures. *Annals of Applied Statistics*, 10(4), 2303–2324. https://doi.org/10.1214/16-AOAS980

UN (2015). A/res/70/1 un general assembly transforming our world: The 2030 agenda for sustainable develop- ment. *Seventieth session of the General Assembly on 25 Sept 2015*.

van de Schoot R., Depaoli S., King R., Kramer B., Märtens K., Tadesse M. G., Vannucci M., Gelman A., Veen D., Willemsen J., & Yau C. (2021). Bayesian statistics and modelling. *Nature Reviews Methods Primers*, 1(1), 1. https://doi.org/10.1038/s43586-020-00001-2

Vecchia A. V. (1988). Estimation and model identification for continuous spatial processes. *Journal of the Royal Statistical Society. Series B (Methodological)*, 50(2), 297–312. https://doi.org/10.1111/j.2517-6161.1988.tb01729.x

Wijeyakulasuriya D. A., Hanks E. M., & Shaby B. A. (2024). Modeling first arrival of migratory birds using a hierarchical max-infinitely divisible process. *Journal of Agricultural, Biological and Environmental Statistics*, https://doi.org/10.1007/s13253-024-00624-y

Wijeyakulasuriya D. A., Hanks E. M., Shaby B. A., & Cross P. C. (2019). Extreme value-based methods for mod- eling elk yearly movements. *Journal of Agricultural, Biological and Environmental Statistics*, 24(1), 73–91. https://doi.org/10.1007/s13253-018-00342-2

Wynn J. (2017). *Citizen science in the digital age: Rhetoric, science, and public engagement*. Rhetoric, science, and public engagement. Rhetoric, science, and public engagement. The University of Alabama Press Tuscaloosa.

Yadav R., Huser R., Opitz T., & Lombardo L. (2023). Joint modelling of landslide counts and sizes using spatial marked point processes with sub-asymptotic mark distributions. *Journal of the Royal Statistical Society Series C: Applied Statistics*, 72(5), 1139–1161. https://doi.org/10.1093/jrssc/qlad077

Youngflesh C., Socolar J., Amaral B. R., Arab A., Guralnick R. P., Hurlbert A. H., LaFrance R., Mayor S. J., Miller D. A., & Tingley M. W. (2021). Migratory strategy drives species-level variation in bird sensitivity to vegetation green-up. *Nature Ecology & Evolution*, 5(7), 987–994. https://doi.org/10.1038/s41559-021-01442-y
# Proposer of the vote of thanks and contribution to the Discussion of the ‘Discussion Meeting on the Analysis of citizen science data’

Ben Swallow<fnref n="1" />,<fnref n="2" /> ID

*Address for correspondence:* Ben Swallow, School of Mathematics and Statistics, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK; Centre for Research into Ecological and Environmental Modelling, University of St Andrews, Buchanan Gardens, St Andrews KY16 9LZ, UK. Email: bts3@st-andrews.ac.uk

## 1 Introduction

Data are a vital part of statistical inference, grounding and calibrating statistical models, and supporting prediction and forecasting. Whilst data are routinely collected for analyses, they can be very expensive and/or resource intensive to collect. This can lead to prohibitively small-scale sampling, making it hard to draw conclusions at a population level. Citizen science (CS) data, however, can provide an alternative that has the potential to provide data at spatial and/or temporal resolutions that would not otherwise be feasible. These data utilize a network of highly skilled amateurs who routinely observe and make record of personal sightings.

Due to the often nonstructured nature of data collection, it is particularly important with CS data to model the observation process in which these data are observed and recorded to help account for the lack of randomized sampling and potential structural biases. This is particularly important given that they often feed into important national and international monitoring schemes, underpinning legal obligations of conservation policy, for example, Newson et al. (2017).

A particular challenge with CS data is the potential for a large range of biases in them, including from oversampling (that is rare or noteworthy events being more likely to be reported, spatial bias in high-population regions, visiting ‘good’ locations with a higher chance of successful observations); miss-identification; partial observations or checklists being submitted; and multiple observations of the same individuals. In order to alleviate these potential biases, rich metadata needs to be collected to help explain the process in which the data were collected, recorded, and submitted.

Other challenges include high data resolution leading to ‘big data’ computational challenges, including hardware demands, a need for efficient algorithms and multiple heterogeneities across the population and data sources. In addition, communication and interaction with domain experts or nonstatisticians are common due to the nature of both data collectors and the end users. A lack of statistical involvement in the design and implementation of data collection can make many of these issues further challenging to overcome.

It is, therefore, particularly welcome to comment on the three novel approaches presented here.

<footnote n="1">¹School of Mathematics and Statistics, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK</footnote>

<footnote n="2">²Centre for Research into Ecological and Environmental Modelling, University of St Andrews, Buchanan Gardens, St Andrews KY16 9LZ, UK</footnote>
## 2 Efficient statistical inference methods for assessing changes in species' populations using citizen science data

The two main developments of this paper are, firstly, to incorporate time and site effects through an annual model integrated within the generalized abundance index. A concentrated likelihood framework is proposed to significantly reduce the dimension of the parameter space. The approach is fitted to both data from the UK Butterfly Monitoring Scheme and the Big Butterfly Count. Secondly, the authors propose the use of variational inference to speed up the efficiency of fitting the occupancy model framework. Using the observed data likelihood, required assumptions of independence between model parameters are avoided and a reparameterization provides an approach for efficient stochastic gradient descent. A small simulation study is implemented to study coverage probabilities of credible intervals for the parameters, with the intervals having nominal coverage. The approach is also applied to a common species of butterfly from the Butterflies for the New Millennium recording scheme.

### 2.1 Comments and questions

It seems the new development of a two-stage model offers differing estimates of uncertainties, however, which of these is correct is uncertain. Some simulations would be advantageous to determine which of these is more realistic. Wald-type standard errors from the Hessian will likely be unreliable in rarer species and/or small studies.

An advantage of the approach is that some very simple adjustments to the model, such as changing the discrete family distribution, can be made to determine impacts on uncertainty propagation. For example, a negative binomial or quasi-Poisson distribution could be compared using simple (quasi-)information criteria to select between these. There is an advantage in allowing different levels of variation across different levels of heterogeneity in the population. Global dispersion parameters from these alternative families can determine species-level variation across spatial and temporal regions, with additional variation allocated to the random effects in a hierarchical manner.

Furthermore, species populations and associated distribution are inherently spatial phenomena. This paper, amongst others, paves the way for accounting for the spatial correlation beyond independent site effects, combined with efficient inference for binomial models (e.g. Belmont et al., 2024; Held & Holmes, 2006). An alternative without quite as much computational overhead would be to incorporate spatial covariates that may account for much of the correlation.

Finally, as there are many species, some very rare, some discussion of multispecies approaches to learn across similar species would be fruitful, perhaps in conjunction with integrated population models that are already discussed in the paper. A common species with a similar ecological niche could be highly informative of the distributions of an undersampled species.

## 3 Frequentist prediction sets for species abundance using indirect information

The authors develop a frequentist-valid prediction sets framework for eBird data in North Carolina counties, assuming a nonparametric framework to obtain sets for a multinomial random sample. The approach uses data from more sampled counties to better inform undersampled counties through an integrated approach. An empirical Bayes framework for estimating prediction sets uses the maximum marginal likelihood calculated using the Newton-Raphson method to determine parameter values. The approach utilizes neighbouring counties to indirectly inform estimates, where more data diminishes the impact of neighbouring counties and data-scarce counties benefit from updated posteriors from informed neighbours. This has benefits when sampling is comparable and unbiased by reducing the size of the prediction set in most of the species tested.

### 3.1 Comments and questions

The authors show that assuming independence of observations across county, frequentist coverage probabilities are theoretically guaranteed. As discussed above, the assumption is perhaps not valid with CS data due to preferential sampling of observers, incomplete lists and observer
heterogeneity. The coverage probabilities are above expected for moderately large *K* (see Figure 2), so some discussion of this in relation to potential over-confident standard errors would be welcome. Further attempts to account for the sampling process could alleviate this (e.g. Bird et al., 2014). Further work could develop a simulation incorporating this heterogeneity to determine the impact of preferential sampling on the prediction sets and determine sensitivity to this assumption or how to account for it if not realistic.

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

The authors study the first arrival dates of migrating birds from the global eBird CS project, and the more structured Breeding Bird Survey. Focusing on the northeastern US states, they treat the first arrival date as an extreme event, and develop an extreme value theory approach to predicting first arrival dates, accounting for survey effort and environmental niche. A complex hierarchical Bayesian model incorporating four latent Gaussian process priors to explain the spatial covariance in arrival dates across regions, sampling effort, and environmental niches is proposed. The data from two sources are fused together with common model structures to alleviate sampling bias in the CS data.

## 4.1 Comments and questions

There is evidently difficulty with validating latent variable models such as these. Variation in goodness-of-fit tests and information criteria for latent variable models are made more difficult in the complex observation process of CS data. More work is definitely needed here, particularly in estimating effective numbers of parameters. The suggestion of combining CS with more structured data definitely shows promise (e.g. Kelling et al., 2019; Robinson et al., 2020) and needs further work to determine when and how this should be done and how valuable it is in a specific scenario.

Multivariate extreme approaches are in their infancy but determining marginal dependence, e.g. through copulas, multivariate Generalized Pareto Distributions or tensor products of Gaussian Processes, would be a worthwhile extension to this approach. This may allow for sharing of information across species, where early arrivals of one species may be able to suggest similar species may also arrive early.

The Markov chain Monte Carlo (MCMC) chain mixing shown in Figure 12 shows there is a strong correlation between model components that may not be fully alleviated by the Metropolis-adjusted Langevin algorithm approach used. The authors highlight that a Gauss Markov random field may not be feasible due to nonlinearity of the model, but perhaps some components could be embedded within the INLA-SPDE approach with others estimated by MCMC (Gómez-Rubio & Rue, 2018)?

The accuracy of the Vecchia approximation used will also depend on the validity of the sparsity assumption, which is linked to the smoothness of the surface. This may work well for some species but less well for others with more patchy distributions.

Finally, all species chosen had almost identical preference parameter estimates for $\hat{\theta}^{\text{pref}}$ (Table 1). It would be interesting to choose some species that clearly show a very strong or weak observer preference, to show these are estimable and not prior-dependent.

# 5 Summary and looking forward

The three papers presented here show a wide computational variability of statistical approaches for tackling applied ecological questions using CS data. These include novelties covering a mix of computational, methodological, and philosophical approaches to model these complex data. Harnessing the power of information contained within CS data, they aim to account for (some) of the potential biases or sub-optimality of sampling inherent in their chosen schemes.

Ecological sciences have made CS data a cornerstone of the domain, giving much opportunity for methodological development; however, there remains many other areas of applied sciences where the power of these data could be better realized (e.g. epidemiology/infectious disease dynamics; environmental surveillance; public health). Further computational demands (including methods, hardware integration, and software) are needed to analyse the vast quantities of data that can be obtained
from these schemes. As shown in these papers, the combination of differing data sources and their associated heterogeneities can provide significant benefits over single sources of data. With those benefits comes a greater need for improved methods for relative and absolute goodness-of-fit for these data and models, however, to allow for differing resolutions and collection procedures. The three contributions here show the potential for further development in analysing CS data and highlight that further work is still needed from statisticians working in this area.

*It therefore gives me great pleasure to propose the vote of thanks.*

*Conflicts of interest:* None declared.

## References

Belmont J., Martino S., Illian J., & Rue H. (2024). 'Spatio-temporal occupancy models with INLA.' *Methods in Ecology and Evolution*, 15(11), 2087-2100. https://doi.org/10.1111/2041-210X.14422.

Bird T. J., Bates A. E., Lefcheck J. S., Hill N. A., Thomson R. J., Edgar G. J., Stuart-Smith R. D., Wotherspoon S., Krkosek M., Stuart-Smith J. F., Pecl G. T., Barrett N., & Frusher S. (2014). Statistical solutions for error and bias in global citizen science datasets. *Biological Conservation*, 173(6), 144-154. https://doi.org/10.1016/j.biocon.2013.07.037

Gómez-Rubio V., & Rue H. (2018). Markov chain Monte Carlo with the integrated nested Laplace approximation. *Statistics and Computing*, 28(5), 1033-1051. https://doi.org/10.1007/s11222-017-9778-y

Held L., & Holmes C. C. (2006). Bayesian auxiliary variable models for binary and multinomial regression. *Bayesian Analysis*, 1(1), 145-168. https://doi.org/10.1214/06-BA105

Kelling S., Johnston A., Bonn A., Fink D., Ruiz-Gutierrez V., Bonney R., Fernandez M., Hochachka W. M., Julliard R., Kraemer R., & Guralnick R. (2019). Using semistructured surveys to improve citizen science data for monitoring biodiversity. *Bioscience*, 69(3), 170-179. https://doi.org/10.1093/biosci/biz010

Newson S. E., Evans H. E., Gillings S., Jarrett D., Raynor R., & Wilson M. W. (2017). Large-scale citizen science improves assessment of risk posed by wind farms to bats in southern Scotland. *Biological Conservation*, 215(1), 61-71. https://doi.org/10.1016/j.biocon.2017.09.004

Robinson O. J., Ruiz-Gutierrez V., Reynolds M. D., Golet G. H., Strimas-Mackey M., & Fink D. (2020). Integrating citizen science data with expert surveys increases accuracy and spatial extent of species distribution models. *Diversity & Distributions*, 26(8), 976-986. https://doi.org/10.1111/ddi.v26.8

https://doi.org/10.1093/jrsssa/qnaf009
Advance access publication 12 February 2025

# Seconder of the vote of thanks and contribution to the Discussion of 'the Discussion Meeting on the Analysis of citizen science data'

## Kerrie Mengersen

QUT Centre for Data Science, Queensland University of Technology (QUT), Brisbane, Australia

*Address for correspondence:* Queensland University of Technology (QUT), Centre for Data Science.

Email: k.mengersen@qut.edu.au

Thank you for the opportunity to second the vote of thanks at this meeting. I want to start by saying that after reading the three papers presented in this session, it's such a wonderful time to be a statistician! There are three reasons for saying this.
The first is that, as we all know, we are in the best of professions. Statistics is currently one of the most in-demand professions globally, in no small part due to the need to 'do something' with the wealth of data that is now available. This creates an ideal space for us to flourish as researchers and practitioners. The high profile of statistics as a trusted 'data-whisperer' is in no small part due to the great efforts of our professional societies, and the RSS is one of the global leaders in this endeavour, so thank you.

The second reason is that we not only have a fantastic new range of tools thanks to the rise of statistical machine learning, but we also have the compute power to implement our methods and shift the horizon. Models and algorithms that were 'nice but not practical' a decade ago are now common tools and established research directions. Emerging algorithmic enablers—such as foundation models and quantum computing—promise to open even more doors, creating again a new round of challenges and opportunities for statisticians as we formulate, validate and integrate solutions.

The third reason is that we now have a huge range of new data sources available to us, to wrangle, model and understand. Data from digital sensors, cameras, satellites, social media, and—yes—citizen science provide rich information in their own right and as a complement to more traditional sources of data. As statisticians, we have an obligation to embrace these new data sources and be the trusted voice for researchers, practitioners, managers, and the public regarding their benefits and deficits. If the data are rubbish, then we can—and should—say so. However, if there is merit, we have an opportunity—indeed a responsibility—to contribute to research, tools, and applications that enhance the utility of these data, extract information in a rigorous manner and merge the resultant insights with other sources of knowledge, thereby creating a true data-informed learning and decision-making environment.

Which brings me to the three papers under discussion at this meeting. In accord with the theme of the meeting, the common ambition of these papers is to enhance the utility of citizen science data. Although all three papers focus on applications in ecology, the methodological innovations are generalizable to a wide range of other contexts.

Within this broad remit, the three papers take slightly different directions. Two papers focus primarily on improved models for estimation and inference which correct for sampling bias. Elizabeth and Peter articulate a two-stage approach, with neighbourhood information entering as 'indirect information' through the concentration parameter of the Bayesian prior in one of the stages. In contrast, Jonathan and Thomas propose a strictly Bayesian hierarchical spatio-temporal model that similarly adjusts for heterogeneous sampling effort. Both papers comment on computational efficiency. Emily and coauthors focus primarily on this aspect and propose variational inference as a solution.

I have so many comments and questions about each of these papers, but in this forum I would like to make just three points.

First, all three authors deal with what I would call 'well-behaved' citizen science data. That is, the data are sourced from large, public studies with an underlying design, albeit not to the standard of a randomized controlled study. This is great, and the contributions are very worthwhile. However, there is more to be done: there are other biases in these data that need to be addressed to increase their trust and utility, and there are citizen science data from more 'wild' studies that also have great value. So keep going!

Second, the papers advance the practice of data fusion. Elizabeth and Peter focus on one study—the eBird survey—and fuse information across the geographic sites; Jonathan and Thomas 'see this and raise it' to the merger of two separate studies with different protocols, the eBird study and the North American Breeding Bird Survey. Emily and coauthors raise the ante again to the merger of three butterfly datasets. This work touches on a large literature that rectifies biases and gaps in this way—I call out earlier work of Sylvia Richardson in epidemiology, and the United Nations Global Platform for enabling the use of different types of remotely sensed data to monitor the Sustainable Development Goals, as examples.

The papers also advance the practice of model fusion. In particular, Elizabeth and Peter suggest a merger of frequentist and Bayesian methods to achieve the desired improvement in the size and precision of the area-specific species prediction sets. This 'smooshing' of models and philosophies is now commonplace and potentially very effective, but again as statisticians we need to ensure that this is undertaken in a rigorous manner, that all of the properties of the resultant approach
are understood, that comparisons with other approaches are thorough and fair, and that ‘edge
cases’ are carefully considered. Furthermore, we should be actively stepping into the space of fus-
ing machine learning and statistical models, particularly for the analysis of citizen science data. I
think that there is much more work to be done in this area.

My third point is that all of the papers work hard to ‘respect’ the citizen science data. That is,
their new tools aim to address the characteristics of the data rather than ‘make the shoe fit’ with
existing tools. This is important, and pushes us to consider the relative benefits and drawbacks of
employing models at all in this context. This is a longer discussion for another time.

It would be remiss of me to not mention the efforts of colleagues in my part of the world to capitalize
on the enormous potential of citizen science data. A special call-out to early work by Mark Burgman,
Hugh Possingham, and Sama Low Choy in eliciting expert information to enhance ecological model-
ling, and more recently Edgar Santos-Fernandez and Julie Vercelloni for their leadership in using
citizen-contributed classifications of coral cover in underwater images to improve monitoring of the
health of the Great Barrier Reef. Through Julie, this effort has expanded to the Pacific Islands, which
are particularly vulnerable to socio-industrial development and climate change.

In summary, I applaud all of the authors in this session for their very welcome, substantive
contributions to the three areas that I touched on in my introduction: to raising the profile of sta-
tisticians as ‘data-whisperers’ in the development and application of trusted methods, to extending
our computational capability, and to advancing the analysis of new data sources—in particular,
citizen science data. It is with great pleasure that I second the vote of thanks for these stimulating,
cutting-edge papers.

*Conflicts of interest:* none declared.

https://doi.org/10.1093/jrsssa/qnaf013
Advance access publication 11 February 2025

Andrej Srakar's contribution to the
# Discussion of 'the Discussion Meeting
on the Analysis of citizen science data'

Andrej Srakar ID

Institute for Economic Research and University of Ljubljana, Ljubljana, Slovenia

Address for correspondence: Andrej Srakar, Institute for Economic Research, Kardeljeva pl. 17, 1000 Ljubljana, Slovenia.

Email: srakara@ier.si

This contribution discusses the papers of Dennis et al., Bersson and Hoff, and Koh and Opitz.

Paper by Dennis et al. (2024) mentions reducing the functional model parameter space. As a mathematician, I wonder which mathematical spaces (such as metric spaces of Banach or Hilbert space variety) to use to study citizen science (CS) data. Could inclusion of tensor data approaches and tensorial ‘spaces’ be of help? Namely studying CS data leads to high-dimensional problems which gives this question a special importance. As the paper uses concentrated likelihood functions, I wonder whether using approaches which rely on intractable likelihoods (such as indirect inference and in particular approximate Bayesian computation) would provide responses to many open problems addressing analysis of CS data. In general,

This is a discussion contribution to the Discussion Meeting on the ‘Analysis of citizen science data’ at the 2024 Royal Statistical Society conference in Brighton, United Kingdom.
do you see any special properties of CS data which would give priority to Bayesian over frequentist approaches?

Paper by Bersson and Hoff (2024) uses borrowing of information, and approach related to empirical Bayes procedures. Interesting could be to refer to and use recent thread of research from Ascolani, Prünster, Lijoi, and coauthors (e.g. Ascolani, Franzolini et al. 2024) on nonparametric priors with full-range borrowing of information. Also, addressing geographic adjacency structure seems an extremely interesting aspect of work with CS data. This immediately leads to spatiotemporal but also random matrix perspectives, and it would be great to hear more about future possibilities here.

Paper by Koh and Opitz (2024) mentions development of bird and ecology related applications of CS data. Could they be found in economics, in medicine and healthcare, in robotics or in law? Could they be fruitfully applied and found in any scientific and sector field? Usage of Gaussian latent random fields bring questions how to study asymptotics of approaches studying CS data—is there any aspect to be careful about? Using MCMC and MALA algorithms questions how to approach computational efficiency for modelling CS data in the future? One also wonders whether modelling species behaviour using microsimulation or even agent-based simulation approaches be used combined with CS data and how.

Papers providing mathematical and statistical approaches to study CS data are highly important based on the present spread of such initiatives as well as the fact that leading journals reject paper using CS data on the ground of their questionable validity. It would be great for all three papers to explain how their findings can be used to mitigate such problems in the future.

*Conflicts of interest:* None declared.

## References

Bersson E., & Hoff P. D. (2025). Frequentist prediction sets for species abundance using indirect information. *Journal of the Royal Statistical Society: Series A*, **188**(3), 658–673. https://doi.org/10.1093/jrsssa/qnae096

Dennis E. B., Diana A., Matechou E., & Morgan B. J. T. (2025). Efficient statistical inference methods for assessing changes in species' populations using citizen science data. *Journal of the Royal Statistical Society: Series A*, **188**(3), 641–657. https://doi.org/10.1093/jrsssa/qnae105

Koh J., & Opitz T. (2025). Extreme-value modelling of migratory bird arrival dates: Insights from citizen science data. *Journal of the Royal Statistical Society: Series A*, **188**(3), 674–699. https://doi.org/10.1093/jrsssa/qnae108

https://doi.org/10.1093/jrsssa/qnaf017

Advance access publication 8 February 2025

# Benjamin R. Baer's contribution to the Discussion of 'the Discussion Meeting on Analysis of citizen science data'

Benjamin R. Baer ID

Centre for Research into Ecological and Environmental Modelling, School of Mathematics and Statistics, University of St Andrews, St Andrews KY16 9SS, Scotland

Address for correspondence: Benjamin R. Baer, Centre for Research into Ecological and Environmental Modelling, School of Mathematics and Statistics, University of St Andrews, St Andrews KY16 9SS, Scotland. Email: benjamin.baer@st-andrews.ac.uk

We congratulate Dennis, Diana, Matechou, and Morgan, the authors of 'Efficient statistical inference methods for assessing species' populations using citizen science data', on an interesting
contribution to the literature on computationally efficient algorithms for statistical problems involving citizen science data and thank them for an interesting conversation about their work.

We start by briefly reviewing the problem set up of the first example. Consider counts $Y_{svr}$ of butterflies where $s = 1, \dots, 4,328$ denotes sites, $r = 1, \dots, 47$ denotes years, and $v = 1, \dots, 26$ denotes weeks/visits. Dennis et al. (2024) define the model

$$ 
\begin{aligned}
Y_{svr} &\sim \text{Pois}(n_{sr}a_{vr}), \\
\log n_{sr} &= a_s + \beta_r, \\
a_{vr} &= \varphi(v; \mu_r, \sigma_r),
\end{aligned}
$$

where $\varphi$ denotes the normal density. With independent observations, the log likelihood up to an additive constant is

$$ \ell(\alpha, \beta, \mu, \sigma) := \sum_{s,v,r} -e^{a_s + \beta_r} \varphi(v; \mu_r, \sigma_r) + Y_{svr}\{\alpha_s + \beta_r + \log \varphi(v; \mu_r, \sigma_r)\}. $$

Let $(\hat{\alpha}, \hat{\beta}, \hat{\mu}, \hat{\sigma})$ be maximizers of $\ell$. Throughout a ‘•’ in a subscript denotes aggregation over that index. Applying the first order condition, Dennis et al. (2024, sec. S1) show that $e^{\hat{\alpha}_s} = Y_{s\bullet\bullet} / \sum_{v,r} e^{\hat{\beta}_r} \varphi(v; \hat{\mu}_r, \hat{\sigma}_r)$ for all $s$. Thus the log profile likelihood $\ell_p(\beta, \mu, \sigma) := \ell(\hat{\alpha}, \hat{\beta}, \mu, \sigma)$ admits $(\hat{\beta}, \hat{\mu}, \hat{\sigma})$ as maximizers. In their data example with Gatekeepers, the sample size is $n = 527,588$, and the (profiled) parameter dimension is $p = 141 = 3.47$ so that $n \gg p$. Without the profiling step, there would be an additional 4,328 parameters, so profiling results in a large reduction.

Now, we consider a further computational speedup. The parameter $\alpha$ is the only one that varies with site $s$: once it is profiled out to create $\ell_p$, the remaining parameters do not depend on site. The following calculation shows that the observations may, therefore, be aggregated over site.

**Proposition 1** Up to an additive constant,

$$ \ell_p(\beta, \mu, \sigma) = \sum_{v,r} Y_{\bullet vr}\{\beta_r + \log \varphi(v; \mu_r, \sigma_r)\} - Y_{\bullet \bullet} \log \left\{ \sum_{v,r} e^{\beta_r} \varphi(v; \mu_r, \sigma_r) \right\}. $$

**Proof.** Denote ‘$=$’ as equality up to additive constants. The log profile likelihood satisfies

$$ 
\begin{aligned}
& \sum_{s,v,r} -e^{\hat{\alpha}_s + \beta_r} \varphi(v; \mu_r, \sigma_r) + Y_{svr}\{\hat{\alpha}_s + \beta_r + \log \varphi(v; \mu_r, \sigma_r)\} \\
&= \sum_{s,v,r} -e^{\hat{\beta}_r} \frac{Y_{s\bullet\bullet}}{\sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'})} \varphi(v; \mu_r, \sigma_r) + Y_{svr} \left[ \log \left\{ \frac{Y_{s\bullet\bullet}}{\sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'})} \right\} + \beta_r + \log \varphi(v; \mu_r, \sigma_r) \right] \\
&\stackrel{\cdot}{=} \sum_{s,v,r} -Y_{s\bullet\bullet} \frac{e^{\hat{\beta}_r} \varphi(v; \mu_r, \sigma_r)}{\sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'})} + Y_{svr}\{\beta_r + \log \varphi(v; \mu_r, \sigma_r)\} - Y_{svr} \log \left\{ \sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'}) \right\} \\
&= \sum_{v,r} -Y_{\bullet\bullet\bullet} \frac{e^{\hat{\beta}_r} \varphi(v; \mu_r, \sigma_r)}{\sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'})} + \sum_{v,r} Y_{\bullet vr}\{\beta_r + \log \varphi(v; \mu_r, \sigma_r)\} - Y_{\bullet\bullet\bullet} \log \left\{ \sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'}) \right\} \\
&\stackrel{\cdot}{=} \sum_{v,r} Y_{\bullet vr}\{\beta_r + \log \varphi(v; \mu_r, \sigma_r)\} - Y_{\bullet\bullet\bullet} \log \left\{ \sum_{v',r'} e^{\hat{\beta}_{v'}} \varphi(v'; \mu_{v'}, \sigma_{v'}) \right\}.
\end{aligned}
\quad \square $$

The form of the profile likelihood involves a sum over only visits (denoted by $v$) and years (denoted by $r$). This reduces the size of the problem from $n = 527,588$ to $1,222 = 26.47$, which
reduces the run time to evaluate the profile likelihood and, hence, the run time of an algorithm to compute associated maximum likelihood estimators.

This note illustrates that inspecting the form of objective functions can lead to simplifications and computational speedups. Separately, if the simplification had not occurred and a sum over $n = 527,588$ observations remained in the objective function, then an algorithm that scales well with large $n$, such as stochastic gradient descent, may be preferred to limited-memory BFGS, the algorithm employed by Dennis et al. (2024).

*Conflicts of interest:* The author has no conflicts of interest.

## Reference

Dennis E. B., Diana A., Matechou E., & Morgan B. J. (2024). Efficient statistical inference methods for assessing changes in species' populations using citizen science data. *Journal of the Royal Statistical Society: Series A*, 188(3), 641–657. http://dx.doi.org/10.1093/jrsssa/qnae105

https://doi.org/10.1093/jrsssa/qnaf011

Advance access publication 11 February 2025

# Jafet Belmont, Sara Martino, Janine Illian, and Håvard Rue's contribution to the Discussion of the 'Discussion meeting on the analysis of citizen science data'

Jafet Belmont<fnref n="1" /> ID, Sara Martino<fnref n="2" />, Janine Illian<fnref n="1" /> and Håvard Rue<fnref n="3" />

*Address for correspondence:* Jafet Belmont, School of Mathematics and Statistics, University of Glasgow, 132 University Pl, Glasgow G12 8QA, Glasgow, G12 8QQ, UK. Email: Jafet.BelmontOsuna@glasgow.ac.uk

The authors highlight the prominent role of occupancy modelling in the context of citizen science (CS) data and the need to account for imperfect detection. They comment on the computational cost of fitting these models and, consequently, the importance of computational efficiency, which they address with a variational inference (VI) approach.

We suggest considering integrated-nested Laplace approximation (INLA) as an alternative, more efficient fitting approach that can be used in the current context since occupancy models can now be fitted using the library *R-INLA* (Belmont et al., 2024), provided detection depends on fixed effects only. To illustrate this, we have used INLA to fit the model discussed in the paper to the Ringlet butterfly dataset (R code accessible through https://github.com/Ecol-Stats/CS_data_analysis). This took about 30 s on a standard Apple M2 Pro machine, i.e. about 120 times faster than the VI approach with very similar estimates, see Figure 1a,b.

As an additional benefit, it is straightforward to include other complex model structures, such as smooth terms or, specifically, spatio-temporal effects as these are readily available in *R-INLA*. The latter are particularly relevant in the context of ecological data as spatial random effects can account for spatial dependencies not captured by covariates avoiding wrong inference (Wright

<footnote n="1">¹School of Mathematics and Statistics, University of Glasgow, 132 University Pl, Glasgow G12 8QA, UK</footnote>

<footnote n="2">²Department of Mathematics, Norwegian University of Science and Technology, Trondheim, Norway</footnote>

<footnote n="3">³CEMSE Division, King Abdullah University of Science and Technology, KAUST, Thuwal, Saudi Arabia</footnote>
<figure><img src="image_24.png" /><figcaption>**Figure 1.** (a) Posterior median and 96% credible interval for yearly probability of occurrence. (b) Posterior median for the probability of detection. Both results refer to the model in Dennis et al. (2024).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 2.** Posterior median of the probability of occurrence of the Ringlet butterfly (*Aphantopus hyperantus*) in space and across time, for the model including a spatial term.</figcaption></figure>

et al., 2019). In an extended model, we have included a spatially structured random effect using the SPDE approach (Lindgren et al., 2011) as implemented in R-INLA, similar to the model in Diana et al. (2023), but omitting yearly random effects for the detection probabilities (computation time approximately 5 min). Figure 2 shows the estimated probability of occurrence over the whole domain for four different years; the results resemble those in Diana et al. (2023).

We agree with the authors on the need for developing accurate, computationally efficient—yet still practically relevant and accessible—methods to address the computational challenges resulting from large CS data. The INLA framework, along with the R-INLA package, addresses exactly this. Therefore, we believe that the growing community of R-INLA users can greatly benefit from its relevance to CS data as well as its general versatility. For instance, the package provides goodness of fit and prediction accuracy measures (Adin et al., 2024; Van Niekerk et al., 2023), and a straight forward way of combining different data structures into a single global model (Illian et al., 2013; Martino et al., 2021; Panunzi et al., 2024).
However, we also agree that for some complex data structures, e.g. where detection may have to depend on random effects, Markov chain Monte Carlo as a less computationally efficient approach remains the method of choice for Bayesian inference. Striking the right balance here requires a continuous dialogue between practitioners and method developers to adequately analyse CS data and to encourage innovation in tackling complex ecological questions (Illian & Burslem, 2017).

*Conflicts of interest:* All the authors listed made substantial contributions to the manuscript and qualify for authorship, and no authors have been omitted. We warrant that none of the authors has any conflict of interest in regard to this manuscript.

## References

Adin A., Krainski E. T., Lenzi A., Liu Z., Martínez-Minaya J., & Rue H. (2024). Automatic cross-validation in structured models: Is it time to leave out leave-one-out? *Spatial Statistics*, 62(1), 100843. https://doi.org/10.1016/j.spasta.2024.100843

Belmont J., Martino S., Illian J., & Rue H. (2024). Spatio-temporal occupancy models with INLA. *Methods in Ecology and Evolution*, 15(11), 2087–2100. http://dx.doi.org/10.1111/2041-210X.14422

Dennis E. B., Diana A., Matechou E., & Morgan B. J. (2025). Efficient statistical inference methods for assessing changes in species’ populations using citizen science data. *Journal of the Royal Statistical Society: Series A*, 188(3), 641–657. https://doi.org/http://dx.doi.org/10.1093/jrsssa/qnae105

Diana A., Dennis E. B., Matechou E., & Morgan B. J. T. (2023). Fast Bayesian inference for large occupancy datasets. *Biometrics*, 79(3), 2503–2515. https://doi.org/10.1111/biom.13816

Illian J. B., & Burslem D. F. (2017). Improving the usability of spatial point process methodology: An interdisciplinary dialogue between statistics and ecology. *Advances in Statistical Analysis: AStA: A Journal of the German Statistical Society*, 101(4), 495–520. https://doi.org/10.1007/s10182-017-0301-8

Illian J. B., Martino S., Sørbye S. H., Gallego-Fernández J. B., Zunzunegui M., Esquivias M. P., & Travis J. M. (2013). Fitting complex ecological point process models with integrated nested Laplace approximation. *Methods in Ecology and Evolution*, 4(4), 305–315. https://doi.org/10.1111/mee3.2013.4 issue-4

Lindgren F., Rue H., & Lindström J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 73(4), 423–498. https://doi.org/10.1111/j.1467-9868.2011.00777.x

Martino S., Pace D. S., Moro S., Casoli E., Ventura D., Frachea A., Silvestri M., Arcangeli A., Giacomini G., Ardizzone G., & Jona Lasinio G. (2021). Integration of presence-only data from several sources: A case study on dolphins’ spatial distribution. *Ecography*, 44(10), 1533–1543. https://doi.org/10.1111/ecog.2021.v44.i10

Panunzi G., Moro, S., Marques I., Martino S., Colloca F., Ferretti F., & Jona Lasinio G. (2024). Estimating the spatial distribution of the white shark in the Mediterranean Sea via an integrated species distribution model accounting for physical barriers. *Environmetrics*, 36(1), e2876. http://dx.doi.org/10.1002/env.2876

Van Niekerk J., Krainski E., Rustand D., & Rue H. (2023). A new avenue for Bayesian inference with INLA. *Computational Statistics & Data Analysis*, 181(1), 107692. https://doi.org/10.1016/j.csda.2023.107692

Wright W. J., Irvine K. M., & Higgs M. D. (2019). Identifying occupancy model inadequacies: Can residuals separately assess detection and presence? *Ecology*, 100(6), e02703. https://doi.org/10.1002/ecy.2019.100 issue-6
# Léo R. Belzile and Rishikesh Yadav's
contribution to the Discussion of 'the
Discussion Meeting on the Analysis of
citizen science data'

Léo R. Belzile ID and Rishikesh Yadav ID

Department of Decision Sciences, HEC Montréal, Montréal, Canada

Address for correspondence: Léo R. Belzile, Department of Decision Sciences, HEC Montréal, 3000, ch. Côte-Sainte-Catherine, Montréal (Québec), Canada H3T 2A7. Email: leo.belzile@hec.ca

We congratulate Koh and Opitz for this stimulating piece of work. The complex framework adopted by the authors offers crucial insights into migratory patterns and species behaviours that simpler models may overlook. Much effort has gone into the model specification to account for sampling biases and to ensure interpretability, and it shows.

Hierarchical Bayesian models are well suited for this type of modelling exercise, but the inference is complicated and preprocessing tedious. Assuming separability of space and time effects simplifies the problem, and offers the possibility of using leave-one-year-out cross-validation, although this is computationally intensive. The sharing of random effects allows one to borrow strength across data sources, but may lead to model misspecification without great care. Adding covariates, such as land cover, could reduce the residual variability, but we acknowledge that their effect may be nonlinear, and suitable smooths would add multiple fixed effect parameters.

One concern is the slow convergence and poor mixing observed in Figure 12, even after thinning. We wonder what the effective sample size is after burn-in. Although standard methods, such as adaptive MCMC (Andrieu & Thoms, 2008; Rosenthal, 2011), could help, hyperparameters may be strongly correlated due to shared components, and joint updates may be necessary to increase the efficiency of the sampler. The performance of Metropolis adjusted Langevin algorithm (MALA) is highly sensitive to the global tuning parameter or prewhitening matrix; locally adaptive schemes could fare better (e.g. Girolami & Calderhead, 2011; Rue & Held, 2005, Section 4.4.1).

Although the Markov chains seem to stabilize eventually, running multiple chains could be used to check whether they reach a unique stationary distribution. Model predictions at the data level (e.g. Figure 9) look sensible, but it is unclear whether individual model components are identifiable. For example, consider the function relating the generalized extreme value distribution location parameter $\mu$ with the sampling effort, $g(x_{bound}, x_{effort}) = \exp(x_{bound})/(1 + \exp(-x_{effort}))$. The functional form of eq. (4) implies that we cannot distinguish between parameters for $x_{bound}$ when $x_{effort}$ is low. Increases in $g$ (and thus in $\mu$) lead to an earlier minimum arrival rate.

We believe data fusion of related databases observed at different locations or resolutions has great potential in the field of spatial extreme value analysis as there is limited information available and pooling can help partly alleviate this.

*Conflicts of interest:* The authors have no conflict of interest to declare.

## References

Andrieu C., & Thoms J. (2008). A tutorial on adaptive MCMC. *Statistics and Computing*, 18(4), 343–373.
https://doi.org/10.1007/s11222-008-9110-y

Girolami M., & Calderhead B. (2011, March). Riemann manifold Langevin and Hamiltonian Monte Carlo methods. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 73(2), 123–214. https://doi.org/10.1111/j.1467-9868.2010.00765.x. ISSN 1369-7412
Rosenthal J. S. (2011). Optimal proposal distributions and adaptive MCMC. In S. Brooks, A. Gelman, G. Jones, & X. Meng, (Eds.), *Handbook of Markov Chain Monte Carlo* (pp. 93–111). CRC Press. https://doi.org/10.1201/b10905-5

Rue H., & Held L. (2005). *Gaussian Markov random fields: Theory and applications*. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. CRC Press. https://doi.org/10.1201/9780203492024. ISBN 9780203492024.

https://doi.org/10.1093/jrsssa/qnaf016
Advance access publication 8 February 2025

# Abdelaati Daouia and Gilles Stupfler's contribution to the Discussion of the 'Discussion Meeting on the Analysis of citizen science data'

Abdelaati Daouia<fnref n="1" /> and Gilles Stupfler<fnref n="2" />

*Address for correspondence:* Gilles Stupfler, Laboratoire Angevin de Recherche en Mathématiques, 2 Boulevard de Lavoisier, Angers Cedex 49045, France. Email: gilles.stupfler@univ-angers.fr

We congratulate Jonathan Koh and Thomas Opitz for developing a sophisticated but fully oper-
ational spatiotemporal Bayesian hierarchical model for an explicit bias-corrected estimation of
arrival dates from combined heterogeneous datasets. An econometrician may have thought, how-
ever, that this problem also belongs to the area of (frequentist) frontier analysis, which provides a
rather different, and, we hope, useful perspective. More precisely, for a given species of birds and
each pixel of a spatial mesh covering the study area, let Y represent the arrival date reported by an
observer and T be the corresponding year of observation. As highlighted in the paper, the condi-
tional distribution of Y given T obviously has a tail bounded to the left. The finite lower extremity
of the support of Y given T = t defines the frontier point that corresponds to the population left-
endpoint of all dates of occurrence of the species during the year t. This frontier function of t can
be estimated by a global envelopment spline smoother under/without shape (e.g. concavity) con-
straints with automatic selection information criteria of the number and location of knots (Daouia
et al., 2016), or by a local polynomial approach based on local extreme value statistics with adaptive
selection of the tuning parameters (Jirak et al., 2014). Both approaches are based on the char-
acterization of the yearly first arrival dates as a regression function in a nonparametric regression
model with one-sided errors. The resulting spline or local polynomial boundary would then de-
scribe the yearly evolution of the desired first arrival date at the local zone considered. This avoids
setting up a complex Bayesian hierarchical structure, as well as the Generalized Extreme Value
approach whose goodness-of-fit to the data is not always easy to assess. Furthermore, instead
of restricting the analysis to each pixel, thus resulting in a discrete geographical differentiation
of the first arrival date, we suggest using all the data covering the study area and incorporating
the observation location (i.e. longitude and latitude), as well as exogenous climate, land cover
and other ecological factors into the estimation procedure as covariates (Wang et al., 2020), which
would provide spatially smoother predictions of the true first arrivals. Proceeding in this way
would also implicitly integrate the sampling effort, reflected by both preference and activity di-
mensions, as this would naturally show up in the confidence intervals, whether produced by

<footnote n="1">¹Toulouse School of Economics, University of Toulouse Capitole, Toulouse, France</footnote>

<footnote n="2">²Department of Mathematics & CNRS, University of Angers, Angers, France</footnote>
asymptotics or bootstrapping. We think that this would be a useful complement to the very interesting approach constructed here.

*Conflicts of interest:* None declared.

## Funding

This research was supported by the French National Research Agency under the grants ANR-19-CE40-0013 (ExtremReg project), ANR-23-CE40-0009 (EXSTA project), ANR-17-EURE-0010 (EUR CHESS) and ANR-11-LABX-0020-01 (Centre Henri Lebesgue). A. Daouia and G. Stupfler acknowledge financial support from the TSE-HEC ACPR Chair “Regulation and systemic risks”. G. Stupfler acknowledges further support from the Chair Stress Test, RISK Management and Financial Steering of the Foundation Ecole Polytechnique.

## References

* Daouia A., Noh H., & Park B. U. (2016). Data envelope fitting with constrained polynomial splines. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 78(1), 3–30. https://doi.org/10.1111/rssb.12098

* Jirak M., Meister A., & Reiss M. (2014). Adaptive function estimation in nonparametric regression with one-sided errors. *Annals of Statistics*, 42(5), 1970–2002. https://doi.org/10.1214/14-AOS1248

* Wang L., Xue L., & Yang L. (2020). Estimation of additive frontier functions with shape constraints. *Journal of Nonparametric Statistics*, 32(2), 262–293. https://doi.org/10.1080/10485252.2020.1721494

https://doi.org/10.1093/jrsssa/qnaf018
Advance access publication 10 February 2025

# Dani Gamerman's contribution to the Discussion of 'the Discussion Meeting on the Analysis of citizen science data'

Dani Gamerman Sr. ID

Department of Statistical Methods, Institute of Mathematics, Universidade Federal do Rio de Janeiro,
Rio de Janeiro, RJ, Brazil

*Address for correspondence:* Dani Gamerman, Department of Statistical Methods, Institute of Mathematics, Universidade Federal do Rio de Janeiro, Av. Athos da Silveira Ramos, 149 Cidade Universitária, Rio de Janeiro, RJ 21941-909, Brazil. Email: dani@im.ufrj.br

This contribution discusses the paper by Dennis et al. (2025). The RSS must be congratulated for this timely discussion session, with three great papers.

The model-based approach adopted in this paper is welcomed. It was particularly pleasing to see (observed and unobserved) temporal components presented and discussed in this paper. This paper will certainly be a useful reference for the analysis of archaeologically dated ancient remains we are currently contemplating. My contribution to the discussion addresses the possible extension of the ideas of this paper to temporal sequences of infinite-dimensional observational processes.

Locations of occurrences are as relevant as their counts for point pattern analyses. Inhomogeneous Poisson processes (IPP) are frequently used as their observational distribution. Their intensity functions are based on two components, as in the paper: one for explaining occurrences and one for explaining detection of occurrences. This was implemented for presence-only analyses in Moreira and Gamerman (2022), also using Bayes, but with no model approximation. Prediction of the locations
and the total number of unobserved occurrences are easily provided. This prediction was obtained in the multidisciplinary study of the impact of pre-Columbian occupation of Amazonia (Peripato et al., 2023), and it attracted a sizeable media attention (see Gamerman, 2024). Could the approach of this paper be adapted to such scenario or to its temporal generalization when marks (e.g. dates) are associated with sequences of point patterns (e.g. archaeological remains)?

I will not dwell much over computational aspects other than sharing the concern of the authors about possible errors associated with approximations to the likelihood and models. I appreciate the efforts reported in Table 1. Could the results be expanded to include comparisons against the respective true posteriors (e.g. using Markov Chain Monte Carlo as a gold standard), other parameters, and their (co)variances? The analyses of Amazonia reported results about the true posterior despite the IPP likelihood intractability and despite having to work with databanks with millions of rows.

There are many nice features in the paper by Dennis et al. (2025) and I congratulate the authors for all their comprehensive and hard work.

*Conflicts of interest:* none declared.

## References

Dennis, E. B., Diana, A., Matechou, E. & Morgan, B. J. T. (2025). Efficient statistical inference methods for assessing changes in species' populations using citizen science data. *Journal of the Royal Statistical Society Series A: Statistics in Society, 188*(3), 641–657. https://doi.org/10.1093/jrsssa/qnae105

Gamerman, D. (2024). Uncovering Amazonia: How statistics dismantled a myth. *Significance*, 21(5): 6–10. https://doi.org/10.1093/jrssig/qmae071

Moreira, G. A., & Gamerman, D. (2022). Analysis of presence-only data via exact Bayes, with model and effects identification. *The Annals of Applied Statistics*, 16(3), 1848–1867. https://doi.org/10.1214/21-AOAS1569

Peripato, V., Levis, C., Moreira, G. A., Gamerman, D., Aragao, L. E. O. C., et al. (2023). More than 10,000 pre-Columbian earthworks are still hidden throughout Amazonia. *Science*, 382(6666), 103–109. https://doi.org/10.1126/science.ade2541

https://doi.org/10.1093/jrsssa/qnaf010

Advance access publication 11 February 2025

# Raphaël Huser and Andrew Zammit-Mangion's contribution to the Discussion of the 'Discussion Meeting on the Analysis of citizen science data'

Raphaël Huser<fnref n="1" /> and Andrew Zammit-Mangion<fnref n="2" />

We are grateful for the opportunity to discuss the paper entitled ‘Extreme-value modelling of migratory bird arrival dates: Insights from citizen science data’ (Koh & Opitz, 2025). We congratulate the authors for conducting such an elaborate data analysis of bird arrivals showing

<footnote n="1">¹Statistics Program Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia</footnote>

<footnote n="2">²School of Mathematics and Applied Statistics, University of Wollongong, Wollongong, NSW 2522, Australia</footnote>

<footnote>*Address for correspondence:* Raphaël Huser, Statistics Program Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia. Email: raphael.huser@kaust.edu.sa</footnote>
how multiple citizen-science data products of varying quality, completeness, and space-time res-
olutions, can be fused together and modelled jointly using a complex Bayesian hierarchical model.
The model implements different likelihood functions over a shared latent structure, and accounts
for the sampling effort to remove systematic biases.

Our comment focuses on the methods for making inference, and how alternative approaches
can be used to substantially reduce the computational cost. The authors perform Bayesian infer-
ence using a customized Metropolis-Hastings algorithm, where hyperparameters are updated us-
ing Gibbs sampling, and where the latent Gaussian components are updated using block proposals
based on the Metropolis-adjusted Langevin algorithm. They report that each iteration takes ap-
proximately 6 sec., so a single model fit obtained by drawing 80 000 posterior samples takes about
5.5 days. While such an inference time may be acceptable if the model needs to be fitted only a few
times, it is clearly too expensive if repeated fitting is necessary, e.g., to conduct an extensive cross-
validation study, to analyse different bird species in different study regions, or to incorporate new
data, as is often required in citizen science. Recent years have seen the emergence of neural net-
works being used for amortized inference (see Zammit-Mangion et al., 2025, for a recent review),
both in the context of point estimation (e.g. Sainsbury-Dale et al., 2023, 2024) and full posterior
inference (e.g., Radev et al., 2022; Radev, Schmitt, Pratz, et al., 2023; Radev, Schmitt,
Schumacher, et al., 2023) for models with intractable or unavailable likelihood functions. With
such neural ‘amortized’ approaches, inference with new data can be performed repeatedly at a
fraction of the time needed using standard Markov chain Monte Carlo methods, after an initial
computational cost is incurred to train a neural network using training data simulated from the
model of interest. Amortized methods for complex spatio-temporal Bayesian hierarchical models
are still in their infancy, but see Zammit-Mangion and Wikle (2020) for an early example.
Although neural inference methods can be used to estimate model parameters, they will likely
struggle to make inference on latent variables, which are often numerous and highly correlated
a posteriori. One solution is to adopt an empirical Bayes approach, where hyperparameters are
estimated in a first step (e.g., using neural Bayes estimators; see Sainsbury-Dale et al., 2024),
and latent variables are then inferred in a second step conditional on the estimated hyperpara-
meters, perhaps by using a learning network such as that proposed by Liu and Liu (2020). This
strategy remains to be tested, especially in cases where the model has a relatively large number
of hyperparameters (i.e., more than 20, say) and when the trained neural networks need to be
adaptable to small model changes (e.g., to different covariates values). In principle, however,
this strategy would allow for fast inference for a wide range of Bayesian hierarchical models, in-
cluding models for data fusion, thus enabling the analysis of citizen-science data at unprecedented
scale.

***Conflicts of interest:*** None declared.

References

Koh J., & Opitz T. (2025). Extreme-value modelling of migratory bird arrival dates: Insights from citizen science data. *Journal of the Royal Statistical Society: Series A*, **188**(3), 674–699. https://doi.org/10.1093/jrsssa/qnaf012

Liu L., & Liu L. (2020). Localizing and amortizing: Efficient inference for Gaussian processes. In S. J. Pan & M. Sugiyama (Eds.), *Proceedings of The 12th Asian Conference on Machine Learning* (Vol. 129, pp. 823–836). PMLR.

Radev S. T., Mertens U. K., Voss A., Ardizzone L., & Köthe U. (2022). BayesFlow: Learning complex stochastic models with invertible neural networks. *IEEE Transactions on Neural Networks and Learning Systems*, **33**(4), 1452–1466. https://doi.org/10.1109/TNNLS.2020.3042395

Radev S. T., Schmitt M., Pratz V., Picchini U., Köethe U., & Buerkner P.-C. (2023). JANA: Jointly Amortized Neural Approximation of complex Bayesian models. In *Proceedings of the 39th Conference on Uncertainty in Artificial Intelligence (UAI 2023)* (Vol. 216, pp. 1695–1706). PMLR.

Radev S. T., Schmitt M., Schumacher L., Elsemüller L., Pratz V., Schälte Y., Köthe U., & Bürkner P.-C. (2023). BayesFlow: Amortized Bayesian workflows with neural networks. *Journal of Open Source Software*, **8**(89), 5702. https://doi.org/10.21105/joss

Sainsbury-Dale M., Richards J., Zammit-Mangion A., & Huser R. (2024). Neural Bayes estimators for irregular spatial data using graph neural networks. *Journal of Computational and Graphical Statistics*.
Sainsbury-Dale M., Zammit-Mangion A., & Huser R. (2024). Likelihood-free parameter estimation with neural Bayes estimators. *The American Statistician*, 78(1), 1-14. https://doi.org/10.1080/00031305.2023.2249522

Zammit-Mangion A., Sainsbury-Dale M., & Huser R. (2025). Neural methods for amortized inference. *Annual Reviews of Statistics and Its Application*, 12. https://doi.org/10.1146/annurev-statistics-112723-034123

Zammit-Mangion A., & Wikle C. K. (2020). Deep integro-difference equation models for spatio-temporal forecasting. *Spatial Statistics*, 37, 100408. https://doi.org/10.1016/j.spasta.2020.100408

https://doi.org/10.1093/jrsssa/qnaf012
Advance access publication 11 February 2025

# Kuldeep Kumar's contribution to the Discussion of 'the Discussion Meeting on the Analysis of citizen science data'

Kuldeep Kumar

Bond University, Australia

Address for correspondence: Kuldeep Kumar, Bond University, Australia. Email: kkumar@bond.edu.au

This contribution discusses the papers of 'Extreme-value modelling of migratory bird arrival dates: Insights from citizen science data' by Jonathan Koh and Thomas Opitz.

First of all, I would like to congratulate the authors for modelling such a complex issue as migratory bird arrival using citizen science data. As mentioned by the authors, citizen science data often suffers from inconsistent quality, including a lack of standardization. The modelling becomes more complex because the response variables include Generalized Extreme Value distributions, as well as values distributed according to Binomial and Poisson distributions. Migratory bird arrival dates can be influenced by the interaction of various factors, reflecting the complex interplay between environmental conditions, biological processes, and ecological interactions. I am unsure how many variables the authors have considered; naturally, the model will become more complex if additional variables are included. There are several cutting-edge machine learning models, such as Artificial Neural Networks and Stochastic Gradient Boosting, which could be applied to extreme value modelling. The validity of the data in these models can be tested by dividing the data into training and testing sets, which is particularly crucial in the context of citizen science data. I am curious whether the authors have compared the validity of their predictions using Bayesian hierarchical models with these machine learning models.

Conflicts of interest: none declared.

https://doi.org/10.1093/jrsssa/qnaf014
Advance access publication 11 February 2025
# Subhash Lele's contribution to the Discussion of 'the Discussion Meeting on the Analysis of citizen science data'

## Subhash Lele

Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, AB, Canada T6G 2R3

*Address for correspondence:* Subhash Lele, Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, AB, Canada T6G 2R3. Email: slele@ualberta.ca

I applaud Drs. Bersson and Hoff for the paper 'Frequentist Prediction Sets for Species Abundance using Indirect Information' for bringing to attention data collected by citizen scientists. These data are particularly important because they increase the spatial and temporal extent for many bird monitoring surveys conducted by trained researchers. I have a couple of concerns and a suggestion.

1. Detection error occurs when a bird species that is present is not detected because of various factors such as season, time of the day, rainy day, and such. Correcting for detection error is a major issue in ecological studies. How do authors take into account detection error?

2. Structural zeros occur when a species may not be present at a location due to specific habitat characteristics. How do authors take into account such structural zeros? In their presentation, the number of species K is assumed to be the same in every location. This assumption does not seem realistic.

3. I would like to draw attention to a hierarchical model approach described in Lele and Allen (2006) that calibrates expert opinion or soft data such as collected by E-bird with hard observations collected by trained researchers. Florio et al. (2004) combined satellite data with large spatial extent and ground truth data with smaller spatial extent using a similar approach. In Lele and Allen (2006) formulation, the weight of the soft data in the analysis is automatically determined. As was illustrated, sometimes soft data can, in fact, deteriorate the quality of inference. The moral was that not all soft data are informative and hence should not be used without checking its usefulness and relevance. Could one conduct similar analysis to get prediction sets that combine not just the citizen science data from neighbouring counties but also data systematically collected by bird researchers?

*Conflicts of interest:* none declared.

## References

Florio, E. N., Lele, S. R., Chi Chang, Y., Sterner, R., & Glass, G. E., (2004). Integrating AVHRR satellite data and NOAA ground observations to predict surface air temperature: A statistical approach. *International Journal of Remote Sensing, 25*(15), 2979–2994. https://doi.org/10.1080/01431160310001624593

Lele, S. R., & Allen, K. L. (2006). On using expert opinion in ecological analyses: A frequentist approach. *Environmetrics, 17*(7), 683–704. https://doi.org/10.1002/env.786
# Allan Reese's contribution to the Discussion of the 'Discussion Meeting on the Analysis of citizen science data'

## R. Allan Reese

Forston, Dorset, UK

Address for correspondence: R. Allan Reese, Forston, Dorset, UK. Email: allan@eurolyme.uk

Koh and Opitz present an impressive piece of modelling that covers a far larger geographical area than our own reporting of first-arrival dates (Reese & Tucker, 2019), which by comparison equates to some 9 by 9 20 km pixels. Our data, however, covered more than 130 years (1880s–2010s). We speculated on sources of bias, some of which are corroborated in the new work, and how these may have changed over time. Koh and Opitz’s work would usefully be made more accessible to ornithologists, with the mathematics simplified and more biological interpretation.

One minor point is the coding of dates. The crude day number within a year introduces a small anomaly: each date after February is larger by one in leap years. This effect is easily removed by coding March 1 as Day 1, using negative values for January and February—which do not anyway occur in these data. The authors could also check for ‘weekend bias’, which we found in more recent (post-1950) reports.

Most studies, including our own, consider first arrival dates in a strict sense: ‘the first Cuckoo of spring’. The definition here, ‘first arrivals of migratory birds at their destination’ is more obviously a latent variable. Many migrants may stop in a suitable habitat and, for example, sing before continuing their journey. The second definition, ‘a date that can be viewed as representative of the arrival of birds at their breeding site’, smacks more of a population average than the earliest pioneer. The span of dates in Figure 1 may reflect sampling variation, but also the large geographical spread. Surely one should be modelling the contours of the northward progress? The pixels appear to be coded with a single index ($i = 1,...,1,268$) rather than in rows successively northward.

Our paper looked briefly at ‘effort’ based on the estimated number of observers, with equivocal results. For me, Figure 1 needed some decoding, and I would advise revising the labels and the caption. Values of $log(N)$—I assume base e—do not easily convey a magnitude: the value of 10 equates to 22,000 observations in a pixel-year. I zoomed in on the PDF, and the observations appear to extend to at least 9 (8,000). If these were spread randomly across the year (e.g. not just wildfowl counts in winter), there is an extremely small Poisson probability of a day without data. Labelling with natural numbers, including the minimum and maximum number of observations, would be more understandable. Nor is it obvious why the Y-axis is inverted. Plotting the day number ascendant leads to the more natural description: more observations in an area lead to earlier dates with smaller variance. I would try a curve with a single inflexion and an asymptote.

A final point is that the interest in first dates is generally in relation to change. If the sources of bias and the direction of effect are largely unchanged, then changes in arrival dates have importance and need interpretation.

*Conflicts of interest:* none declared.

## Reference

Reese R. A., & Tucker J. J. (2019). Descriptive models for first arrival dates of migrants in an inland UK county. *Ringing & Migration*, 34(1), 63–69. https://doi.org/10.1080/03078698.2019.1759896
# Stefano Rizzelli's contribution to the Discussion of the 'Discussion Meeting on the Analysis of citizen science data'

Stefano Rizzelli ID

Department of Statistics, University of Padova, Padova, Italy

Address for correspondence: Stefano Rizzelli, Department of Statistics, University of Padova, Padova, Italy.

Email: stefano.rizzelli@unipd.it

This contribution discusses the paper 'Extreme-value modelling of migratory bird arrival dates: Insights from citizen science data' by Koh and Opitz, whom I congratulate on this work. Of particular interest is their use of Generalized Extreme Value distributions for modelling maxima $Z_i$ of (transformed) birds arrival times, performing parameter inference, and point prediction through a Bayesian approach. One of its advantages is the possibility of incorporating field knowledge through prior distributions. Since the variables $Z_i$ plausibly have an upper bound, it would be interesting to incorporate this information by assigning $\xi$ a prior supported on negative values, and to assess the sensitivity of the proposed method in moving from a vague to an informative prior specification. Another advantage is the possibility of issuing forecasts of maxima through predictive distributions, as remarked in Coles and Pericchi (2003) and Padoan and Rizzelli (2022, 2024). Expanding the analysis in §5.3 to include predictive median (or mode), checking how much predictive distributions concentrate around it and whether predictive intervals cover truly observed arrival times, would also be interesting.

As the authors remark, fitting models grounded on extreme value theory (EVT) to extremes of biotic processes is a rare practice, which may seem surprising given the relevance of extreme events in life sciences. Besides EVT primary focus on continuous variables, this may be due to a relative inflexibility of the more established extreme value models, arising from EVT fundamental theorems (de Haan & Ferreira, 2006). While they allow to successfully account e.g. for spatial dependence (Davison et al., 2012), there are as yet no commonly agreed-upon assumptions to derive extreme value models in the presence of covariates. The need for extensions of classical EVT in this sense is evidenced by the plethora of applications for which traditional extreme value models have been expanded in different ways, according to problem-specific desiderata, without, however, dwelling with probabilistic foundations of such extensions and mathematical guarantees on the proposed methodologies, which is far from trivial but crucial to establish how reliably they can be used in other areas. Some efforts in this direction have been made by Bobbia et al. (2021), Chernozhukov (2005), Daouia et al. (2022), and Dombry et al. (2023), among others. Rigorous attempts to also account for time heterogeneity can be found in Bücher Jennessen (2024) and Einmahl et al. (2016, 2022).

Recent progress is remarkable, but there is still a long way to go to develop and consolidate rich, interpretable and mathematically sound models. Bridging the gap between theory and practice remains a key challenge for the advancement of extreme value analysis, making its toolbox increasingly usable and trustworthy for practitioners among different fields.

*Conflicts of interest:* The author declares to have no conflict of interest.

## References

Bobbia B., Dombry C., & Varron D. (2021). Extreme quantile regression in a proportional tail framework. *Transactions of A. Razmadze Mathematical Institute*, 175, 13–32. https://rmi.tsu.ge/transactions/TRMI-volumes/175-1/v175(1)-2.pdf

Bücher A., & Jennessen T. (2024). Statistics for heteroscedastic time series extremes. *Bernoulli*, 30(1), 46–71. https://doi.org/10.3150/22-BEJ1560
Chernozhukov V. (2005). Extremal quantile regression. *Annals of Statistics*, 33(2), 806-839. https://doi.org/10.1214/009053604000001165

Coles S., & Pericchi L. (2003). Anticipating catastrophes through extreme value modelling. *Journal of the Royal Statistical Society: Series C, Applied Statistics*, 52(4), 405-416. https://doi.org/10.1111/1467-9876.00413

Daouia A., Gijbels I., & Stupfler G. (2022). Extremile regression. *Journal of the American Statistical Association*, 117(539), 1579-1586. https://doi.org/10.1080/01621459.2021.1875837

Davison A. C., Padoan S. A., & Ribatet M. (2012). Statistical modeling of spatial extremes. *Statistical Science*, 27(2), 161-186. https://doi.org/10.1214/11-STS376

de Haan L., & Ferreira A. (2006). *Extreme value theory: An introduction*. Springer.

Dombry C., Padoan S. A., & Rizzelli S. (2023). *Asymptotic theory for Bayesian inference and prediction: From the ordinary to a conditional peaks-over-threshold method*. arXiv eprint 2310.06720, available at: https://arxiv.org/abs/2310.06720.

Einmahl J. H. J., de Haan L., & Zhou C. (2016). Statistics of heteroscedastic extremes. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 78(1), 31-51. https://doi.org/10.1111/rssb.12099

Einmahl J. H. J., Ferreira A., de Haan L., Neves C., & Zhou C. (2022). Spatial dependence and space-time trend in extreme events. *Annals of Statistics*, 50(1), 30-52. https://doi.org/10.1214/21-AOS2067

Padoan S. A., & Rizzelli S. (2022). Consistency of Bayesian inference for multivariate max-stable distributions. *Annals of Statistics*, 50(3), 1490-1518. https://doi.org/10.1214/21-AOS2160

Padoan S. A., & Rizzelli S. (2024). Empirical Bayes inference for the block maxima method. *Bernoulli*, 30(3), 2154-2184. https://doi.org/10.3150/23-BEJ1668

https://doi.org/10.1093/jrsssa/qnaf020

Advance access publication 11 February 2025

# Maozai Tian's contribution to the Discussion of the 'Discussion meeting on the Analysis of citizen science data'

Maozai Tian

Center for Applied Statistics, School of Statistics, Renmin University of China, Beijing 100872, China

Address for correspondence: Maozai Tian, Center for Applied Statistics, School of Statistics, Renmin University of China, Beijing 100872, China. Email: mztian@ruc.edu.cn

## 1 Comments and suggestions

I enthusiastically congratulate the authors on a stimulating article that provides a cutting-edge technology for implementing data fusion of heterogenous datasets and may inspire many follow-up works.

The paper develops a spatiotemporal Bayesian hierarchical model for bias-corrected estimation of arrival dates of the first migratory bird individuals at a breeding site, devises an MCMC scheme, and checks by simulation that the latent process components are identifiable. The proposed model is used to several migratory bird species and resulted in an interesting finding that the sampling effort significantly modulates the observed first arrival date. The relationship is effective in making bias-correct predictions of the true first arrivals.

Here are three comments:

First, the article develops a spatiotemporal Bayesian hierarchical model with two levels. The level-1 model is with the hierarchical system of the four regression equations in Section 3.3 on
pages 7 and 8,

$$ 
\left\{ 
\begin{array}{l} 
N_j^{\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}), \xi \} 
\end{array} 
\right. 
\quad (1) 
$$

At level-2 model, the level-1 hyperparameters $\lambda^{\text{BBS}}$, $\lambda^{\text{ckl}}$, $p^{\text{spc}}$, $\mu$, and $\sigma$ become outcomes:

$$ 
\left\{ 
\begin{array}{l} 
\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_0^{\text{act}}}{d_{s_i, t_i}} \\ 
\mu(s_i, t_i) = g \left\{ \beta_0^{\text{GEV-}\mu} + X^{\text{GEV-}\mu}(s_i) + \beta_1^{\text{GEV-}\mu} \text{NAO}_{t_i} \right. \\ 
\left. + \theta^{\text{niche-GEV}} X^{\text{niche}}(s_i), x_{\text{effort}}(s_i, t_i) \right\} \\ 
\log \sigma(s_i) = \beta_0^{\text{GEV-}\sigma} + X^{\text{GEV-}\sigma}(s_i). 
\end{array} 
\right. 
\quad (2) 
$$

In practice, the existence of such citizen science data with hierarchies is neither accidental nor ignorable, it is a common phenomenon. To ignore this, hierarchical data structure risks overlooking the importance of group effects and may also render invalid many of the traditional statistical analysis techniques used for studying data relationships. On the other hand, hierarchical models take hierarchical data structure into account and have also many applications in statistics. However, the aforementioned level-2 model fails in taking into account the inherent association among the four regression equations at level-1. So, it needs to impose, for example, some strict constraints on the joint distribution of the latent effect vector ($X^{\text{niche}}(s_i)$, $X^{\text{year}}(t_i)$, $X^{\text{pref}}(s_i)$, $X^{\text{niche}}(s_i)$, $X^{\text{GEV-}\mu}(s_i)$, $X^{\text{GEV-}\sigma}(s_i)$). Obviously, it is not suitable to assume that these latent effects be assigned independent and identical distribution, see, for example, Amini et al. (2021), Li et al. (2014), Ma and Tian (2024), Rotejanaprasert et al. (2023), Yu et al. (2022), and Tian and Yu (2024).

Second, what are the conditions for the identifiability of the model? There are seven hyperparameters $\beta_0^{\text{BBS}}$, $\beta_0^{\text{ckl}}$, $\beta_0^{\text{spc}}$, $\beta_0^{\text{act}}$, $\beta_0^{\text{GEV-}\mu}$, $\beta_1^{\text{GEV-}\mu}$, and $\beta_0^{\text{GEV-}\sigma}$ as the fixed effects to estimate. Whether so many hyperparameter estimators will lead to the superposition of errors and the failure of the method. Refer to Yu et al. (2022) and Tian and Yu (2024).

Third, how to monitor the convergence of Markov Chain Monte Carlo algorithm for the proposed spatiotemporal Bayesian hierarchical model at each step? See, for example, Gnecco et al. (2024).

## Acknowledgments

The authors thank the RSS Journals Manager Judith Shorten for her kind invitation to make presentation with discussion of the paper. Prof. Maozai Tian's work was supported by the Natural Science Foundation of Beijing (No. 1242005).

*Conflicts of interest:* None declared.

## References

Amini P., Moghimbeigi A., Zayeri F., Tapak L., Maroufizadeh S., & Verbeke G. (2021). Longitudinal joint modelling of ordinal and overdispersed count outcomes: A bridge distribution for the ordinal random intercept. *Computational and Mathematical Methods in Medicine*, 2021(1), 5521881. https://doi.org/10.1155/2021/5521881
Gnecco N., Terefe E. M., & Engelke S. (2024). Extremal random forests. *Journal of the American Statistical Association*, 119(548), 3059-3072. https://doi.org/10.1080/01621459.2023.2300522

Li Q., Pan J., & Belcher J. (2014). Bayesian inference for joint modelling of longitudinal continuous, binary and ordinal events. *Statistical Methods in Medical Research*, 25(6), 2521-2540. https://doi.org/10.1177/0962280214526199

Ma S. P., & Tian M. Z. (2024). A censored quantile transformation model for Alzheimer's disease data with multiple functional covariates. *Journal of the Royal Statistical Society Series A: Statistics in Society*, 188(2), 515-538. https://doi.org/10.1093/jrsssa/qnae061

Rotejanaprasert C., Lawpoolsri S., & Sa-angchai P. (2023). Projecting malaria elimination in Thailand using Bayesian hierarchical spatiotemporal models. *Scientific Reports*, 13(1), 7799. https://doi.org/10.1038/s41598-023-35007-9

Tian M. Z., & Yu K. M. (2024). Maozai Tian and Keming Yu's contributions to the Discussion of “Safe Testing” by Peter Grüunwald, Rianne de Heide, and Wouter Koolen. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 86(5), 1160. https://doi.org/10.1093/jrsssb/qkae076

Yu Z., Yu K. M., Härdle W. K., Zhang X. L., Wang K., & Tian M. Z. (2022). Bayesian spatio-temporal modeling for inpatient hospital costs of alcohol-related disorder. *Journal of the Royal Statistical Society Series A: Statistics in Society*, 185(2), 644-667. https://doi.org/10.1111/rssa.12963

https://doi.org/10.1093/jrsssa/qnaf021

Advance access publication 10 February 2025

# Authors' reply to the Discussion of 'Efficient statistical inference methods for assessing changes in species' populations using citizen science data' at the 'Discussion meeting on the analysis of 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, CT2 7NF Kent, UK. Email: B.J.T.Morgan@kent.ac.uk

Emily Dennis is unable to comment on the issues raised by the Discussion due to maternity leave.

We are grateful for the wide-ranging comments on our paper, which have given us encouragement and provide pointers to future research.

We have focused on providing analyses of data provided by *Butterfly Conservation* in the UK. For both abundance and distribution data, analyses are repeated regularly on sets of 59 species of butterflies and a subset of the approximately 900 species of larger moths. What we need therefore are efficient, flexible methods, where citizen-science records typically contain minimal covariate information.

## Ben Swallow

We thank Ben Swallow for his wide-ranging comments covering citizen-science data, and for the particular comments relating to our paper.

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

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

<footnote n="3">³School of Mathematics, Statistics and Actuarial Science, University of Essex, Wivenhoe Park, Colchester, CO4 3SQ Essex, UK</footnote>
He mentions the importance of modelling the observation process, which resonates with the paper by Moreira and Gamerman (2022), mentioned below. We note here one way in which this is not possible in comparison with planned surveys conducted by known individuals. In the latter case, there is often information available regarding the observers themselves, for instance regarding age and experience, and these are typically not available with citizen-science data.

As we say, the extended Generalised Abundance Index (GAI), incorporating the annual model, is at a preliminary stage of development. However, the examples that we have shown are promising, and we are in the process of undertaking a wider study. The annual model is much used in ecology and we believe that our new approach should have wide-ranging applications. This ties in also with the comments of Ben Baer and our responses to these given below. It is not uncommon for fitted models in ecology to result in parameter estimates which are then analysed further, for example using Generalised Linear Models (GLMs). In order then to account for error propagation some form of time-consuming bootstrapping is used. Differences in the estimation of uncertainty are mentioned, requiring more research, however incorporating an annual model within the GAI avoids the arbitrariness of a post hoc analysis which might, for example, depend upon the use of a particular GLM.

The attraction of the GAI family of models is that it allows for a range of different discrete distributions for the analysis of count data, and for comparisons to be made using information criteria—see Dennis et al. (2016). In addition relevant covariates, such as northing, easting and relating to land cover, are easily included. Northing and easting can reveal spatial features; an illustration is given in Dennis et al. (2022), showing how the Common Blue butterfly, *Polyommatus icarus*, moves from having two broods in the south of Britain to having only one in the north, where the colder climate only allows one brood to be completed each year. More complex inclusion of spatial covariates results in the occupancy maps in Diana et al. (2023), which are reproduced in the Discussion contribution by Jafet and coauthors using Integrated nested Laplace approximations (INLA), which also incorporate spatial correlation.

The possibility of multi-species modelling is a good one, and applications are mentioned in McCrea et al. (2023) and Moreira and Gamerman (2022). The benchmarking used in the paper to produce presence-absence data uses observations on species other than target species to induce absences, and so is a particular type of multi-species analysis. UK butterflies are generally classified as habitat specialists, migrants, and wider countryside; specialists are only found where the habitat is favourable to them. This point is relevant to the utility of multi-species analyses, related to the construction of indicators which describe how groups of individual species behave. Indicators are widely used to summarize complex, multi-species analyses. However, indicators are only useful if they describe relatively homogeneous groups of species. Methods of cluster analysis and ordination are used to check this, and illustrations using functional principal component analysis are in Dennis et al. (2019), which illustrate in a two-dimensional plot the relationships between butterfly occupancy indices. It is an interesting point that modelling of data on specialists might benefit from appropriate analyses of data from wider countryside species with overlapping ranges.

## Kerrie Mengersen

We welcome Kerrie’s mention of relevant work in other areas and also the citizen-science research taking place in the antipodes. We very much endorse the need for rigor in model fusion, emphasized also in McCrea et al. (2023). As a comment on integrated modelling, in our paper, we analyse both count and distribution data on British butterflies, but separately, using different models in each case as abundance and distribution are inherently different. Typically data in the first case relate to individuals while in the second case they refer to species. As we mention in the Discussion to our paper, we are developing a new urban butterfly indicator based on an integrated analysis of UKBMS and BBC data, building on the material of Section 3 of the paper.

Citizen-science data are increasing in size, spatio-temporal coverage and species monitored, especially with the emergence of new technologies, such as cameras, acoustic, and DNA-based approaches. Therefore, they provide considerable opportunities for species’ monitoring and timely intervention. At the same time, the unstructured, often opportunistic and biased process by which such data are collected presents modelling challenges, many of which remain open. For the data sets considered in our paper, reporting preferences, observer differences and spatial bias
(Johnston et al., 2023) are just some of the known and unavoidable biases in the data generating process. To maximize the benefit of citizen-science data for reliable and timely monitoring, these biases have to be addressed within a modelling framework, and we agree with Kerrie that more work is needed in the area.

Machine learning methods, especially deep models, are some of the most flexible and scalable tools currently available to us, and therefore it is indeed beneficial to incorporate them within statistical models of the type we present in our paper. For example, deep models can be used in place of logistic regression models to increase the flexibility of our framework. Variational inference approaches, such as the one we developed, provide a natural way forward for this task, thanks to the now widely available automatic differentiation tools.

## Andrej Srakar

We thank Andrej for the interesting and thought-provoking ideas.

Our use of concentrated likelihood is a simple way of reducing the effective size of parameter spaces for particular survey data. More generally, the approach has been found to be useful in time-series analysis for example; see Durbin and Koopman (1997).

As we have shown, complex models in ecology may well have many parameters, and one way to deal with this is to use appropriate random effects. Dennis (2015) compared this approach with using concentrated likelihood and found the latter to be preferable for our models. In his contribution to the paper Discussion, Ben Baer both simplified the expression for a likelihood and computed a dimension reduction in one application from 527588 to 1222; further discussion with him has shown that additional concentration of the likelihood in this case allows us to go even further. Thus for the GAI-based models in particular, and the data being analysed, we can probably achieve the efficiency that we need at present. In terms of efficiency, Ben also suggests using stochastic gradient descent (which is used in the paper), and we need to look into the value of that for the GAI-based models. The potential of using Banach, Hilbert, or tensor spaces for our work is as yet unknown. However data and models are becoming more extensive and complex—McCrea et al. (2023)—and these suggestions may well gain traction.

In some of our work, likelihoods are not intractable, when methods of indirect inference may not be the best approach. Approximate Bayesian Computation, ABC, now has a wide literature—see for example Beaumont (2019) and Gutmann and Corander (2016) as an illustration of efficiency in computing. Each application requires its own simulators and data-matching functions, in both cases with the potential for time-consuming computation. One of the drivers for our work has been the need for methods to analyse data arising from established surveys, which can be used by ecologists with limited computational power. The main issues with ABC are that it relies on identifying summary statistics, which can be challenging with large spatio-temporal data and models with many parameters, of the type we consider here. ABC also does not scale well as the number of parameters increases (the curse of dimensionality, Cranmer et al., 2020). Instead, we think that when the likelihood is tractable, as is the case for the data and corresponding models of this paper then likelihood-based inference and efficient computational tools, such as variational inference, hold considerable promise. However, we expect that Andrej's suggestions are highly relevant for the new data that are regularly appearing, as with the use of autonomous machines for data collection for example.

In our paper, we have used both classical inference and variational inference. The two classical inference illustrations are easily performed and described. In general, the value of Bayesian methods for certain data analysis is evident from examples: for instance the variational inference illustration and also Diana et al. (2023), combining a Pólya-Gamma scheme and spatio-temporal random effects using Gaussian processes, the exact Bayes method of Moreira and Gamerman (2022), the extension of INLA to occupancy problems by Belmont et al. (2024), and the sampling approach of King et al. (2023) combined with importance sampling.

## Dani Gamerman

The approach proposed in Moreira and Gamerman (2023) is particularly useful for presence-only data. However, the data and model considered in our paper occupy an intermediate position between presence-only and presence-absence data. While it is true that observations are initially recorded as presence-only data, observers document all species they encounter. This allows absences
to be inferred whenever an observer records a different species but not the one under consider-
ation. This process, known as benchmarking, enables the identification of absences in addition
to presences, distinguishing the data from typical presence-only data. We note that if observers
were to report on only one species, the data would indeed be presence-only. We do recognize, how-
ever, that the data are not pure presence/absence data either, since that would require observers to
explicitly document absences whenever no species is recorded, which is not the case here. We be-
lieve that in this case it is preferable to infer non-detections for each species using records of other
species (assuming that all species share the same reporting probability, i.e. their detections are
equally likely to be reported) than to ignore that information all together and treat the data as
presence-only.

We note with interest the approach of Moreira and Gamerman (2022), of assuming a latent pro-
cess of unobserved occurrences as a way to mitigate preferential sampling. Our approach does not
account for potential preferential sampling, which is a known issue in citizen-science surveys since
observers might tend to visit sites that are more likely to be occupied, which leads to the sampled
sites being a biased sample of all sites. Preferential sampling is a well-known but open problem
when modelling citizen-science data; see Section 3.1 in Johnston et al. (2023), who comment
that ‘... (preferential sampling is) particularly a challenge if the drivers of observer site selection
are aligned with the ecological process of interest, and if there are no variables that describe
this preference in the model.’ To the best of our knowledge, preferential sampling has not been
considered in models of abundance or occupancy of the type we develop and use in this paper—
though see Conn et al. (2017)—so this is certainly an exciting future research direction.

The approach of Moreira and Gamerman (2022) relies on the presence of relevant covariates, and may not apply to various types of citizen-science data. For example in the UK, a recorder using iRecord is required to add a location to a submitted presence entry. We do not think that in general there is information on other covariates. As also observed in Moreira and Gamerman (2022), observer information in citizen science data can be difficult to obtain. Another feature might relate to sample sizes, as data can be very large; see for example Diana et al. (2023).

Finally, on the question regarding (co)variances, simulation results in our paper show that the marginal posterior distributions of parameters from our variational inference approach and Markov chain Monte Carlo (MCMC), are very similar. In response to the question about the posterior covariances, we have run a new set of simulations to compare the posterior covariances of the parameters ($\beta_1^{\prime\prime}$, ..., $\beta_{p_1}^{\prime\prime}$, $\beta_1^p$, ..., $\beta_{p_2}^p$) obtained by variational inference against those ones obtained by MCMC (which should be the true ones). We have used 2 covariates for both detection and occupancy (including the intercept). Across 50 simulations and the 6 posterior covariances, the average relative error was 38%. It is likely that the bias is still due to the low number of replications or to the lack of convergence of the variational inference procedure, since the estimator of the posterior covariance matrix tends to exhibit slower convergence for the off-diagonal elements.

Benjamin R. Baer

The extended GAI presents a way of including the main-effects, annual, model within the standard
GAI framework. This simple idea allows the GAI to be applied in a true dynamic fashion, for ana-
lysis of count data over time. It is preferable to the current approach in which the static GAI is ap-
plied separately to each year, followed by a Poisson GLM in order to obtain time trends, as
explained. Having a dynamic model avoids the need for time-consuming bootstrap analyses in or-
der to deal with the variance propagation resulting from the GLM applying to estimated values. It
is also desirable to remove reliance on a particular GLM. This new approach is preliminary, and
we plan extensive checking on simulated data to investigate performance and properties, building
on the initial good results presented for real data. We thank Ben for his useful development of the
model. We focused on using concentrated likelihood to estimate the $\alpha$ parameters to reduce the
dimensionality of the parameter space for maximum likelihood. Ben has taken this work further,
for a particular submodel in which the phenological function $a(\cdot, \cdot, \cdot)$ is not a function of site. This
restriction is in fact one that is often used in practice—see Dennis et al. (2013). He then proceeds to
derive an analytical form for the likelihood which will also enhance computational efficiency. We
have discussed this work with Ben which has resulted in an additional simplification due to the
concentrated likelihood for the $\beta$ parameters also resulting in explicit parameter estimates. This
in turn provides a further improvement in efficiency which we look forward to examining in detail
in the near future. The point about stochastic steepest descent is an interesting one, as it has in fact
been used in Section 4 of the paper. We shall certainly explore its wider relevance.

Jafet Belmont, Sara Martino, Janine Illian, and Håvard Rue

The flexibility and speed of R-INLA provides considerable modelling opportunities for spatio-
temporal data, such as those considered in our paper. It is interesting to see how observation error,
in the form of the probability of detection given presence in occupancy models—accounting for
false negative observation error—can now be accounted for within R-INLA, and undoubtedly
this will be widely used by practitioners to model their data in the future.

We agree with the reviewers that in the specific case of an occupancy model, using INLA is
more efficient than the variational inference approach. However, we believe that the advantage
of variational inference is its greater flexibility in considering potentially more complex sam-
pling scenarios that arise in citizen-science data, since the approach only requires us to be
able to express the likelihood in closed analytic form (potentially having to integrate over dis-
crete variables if these are present). Moreover, using variational inference, it is possible to take
advantage of the developments on automatic differentiation and potentially model occupancy
and detection effects using deep models, as we discuss in response to Kerrie's comments.
Although in the specific case study considered there are more efficient approaches, as shown
by the reviewers, we think that in the complex observational processes that give rise to different
types of citizen-science data there would be considerable scope to perform Bayesian (variation-
al) inference.

*Conflicts of interest:* None declared.

Funding

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 contributed 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. B.M. was supported by a Leverhulme Emeritus Fellowship. This work was partly funded by NERC grant NE/T010045/1.

References

Beaumont M. A. (2019). Approximate Bayesian computation. *Annual Review of Statistics and Its Application*, **6**, 379–403. https://doi.org/10.1146/annurev-statistics-030718-105212

Belmont J., Martino S., Illian J., & Rue H. (2024). Spatio-temporal occupancy models with INLA. *Methods in Ecology and Evolution*, *15*(11), 1923–2167. https://doi.org/10.1111/mee3.v15.t11

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

Cranmer K., Brehmer J., & Louppe G. (2020). The frontier of simulation-based inference. *PNAS*, *117*(48), 30055–30062. https://doi.org/10.1073/pnas.1912789117

Dennis E. B. (2015). *Development of statistical methods for monitoring insect abundance* [Ph.D. thesis]. University of Kent. https://kar.kent.ac/49079/

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.v12.8

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*, 637–645. https://doi.org/10.1111/2041-210X.12053

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., Roy D. B., & Brereton T. M. (2019). Functional data analysis of multi-species abundance and occupancy data sets. *Ecological Indicators*, *104*(1–2), 156–165. https://doi.org/10.1016/j.ecolind.2019.04.070
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(3), 2503-2515. https://doi.org/10.1111/biom.13816

Durbin J., & Koopman S. J. (1997). Monte Carlo maximum likelihood estimation for non-Gaussian state space models. *Biometrika*, 84(3), 669-684. https://doi.org/10.1093/biomet/84.3.669

Gutmann M., & Corander J. (2016). Bayesian optimization for likelihood-free inference of simulator-based statistical models. *Journal of Machine Learning Research*, 17, 1-47. http://jmlr.csail.mit.edu/papers/volume17/15-017/15-017.pdf

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(1), 103-116. https://doi.org/10.1111/mee3.v14.1

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(4), 3112-3132. https://doi.org/10.1214/23-AOAS1753

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(1), 4-11. https://doi.org/10.1111/mee3.v14.1

Moreira G. A., & Gamerman D. (2022). Analysis of presence-only data via exact Bayes, with model and effects identification. *The Annals of Applied Statistics*, 16(3), 1848-1867. https://doi.org/10.1214/21-AOAS1569

https://doi.org/10.1093/jrsssa/qnaf041
Advance access publication 11 April 2025

# Authors' reply to the Discussion of
'Frequentist prediction sets for species
abundance using indirect information'
at the 'Discussion meeting on the
analysis of citizen science data'

Elizabeth Bersson<fnref n="1" /> ID and Peter D. Hoff<fnref n="2" />

To begin, we would like to thank all discussants for their insightful comments on our article, and congratulate Dennis, Diana, Matechou, Morgan, Koh, and Opitz on their contributions to the discussion meeting that enabled such an enriching conversation on analysis of citizen science data. Much of the discussion comments pertain to how our proposed framework coincides with unique concerns of citizen science data. In what follows, we group common themes that arose in the discussion and respond.

## 1 Assumptions and applicability to citizen science data

Our work presents methodology to construct a prediction set based on avian species abundance count data from a given county. Moreover, we present an empirical Bayesian framework to incorporate indirect information to improve prediction set precision using auxillary data from

<footnote n="1">¹Laboratory for Information and Decision Systems, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="2">²Department of Statistical Science, Duke University, Durham, NC 27701, USA</footnote>

<footnote>*Address for correspondence:* Elizabeth Bersson, Laboratory for Information and Decision Systems, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. Email: ebersson@mit.edu</footnote>
neighbouring counties. The prediction set for a given county will have guaranteed frequentist coverage if two primary assumptions are upheld: The data are exchangeable within the county and are independent from the auxiliary data. To this end, Swallow, Mengersen, Srakar, and Lele each commented on some unique aspects of citizen science data that require particular attention in any statistical methodology, which opens a discussion on the practicality of the assumptions required for our approach.

## 1.1 Within-county exchangeability

For the frequentist coverage guarantee to hold, our proposed method relies on an assumption of exchangeability of samples. Swallow points out that this may be a concern as sampling bias is a warranted concern with citizen science data. Some concerns of biases with citizen science-collected data may be alleviated by pre-processing data such that only entries collected with some sufficient level of effort are used in an analysis. In the eBird database, data entries may be filtered using some of the ample meta-data available regarding the data collection session. For example, utilizing data from complete checklists (eBird, 2020) and from checklists uploaded by frequent contributors can ensure only data with a sufficient level of skill and effort are being utilized in an analysis (Horns et al., 2018).

More generally, for the coverage guarantee of our approach to hold, the sum of the count vectors obtained from each birding session must follow a multinomial distribution when conditioning on the total number of observed birds. This assumption will be violated, for example, if repeated samples are uploaded from individuals in a group who collected data together. These entries will be correlated, and the sum of the corresponding count vectors will not follow a multinomial distribution. In the eBird database, there is meta-data identifying if separate checklists were uploaded from a group birding session, and duplicate lists may be excluded from an analysis using this information. There is recent work that relaxes the exchangeability assumption in conformal prediction (Barber et al., 2023), and extending our method in this direction may be useful for citizen science datasets with more prevalent biases.

## 1.2 Between-county independence

In this work, we are motivated by utilizing auxiliary data from neighbouring counties to improve prediction set precision for a given county. This is accomplished via a hyperparameter that may be estimated from data that is independent from the sample in the county of interest.

As Srakar highlights, geographic adjacency structure is an inherently interesting aspect of much citizen science data. In our approach, we propose an estimation procedure that estimates the hyperparameter with auxiliary data from the nearest neighboring counties. In this way, we allow for information to be shared among regions that are geographically near one another.

Instead, we may aim to share information among counties that exhibit similar ecological patterns in a more nuanced manner, perhaps allowing for greater geographical distance. For example, in some contexts, it may be more useful to share information across counties with similar county-level environmental covariates. Accounting for more complex dependencies between the counties in this way can be integrated into a more heuristic-based estimation procedure for the indirect information. As an example, the distribution of the transformation of counts of observations of a single species can be approximately modelled with a normal distribution with a spatial dependence structure, such as a simultaneous auto-regressive (SAR) model (Singh et al., 2005) that includes a regression with county-level covariates. Then, the estimated hyperparameter may be taken to be an empirical Bayes estimate of the mean from this distribution. This perspective also pertains to Lele's suggestion of incorporation of alternative expert data.

## 1.3 Additional data concerns

In this work, we analyze eBird data from the state of North Carolina in the United States. For inference on each county, we take the categorical support to include any species seen in the entire state during the time frame we analyze. Lele notes this all-encompassing approach may result in structural zeros such that a species may have near zero probability of occurrence in a given county. While the support is the same in each county, the species' occurrence probability vector varies across state, so this is theoretically coherent.
An additional concern with citizen science data, raised by Lele, is detection error that arises when species are not detected due to various extenuating circumstances that affect data collection, including inclement weather, time of day, and others. Some of these concerns may be alleviated by aggregation, however, this solution is dependent upon the inferential question of interest.

## 2 Auxillary information

Lele discusses the nuance around incorporating auxiliary information in an analysis. Although in some settings it may be undesirable to incorporate indirect information that is uninformative, our approach is adaptable to such a scenario. For example, if many species with small observed counts have large hyperparameter values, then this reflects a general uncertainty of knowledge, and, consequently, will result in a large prediction set. Alternatively, if the indirect information is aligned with the sample, then a smaller prediction set will be obtained.

## 3 Coverage probability

Swallow points out that the coverage probabilities can be above the nominal rate, as seen in Figure 2. Indeed, our method guarantees a lower bound on the coverage rate of a prediction set (Eqn 2). The realized coverage rate of the indirect approach depends on, among others, the sample size *N*, as described in Bersson and Hoff (2024). Exact coverage may be achieved by utilizing a randomized confidence set algorithm, see Vovk et al. (2022) §2.2.6. In contrast to the α-validity guaranteed by the indirect approach, a Bayesian or empirical Bayesian prediction set can have arbitrarily poor coverage, depending on the prior hyperparameter.

*Conflicts of interest:* None declared.

## References

* Barber R. F., Candès E. J., Ramdas A., & Tibshirani R. J. (2023). Conformal prediction beyond exchangeability. *Annals of Statistics*, *51*(2), 816–845. ISSN 21688966. https://doi.org/10.1214/23-AOS2276

* Bersson E., & Hoff P. D. (2024). Optimal conformal prediction for small areas. *Journal of Survey Statistics and Methodology*, *12*(5), 1464–1488. https://doi.org/10.1093/jssam/smae010

* eBird (2020). *Birding as your 'Primary Purpose' and complete checklists* (Technical Report). https://support.ebird.org/en/support/solutions/articles/48000967748-birding-as-your-primary-purpose-and-complete-checklists.

* Horns J. J., Adler F. R., & Şekercioğlu Ç H. (2018). Using opportunistic citizen science data to estimate avian population trends. *Biological Conservation*, *221*(5), 151–159. ISSN 00063207. https://doi.org/10.1016/j.biocon.2018.02.027

* Singh B. B., Shukla G. K., & Kundu D. (2005). Spatio-temporal models in small area estimation. *Survey Methodology*, *31*(2), 183–195. https://www.researchgate.net/publication/237110155

* Vovk V., Gammerman A., & Shafer G. (2022). *Algorithmic learning in a random world* (2nd ed.). Springer US.
# Authors' reply to the Discussion of 'Extreme-value modelling of migratory bird arrival dates: insights from citizen-science data'

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

*Address for correspondence:* Jonathan Koh, Seminar for Statistics, ETH Zurich, Rämistrasse 101, 8092 Zurich, Switzerland. Email: jonathan.koh@stat.math.ethz.ch

## Abstract

We respond to the discussion comments of the Proposer and Seconder of the vote of thanks, and to the eight other contributions discussing our work.

We warmly thank the Proposer and Seconder of the vote of thanks, and all discussants. We appreciate their numerous observations, suggestions, and questions on statistics for citizen science (CS) data in general, and their insightful and encouraging comments on our work in particular.

The Proposer raises perspicacious points on the challenges of working with CS data, highlighting the inherent structural biases and lack of appropriately randomized sampling involved with data collection. These challenges resonate strongly with the difficulties faced in our work, and have guided our considerations when building our model. Besides fostering stronger interactions with data collectors and end users of model output, better communication between statisticians, domain experts, and citizen scientists should continuously be stimulated, not just to ensure its continuity by increasing public awareness of the scientific process (Bonney et al., 2016) but also to improve data and model quality. As observed by Reese, ‘Koh and Opitz’s work would usefully be made more accessible to ornithologists, with the mathematics simplified and more biological interpretation’.

We are inspired by the positive comments of the Seconder. We agree with the sentiment that statisticians are increasingly being trusted as ‘data whisperers’, and this comes at an exciting time with new combinations of models, algorithms, and computing power. Though we acknowledge the need to propose methods appropriate for other ‘wilder’ datasets in this domain, we hope that our work on a ‘well-behaved’ CS dataset will trigger new extensions and collaborations, and can be a reference point for further discussion on this important topic.

Several key aspects are fundamental for statistical methods to extract useful insights from the wealth of data sources and powerful AI tools that should guide new methodological developments. We have to explain the phenomena and achieve good predictive behaviour by separating the ecological signal from the noise and biases, thus allowing for generalization and extrapolation of results. Statistical models representing bias and noise explicitly with latent variables (e.g. through spatio-temporal random effects) provide interesting paths to this goal. Moreover, we should provide decision support through statistical inference (e.g. selecting variables and finding significant effects), and provide probabilistic statements to assess remaining uncertainties. We believe that Bayesian hierarchical models (BHMs) provide an appropriate framework to jointly address these aspects.

Regarding the discussion contributions on our work, we have thematically regrouped their comments along with our responses in the five sections below. Extreme-value aspects are addressed in Section 1; Section 2 discusses validating complex BHMs, especially when weakly structured CS data are involved. Considerations on model identifiability, efficient estimation algorithms, and

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

<footnote n="2">²Seminar for Statistics, ETH Zurich, Rämistrasse 101, 8092 Zurich, Switzerland</footnote>

<footnote n="3">³Biostatistics and Spatial Processes (UR546), INRAE, 228 route de l'Aérodrome, 84914 Avignon, France</footnote>
the potential of machine learning (ML) techniques follow in Section 3. Possible new applications of CS statistics beyond ecology are summarized in Section 4. Section 5 concludes with a discussion on the opportunities of CS data and model fusion.

# 1 Extreme-value distributions and their priors

A key focus of our work is the analysis of extreme phenological events (i.e. life-history events of individuals of a species) using probability distributions motivated by extreme-value theory (EVT). We iterate our hope that this work will spur broader applications of these techniques in ecology and other fields concerned with biotic processes where they are still scarce, and we are appreciative of the insightful comments by the discussants on this theme.

Reese scrutinizes the definition of first arrival dates that we model with the generalized extreme value (GEV). It is indeed difficult for an observer to know whether a spotted bird has arrived at its *destination* (i.e. the breeding site) or has just stopped before continuing its journey. Nevertheless, we were able to implement some basic data filtering operations using the ‘flyover’ flag in eBird, which observers can toggle when they are unsure whether the bird has already arrived at its destination. We also believe that there might be a ‘weekend bias’, though this bias could be evened out across pixel-years in our spatio-temporal model. Reese’s idea to model dates consistently over leap years, for example by modelling the first date as 1st March, is also interesting; this is easy to implement and would only lead to a location shift for the fitted GEV. We agree that there should be more discussion between ecologists and statisticians to better differentiate and estimate both individual and populational behaviour of migratory bird arrivals.

To address Reese’s comments on Figure 1 of our main paper, our reasoning to invert the y-axis there is based on the EVT convention of modelling maxima instead of minima; extreme values are typically associated with the right tail. In our case, the higher our transformed variable (which we later use for modelling), the earlier the date. We also chose to have $log(N_i)$, with $N_i$ the number of checklists in pixel-year $i$, in the x-axes of that figure since we later model the checklist count with a log-link function. We appreciate Reese’s suggestion to estimate contours of northward migration of birds and will note this for follow-up work. Relatedly, Figure 1 shows a strong dependence of the first arrivals on latitude, a feature that has also been detected by our model; see Figure 9 of the main paper.

The Proposer suggests extending our approach to multivariate modelling by using copulas and especially multivariate extreme-value distributions, to share information across species and to enable joint inference of the first arrivals of multiple species. We agree that multivariate modelling would be a natural and worthwhile extension and could provide further insights into bird community structure. Regarding our study, we argue that it would pertain to a different modelling philosophy. It may require sacrificing important benefits of mono-specific modelling and of the conceptually and computationally useful conditional independence assumption in BHMs. Using multivariate probability distributions for the response variable would entail building our model components differently, and care would have to be taken not to lose computational tractability or omit important model components. A major drawback of the joint modelling approach is the inability to simply parallelize computations across tens to thousands of species as we do in our work. Choosing the right multivariate distribution also raises questions about choosing the copula that can capture the appropriate dependence class for the first arrivals of multiple species, and this is an avenue yet to be explored. An alternative approach to share information between species, more in line with our current model, would be to keep the conditional independence assumption inherent in our BHM but to additionally incorporate cross-species dependencies into the latent processes, for instance using tensor products of Gaussian processes as suggested by the Proposer. Multivariate Gaussian processes are already routinely used in joint species distribution models, which have their foundations on Chib and Greenberg (1998), with many recent developments (e.g. Chiquet et al., 2021; Ovaskainen & Abrego, 2020; Tikhonov et al., 2020). We concur that using computationally efficient representations (e.g. Gauss Markov random fields) could indeed make certain multivariate extensions feasible, even for complex models combining multiple response variables for each species.

Rizzelli advocates for more informative prior distributions for the shape parameter $\xi$ of the GEV, which is expected to be negative in our case study due to the finite upper bound of the variable of
<figure><img src="image_26.png" /><figcaption>**Figure 1.** Scatterplots of the latitude (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 21 binned estimates, the black line a smooth fitted curve, and the shaded region the 95% pointwise confidence intervals of the curve.</figcaption></figure>

transformed first arrival days. This is a valid point, and a strategy in this direction is to use penalized complexity priors (Simpson et al., 2017), though the choice of the reference distribution when building this prior remains less clear in our application. Choosing the reference distribution as the GEV with $\xi = 0$, as advocated in Opitz et al. (2018), penalizes relatively heavy or short tails in the GEV, and we could truncate this prior such that only negatives values of $\xi$ are possible. However, this would problematically push the posterior towards earlier arrivals in our application, especially when observations are scarce. A different approach is to model $\xi = \log(-\xi)$ instead, forcing $\xi$ to be negative.

Modelling $\xi$ instead of $\xi$ falls under the general strategy of reparametrization, which could also be motivated by parameter identifiability and computational reasons; we discuss this more in Section 3. One could also directly model the finite upper bound of the GEV as a parameter $B = \mu - \sigma/\xi$ with the usual parameter notation, and in our case incorporate the $g$ function on this bound. Alternatively, one could avoid estimating this upper bound $B$ but rather fix it to a known value, e.g. the value corresponding to 1st January or 1st March, or based on other physical constraints; see Noyelle et al. (2024) who develop this approach in the context of temperature extremes. In a reparametrized model, one of the three GEV parameters could be equal to the upper bound.

Rizzelli also discusses the general inflexibility of extreme-value models, especially for modelling contexts with biotic processes, and relatedly discusses EVT in the presence of covariates. When working with Peaks-Over-Threshold models and the resulting asymptotic Generalized Pareto distribution, common practice in the covariate setting is to assume conditional independence of the response given the covariates, such that the usual asymptotics from the classical 'i.i.d.' setting hold after the conditioning (e.g. see Koh et al., 2025). By contrast, deriving theoretical asymptotics for maxima over non-identically distributed variables is more involved; for example, see the results derived for extreme-value limits of upper-bounded distributions with seasonal patterns in Stein (2017), and for semiparametric tail estimation (i.e. without a limit probability distribution such as the GEV) in Einmahl and He (2023). We agree that rigorous attempts to provide mathematical guarantees in the presence of covariates, as done for example in the temporal heterogeneity setting, should be further explored. Though asymptotically justified distributions confer more robustness in data-scarce contexts, asymptotic properties remain difficult to verify formally in many data applications.

Daouia and Stupfler propose an interesting approach from econometrics based on frequentist frontier analysis (Daouia et al., 2016), with the frontier defined as the lower bound of arrivals, i.e. the first arrival date. They suggest using all data covering the study region and incorporating as covariates the spatial coordinates, exogenous climate, land cover, and other ecological factors into the estimation procedure in a nonparametric regression framework with one-sided errors. This could be applied in continuous space without needing a spatial pixel discretization. Instead, one would use a kernel with a data-driven choice of its bandwidth to pool nearby observations to estimate local parameters. However, as far as we
understand, this approach would not attempt to explicitly debias the boundary estimate, which is the direction taken in our work and others in ecological statistics that explicitly capture the observation process (Gelfand & Shirota, 2019). We acknowledge that the sampling effort from the preference dimension would be captured indirectly by the width of confidence intervals of the estimated frontier, though it is less clear to us how one should systematically incorporate the activity component (i.e. the ‘temporal effort’ of an observer during an observation event). Pooling data locally around each prediction location could certainly provide spatially smoother predictions of the first arrivals. However, we believe that it would also be quite sensitive to the spatio-temporally heterogeneous distribution of sampled locations in our application, and it might be ill-adapted to the sharp spatial changes in land cover configurations and the corresponding changes in ecological processes. Nevertheless, we would be delighted to see future work combining our model with the one proposed.

## 2 Validation of inferences and debiased predictions

A central question when modelling CS data is how models, and the bias corrections for varying and unknown observational effort they provide, can be validated. Attenuating the possible biases in model-based predictions requires either explicitly identifying the effort or predicting only quantities that are not affected by varying effort under realistic assumptions. This is a challenging and fundamental topic of research in CS statistics, and not surprisingly appears as a common theme among many of the discussants.

The Proposer poses questions about validating latent variable models in the presence of complex observational processes. We reiterate that using structured data (i.e. collected by following a known predefined sampling design) to define validation datasets, or constructing integrated models combining structured and unstructured data sources, remain the best options to reliably validate inferences on observation biases. The Proposer also questions how one should estimate the effective number of parameters in our setting to assess model complexity and calculate information criteria that allow comparing different models. This appears challenging, and a reasonable approach could focus on the parameters describing the ecological process by integrating out the observational process before estimating the model’s remaining degrees of freedom.

Our model’s debiasing mechanism is based on appropriately identifying the observation process jointly with the other model components. An approach to better validate this mechanism is to randomly remove observation days to artificially reduce the observational effort, and then verify the observed first arrivals, calculated including the hold-out data, with the model-debiased ones estimated under lower observational effort. Similarly, one could assume that regions with high observational effort should already be at the saturation point of our g function, and verify whether the observed first arrivals in these regions are not significantly different from the model-debiased ones under high observational effort.

Rizzelli highlights the advantages of the Bayesian approach, such as the possibility of using predictive distributions for forecast evaluation, including using summary statistics like the predictive mode and median. We agree that validation schemes holding out recent data to mimic prediction with relatively short lead times, along with coverage checks on predictive distributions, are relevant when models are used to forecast into the near future. These approaches have been considered in previous studies with BHMs (e.g. in wildfire modelling contexts by Koh et al., 2023).

Tian emphasizes building hierarchical frameworks with identifiable components. Though we verify this for our model with a simulation study, we resonate with his call to explore the conditions that ensure identifiability of parameters in complex BHMs more generally.

We also agree with Belzile and Yadav that a possible model assessment approach is to use *k*-fold cross-validation designs appropriately stratified for spatial and temporal dependence. This is generally good practice, but it would be computationally demanding in our setting.

## 3 Conducting computations and inference

We here respond to the numerous comments on how model building, Bayesian computation and inference can be conducted for complex models such as ours.
Achieving computational efficiency while staying flexible with respect to model structure and
inferrable quantities remains challenging despite continual progress in computing power and al-
gorithmic efficiency. A recommended practice for BHMs is to structure their latent variables
and hyperparameters to maximally facilitate their identification, and to keep auto-correlations
low among the posterior samples from hyperparameters and latent variables.

## 3.1 Identifiability of latent spatio-temporal fields

Belzile and Yadav remark that some components of our latent model may be difficult to identify from data. Specifically, they notice that the field $x_{\text{bound}}(s, t)$, which is in the GEV location parameter $\mu(s, t) = \exp\{x_{\text{bound}}(s, t)\}/\{1 - \exp(-x_{\text{effort}}(s, t))\}$, could be hard to identify at space-time locations $(s, t)$ where the field $x_{\text{effort}}(s, t)$ representing observational effort is low. We consider the way our model works in this situation rather as a feature than a problem: if the effort is low at $(s, t)$, then the GEV location parameter will tend to be smaller, which implies that observed first arrivals are expected to occur relatively late in the year $t$, or not at all if observations are missing; this feature is coherent. Moreover, if the effort is low at $(s, t)$, then the posterior of the field $x_{\text{bound}}(s, t)$ will be strongly influenced by the behaviour at surrounding locations and its intercept term, i.e. the model will borrow information from nearby locations where more data are available and from global behaviour. We do not see an identifiability issue here.

Regarding the identifiability of our latent fields, one could explore alternative parameterizations of the GEV that were shown to yield weaker correlations between the posterior samples of these three parameters, as outlined in Section 1. Moins et al. (2023) propose a reparameterization that reduces the dependence between the posterior samples of the three GEV parameters and improves identifiability. Using this alternative parametrization in our BHM could facilitate identification of the latent variables and parameters we incorporated into them. Another idea to aid identifiability is to impose additional modelling assumptions, e.g. as in Section 2, one could assume that data from high observational effort regions should already be the g function's saturation point, and impose this as a constraint in the model.

Tian questions whether additional constraints or prior dependencies should be imposed to further link the different Gaussian random effects in our model and ensure identifiability. We underline that these effects are of quite different nature, characterizing, for example, spatial niche variability, interannual variability, spatial preferential sampling, and spatial variation of true first arrivals. We do not see strong a priori links between these different effects and their sources of variability. Therefore, we believe that independent Gaussian process priors are more appropriate for our model.

## 3.2 MCMC-based inference and alternatives

Several discussants, including the Proposer, have commented on the relative slow convergence and mixing of some of the posterior chains in our model, and have asked (e.g. Tian) how one could efficiently monitor convergence of MCMC algorithms for complex hierarchical space-time model. Suggestions for algorithmic improvement from Belzile and Yadav include closer monitoring of effective sample size (ESS), using adaptive MCMC schemes instead of the global Metropolis-adjusted Langevin algorithm (MALA) we implemented, performing joint hyperparameter updates when posteriors of hyperparameters are correlated, and running multiple chains. We agree that these techniques could help improve our MCMC algorithm and its diagnostics. Another interesting option we see would be to use Hamiltonian MCMC instead of MALA to better explore high-density regions of the posterior distribution through the Hamiltonian dynamics (Neal, 2011). Moreover, restructuring the model form to improve identifiability, as discussed in the Section 3.1, could also aid convergence.

As the Proposer recalls, Gauss-Markov random field (GRMF) structures facilitate fast calculations with latent fields in spatial Bayesian frameworks. The INLA-SPDE approach is a well-known example (Lindgren et al., 2011; Rue et al., 2009); it uses GMRF approximations of spatial Matérn covariance functions based on the stochastic partial differential equation (SPDE) approach, combined with fast and accurate deterministic posterior approximations with the INLA method. The Proposer suggested that simulation-based MCMC inference, as implemented in our work, could be combined with the INLA-SPDE approach depending on the structure of the components of the
BHM. This could leverage INLA for model components where linear combinations of latent variables in GMRF form arise, whereas slower MCMC inference can be used only for parts with strongly nonlinear structures. We agree that such extensions bear promise. Combined inference of this type has been developed by Gómez-Rubio and Rue (2018), where INLA is run within an MCMC algorithm by conditioning on those other variables that would have impeded full INLA-based inference. This approach could indeed lead to performance gains for our model, although its implementation appears complex, and further research into such combined methods would be welcome in the setting of complex spatio-temporal Bayesian inference. Our main algorithmic difficulty preventing us from using faster inference tools comes from the nonstandard, nonlinear and parametric transformation of GMRFs that achieves a saturating effect when the observational effort increases to infinity, as also recognized by the Proposer. In this respect, the inlabru package (Bachl et al., 2019) with its extensions of the classical INLA method could offer help through implementation of methods that handle nonstandard nonlinear transformations of GMRFs, but to our knowledge, it does not yet cover the parametric transformations we would need. In response to Tian, who suspected possible issues with the relatively large number of fixed-effect coefficients (e.g. regression constants) in our model, we reiterate that they are endowed with Gaussian priors and are part of the linear predictors in the model (i.e. of multivariate Gaussian vectors), such that their estimation does not pose specific challenges.

## 3.3 Vecchia approximation

The Proposer questions whether the Vecchia approximation provides a sufficiently accurate representation of the underlying covariance function. Indeed, differences could arise between the approximation and the exact covariance function. However, we emphasize that the principal role of the Vecchia-approximated Gaussian field in our model is to serve as a rich prior with spatial dependence—what matters is the flexibility it provides for posterior inferences. We do not necessarily want to exactly reproduce the underlying covariance function, but rather obtain a prior that has similar and interpretable properties in terms of variance and correlation range. We do not use the Vecchia approximation to simulate data directly from the prior model, where stronger deviations could arise between the exact and approximated simulations.

## 3.4 Likelihood-free amortized inference

Huser and Zammit-Mangion outline the advantages of likelihood-free inference using neural networks (NNs) and recent developments in this area. Neural Bayes estimators are NNs trained on large amounts of simulated data to estimate (i.e. *predict* through the NN) the parameters of a stochastic process. They have emerged as powerful alternatives to likelihood-based estimation (Zammit-Mangion et al., 2025). The approach is still limited to models with a moderate number of parameters, since the parameter configurations used for simulating training data should cover the possible range of parameters to be estimated from real data. While estimating the large number of latent variables in a BHM through this approach appears unattainable, the discussants question whether hyperparameters, usually of moderate number in a BHM, could be estimated with the NN approach. Other parameters, especially the latent variables representing fixed effects and spatio-temporal effects, would be estimated separately. Intercepts could be considered as hyperparameters (i.e. estimated by the NN) or as latent variables (i.e. estimated separately). Given the proven benefits of neural Bayes estimators for estimating parameters of various spatial stochastic process classes, we recognize that they could also facilitate inference for spatial latent-process models. They would be especially beneficial in situations where parameters of the same model have to be estimated for a large number of different datasets (e.g. for different regions of the globe).

We comment on possible challenges for developing neural Bayes estimators in the setting of latent models for discrete responses with covariates. Three estimation steps could be necessary: (i) covariate coefficients (including intercepts) using a pseudo-likelihood regression technique; (ii) hyperparameters using a neural Bayes estimator; (iii) latent variables using a chosen technique (e.g. MCMC or INLA with fixed covariate coefficients and hyperparameters, or new NN-based approaches). Clearly, one would lose the benefits of performing joint inference on all model components, but could in turn gain in robustness and speed.
Neural Bayes estimators may be more difficult to implement when there are covariates. In latent-variable models for discrete response data (e.g. binary or count data), it is not possible to define classical standardized residuals, and marginal distributions and the dependence structure cannot be fully separated as in copula techniques for continuous margins. Therefore, it is more intricate to separate fixed covariate effects acting on marginal distributions from the random spatio-temporal effects governed by hyperparameters. For example, with binary and count data, the residuals take values outside the discrete support of the response. Recall that this is unlike continuous response distributions, such as the normal $N(\mu_i, \sigma_i^2)$, where we could estimate marginal parameters $\mu_i$ and $\sigma_i^2$ in a first step, then generate standardized residuals by normalizing observations $y_i$ to $(y_i - \mu_i)/\sigma_i$, and finally infer dependence hyperparameters from standardized residuals, using methods such as the NN estimator for the last step. This is not possible with discrete data, which also explains the popularity of fully Bayesian inference (i.e. with joint estimation of the posterior distributions of hyperparameters and latent variables) when modelling this data, which are very common in fields such as ecology.

To make the NN approach for estimating hyperparameters feasible in the presence of covariates, one could estimate the fixed-effect part with covariates in a preliminary step and then feed it into a NN as one of its inputs used for predicting the hyperparameters. However, there are essentially no constraints on the values of the fixed effect field, which presents a very high-dimensional input to the NN. This could make it very challenging to simulate a sufficiently large number of NN training instances that ensure good generalization behaviour.

We would be very interested to see how research on this topic could lead to NN estimators for latent-variable models that cope well with the aforementioned challenges.

## 3.5 Simplified estimation of shared latent fields

Instead of estimating certain hyperparameters differently, or separate from the other parameters, as suggested by Huser and Zammit-Mangion, one could also simplify the model by separately estimating some of the latent variable fields, or some of the regression equations.

We emphasize that the motivation for our BHM consists of jointly estimating all model components to allow for optimal information flow between components and to precisely assess posterior uncertainties. However, to alleviate computational cost and improve identifiability, one could instead separately estimate some of the shared fields by using only a subset of the response variables for which they appear in the predictor. The species niche field $X^{\text{niche}}(s)$ could be estimated using only the BBS data (with its $\lambda^{\text{BBS}}(s)$ parameter) and the binary presence-absence indicator for each checklist (with its $p^{\text{spc}}(s, t)$ parameter). We would not use the GEV response of first arrival dates to estimate the species niche field, for which it appears in the GEV location parameter $\mu(s, t)$. Now, assume we have estimated the species niche field separately from the GEV model in a first step. Then, in a second step, posterior means of the separately estimated fields could be included as a simple covariate (i.e. as a fixed effect) in the GEV model. We could proceed similarly for other spatial fields (e.g. the observational effort $\lambda^{\text{ckl}}(s, t)$) that we share towards the GEV response by estimating them separately in a preliminary step. Implementing this simplified model structure should drastically reduce computational cost while maintaining identifiability and reliable predictions.

## 3.6 Using ML algorithms

The main proposed utility of our model is the separation of the ecological signal from heterogeneous observational effort. This is made possible by identifying how the latter modifies outcomes when we consider several data sources with different observation design. Disentangling the roles of these different latent processes cannot be easily achieved with standard ML algorithms, even more so if they have black-box like behaviour. Nevertheless, we believe that ML models can still be useful if the predictor and response variables are adapted appropriately, but we are not aware of any existing general framework proposed to this end in the literature.

We provide a specific example where standard ML algorithms could be leveraged. Assume that only the presences of the individuals of various species were reported, with non exhaustive and varying observational effort, i.e. we have opportunistic presence-only data. Using a point-process framework, we can construct a response variable that can be estimated with essentially any
classification algorithm as follows: The points of the point process are defined as the reported presences of a focal species for which we want to make predictions. A set of background points representing pseudo-absences of this species are constructed from the reported presences of other species for which similar observational effort can be reasonably assumed. Therefore, the background points represent a proxy for the general observational effort. Any classification algorithm trained to classify the points into ‘species’ vs. ‘background’ can then be used to predict the probability of species presence conditional on environmental predictors and location. The predicted probabilities can be interpreted as the relative abundances of the focal species with respect to the group of background species, which can provide valuable ecological insights. This so-called *target-group background* approach can be deployed even for purely opportunistic data without the need for strong assumptions on sampling design (Botella et al., 2020, 2021, Phillips et al., 2009).

In the above example, and certainly in many other CS models, ML approaches could be highly efficient, such as the use of NNs or stochastic gradient boosting techniques advocated by Kumar. Crucially, the model and the input data must be appropriately structured, which requires solid assumptions on how data were collected and how ecological and observational processes interact.

Kumar also encourages using ML models in extreme-value modelling. We concur with this sentiment, though this approach is less common when modelling count data; for an example, Koh (2023) uses gradient boosting to predict the parameters of discrete generalized Pareto distributions.

# 4 Pushing CS beyond ecology

CS data are commonly exploited in statistical ecology, and especially for species distribution modelling. Consequently, this research community has become a major driver for advancing statistical methods for CS purpose. As the Proposer and Srakar suggest, CS-based approaches are promising in various other areas of applied science, including epidemiology and public health, environmental surveillance, and possibly even in fields such as medicine, robotics, and law. We absolutely share this view. Due to its potential for collecting large amounts of diverse data, CS can provide valuable contributions to integrated modelling, for instance of entire ecosystems in the *One Health* perspective (Destoumieux-Garzón et al., 2018). CS also appears highly relevant for filling data gaps through the collection of data using mobile devices, especially in countries with low coverage of standardized datasets (e.g. in the Global South), and can help monitor and assess the sustainable development goals formulated by the United Nations.

It is challenging to establish statistical guarantees for reliable CS-based inferences similar to the asymptotic results in classical statistics, and no general theory has emerged so far. This is certainly due to the wide variety of application-specific data collection processes and the lack of well-known sampling designs, especially with opportunistic data. As Srakar highlights, more work should be invested into theoretical guarantees for CS-based approaches. We second this opinion and would like to incite statisticians to propose appropriate sampling design frameworks that preserve the simplicity and low cost of CS-based data collection while ensuring reliable inferences and predictions. Checklist-based approaches are an interesting example in species monitoring: observers must comply to exhaustively report relevant observations. Alternatively, the software interface used to report the species occurrences could explicitly ask the observers to confirm the nonobservation of certain other species not reported. Datasets collected with this protocol allows modellers to discern the absence of an observer (i.e. missing data) from the absence, or non detection, of any species individuals. One can then define pseudo-absences of individuals, which is crucial for incorporating observational effort and reducing bias in species distribution mapping, as highlighted in our case study.

However, ecological datasets covering past periods often contain only opportunistic presence-only data: the absence of a species is not explicitly reported, and the reporting is not exhaustive, i.e. it is possible that the observer detects an individual of a species but does not report it, so we are not in the convenient checklist setting mentioned before. Then, a strategy is needed to generate useful pseudo-absences, for example by using the target-group background approach outlined in Section 3.6. Given a focal species that one wants to model, the space-time locations where other species (the *target group*) were reported could be declared as pseudo-absences for the focal species. This requires carefully choosing the target group of species and the subsample scheme for its occurrences (Botella et al., 2020, 2021, Phillips et al., 2009; Van Strien et al., 2013). Knowledge of ecological
experts is indispensable for this step, and it highlights the strong need for interdisciplinary work between data scientists and domain experts when developing statistical methods for CS data.

# 5 Opportunities of data and model fusion

By developing sampling approaches and models that integrate unstructured CS data with structured ones, we could benefit from both the often dense spatio-temporal coverage of CS data and the precisely known observation design and effort in structured data. The fundamental question was asked by the Proposer: When and how should we do this? Regarding spatial species distribution modelling, recent studies (e.g. Simmonds et al., 2020) highlight considerable added value from data integration, provided that models are appropriately structured using spatial random effects. Ideally, one would need covariates that strongly correlate with observational effort, but not with the species presence or abundance. However, such covariate data are typically unavailable in practice, except in very specific sampling situations (e.g. distance sampling). Instead, one usually needs both dataset-specific and shared spatial random effects to disentangle the observational effort from ecological effects. The flip side of implementing multiple spatial random effects and regression equations in these approaches is that estimation becomes computationally very intensive, and setting up models and algorithms may require considerable expertise and experience.

The Seconder further encourages statisticians to compare and hybridize statistical and ML approaches to tackle such complex statistical questions. We are curious how ML techniques could be adapted to mimic the estimation of spatial effects; for example, one could try to achieve this by feeding the predictor data, spatial coordinates, and response data into the prediction algorithm in a specific way. The example we outline in Section 3.6, where the modelling question is transformed into a classification problem for ML algorithms, is an example where combining expertise on statistical modelling with ML leverages the simple use of ML.

We also suggest to investigate how existing large CS datasets with limited knowledge about observational effort could be efficiently complemented by the additional collection of relatively small datasets in controlled surveys. The sampling design for newly collected standardized data would have to ensure reliable inferences and identification of observational effort through integrated models combining CS and survey data, while the cost of sampling in terms of required human and material resources must remain manageable. Modern technology for automated data collection (e.g. camera traps and other remote sensing devices) could be used for collecting data. Principles and techniques of optimal design can be applied (Müller, 2005; Pukelsheim, 2006) to develop variants of stratified or adaptive sampling, for example in a model-based perspective. BHMs, such as modifications of our model, could be useful for Bayesian optimal design experiments (Ryan et al., 2016) adapted to geostatistical applications (Diggle & Lophaven, 2006).

Once we have general standardized statistical methods to obtain trustworthy CS-based insights, routine publication of scientific results from CS would be facilitated. We are aware of the long road ahead, but are optimistic about the future.

# Acknowledgments

The authors thank the Royal Statistical Society for facilitating the discussion meeting.

*Conflict of interests:* The authors declare no conflict of interests.

# References

Bachl F. E., Lindgren F., Borchers D. L., & Illian J. B. (2019). inlabru: An R package for Bayesian spatial modelling from ecological survey data. *Methods in Ecology and Evolution, 10*(6), 760–766. https://doi.org/10.1111/mee3.2019.10.issue-6

Bonney R., Phillips T. B., Ballard H. L., & Enck J. W. (2016). Can citizen science enhance public understanding of science? *Public Understanding of Science*, 25(1), 2–16. https://doi.org/10.1177/0963662515607406

Botella C., Joly A., Bonnet P., Munoz F., & Monestiez P. (2021). Jointly estimating spatial sampling effort and habitat suitability for multiple species from opportunistic presence-only data. *Methods in Ecology and Evolution, 12*(5), 933–945. https://doi.org/10.1111/mee3.v12.5
Botella C., Joly A., Monestiez P., Bonnet P., & Munoz F. (2020). Bias in presence-only niche models related to sampling effort and species niches: Lessons for background point selection. *PLOS One*, 15(5), e0232078. https://doi.org/10.1371/journal.pone.0232078

Chib S., & Greenberg E. (1998). Analysis of multivariate probit models. *Biometrika*, 85(2), 347–361. https://doi.org/10.1093/biomet/85.2.347

Chiquet J., Mariadassou M., & Robin S. (2021). The Poisson-lognormal model as a versatile framework for the joint analysis of species abundances. *Frontiers in Ecology and Evolution*, 9, 588292. https://doi.org/10.3389/fevo.2021.588292

Daouia A., Noh H., & Park B. U. (2016). Data envelope fitting with constrained polynomial splines. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 78(1), 3–30. https://doi.org/10.1111/rssb.12098

Destoumieux-Garzón D., Mavingui P., Boetsch G., Boissier J., Darriet F., Duboz P., Fritsch C., Giraudoux P., Le Roux F., & Morand S. (2018). The one health concept: 10 years old and a long road ahead. *Frontiers in Veterinary Science*, 5, 14. https://doi.org/10.3389/fvets.2018.00014

Diggle P., & Lophaven S. (2006). Bayesian geostatistical design. *Scandinavian Journal of Statistics*, 33(1), 53–64. https://doi.org/10.1111/sjos.2006.33 issue-1

Einmahl J. H., & He Y. (2023). Extreme value inference for heterogeneous power law data. *Annals of Statistics*, 51(3), 1331–1356. https://doi.org/10.1214/23-AOS2294

Gelfand A. E., & Shirota S. (2019). Preferential sampling for presence/absence data and for fusion of presence/absence data with presence-only data. *Ecological Monographs*, 89(3), e01372. https://doi.org/10.1002/ecm.2019.89 issue-3

Gómez-Rubio V., & Rue H. (2018). Markov Chain Monte Carlo with the integrated nested Laplace approximation. *Statistics and Computing*, 28(5), 1033–1051. https://doi.org/10.1007/s11222-017-9778-y

Koh J. (2023). Gradient boosting with extreme-value theory for wildfire prediction. *Extremes*, 26(2), 273–299. https://doi.org/10.1007/s10687-022-00454-6

Koh J., Pimont F., Dupuy J.-L., & Opitz T. (2023). Spatiotemporal wildfire modelling through point processes with moderate and extreme marks. *The Annals of Applied Statistics*, 17(1), 560–582. https://doi.org/10.1214/22-AOAS1642

Koh J., Steinfeld D., & Martius O. (2025). Using spatial extreme-value theory with machine learning to model and understand spatially compounding weather extremes. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 481, 20240763. https://doi.org/10.1098/rspa.2024.0763

Lindgren F., Rue H., & Lindström J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 73(4), 423–498. https://doi.org/10.1111/j.1467-9868.2011.00777.x

Moins T., Arbel J., Girard S., & Dutfoy A. (2023). Reparameterization of extreme value framework for improved Bayesian workflow. *Computational Statistics & Data Analysis*, 187, 107807. https://doi.org/10.1016/j.csda.2023.107807

Müller P. (2005). Simulation based optimal design. *Handbook of Statistics*, 25, 509–518. https://doi.org/10.1016/S0169-7161(05)25017-4

Neal R. M. (2011). MCMC using Hamiltonian dynamics. In *Handbook of Markov Chain Monte Carlo* (Vol:2(11), p. 2).

Noyelle R., Robin Y., Naveau P., Yiou P., & Faranda D. (2024). Integration of physical bound constraints to al-leviate shortcomings of statistical models for extreme temperatures. HAL-04479249.

Opitz T., Huser R., Bakka H., & Rue H. (2018). INLA goes extreme: Bayesian tail regression for the estimation of high spatio-temporal quantiles. *Extremes*, 21(3), 441–462. https://doi.org/10.1007/s10687-018-0324-x

Ovaskainen O., & Abrego N. (2020). *Joint species distribution modelling: With applications in R*. Cambridge University Press.

Phillips S. J., Dudík M., Elith J., Graham C. H., Lehmann A., Leathwick J., & Ferrier S. (2009). Sample selection bias and presence-only distribution models: Implications for background and pseudo-absence data. *Ecological Applications*, 19(1), 181–197. https://doi.org/10.1890/07-2153.1

Pukelsheim F. (2006). *Optimal design of experiments*. SIAM.

Rue H., Martino S., & Chopin N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. *Journal of the Royal Statistical Society: Series B, Statistical Methodology*, 71(2), 319–392. https://doi.org/10.1111/j.1467-9868.2008.00700.x

Ryan E. G., Drovandi C. C., McGree J. M., & Pettitt A. N. (2016). A review of modern computational algorithms for Bayesian optimal design. *International Statistical Review*, 84(1), 128–154. https://doi.org/10.1111/insr. v84.1

Simmonds E. G., Jarvis S. G., Henrys P. A., Isaac N. J., & O'Hara R. B. (2020). Is more data always better? A simulation study of benefits and limitations of integrated distribution models. *Ecography*, 43(10), 1413–1422. https://doi.org/10.1111/ecog.2020.v43.i10
Simpson D., Rue H., Riebler A., Martins T. G., & Sørbye S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. *Statistical Science*, 32(1), 1–28. https://doi.org/10.1214/16-STS576

Stein M. (2017). Should annual maximum temperatures follow a generalized extreme value distribution? *Biometrika*, 104(1), 1–16. https://doi.org/10.1093/biomet/asw070

Tikhonov G., Opedal Ø. H., Abrego N., Lehikoinen A., de Jonge M. M., Oksanen J., & Ovaskainen O. (2020). Joint species distribution modelling with the R-package Hmsc. *Methods in Ecology and Evolution*, 11(3), 442–447. https://doi.org/10.1111/mee3.v11.3

Van Strien A. J., Van Swaay C. A., & Termaat T. (2013). Opportunistic citizen science data of animal species produce reliable estimates of distribution trends if analysed with occupancy models. *Journal of Applied Ecology*, 50(6), 1450–1458. https://doi.org/10.1111/jpe.2013.50 issue-6

Zammit-Mangion A., Sainsbury-Dale M., & Huser R. (2025). Neural methods for amortized inference. *Annual Review of Statistics and Its Application*, 12, 311–335. https://doi.org/10.1146/annurev-statistics-112723-034123

https://doi.org/10.1093/jrsssa/qnaf058
Advance access publication 14 May 2025