Bayesian Spatiotemporal Occupancy Model

• The model is a hierarchical framework that estimates species occupancy and dynamic changes, explicitly accounting for detection imperfections and spatial-temporal autocorrelation.
• It uses a combination of Gaussian processes and BYM-2 spatial random effects, with Hamiltonian Monte Carlo and Polya-Gamma augmentation enhancing computational efficiency.
• Applications include analyzing citizen science data, mapping retrospective distributions, and providing robust ecological inferences for conservation and management.

A Bayesian spatiotemporal site-occupancy model is a hierarchical probabilistic framework designed to infer species occupancy and dynamics across space and time, explicitly accommodating detection imperfections and spatiotemporal autocorrelation. This methodology is essential for analyzing opportunistically collected biodiversity data, such as those sourced from citizen science platforms, where the data-generating process is complex, observer effort is highly variable, and spatial/temporal coverage is irregular. Bayesian inference allows integration of structured prior knowledge and propagation of uncertainty throughout model components (Fajgenblat et al., 11 Nov 2025, Diana et al., 2021).

1. Hierarchical Model Formulation

The Bayesian spatiotemporal site-occupancy model introduces a multi-level structure indexed by spatial sites $s=1,\ldots,S$ and temporal units (often years) $t=1,\ldots,T$. Each site-year combination $(s,t)$ is associated with a binary latent state $z_{s,t} \sim \mathrm{Bernoulli}(\psi_{s,t})$, where $\psi_{s,t}$ denotes the probability of occupancy during year $t$ at site $s$ (Fajgenblat et al., 11 Nov 2025).

The site-year occupancy probability is logit-modeled as: $\logit(\psi_{s,t}) = \beta^{\psi}_0 + X^{\psi}_s \beta^{\psi} + \delta_t + \upsilon_s + \varsigma_s t^*,$ where $\beta^{\psi}_0$ is a global intercept, $X^{\psi}_s$ are site-level covariates, $\delta_t$ is a temporal random effect, $\upsilon_s$ is a spatial random effect (BYM-2 mixture of unstructured and B-spline-GP components), and $\varsigma_s t^*$ models local deviations in trend as a site-specific slope on scaled time $t^* \in [-0.5,0.5]$ (Fajgenblat et al., 11 Nov 2025).

Detection is modeled with a conditional Bernoulli likelihood for each visit $v$ during year $t(v)$ at site $s(v)$: $$y_v|z_{s(v),t(v)} \sim \mathrm{Bernoulli}(p_v z_{s(v),t(v)}),$$ and detection probability $p_v$ is logit-linear: $\logit(p_v) = \beta^p_0 + X^p_v\beta^p + b^{obs}_{o(v)} + f_{phen}(w(v)),$ with $b^{obs}_{o(v)}$ a random observer effect and $f_{phen}(w(v))$ a periodic Gaussian process to capture phenology (Fajgenblat et al., 11 Nov 2025).

2. Random Effects and Gaussian Process Priors

Temporal structure is decomposed as: $$\delta_t = \beta^{\delta} t^* + \delta_t^{GP} + \delta_t^{i.i.d.},$$ where $\delta_t^{GP} \sim \mathcal{GP}(0,k_\delta)$ with $$k_\delta(t,t') = \sigma^2_\delta \exp\left(-\frac{(t-t')^2}{2\ell^2_\delta}\right)$$, and $$\delta_t^{i.i.d.} \sim \mathcal{N}(0,\sigma^2_{\delta,i})$$ (Fajgenblat et al., 11 Nov 2025).

Spatial random effects use a BYM-2 style mixture: $$\upsilon_s = (1-p)\upsilon^{unstr}_s + p\upsilon^{str}_s,$$ with $$\upsilon^{unstr}_s \sim \mathcal{N}(0,\sigma^2_{space})$$ and $\upsilon^{str}_s$ a B-spline-projected GP over longitude/latitude knots.

The detection model includes a periodic GP for phenology: $f_{phen} \sim \mathcal{GP}(0, k_{phen}),$ with kernel $$k_{phen}(w, w') = \sigma^2_{phen} \exp\left(-\frac{2\sin^2(\pi|w-w'|/53)}{\ell^2_{phen}}\right)$$ to induce 53-week periodicity (Fajgenblat et al., 11 Nov 2025).

Priors for regression and variance parameters are weakly/mildly informative: $$\beta \sim \mathcal{N}(0, 3^2),\quad \sigma \sim \mathrm{Half\text{-}Normal}(0,3),\quad (\sigma^2, \ell) \sim \mathrm{InvGamma}(5,5),$$ with all coordinates/time scaled to $[-1,1]$ (Fajgenblat et al., 11 Nov 2025).

3. Computational Inference and Scalability

