A Bayesian Joint Model for Multiple Point Processes with Application to Presence-Only Data

Abstract: Joint modeling of multiple point processes is relevant in applications where relationships among processes are of interest, such as in ecological and archaeological studies. Statistical inference becomes particularly challenging when multiple processes are analyzed jointly and the observed data correspond to presence-only patterns, which are subject to preferential sampling and partial observability. This paper proposes a Bayesian joint model for multiple point processes, with application to the presence-only setting. The dependence between processes is explicitly incorporated into the probabilistic specification of the model using Bayesian networks. Direct use of the likelihood leads to intractable likelihood functions. Latent data processes are then introduced so that the augmented likelihood function becomes tractable and can be exactly evaluated. This formulation also enables direct inference on the number and the spatial distribution of unobserved occurrences of any of the point patterns. Inference is carried out using Markov chain Monte Carlo with blocked Gibbs sampling. Simulation studies demonstrate that the proposed inferential scheme is able to recover the true model parameters. The proposed model is applied to real presence-only data of archaeological sites and tree species from Amazonia, as part of the study of the effect that pre-Columbian Indigenous presence might have on the occurrences of relevant tree species. The results are consistent with the findings reported in the literature. They also illustrate how the proposed model enables inference on the existence and on the magnitude of the relation between processes, in addition to their association with environmental covariates.

• The paper presents a Bayesian joint model that explicitly encodes dependency among spatial point processes using Bayesian networks.
• It employs data augmentation and blocked Gibbs sampling to overcome intractable Poisson likelihood integrals, ensuring exact Bayesian inference.
• The model quantifies ecological and archaeological dependencies, revealing sharp occurrence variations near Amazonian features like ADEs and earthworks.

Bayesian Joint Modeling of Multivariate Point Processes: Methodology and Application to Amazonian Presence-Only Data

Introduction and Motivation

This paper introduces a comprehensive Bayesian framework for the joint modeling of multiple spatial point processes, motivated by the analytical demands encountered in ecological and archaeological studies. The primary technical focus is on presence-only datasets, where only observed occurrence locations are available, resulting in inherent preferential sampling and observation bias. Traditional methods for multivariate point process modeling—especially in ecology through joint species distribution models (JSDMs)—predominantly treat inter-process dependence via shared latent variables or sampling efforts, rather than encoding explicit, parametric dependencies.

The proposed methodology instead operationalizes dependencies using Bayesian networks: process realizations themselves function as covariates in the intensity functions of others, thus enabling identification and quantification of direct and conditional relationships in a statistically rigorous manner. Frequent intractability of the Poisson process likelihood due to integrals over complex intensity functions is addressed by a data augmentation scheme, rendering the full model amenable to exact Bayesian inference via MCMC.

Model Specification and Inference

The joint Bayesian model is predicated on a DAG (Bayesian network) that encodes conditional independence structure among $N$ point processes $X_1, \ldots, X_N$, each defined on a spatial domain $D \subset \mathbb{R}^p$. Each process is modeled as an inhomogeneous Poisson process with intensity $\lambda_i(s,\mathrm{Pa}(X_i))$, where $\mathrm{Pa}(X_i)$ denotes the set of parent processes in the network for $X_i$.

Likelihood Augmentation

To circumvent intractable integrals in the Poisson likelihood, latent processes $X_i'$ (unobserved/undetected events) and $U_i$ (auxiliary) are introduced for each node, yielding tractable hierarchical relationships. The augmentation leads to a likelihood where the intractable integrals collapse to terms involving only the domain volume and upper-bound intensities $\lambda_i^*$, critically enabling exact MCMC for the full posterior. The conditional utilization of both observed $X_j$ and latent $X_1, \ldots, X_N$0 as predictors for a child $X_1, \ldots, X_N$1 is a central feature, supporting dependency modeling via functions $X_1, \ldots, X_N$2—such as minimum interpoint distances.

Covariate and Dependency Structure

Intensity and detectability are modeled using environmental covariates (e.g., bioclimatic or topographic variables) and detection-related covariates (e.g., proximity to roads, canopy cover), linked via $X_1, \ldots, X_N$3 for intensity and $X_1, \ldots, X_N$4 for detectability. The explicit representation of observed and latent events from parent processes in these predictors allows for expressive and interpretable estimates of inter-process influences.

Bayesian Computation

Full conditionals of model parameters and latent processes are derived to exploit prior-factorization and data augmentation. Blocked Gibbs sampling, with steps for regression coefficients (using Polya–Gamma or truncated normals for logistic/probit links), latent process locations, process-specific $X_1, \ldots, X_N$5, and compliance with topological orderings in $X_1, \ldots, X_N$6, ensures efficient and valid posterior exploration.

