ABSTRACT Multiple seasonality have been widely studied in continuous time series using models such as TBATS (Trigonometric, Box‐Cox transform, ARMA errors, Trend and Seasonal components) and MSTL (Multiple Seasonal‐Trend decomposition using LOESS). However, their treatment in categorical time series, such as air quality index (AQI) data, remains limited. Categorical AQI often exhibits distinct seasonal patterns at multiple frequencies, which are not captured by standard models. In this paper, we propose a framework that models multiple seasonality using Fourier series and indicator functions. The approach accommodates the ordinal nature of AQI categories while explicitly capturing weekly, monthly and yearly seasonal cycles. Simulation studies demonstrate the empirical consistency of parameter estimates and the forecasting performance of the proposed models under different data generating processes. We further illustrate its applicability using real categorical AQI data from Kolkata and Mumbai and compare forecasting performance with different methods.
ABSTRACT Understanding temporal dynamics in complex systems often requires identifying abrupt structural changes, known as change points in multivariate time series. Traditional vector autoregressive (VAR) models have been widely used for modeling dependencies across time, yet their parameter space grows quadratically with the number of variables, leading to computational and estimation challenges in high‐dimensional settings. The recently proposed network autoregressive (NAR) modeling framework offers a computationally efficient alternative by reducing parameter complexity through a network‐based representation. However, existing NAR models either assume homogeneous temporal behavior across all nodes, overlooking node‐specific dynamics that frequently arise in environmental and socio‐economic systems, or do not allow for structural breaks. In this work, we propose NLDNAR‐CP , a novel change point detection method within a node‐specific NAR framework that accommodates heterogeneous temporal dependencies across variables. The proposed approach efficiently detects multiple structural breaks while preserving scalability to high‐dimensional networks. We demonstrate the method's superior empirical performance through extensive simulations and a real‐world environmental application.
ABSTRACT Health benefits assessment from current exposure to fine particulate sources involves calculation of the number of Attributable adverse health events, while hypothetical future reductions in source‐specific exposure for the same population involve calculation of the number of Avoidable events. These two estimates can be very different depending on the shape of the relative risk function. We develop the mathematical justification for the calculation of Attributable and Avoidable events and characterize how they differ depending on the size of the source contribution and where on the relative risk curve exposure changes are being evaluated. The ratio of Attributable to Avoidable events is modeled as the ratio of the average derivative of the relative risk function over the total particulate mass range to the average derivative over the source‐specific mass range. This ratio can be greater or less than unity depending on the proportional size of the source to total mass and the curvature of the derivative. Relative risk functions at high global particulate concentrations can display little change in magnitude, resulting in ratios reaching into the hundreds. The difference between these estimates can be diminished if multiple particle sources are examined simultaneously. The use of highly non‐linear relative risk functions can pose risk communication challenges, with current estimates of health burden being very different than proposed improvements in air quality due to a reduction of a specific source for the same population. Widening the scope of health assessments to include multiple particulate sources can address these challenges.
ABSTRACT Urban traffic is widely recognized as a major contributor to air pollution in densely populated areas. To address this issue, many cities have implemented traffic restriction zones to reduce pollutant concentrations. Area C , a congestion charge and traffic restriction zone in central Milan, introduced in January 2012 and still in operation today, is a prominent example of such an initiative. This study investigates the causal impact of Area C on air quality by applying advanced statistical learning techniques, specifically Matrix Completion. These methods enable robust counterfactual analysis while relaxing traditional econometric assumptions, such as parallel trends. Using monthly pollution data from 2008 to 2019 in Lombardy and incorporating meteorological variables to control for confounding influences, we find a statistically significant reduction in concentrations within Area C following the policy's implementation. However, no consistent effect is observed for nitrogen oxides (), suggesting that additional or alternative interventions may be required to address gaseous pollutants. Our findings underscore the effectiveness of targeted traffic restrictions in reducing particulate pollution and highlight the value of statistical learning methods for the evaluation of environmental policy.
ABSTRACT Accurate wind speed nowcasting is crucial for optimizing wind energy production and grid stability, especially for countries such as Saudi Arabia with ambitious renewable energy targets. This article introduces a novel deep learning framework combining Deep Echo State Networks (DESNs) with spatial information through ‐nearest neighbors for high‐resolution wind speed prediction. Specifically, we develop a Spatio‐Temporal Attention Graph Autoencoder (STAGA) that effectively reduces spatial dimensionality while preserving critical spatio‐temporal patterns, enabling efficient processing of large‐scale meteorological data. Experiments using high‐resolution simulated wind data from Saudi Arabia demonstrate that our model consistently outperforms traditional methods. In addition, we provide uncertainty quantification through conformal prediction and demonstrate its practical value by assessing wind power estimates. The proposed methodology offers significant improvements for operational wind forecasting systems, supporting the efficient integration of wind energy into power grids.
ABSTRACT Understanding how environmental exposures influence health outcomes through biological pathways is critical for environmental risk assessment. We propose a Bayesian mediation analysis framework that integrates Bayesian Additive Regression Trees (BART) to model complex, nonlinear, and high‐dimensional relationships among exposures, mediators, and outcomes. The method accommodates multiple mediators, interaction effects, and uncertainty quantification through posterior inference. We introduce two estimation strategies for direct and indirect effects and develop an algorithm for computing the Deviance Information Criterion (DIC) for model evaluation. An open‐source R package, bmabart , implements the approach with visualization tools. Simulation studies demonstrate robustness in identifying significant mediators under challenging correlation structures. We apply the method to Louisiana Triple Negative Breast Cancer (TNBC) data, investigating microRNAs that mediate the association between air pollution burden, measured by the Environmental Justice Index, and cancer stage at diagnosis. Results highlight key biomarkers linking environmental exposures to disease progression, offering insights into mechanisms underlying health disparities. This framework provides a flexible and reproducible tool for environmental health research where complex mediation structures are common.
Fractional Ornstein–Uhlenbeck (fOU) processes model temporal dependence and memory, including long-range dependence, while retaining the classical Ornstein–Uhlenbeck process as a special case. We extend the integral fractional Ornstein–Uhlenbeck (ifOU) process to a multidimensional setting for animal telemetry. Longitude, Latitude, and Altitude velocities are represented by coordinate-specific fOU processes driven by a multivariate fractional Brownian motion, allowing each coordinate to retain its own damping, scale, and Hurst parameters. We establish covariance validity, characterize the admissible cross-correlation region, and derive the asymptotic behavior of cross-covariances and separated increments. We develop procedures for finite-dimensional simulation, Gaussian likelihood inference, and conditional velocity reconstruction. Replicated simulations examine estimation of cross-coordinate correlations, while joint estimation of the complete parameter vector is illustrated for one three-dimensional trajectory. The proposed model is applied to telemetry records from five common noctule bats migrating in Germany, including three trajectories with Altitude measurements.
ABSTRACT Gaussian random fields (GRFs) are often used to model spatial dependence in many geostatistical applications. However, environmental data, such as water quality measurements, often exhibit substantial skewness which violates the Gaussian assumption. Ignoring this skewness can result in biased parameter estimates and unreliable spatial predictions. To address this limitation, we propose a spatial extension of the multivariate mean‐mixture Gaussian (MMG) distribution that accommodates unbounded skewness while preserving a coherent spatial dependence structure. Our approach incorporates a global latent mixing variable to induce marginal skewness and utilizes a Gaussian process to capture spatial correlation. Identifiability is ensured by a centering transformation and replicated spatial observations, which separate mixture variance from spatial covariance. Model estimation is performed using a generalized expectation‐maximization (GEM) and quasi‐Newton algorithm, providing stable, closed‐form updates for both regression and skewness parameters. We evaluate the model through simulation studies and apply it to highly skewed geo‐referenced electrical conductivity (EC) data from groundwater monitoring stations in Golestan Province, Iran, with the primary aim of predicting values at unknown spatial locations. The results demonstrate substantial improvements over standard Gaussian and transformed Gaussian models in both estimation and prediction accuracy.
ABSTRACT Murphy et al. present a methodology to provide updated and improved estimates of the upper tail of induced earthquake magnitude distributions for the Groningen gas field, the Netherlands. In particular, they propose regression models to characterize the evolving nature of both the magnitude of completion (which corresponds to the threshold in the classical peaks over threshold extreme value model) and the exceedance magnitude. These regression approaches are routinely used to analyze environmental extremes. Our discussion focuses on some aspects of their suitability and applicability to investigate changes in extremes. We consider the risk of confounding physical and observational non‐stationarity, the difficulty of attributing changes in high quantiles to the threshold model or the exceedance model when both vary, and the role of threshold stability in estimating and interpreting the upper endpoint of the Generalized Pareto distribution.
ABSTRACT On April 6, 2009, central Italy was hit by a earthquake that caused 308 victims in the city and province of L'Aquila; subsequently, in 2016, two shocks were recorded in an area located a few dozen kilometers further north, respectively in Amatrice and Norcia. Since, like many of the physical phenomena we observe on the Earth, seismic generation processes are characterized by long‐term dependence and governed by power laws, the magnitudes of the two seismic sequences associated with these strong events were analyzed separately in the framework of nonextensive statistical mechanics to examine the connection between the variations of the magnitude probability distribution and the phases of a seismic crisis. In this work, instead, we consider all the events recorded in the same area as a whole, starting from the beginning of the most complete part of the Italian catalog ISIDe, that is, from 2005 to 2024. The aim is to verify whether the variations observed in the Tsallis entropy and in its parameters before both L'Aquila and Amatrice‐Norcia earthquakes are sufficient, as well as necessary, conditions for the occurrence of strong shocks, so that they can be considered as reliable seismic precursors. The same analysis is repeated for the corresponding data set drawn from the most recent HOmogenized instRUmental Seismic catalog (HORUS) released since July 2020. It turns out that the predictive capacity of such precursors increases when joint variations of more parameters are taken into account.
ABSTRACT In many industrial applications, such as healthcare, economics, and environmental science, unsupervised learning tasks are often associated with multivariate time series in multi‐dimensional datasets. This type of data has unique challenges and needs robust machine learning models. However, the existing literature has shown that the Euclidean distance‐based clustering techniques are ineffective. To this end, we propose two hybrid unsupervised learning models for secondary cluster analysis of multi‐dimensional multivariate time series datasets using multiple dendrograms of data matrices using fuzzy learning, traditional Dynamic Time Warping distance metric and its robust variant. A fuzzy least squares version is minimized with respect to ultrametric distance matrices and fuzzy membership degrees to obtain an optimal secondary fuzzy partition of the set of primary dendrograms. Furthermore, we also modify the previously proposed evaluation metric to assess new methods on real‐world data. Because the ground truth labels for the real‐world data are not available, particularly for the secondary partition, and the existing evaluation metric can only be used to measure similarity between exactly two dendrograms. The experimental results demonstrated the best performance of the new approaches compared to the existing Euclidean distance‐based approach on simulated and real‐world data.
ABSTRACT The rapid expansion of solar PV and wind generation intensifies the challenge of day‐ahead grid planning under stochastic meteorological conditions. While the literature is dominated by sub‐day forecasting studies, day‐ahead predictions remain critical for unit commitment and grid scheduling. This study proposes a stochastic stacked Long Short‐Term Memory (LSTM) framework, tuned via Bayesian optimization, for simultaneous day‐ahead forecasting (96 step‐ahead) of wind speed and solar irradiation using high‐resolution 15‐min measurements from the Adrar 20 MW solar PV facility in Algeria. Four input strategies were rigorously evaluated: full autoregressive lags (S‐LSTM), Random Forest importance‐based lag selection (RF‐LSTM), engineered auxiliary features (AUX‐RF‐LSTM), and exogenous meteorological variables (EXG‐RF‐LSTM), benchmarked against Holt‐Winters and naïve persistence baselines. The parsimonious RF‐LSTM architecture demonstrated superior performance across all evaluation frameworks. In deterministic point forecasting, it achieved the highest accuracy for both variables ( for wind speed and for solar irradiation) while substantially reducing residual serial dependence relative to the Holt‐Winters benchmark. These gains were confirmed as statistically significant via the Harvey‐Leybourne‐Newbold corrected Diebold‐Mariano test across all pairwise comparisons. In probabilistic forecasting via Monte Carlo Dropout, the RF‐LSTM achieved near‐perfect calibration for solar irradiation () and acceptable risk‐aware uncertainty boundaries for wind profiles. Interestingly, expanding the input space with engineered or exogenous features systematically degraded both point and probabilistic accuracy by introducing multicollinearity and optimization complexity.
Classic Bayesian methods with complex environmental models are frequently infeasible due to an intractable likelihood. Simulation‐based inference methods, such as neural posterior estimation, calculate posteriors without accessing a likelihood function by leveraging the fact that data can be quickly simulated from the model, but converge slowly and/or poorly in high‐dimensional settings. In this paper, we suggest that imposing strict variational assumptions on the form of the posterior can often combat these computational issues. Posterior distributions of model parameters are efficiently obtained by assuming a parametric form for the posterior, parametrized by the machine learning model, which is trained with the simulated data as inputs and the associated parameters as outputs. We show theoretically that if the parametric family of the variational posterior is correct, our posteriors converge to the true posteriors in Kullback–Leibler divergence. We also provide tools to help us identify if our parametric assumption is close to the true posterior, and modeling options if that is not the case. Comprehensive simulation studies using environmental models highlight our method's robustness and versatility. An analysis of the Zika virus in Brazil provides a thorough case study.
Earthquake prediction remains one of the most challenging tasks in natural hazard research due to the complexity and heterogeneity of seismic processes in space and time. To address this, we propose a deep learning (DL) framework based on graph convolutional and recurrent neural networks (RNNs) to forecast the maximum seismic magnitude in different regions of Chile. Using a cleaned and spatially segmented catalog of seismic events, we construct a graph where each node represents a seismic cluster derived from K-means clustering, with edges reflecting spatial proximity. Two models are evaluated: a standard Long Short-Term Memory (LSTM) network and a hybrid Graph Convolutional Network-LSTM (GCN-LSTM), which incorporates both temporal dynamics and spatial dependencies. Our results show that the GCN-LSTM model significantly outperforms the simple LSTM in terms of F1-score and recall, especially in regions with complex seismic activity. This demonstrates the advantage of graph-based neural models in capturing spatial correlations and improving earthquake magnitude prediction at a regional scale.
Outlier detection in functional time series is challenging due to temporal dependence and the simultaneous presence of magnitude, shape, and partial anomalies. Existing methods often assume independence or rely on model based approaches, such as the Standard Smoothed Bootstrap on Residuals (SmBoR), which may not work well if the model is misspecified. Model free alternatives, based on the moving block bootstrap, improve robustness but may detect only a limited number of magnitude anomalies. This work proposes a fully model free pipeline with two components. First, the Directional Outlyingness (DirOut) framework is extended by recalibrating its cutoff via an outlier detection procedure based on the moving block bootstratp (MBBo), improving the detection of shape and partial outliers while controlling false positives. Second, a Sliding Window Functional Boxplot (SWOD) is used to focus on local temporal neighborhoods and detect magnitude anomalies that other methods may miss. Simulations show that SWOD has high detection rates for magnitude outliers, while MBBo calibrated DirOut achieves almost perfect detection for shape and partial anomalies, outperforming SmBoR. The method is also tested on a real temperature dataset, showing its practical usefulness.
Multivariate count data are central in community ecology and related fields, where interest lies in how environmental gradients and management actions jointly shape the abundances of many taxa. The Poisson-lognormal (PLN) model is a natural workhorse in this setting, accommodating overdispersion and cross-taxon dependence via a latent Gaussian layer. Standard PLN regressions, however, treat species-specific slopes as unconstrained, which obscures situations where a covariate represents a finite resource or "budget" that is implicitly shared across taxa. We introduce a Bayesian multivariate count regression model in which the effects of selected covariates are modeled via an estimated overall effect magnitude together with compositional effect shares on a simplex. For a focal covariate, the vector of slopes across taxa is parameterized through a signed magnitude-share decomposition and endowed with a logistic-normal prior on additive logratio coordinates, so that each component can be interpreted as a species' share of the total increasing and/or decreasing effect. This construction is embedded in a multigroup PLN framework with a hierarchical prior, allowing effect shares to vary across groups while borrowing strength. Simulation studies show that the constrained Bayesian estimator improves on constrained maximum likelihood in terms of estimation error and predictive Kullback-Leibler risk, and produces more stable, interpretable effect allocations, particularly at small sample sizes. In an ecological application to dune meadow vegetation, the proposed model quantifies how a moisture gradient redistributes abundance among native grasses by allowing taxa to increase or decrease along the gradient within the same community and yields posterior predictive fits consistent with observed count patterns. The approach provides a principled way to encode and interpret effect-sharing structure in multivariate count regression, and is readily extensible to additional covariates, priors, and likelihoods.
Quantifying the temporally lagged effects of environmental hazards on health risk is key in understanding the impacts of climate change on human health. We present a simple but flexible and practical statistical modelling framework to capture the temporally distributed effects on environmental factors on aggregated health count data. The framework was designed to specifically allow use temporally aggregated health data available on a coarse resolution (such as weekly), to uncover the distributed effects of hazards on a higher temporal resolution (such as daily). The aim is to enable researchers to make use of open-access but low temporal-resolution health data to uncover health risks that would otherwise require higher resolution but inaccessible data. We focus on the example of ambient temperature and how it impacts human mortality. We use simulation experiments and implementation to 15 open source mortality data from various cities and countries, to illustrate the functionality and limitations of the approach but to also assess the efficacy of the approach on real-life data containing outliers and confounding factors. We illustrate how practical implementation in the R package mgcv enables very flexible model configurations.
High-speed winds pose several dangers to both the built environment and the natural environment, affecting human lives and livelihoods, as well as endangering wildlife. On the other hand, high-speed winds can be beneficially harnessed through the optimal orientation of turbines to maximize wind energy production. The management response to high-speed winds is a complex process because not only is their speed of interest, but even more so, their direction. High-speed winds are generally seasonal. In this paper, we propose a Bayesian prediction model for a circular variable (such as wind direction) conditional on a linear variable (such as high-speed wind), model that can account for any sparsity in the (circular, linear) pair of data. The proposed Bayesian method comprises a mixture of conditional cylindrical models to capture potential multimodality in the distribution of the circular component and compares two different cylindrical distributions, namely the Abe-Ley distribution and the Kalaylioglu distribution. Monte Carlo simulation studies conducted for both distributions establish the benefit of our proposed method. Finally, we demonstrate the developed posterior predictive distribution in a real application, where we predict the direction of high-speed winds using hourly in situ wind data from a location in South Africa.
Atlantic Meridional Overturning Circulation (AMOC) is an important component of the climate system that profoundly affects the climate of surrounding areas. AMOC is expected to slowdown due to anthropogenic climate change. The scientific community has long been at odds regarding how much AMOC has already slowed down, and whether there was any weakening at all. Key reasons for this conundrum are a limited direct observational record and the presence of substantial natural variability. The existing reconstruction methods using simple statistical approaches are not adequate to address these challenges. Moreover, they typically provide results without any uncertainty estimates. Here, we utilize a climate model ensemble experiment and a deep learning-based approach to reconstruct the long-term historical AMOC trend based on direct and indirect AMOC observations, while filtering out natural variability and providing adequate uncertainty quantification. Our analysis rules out the possibility of no decline with greater than 95% probability.
Researchers often use offsets to control for variation in exposure or effort in GLMM-style models for count data. In this paper, we study the case where effort is dependent on the underlying process of interest (e.g., Poisson intensity) and where inference focuses on prediction. We first examine a simple Poisson GLM, where we use simulation to show that dependence can result in biased GLM predictions, presumably because observations with greater offsets receive more weight within the fitting process. Possible solutions in this case include (i) including inverse offset values as weights within the fitting process, or (ii) dividing counts by the observed level of effort prior to analysis (we term these "Horvitz-Thompson-like responses"). We then consider an ecological application involving the estimation of animal abundance using transect sampling. In this case, a GLMM fitted to animal counts is used to predict abundance over a gridded study area, where environmental covariates account for variation in animal density, and the offset is a product of area surveyed and detection probability. We used simulation to assess the performance of several approaches for estimating abundance and variance when there is possible dependence between detection probability and abundance, showing that models with Horvitz-Thompson-like responses can potentially outperform other alternatives when such dependence exists. However, when applied to a beluga whale data set where turbidity was related to both detection and abundance, there seemed to be little evidence of bias. Nevertheless, we suggest that analysts first test for dependence between offsets and the underlying count process, and consider remedies such as Horvitz-Thompson responses if such dependence exists. We note connections between our research and the topic of preferential sampling in spatial statistics literature.
