A multi-stage Bayesian approach to fit spatial point process models

Spatial point process (SPP) models are commonly used to analyze point pattern data in many fields, including presence-only data in ecology. Existing exact Bayesian methods for fitting these models are computationally expensive because they require approximating an intractable integral each time parameters are updated and often involve algorithm supervision (i.e., tuning in the Bayesian setting). We propose a flexible, efficient, and exact multi-stage recursive Bayesian approach to fitting SPP models that leverages parallel computing resources to obtain realizations from the joint posterior, which can then be used to obtain inference on derived quantities. We outline potential extensions, including a framework for analyzing study designs with compact observation windows and a neural network basis expansion for increased model flexibility. We demonstrate this approach and its extensions using a simulation study and analyze data from aerial imagery surveys to improve our understanding of spatially explicit abundance of harbor seal (Phoca vitulina) pups in Johns Hopkins Inlet, a protected tidewater glacial fjord in Glacier Bay National Park, Alaska.

💡 Research Summary

This paper introduces a novel multi‑stage recursive Bayesian algorithm for fitting spatial point process (SPP) models, focusing on the inhomogeneous Poisson process (IPP) that allows the intensity function λ(s) to vary across space. Traditional exact Bayesian inference for such models is computationally prohibitive because each MCMC iteration requires re‑evaluating an intractable integral Λ(S)=∫_S λ(s)ds, and the lack of conjugate priors forces extensive tuning of proposal distributions. Approximate methods such as INLA provide marginal posteriors but cannot deliver exact joint posterior samples needed for derived quantities like total abundance.

The authors adopt the prior‑proposal recursive Bayesian (PPRB) framework and extend it to a three‑stage procedure that exploits parallel computing. In the first stage they condition on the observed point locations {s_i} and the observed count n, using the conditional likelihood L_cond = ∏_{i=1}^n λ(s_i)·exp(−Λ(S))^n. Because β₀ (the intercept) drops out of this likelihood, they only need to sample the regression coefficients β. This can be done via a logistic‑regression approximation (adding a large set of background points) or other sampling schemes, yielding an “interim” posterior p(β|{s_i},n). The second, intermediate stage computes Λ(S) for each draw of β in parallel, storing these values for later use. The third stage treats the interim posterior as both prior and proposal distribution for a Metropolis–Hastings update of the full parameter vector (β₀,β). The MH acceptance ratio simplifies dramatically: the point‑pattern part cancels, leaving only the Poisson likelihood for n given Λ(S) and the prior terms. Consequently, no additional evaluation of the point data is required, and the expensive integral is evaluated only once per β draw (already done in parallel).

Key advantages of this design are: (1) elimination of algorithmic tuning because the interim posterior serves as an efficient proposal; (2) massive speed‑ups through parallel evaluation of Λ(S); (3) exact joint posterior samples that enable unbiased inference on derived quantities such as total abundance N, spatially explicit intensity maps, and predictive distributions in unsampled regions. The authors also discuss extensions: (i) handling compact observation windows where points are only observed in a subset of the study region, allowing inference on the total count across the whole domain; (ii) incorporating neural‑network basis expansions to model complex, non‑linear covariate effects, thereby increasing model flexibility.

A simulation study evaluates the method under varying numbers of observed points, background points, and window sizes. Compared with a conventional single‑stage MCMC that recomputes Λ(S) at every iteration, the multi‑stage algorithm achieves 5–10× reductions in runtime while delivering comparable or superior posterior accuracy. In particular, when observation windows are small, the method markedly improves estimation of the total count N.

The approach is applied to real aerial‑imagery data from Johns Hopkins Inlet in Glacier Bay National Park, Alaska, where harbor seal (Phoca vitulina) pup locations were recorded across several disjoint survey windows. Using the multi‑stage Bayesian IPP, the authors estimate a spatial intensity surface, total pup abundance, and associated uncertainties. By augmenting the model with a neural‑network basis for covariates such as water depth and substrate type, they capture non‑linear habitat relationships that standard log‑linear models miss. The resulting abundance estimates have tighter credible intervals and more realistic spatial patterns than those obtained from traditional logistic‑regression‑based species distribution models.

Overall, the paper demonstrates that recursive Bayesian filtering, when combined with parallel computation and clever decomposition of the likelihood, can make exact Bayesian inference for spatial point processes tractable for large ecological datasets. The framework is generalizable to multivariate point processes, spatio‑temporal extensions, and more sophisticated observation error models, opening avenues for future research in ecological monitoring, seismology, and other fields where point‑pattern data are prevalent.