In this paper, we propose a Bayesian matrix-variate spatiotemporal modeling framework for jointly analyzing multiple response variables observed at spatial locations over time. The approach relaxes the standard assumption of spatial isotropy by incorporating a deformation-based mechanism, allowing the covariance structure to capture directional effects and nonstationary spatial dependence. Temporal dynamics are modeled through dynamic linear models, enabling coherent uncertainty propagation within a state-space formulation. Missing observations are handled via a data augmentation strategy that preserves the joint structure of the multivariate responses. The proposed methodology is evaluated through simulation studies and an application to air quality data. Results indicate that accounting for spatial deformation leads to substantial gains in predictive performance in anisotropic settings, while cross-variable dependence plays a secondary role in improving overall fit. The framework is computationally tractable for moderate numbers of spatial locations and responses, and provides a flexible basis for modeling multivariate spatiotemporal processes under incomplete data.

Preferential sampling is a common feature in geostatistics and occurs when the locations are sampled based on information about the phenomena under study. In this case, point pattern models are commonly used as the probability law for the distribution of the locations. However, analytic intractability of the point process likelihood prevents its direct calculation. Many Bayesian (and non-Bayesian) approaches in non-parametric model specifications handle this difficulty with approximations, both to the model and to the computations required for drawing inference. Procedures to approximate the model lead to errors that are sometimes difficult to quantify and can lead to biased inference. This paper presents an approach for performing exact Bayesian inference for this setting without the need for model approximation. A qualitatively minor change on the traditional model is proposed to circumvent the likelihood intractability. This change enables the use of an augmented model strategy. Recent work on Bayesian inference for point pattern models can be adapted to the geostatistics setting and renders computational tractability for exact inference for the proposed methodology. Estimation of model parameters and prediction of the response at unsampled locations can then be obtained from the joint posterior distribution. Simulated studies showed good quality of the proposed model for estimation and prediction in a variety of scenarios. The performance of our approach is illustrated in the analysis of simulated and real datasets and also compares favourably against approximation-based approaches. The paper is concluded with comments regarding extensions and improvements to the proposed methodology.

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.

Journal Article Accepted manuscript Dani Gamerman's contribution to the Discussion of the "Discussion Meeting on the Analysis of citizen science data" Get access Dani Gamerman Dani Gamerman UFRJ, Brazil Email: [email protected] https://orcid.org/0000-0003-0697-4589 Search for other works by this author on: Oxford Academic Google Scholar Journal of the Royal Statistical Society Series A: Statistics in Society, qnaf010, https://doi.org/10.1093/jrsssa/qnaf010 Published: 11 February 2025 Article history Received: 10 September 2024 Accepted: 16 December 2024 Published: 11 February 2025

This work considers the joint analysis of time series for epidemiological count data of neighboring regions. The joint analysis involves parameter estimation and prediction of future outcomes. The literature concentrated on imposing similarities on components of the linear predictor for the mean. However, some hierarchical model specifications for the mean contain non-linear components with similar behavior over neighboring regions. This paper proposes the use of spatial specification for these components. Parametric forms based on a data-driven approach are assumed for the waves of epidemic counts, and multiple waves are considered. The resulting model is tested in simulation studies and applied to real data. Model evaluation is based on the fitting and prediction capabilities. An illustration is provided by the analysis of counts of COVID19 cases, and it compares favorably against alternative models. Finally, the paper concludes with a discussion of the proposed methodology.

Point pattern data often exhibit features such as abrupt changes, hotspots and spatially varying dependence in local intensity. Under a Poisson process framework, these correspond to discontinuities and nonstationarity in the underlying intensity function – features that are difficult to capture with standard modeling approaches. This paper proposes a spatial Cox process model in which nonstationarity is induced through a random partition of the spatial domain, with conditionally independent Gaussian process priors specified across the resulting regions. This construction allows for heterogeneous spatial behavior, including sharp transitions in intensity. To ensure exact inference, a discretization-free MCMC algorithm is developed to target the infinite-dimensional posterior distribution without approximation. The random partition framework also reduces the computational burden typically associated with Gaussian process models. Spatial covariates can be incorporated to account for structured variation in intensity. The proposed methodology is evaluated through synthetic examples and real-world applications, demonstrating its ability to flexibly capture complex spatial structures. The paper concludes with a discussion of potential extensions and directions for future work.

Many techniques have been proposed to model space-varying observation processes with a nonstationary spatial covariance structure and/or anisotropy, usually on a geostatistical framework. Nevertheless, there is an increasing interest in point process applications, and methodologies that take nonstationarity into account are welcomed. In this sense, this work proposes an extension of a class of spatial Cox process using spatial deformation. The proposed method enables the deformation behavior to be data-driven, through a multivariate latent Gaussian process. Inference leads to intractable posterior distributions that are approximated via MCMC. The convergence of algorithms based on the Metropolis–Hastings steps proved to be slow, and the computational efficiency of the Bayesian updating scheme was improved by adopting Hamiltonian Monte Carlo (HMC) methods. Our proposal was also compared against an alternative anisotropic formulation. Studies based on synthetic data provided empirical evidence of the benefit brought by the adoption of nonstationarity through our anisotropic structure. A real data application was conducted on the spatial spread of the Spodoptera frugiperda pest in a corn-producing agricultural area in southern Brazil. Once again, the proposed method demonstrated its benefit over alternatives.

