arXiv Statistics @arxiv_stats@qoto.org

In this paper, we present Spatialize, an open-source library that implements ensemble spatial interpolation, a novel method that combines the simplicity of basic interpolation methods with the power of classical geostatistical tools, like Kriging. It leverages the richness of stochastic modelling and ensemble learning, making it robust, scalable and suitable for large datasets. In addition, Spatialize provides a powerful framework for uncertainty quantification, offering both point estimates and empirical posterior distributions. It is implemented in Python 3.x, with a C++ core for improved performance, and is designed to be easy to use, requiring minimal user intervention. This library aims to bridge the gap between expert and non-expert users of geostatistics by providing automated tools that rival traditional geostatistical methods. Here, we present a detailed description of Spatialize along with a wealth of examples of its use.

The widespread collection of data from mobile and wearable devices has created unprecedented opportunities to study human behavior in fine temporal resolution. One common structure for such data is categorical sequences: ordered, multinomial observations across many time points. These sequences present unique statistical challenges due to their high dimensionality and complex temporal dependence, including both short- and long-term correlations. Yet, there has been relatively little methodological development focusing on principled dimension reduction specifically tailored to this type of data. In this paper, we develop and evaluate approaches to identifying "key" sequence positions which distinguish sequence types. We frame this challenge as a regression problem, introduce a variety of regularization techniques that could be applied to achieve position-based dimension reduction, and evaluate them on the motivating dataset that reflects daily time use patterns collected via a smartphone application. Results show that the thresholded LASSO, a relatively underused technique, performs better than more established methods for data with complex sequential structure.

Canonical correlation analysis (CCA) is a technique for finding correlated sets of features between two datasets. In this paper, we propose a novel extension of CCA to the online, streaming data setting: Sliding Window Informative Canonical Correlation Analysis (SWICCA). Our method uses a streaming principal component analysis (PCA) algorithm as a backend and uses these outputs combined with a small sliding window of samples to estimate the CCA components in real time. We motivate and describe our algorithm, provide numerical simulations to characterize its performance, and provide a theoretical performance guarantee. The SWICCA method is applicable and scalable to extremely high dimensions, and we provide a real-data example that demonstrates this capability.

We consider the problem of variable selection in Bayesian multivariate linear regression models, involving multiple response and predictor variables, under multivariate normal errors. In the absence of a known covariance structure, specifying a model with a non-diagonal covariance matrix is appealing. Modeling dependency in the random errors through a non-diagonal covariance matrix is generally expected to lead to improved estimation of the regression coefficients. In this article, we highlight an interesting exception: modeling the dependency in errors can significantly worsen both estimation and prediction. We demonstrate that Bayesian multi-outcome regression models using several popular variable selection priors can suffer from poor estimation properties in low-information settings--such as scenarios with weak signals, high correlation among predictors and responses, and small sample sizes. In such cases, the simultaneous estimation of all unknown parameters in the model becomes difficult when using a non-diagonal covariance matrix. Through simulation studies and a dataset with measurements from NIR spectroscopy, we illustrate that a two-step procedure--estimating the mean and the covariance matrix separately--can provide more accurate estimates in such cases. Thus, a potential solution to avoid the problem altogether is to routinely perform an additional analysis with a diagonal covariance matrix, even if the errors are expected to be correlated.

We study the zeroth-order query complexity of log-concave sampling, specifically uniform sampling from convex bodies using membership oracles. We propose a simple variant of the proximal sampler that achieves the query complexity with matched Rényi orders between the initial warmness and output guarantee. Specifically, for any $\varepsilon>0$ and $q\geq2$, the sampler, initialized at $π_{0}$, outputs a sample whose law is $\varepsilon$-close in $q$-Rényi divergence to $π$, the uniform distribution over a convex body in $\mathbb{R}^{d}$, using $\widetilde{O}(qM_{q}^{q/(q-1)}d^{2}\,\lVert\operatorname{cov}π\rVert\log\frac{1}{\varepsilon})$ membership queries, where $M_{q}=\lVert\text{d}π_{0}/\text{d}π\rVert_{L^{q}(π)}$. We further introduce a simple annealing scheme that produces a warm start in $q$-Rényi divergence (i.e., $M_{q}=O(1)$) using $\widetilde{O}(qd^{2}R^{3/2}\,\lVert\operatorname{cov}π\rVert^{1/4})$ queries, where $R^{2}=\mathbb{E}_π[|\cdot|^{2}]$. This interpolates between known complexities for warm-start generation in total variation and Rényi-infinity divergence. To relay a Rényi warmness across the annealing scheme, we establish hypercontractivity under simultaneous heat flow and translate it into an improved mixing guarantee for the proximal sampler under a logarithmic Sobolev inequality. These results extend naturally to general log-concave distributions accessible via evaluation oracles, incurring additional quadratic queries.