Model inference is performed via Bayesian posterior sampling using dynamic Hamiltonian Monte Carlo (as implemented in Stan v2.36), with careful diagnostics (e.g., traceplots, $\hat{R}<1.1$, effective sample sizes, absence of divergent transitions) (Fajgenblat et al., 11 Nov 2025). The pipeline uses 8 chains, each with 1,000 iterations and 500 warm-up steps, yielding 4,000 posterior draws.

Alternative scalable inference methods exploit Polya-Gamma data augmentation for logistic models, recasting both occupancy and detection submodels as Gaussian-conjugate updates in the Gibbs sampler. This enables efficient inference even for $S \sim 10^5-10^6$ sites and $N \sim 10^6$ visits in minutes to hours, using design-matrix sparsity and subset-of-data (SoD) approximations for spatial GPs (Diana et al., 2021).

Model Component	Bayesian MCMC (Fajgenblat et al., 11 Nov 2025)	PG-augmented Gibbs Sampler (Diana et al., 2021)
Algorithm	Stan HMC	Gibbs with Polya-Gamma latent variables
Typical Data Size	3M visits × 3k sites	$S\sim10^5-10^6$, $N\sim10^6$
Scalability Features	N/A	SoD for GP, design-matrix sparsity

4. Handling Opportunistic Citizen-Science Data

Site-occupancy models must address unique challenges in citizen science data:

• Pseudo-visits: Defined as any day-site-observer combination with at least one focal taxon record (Fajgenblat et al., 11 Nov 2025).
• Search effort: List-length (butterfly species reported) serves as a search effort proxy and enters detection-level covariates $X^p_v$.
• Observer heterogeneity: Modeled via observer-specific random intercepts $b^{obs}_o$.
• Phenology: Seasonality in detectability is accommodated by a periodic GP over calendar weeks.
• False positives: “Proficient” observers (≥500 historical records) are treated as nearly error-free, so their detections/non-detections directly inform the model likelihood, while for less proficient observers only validated records impact the likelihood. Marginalization over $z_{s,t}$ and custom likelihood formulae handle ambiguous absence/presence information (Fajgenblat et al., 11 Nov 2025).

5. Model Outputs and Retrospective Prediction

The full posterior distribution enables high-resolution spatiotemporal “backcasting” regarding occupancy ($\hat\psi_{s,t}$), spatial trends (site-specific slopes $\hat\zeta_s$), and national/provincial summaries ($\frac{1}{S} \sum_s \psi_{s,t}$) (Fajgenblat et al., 11 Nov 2025). Habitat-preference curves and marginal effect plots of covariates ($X^\psi$) elucidate drivers of occupancy. Detection model components yield observer-specific detection probabilities and phenology curves.

Retrospective predictions have been demonstrated on Belgian butterfly data (3M records), producing annual 1 km$^2$ maps for 2009–2024 and revealing, for instance, dramatic northward/range expansion of Apatura iris after 2016, strong observer effect heterogeneity (mean $p^\wedge_o$ from 1.6%–46%), and strong posterior support for positive effects of broadleaved tree and water cover on occupancy (Fajgenblat et al., 11 Nov 2025).

6. Limitations and Prospective Extensions

Limitations of the current framework include:

• Imperfect closure: Within-year movements or multivoltine generations violate the closure assumption, necessitating interpretation of $\psi_{s,t}$ as a “use probability.”
• Detection–abundance confounding: Not explicitly modeled, potentially upwardly biasing occupancy in high-abundance areas.
• Preferential sampling: Remains a source of bias if observer effort is correlated with true occupancy; attempting to mitigate via visit-level covariates.
• Effort proxies: List-length is a coarse measure.
• “This suggests …” that integrating richer effort metadata (e.g., track length, duration) and accounting for more complex spatiotemporal covariance could further improve inference (Fajgenblat et al., 11 Nov 2025).

Potential expansions comprise:

• Integration of semi-structured visit metadata.
• Multi-species occupancy modeling (e.g., via spOccupancy or flocker frameworks).
• More sophisticated or anisotropic spatiotemporal random fields.
• Formal false-positive models when computationally feasible (Fajgenblat et al., 11 Nov 2025).

7. Application and Impact

Bayesian spatiotemporal site-occupancy models substantially increase the analytical value of large-scale opportunistic biodiversity datasets by explicitly modeling imperfect detection, temporal and spatial autocorrelation, observer variation, and confounding sources of error. They provide retrospective “backcasts” of species distributions at high resolution, with quantifiable uncertainty, facilitating robust inference regarding spatiotemporal trends, habitat associations, and ecological processes in data-poor or unstructured monitoring contexts (Fajgenblat et al., 11 Nov 2025, Diana et al., 2021). These models are adaptable to a wide array of taxa, regions, and data sources.