Preferential sampling is a common feature in geostatistics and occurs when the locations to be sampled are chosen based on information about the phenomena under study. In this case, point pattern models are commonly used as the probability law for the distribution of the locations. However, analytic intractability of the point process likelihood prevents its direct calculation. Many Bayesian (and non-Bayesian) approaches in non-parametric model specifications handle this difficulty with approximation-based methods. These approximations involve errors that are difficult to quantify and can lead to biased inference. This paper presents an approach for performing exact Bayesian inference for this setting without the need for model approximation. A qualitatively minor change on the traditional model is proposed to circumvent the likelihood intractability. This change enables the use of an augmented model strategy. Recent work on Bayesian inference for point pattern models can be adapted to the geostatistics setting and renders computational tractability for exact inference for the proposed methodology. Estimation of model parameters and prediction of the response at unsampled locations can then be obtained from the joint posterior distribution of the augmented model. Simulated studies showed good quality of the proposed model for estimation and prediction in a variety of preferentiality scenarios. The performance of our approach is illustrated in the analysis of real datasets and compares favourably against approximation-based approaches. The paper is concluded with comments regarding extensions of and improvements to the proposed methodology.

Abstract The Amazon rainforest is not the untouched wilderness we think it is, but contains thousands of ancient and mysterious earthworks left behind by indigenous populations. Dani Gamerman describes how statisticians and statistical analysis brought numerical evidence to an interdisciplinary paper that made headlines around the world.

Summary We present a novel inference methodology to perform Bayesian inference for spatiotemporal Cox processes where the intensity function depends on a multivariate Gaussian process. Dynamic Gaussian processes are introduced to enable evolution of the intensity function over discrete time. The novelty of the method lies on the fact that no discretization error is involved despite the non-tractability of the likelihood function and infinite dimensionality of the problem. The method is based on a Markov chain Monte Carlo algorithm that samples from the joint posterior distribution of the parameters and latent variables of the model. A particular choice of the dominating measure to obtain the likelihood function is shown to be crucial to devise a valid Markov chain Monte Carlo algorithm. The models are defined in a general and flexible way but they are amenable to direct sampling from the relevant distributions because of careful characterization of its components. The models also enable the inclusion of regression covariates and/or temporal components to explain the variability of the intensity function. These components may be subject to relevant interaction with space and/or time. Real and simulated examples illustrate the methodology, followed by concluding remarks.

Indigenous societies are known to have occupied the Amazon basin for more than 12,000 years, but the scale of their influence on Amazonian forests remains uncertain. We report the discovery, using LIDAR (light detection and ranging) information from across the basin, of 24 previously undetected pre-Columbian earthworks beneath the forest canopy. Modeled distribution and abundance of large-scale archaeological sites across Amazonia suggest that between 10,272 and 23,648 sites remain to be discovered and that most will be found in the southwest. We also identified 53 domesticated tree species significantly associated with earthwork occurrence probability, likely suggesting past management practices. Closed-canopy forests across Amazonia are likely to contain thousands of undiscovered archaeological sites around which pre-Columbian societies actively modified forests, a discovery that opens opportunities for better understanding the magnitude of ancient human influence on Amazonia and its current state.

This paper provides an exact modeling approach for the analysis of presence-only ecological data. Our proposal is also based on frequently used inhomogeneous Poisson processes but does not rely on model approximations, unlike other approaches. Exactness is achieved via a data augmentation scheme. One of the augmented processes can be interpreted as the unobserved occurrences of the relevant species, and its posterior distribution can be used to make predictions of the species over the region of study beyond the observer bias. The data augmentation also leads to a natural Gibbs sampler to make Bayesian inference through MCMC. The proposal shows better performance than the currently standard method based on Poisson process with intensity function depending log-linearly on the covariates. Additionally, an identification problem that arises in the traditional model does not seem to affect our proposal in the analyses of real ecological data.

Detailed knowledge on the effects of air pollutants on human health is a prerequisite for the development of effective policies to reduce the adverse impact of ambient air pollution. However, measuring the effect of exposure on health outcomes is an extremely difficult task as the health impact of air pollution is known to vary over space and over different exposure periods. In general, standard approaches aggregate the information over space or time to simplify the study but this strategy fails to recognize important regional differences and runs into the well-known risk of confounding the effects. However, modelling directly with the original, disaggregated data requires a highly dimensional model with the curse of dimensionality making inferences unstable; in these cases, the models tend to retain many irrelevant components and most relevant effects tend to be attenuated. The situation clearly calls for an intermediate solution that does not blindly aggregate data while preserving important regional features. We propose a dimension-reduction approach based on latent factors driven by the data. These factors naturally absorb the relevant features provided by the data and establish the link between pollutants and health outcomes, instead of forcing a necessarily high-dimensional link at the observational level. The dynamic structural equation approach is particularly suited for this task. The latent factor approach also provides a simple solution to the spatial misalignment caused by using variables with different spatial resolutions and the state-space representation of the model favours the application of impulse response analysis. Our approach is discussed through the analysis of the short-term effects of air pollution on hospitalization data from Lombardia and Piemonte regions (Italy).