The skew normal law has been introduced in Azzalin (1985) as an alternative to adjusting asymmetric data that share important patterns with the normal law. It has been extensively studied. However, there is so much to do in order to catch the diversity and the richness of the investigation of its normal counterpart. The General Jarque-Berra Test (GJBT) has been devised by Lo et al. (2015), Da et al. (2023) for arbitrary laws with at least finite first eight moments, as a generalization of the Jarque-Bera (1987) test that was specially set up for normal data. Here, we particularize it to skew normal data. When particularized in the skew normal law, this test is proven to be extremely powerful in detecting the true model for any $α\neq 0$ and rejected the normal law ($α=0$) whatever be the size of the data. We introduced the use of the samples duplication method to reach a high level of efficiency for the test.

We introduce the microclustering Ewens--Pitman model for random partitions, obtained by scaling the strength parameter of the Ewens--Pitman model linearly with the sample size. The resulting random partition is shown to have the microclustering property, namely: the size of the largest cluster grows sub-linearly with the sample size, while the number of clusters grows linearly. By leveraging the interplay between the Ewens--Pitman random partition with the Pitman--Yor process, we develop efficient variational inference schemes for posterior computation in entity resolution. Our approach achieves a speed-up of three orders of magnitude over existing Bayesian methods for entity resolution, while maintaining competitive empirical performance.

A/B testing is widely used in modern technology companies for policy evaluation and product deployment, with the goal of comparing the outcomes under a newly-developed policy against a standard control. Various causal inference and reinforcement learning methods developed in the literature are applicable to A/B testing. This paper introduces a two-armed bandit framework designed to improve the power of existing approaches. The proposed procedure consists of three main steps: (i) employing doubly robust estimation to generate pseudo-outcomes, (ii) utilizing a two-armed bandit framework to construct the test statistic, and (iii) applying a permutation-based method to compute the $p$-value. We demonstrate the efficacy of the proposed method through asymptotic theories, numerical experiments and real-world data from a ridesharing company, showing its superior performance in comparison to existing methods.

In plant breeding and variety testing, there is an increasing interest in making use of environmental information to enhance predictions for new environments. Here, we will review linear mixed models that have been proposed for this purpose. The emphasis will be on predictions and on methods to assess the uncertainty of predictions for new environments. Our point of departure is straight-line regression, which may be extended to multiple environmental covariates and genotype-specific responses. When observable environmental covariates are used, this is also known as factorial regression. Early work along these lines can be traced back to Stringfield & Salter (1934) and Yates & Cochran (1938), who proposed a method nowadays best known as Finlay-Wilkinson regression. This method, in turn, has close ties with regression on latent environmental covariates and factor-analytic variance-covariance structures for genotype-environment interaction. Extensions of these approaches - reduced rank regression, kernel- or kinship-based approaches, random coefficient regression, and extended Finlay-Wilkinson regression - will be the focus of this paper. Our objective is to demonstrate how seemingly disparate methods are very closely linked and fall within a common model-based prediction framework. The framework considers environments as random throughout, with genotypes also modelled as random in most cases. We will discuss options for assessing uncertainty of predictions, including cross validation and model-based estimates of uncertainty. The methods are illustrated using a long-term rice variety trial dataset from Bangladesh.

In this paper, we develop a theoretical framework for bounding the CVaR of a random variable $X$ using another related random variable $Y$, under assumptions on their cumulative and density functions. Our results yield practical tools for approximating $\operatorname{CVaR}_α(X)$ when direct information about $X$ is limited or sampling is computationally expensive, by exploiting a more tractable or observable random variable $Y$. Moreover, the derived bounds provide interpretable concentration inequalities that quantify how the tail risk of $X$ can be controlled via $Y$.