Simulation Performance

Two simulation studies verify identifiability and inference quality for both dependence and independence cases under various proportions of observation-vs-latency and graph complexity. Across scenarios, empirical coverage rates for 90% credible intervals are generally at or above nominal levels for both regression and dependency ($X_1, \ldots, X_N$7) parameters. Notably, strong posterior support for nonzero $X_1, \ldots, X_N$8 is seen when true dependence exists, while credible intervals capture zero when dependence is absent. Performance degrades, as expected, with extreme observation sparsity, particularly for intercept estimation, but inferred dependencies remain robust.

Application: Joint Analysis of Archaeological and Botanical Data in Amazonia

The methodology is applied to real presence-only datasets comprising earthworks, Amazonian Dark Earths (ADEs), and three tree species: Handroanthus serratifolius, Bertholletia excelsa, and Dipteryx odorata. The goal is to quantify to what extent current distributions of these economically and ecologically important trees are related to pre-Columbian indigenous land-use indicators.

Figure 2: Presence-only occurrences across the Amazonian region for three tree species, earthworks, and ADE sites.

The Bayesian network is specified such that each tree species process is influenced by the union of observed and latent (i.e., unobserved) occurrences of earthworks and ADEs. The dependency term is operationalized through truncated inverse-distance functions, with parameterization reflecting the effect of proximity to each archaeological site type up to a threshold (e.g., 25 km based on ecological rationale).

Results: Inference on Ecological-Archaeological Dependencies

Posterior summaries reveal several salient patterns:

• For B. excelsa and D. odorata, the 90% credible intervals for the ADE proximity parameter ($X_1, \ldots, X_N$9) exclude zero, suggesting substantially elevated occurrence probabilities near ADEs ($D \subset \mathbb{R}^p$0: 2.01 for B. excelsa, 2.59 for D. odorata).
• Negative estimates for earthworks proximity are observed for both H. serratifolius and D. odorata ($D \subset \mathbb{R}^p$1: –0.41 and –0.20, respectively), with credible intervals excluding zero, indicating suppressed occurrence probabilities with increasing proximity to earthworks.
• B. excelsa shows only weak evidence for an earthwork effect: credible intervals are wide, including zero.
• Covariate effects for environmental variables and detectability are consistent with prior literature, supporting model adequacy.

Figure 4: Posterior response curves for presence probability in relation to distance from archaeological sites, for each tree species and each site type.

The estimated response curves for species occurrence as a function of distance to archaeological sites are steepest for ADEs, confirming strong spatial correlation at short distances. This supports the hypothesis that anthropogenic soil modifications and settlement legacies are reflected in the modern distributions of these taxa.

Estimation of Unobserved Archaeological Sites

The model further enables estimation of regions with high probability of unobserved earthworks and ADEs—quantities critical for archaeological survey planning. Posterior predictive probabilities show extensive spatial heterogeneity, reflecting both archaeological sampling bias and true heterogeneity in site occurrence.

Figure 1: Posterior predictive probabilities (by 1 km$D \subset \mathbb{R}^p$2 cell) for unobserved occurrences of earthworks and ADEs across the Amazon.

Additionally, the credible intervals for the total number of unobserved sites provide uncertainty quantification unavailable via non-joint approaches.

Implications and Future Directions

The explicit joint specification allows not only for improved estimation and prediction of presence-only point patterns but also for mechanistic interrogation of process dependencies—crucial for hypothesis-driven research in archaeology, biogeography, and conservation science. The approach is flexible with respect to covariate specification, graph topology (with variable or fixed structure), and can, in principle, be extended to semi-supervised or multi-modal settings.

Practically, this enables the integration of heterogeneous occurrence datasets—subject to differential detection and spatial bias—within a single coherent inferential pipeline. Theoretically, the parameterization directly supports causal inference (to the extent permitted by the DAG structure and domain knowledge) and the estimation of interaction strengths can facilitate structure learning or model averaging.

Future research can address computational scalability (exploiting conditional independence for parallel MCMC), incorporation of spatial autocorrelation beyond explanatory covariates, and mixed directed/undirected graphical frameworks for more general network configurations.

Conclusion

This work provides a theoretically rigorous and computationally tractable approach for multivariate point process analysis with presence-only data. By employing Bayesian networks for explicit dependence modeling, utilizing data augmentation for exact inference, and demonstrating practical value in an Amazon-scale ecological-archaeological application, this methodology establishes a robust paradigm for spatial data integration and structured discovery in the presence of latent sampling biases and interaction effects.