The number of packages/software for Gaussian State Space models has increased over recent decades. However, there are very few codes available for non-Gaussian State Space (NGSS) models due to analytical intractability that prevents exact calculations. One of the few tractable exceptions is the family of NGSS with exact marginal likelihood, named NGSSEML. In this work, we present the wide range of data formats and distributions handled by NGSSEML and a package in the R language to perform classical and Bayesian inference for them. Special functions for filtering, forecasting, and smoothing procedures and the exact calculation of the marginal likelihood function are provided. The methods implemented in the package are illustrated for count and volatility time series and some reliability/survival models, showing that the codes are easy to handle. Therefore, the NGSSEML family emerges as a simple and interesting option/alternative for modeling non-Gaussian time-varying structures commonly encountered in time series and reliability/survival studies.

Inference over multivariate tails often requires a number of assumptions which may affect the assessment of the extreme dependence structure. Models are usually constructed in such a way that extreme components can either be asymptotically dependent or be independent of each other. Recently, there has been an increasing interest on modelling multivariate extremes more flexibly, by allowing models to bridge both asymptotic dependence regimes. Here we propose a novel semiparametric approach which allows for a variety of dependence patterns, be them extremal or not, by using in a model-based fashion the full dataset. We build on previous work for inference on marginal exceedances over a high, unknown threshold, by combining it with flexible, semiparametric copula specifications to investigate extreme dependence, thus separately modelling marginals and dependence structure. Because of the generality of our approach, bivariate problems are investigated here due to computational challenges, but multivariate extensions are readily available. Empirical results suggest that our approach can provide sound uncertainty statements about the possibility of asymptotic independence, and we propose a criterion to quantify the presence of either extreme regime which performs well in our applications when compared to others available. Estimation of functions of interest for extremes is performed via MCMC algorithms. Attention is also devoted to the prediction of new extreme observations. Our approach is evaluated through simulations, applied to real data and assessed against competing approaches. Evidence demonstrates that the bulk of the data do not bias and improve the inferential process for extremal dependence in our applications.

Point processes are one of the most commonly encountered observation processes in Spatial Statistics. Model-based inference for them depends on the likelihood function. In the most standard setting of Poisson processes, the likelihood depends on the intensity function, and can not be computed analytically. A number of approximating techniques have been proposed to handle this difficulty. In this paper, we review recent work on exact solutions that solve this problem without resorting to approximations. The presentation concentrates more heavily on discrete time but also considers continuous time. The solutions are based on model specifications that impose smoothness constraints on the intensity function. We also review approaches to include a regression component and different ways to accommodate it while accounting for additional heterogeneity. Applications are provided to illustrate the results. Finally, we discuss possible extensions to account for discontinuities and/or jumps in the intensity function.

This article discusses the use of a Bayesian model that incorporates differential item functioning (DIF) in analysing whether cultural differences may affect the performance of students from different countries in the various test items which make up the OECD’s Programme for International Student Assessment (PISA) test of mathematics ability. The PISA tests in mathematics and other subjects are used to compare the educational attainment of fifteen-year old students in different countries. The article first provides a background on PISA, DIF and item response theory (IRT) before describing a hierarchical three-parameter logistic model for the probability of a correct response on an individual item to determine the extent of DIF remaining in the mathematics test of 2003. The results of Bayesian analysis illustrate the importance of appropriately accounting for all sources of heterogeneity present in educational testing and highlight the advantages of the Bayesian paradigm when applied to large-scale educational assessment.

Inference over tails is performed by applying only the results of extreme value theory. Whilst such theory is well defined and flexible enough in the univariate case, multivariate inferential methods often require the imposition of arbitrary constraints not fully justifed by the underlying theory. In contrast, our approach uses only the constraints imposed by theory. We build on previous, theoretically justified work for marginal exceedances over a high, unknown threshold, by combining it with flexible, semiparametric copulae specifications to investigate extreme dependence. Whilst giving probabilistic judgements about the extreme regime of all marginal variables, our approach formally uses the full dataset and allows for a variety of patterns of dependence, be them extremal or not. A new probabilistic criterion quantifying the possibility that the data exhibits asymptotic independence is introduced and its robustness empirically studied. Estimation of functions of interest in extreme value analyses is performed via MCMC algorithms. Attention is also devoted to the prediction of new extreme observations. Our approach is evaluated through a series of simulations, applied to real data sets and assessed against competing approaches. Evidence demonstrates that the bulk of the data does not bias and improves the inferential process for the extremal dependence.