This study integrates causal inference, graph analysis, temporal complexity measures, and machine learning to examine whether individual symptom trajectories can reveal meaningful diagnostic patterns. Testing on a longitudinal dataset of N=45 individuals affected by General Anxiety Disorder (GAD) and/or Major Depressive Disorder (MDD) derived from Fisher et al. 2017, we propose a novel pipeline for the analysis of the temporal dynamics of psychopathological symptoms. First, we employ the PCMCI+ algorithm with nonparametric independence test to determine the causal network of nonlinear dependencies between symptoms in individuals with different mental disorders. We found that the PCMCI+ effectively highlights the individual peculiarities of each symptom network, which could be leveraged towards personalized therapies. At the same time, aggregating the networks by diagnosis sheds light to disorder-specific causal mechanisms, in agreement with previous psychopathological literature. Then, we enrich the dataset by computing complexity-based measures (e.g. entropy, fractal dimension, recurrence) from the symptom time series, and feed it to a suitably selected machine learning algorithm to aid the diagnosis of each individual. The new dataset yields 91% accuracy in the classification of the symptom dynamics, proving to be an effective diagnostic support tool. Overall, these findings highlight how integrating causal modeling and temporal complexity can enhance diagnostic differentiation, offering a principled, data-driven foundation for both personalized assessment in clinical psychology and structural advances in psychological research.

The Complete Mediation Test (CMT) serves as a specialized approach of mediation analysis to assess whether an independent variable A, influences an outcome variable Y exclusively through a mediator M, without any direct effect. An application of CMT lies in Mendelian Randomization (MR) studies, where it can be used to investigate non-pleiotropy, that is, to test whether genetic variants impact a disease outcome solely through their effect on a target exposure variable. Traditionally, CMT has relied on two significance-based criteria and a proportion-based criterion with a heuristic threshold that has not been rigorously evaluated. In this paper, we explored the theoretical properties of conventional CMT, and proposed using standardized absolute proportion of mediation (SAPM) as a criterion for CMT. We, systematically assess the performance of various CMT criteria via simulation, and demonstrate their practical utility in the context of MR studies. Our results indicate that the offers the best performance. We also propose using different optimal thresholds depending on whether the mediator and outcome are continuous or binary. The SAPM with proper thresholds ensures that the indirect pathway meaningfully accounts for the effect of the exposure on the outcome, thereby strengthening the case for complete mediation.

Testing the existence of a quantitative trait locus (QTL) effect is an important task in QTL mapping studies. Most studies concentrate on the case where the phenotype distributions of different QTL groups follow normal distributions with the same unknown variance. In this paper we make a more general assumption that the phenotype distributions come from a location-scale distribution family. We derive the limiting distribution of the LRT for the existence of the QTL effect in both location and scale in genetic backcross studies. We further identify an explicit representation for this limiting distribution. As a complement, we study the limiting distribution of the LRT and its explicit representation for the existence of the QTL effect in the location only. The asymptotic properties of the LRTs under a local alternative are also investigated. Simulation studies are used to evaluate the asymptotic results, and a real-data example is included for illustration.

The recent success of deep neural network models with physical constraints (so-called, Physics-Informed Neural Networks, PINNs) has led to renewed interest in the incorporation of mechanistic information in predictive models. Statisticians and others have long been interested in this problem, which has led to several practical and innovative solutions dating back decades. In this overview, we focus on the problem of data-driven prediction and inference of dynamic spatio-temporal processes that include mechanistic information, such as would be available from partial differential equations, with a strong focus on the quantification of uncertainty associated with data, process, and parameters. We give a brief review of several paradigms and focus our attention on Bayesian implementations given they naturally accommodate uncertainty quantification. We then specify a general Bayesian hierarchical modeling framework for spatio-temporal data that can accommodate these mechanistic-informed process approaches. The advantage of this framework is its flexibility and generalizability. We illustrate the methodology via a simulation study in which a nonlinear Burgers' equation PDE is embedded via a Bayesian PINN.

During the COVID-19 pandemic, regulatory decision-making was hampered by a lack of timely and high-quality data on rare outcomes. Studying rare outcomes following infection and vaccination requires conducting multi-center observational studies, where sharing individual-level data is a privacy concern. In this paper, we conduct a multi-center observational study of thromboembolic events following COVID-19 and COVID-19 vaccination without sharing individual-level data. We accomplish this by developing a novel distributed learning method for constructing Kaplan-Meier (KM) curves and inverse propensity weighted KM curves with statistical inference. We sequentially update curves site-by-site using the KM influence function, which is a measure of the direction in which an observation should shift our estimate and so can be used to incorporate new observations without access to previous data. We show in simulations that our distributed estimator is unbiased and achieves equal efficiency to the combined data estimator. Applying our method to Beaumont Health, Spectrum Health, and Michigan Medicine data, we find a much higher covariate-adjusted incidence of blood clots after SARS-CoV-2 infection (3.13%, 95% CI: [2.93, 3.35]) compared to first COVID-19 vaccine (0.08%, 95% CI: [0.08, 0.09]). This suggests that the protection vaccines provide against COVID-19-related clots outweighs the risk of vaccine-related adverse events, and shows the potential of distributed survival analysis to provide actionable evidence for time-sensitive decision making.

Model-based clustering is widely used for identifying and distinguishing types of diseases. However, modern biomedical data coming with high dimensions make it challenging to perform the model estimation in traditional cluster analysis. The incorporation of factor analyzer into the mixture model provides a way to characterize the large set of data features, but the current estimation method is computationally impractical for massive data due to the intrinsic slow convergence of the embedded algorithms, and the incapability to vary the size of the factor analyzers, preventing the implementation of a generalized mixture of factor analyzers and further characterization of the data clusters. We propose a hybrid matrix-free computational scheme to efficiently estimate the clusters and model parameters based on a Gaussian mixture along with generalized factor analyzers to summarize the large number of variables using a small set of underlying factors. Our approach outperforms the existing method with faster convergence while maintaining high clustering accuracy. Our algorithms are applied to accurately identify and distinguish types of breast cancer based on large tumor samples, and to provide a generalized characterization for subtypes of lymphoma using massive gene records.

Compositional data, such as regional shares of economic sectors or property transactions, are central to understanding structural change in economic systems across space and time. This paper introduces a spatiotemporal multivariate autoregressive model tailored for panel data with composition-valued responses at each areal unit and time point. The proposed framework enables the joint modelling of temporal dynamics and spatial dependence under compositional constraints and is estimated via a quasi maximum likelihood approach. We build on recent theoretical advances to establish identifiability and asymptotic properties of the estimator when both the number of regions and time points grow. The utility and flexibility of the model are demonstrated through two applications: analysing property transaction compositions in an intra-city housing market (Berlin), and regional sectoral compositions in Spain's economy. These case studies highlight how the proposed framework captures key features of spatiotemporal economic processes that are often missed by conventional methods.

In settings where units' outcomes are affected by others' treatments, there has been a proliferation of ways to quantify effects of treatments on outcomes. Here we describe how many proposed estimands can be represented as involving one of two ways of averaging over units and treatment assignments. The more common representation often results in quantities that are irrelevant, or at least insufficient, for optimal choice of policies governing treatment assignment. The other representation often yields quantities that lack an interpretation as summaries of unit-level effects, but that we argue may still be relevant to policy choice. Among various estimands, the expected average outcome -- or its contrast between two different policies -- can be represented both ways and, we argue, merits further attention.

This paper introduces a novel approach, the bivariate generalized autoregressive (BGAR) model, for modeling and forecasting bivariate time series data. The BGAR model generalizes the bivariate vector autoregressive (VAR) models by allowing data that does not necessarily follow a normal distribution. We consider a random vector of two time series and assume each belongs to the canonical exponential family, similarly to the univariate generalized autoregressive moving average (GARMA) model. We include autoregressive terms of one series into the dynamical structure of the other and vice versa. The model parameters are estimated using the conditional maximum likelihood (CML) method. We provide general closed-form expressions for the conditional score vector and conditional Fisher information matrix, encompassing all canonical exponential family distributions. We develop asymptotic confidence intervals and hypothesis tests. We discuss techniques for model selection, residual diagnostic analysis, and forecasting. We carry out Monte Carlo simulation studies to evaluate the performance of the finite sample CML inferences, including point and interval estimation. An application to real data analyzes the number of leptospirosis cases on hospitalizations due to leptospirosis in São Paulo state, Brazil. Competing models such as GARMA, autoregressive integrated moving average (ARIMA), and VAR models are considered for comparison purposes. The new model outperforms the competing models by providing more accurate out-of-sample forecasting and allowing quantification of the lagged effect of the case count series on hospitalizations due to leptospirosis.

As a paradigm for sequential decision making in unknown environments, reinforcement learning (RL) has received a flurry of attention in recent years. However, the explosion of model complexity in emerging applications and the presence of nonconvexity exacerbate the challenge of achieving efficient RL in sample-starved situations, where data collection is expensive, time-consuming, or even high-stakes (e.g., in clinical trials, autonomous systems, and online advertising). How to understand and enhance the sample and computational efficacies of RL algorithms is thus of great interest. In this tutorial, we aim to introduce several important algorithmic and theoretical developments in RL, highlighting the connections between new ideas and classical topics. Employing Markov Decision Processes as the central mathematical model, we cover several distinctive RL scenarios (i.e., RL with a simulator, online RL, offline RL, robust RL, and RL with human feedback), and present several mainstream RL approaches (i.e., model-based approach, value-based approach, and policy optimization). Our discussions gravitate around the issues of sample complexity, computational efficiency, as well as algorithm-dependent and information-theoretic lower bounds from a non-asymptotic viewpoint.