Papers

We describe a new family of experimental designs that extends the principle of stratified randomization to settings with continuous, constrained multivariate, and other irregular treatment spaces. Our approach is to first match units into homogeneous groups, then use Monte Carlo couplings to assign within-group treatments to be highly dispersed over the treatment space. We show that ensuring similar units receive dissimilar treatments improves estimation efficiency. The efficiency gains are proportional to the product of dispersion and match quality, where dispersion measures how spread out the assignments are relative to independent randomization. We develop a new spectral analysis showing how efficiency depends on alignment between the smoothness and shape of the estimator's influence function and the coupling's principal directions. We illustrate these designs with examples from development, behavioral, and labor economics. In particular, our empirical application uses data from a real experiment allocating savings monitors using their position within village social networks.

Hamiltonian Monte Carlo (HMC) and its dynamic extensions, such as the No-U-Turn Sampler (NUTS), are powerful Markov chain Monte Carlo methods for sampling from complex, high-dimensional probability distributions. Riemannian manifold Hamiltonian Monte Carlo (RMHMC) extends HMC by allowing the mass matrix to depend on position, which can substantially improve mixing but also makes implementation considerably more challenging. In this paper, we study an adaptive hierarchical version of RMHMC that is well suited to many hierarchical sampling problems. A key feature of hierarchical RMHMC is that, unlike general RMHMC, it admits a closed-form explicit leapfrog integrator, enabling efficient implementation and direct use within dynamic HMC methods such as NUTS. We introduce an adaptive scheme that automatically tunes the parameters of the hierarchical mass matrix during simulation. Importantly, the target density need not exhibit any hierarchical or block structure; the hierarchy is instead imposed on the mass matrix as a modeling device to capture the local geometry of the target distribution. Numerical experiments demonstrate appealing empirical performance in high-dimensional Bayesian inference problems.

The sequential nature of autoregressive next-token prediction imposes a fundamental speed limit on large language models. While continuous flow models offer a path to parallel generation, they traditionally demand expensive iterative integration. Flow Maps bypass this bottleneck by compressing generative trajectories into single-step mappings, theoretically enabling the generation of full text sequences from noise in a single forward pass. However, standard formulations rely on Euclidean regression losses that are geometrically ill-suited for discrete data. In this work, we resolve this conflict with Discrete Flow Maps, a framework that reconciles trajectory compression with the geometry of the probability simplex. We recast standard flow map training for the discrete domain, aligning the training dynamics with the discrete nature of language. Empirically, this strict geometric alignment allows our method to surpass previous state-of-the-art results in discrete flow modeling.

In this article, we review selective inference, a set of techniques for inference when the statistical question asked is a function of the data. This setting often arises in contemporary scientific workflows, where hypotheses and parameters may be selected from the data, rather than specified in advance. In this setting, classical inferential techniques do not achieve "classical" guarantees, such as nominal coverage of confidence intervals. We consider three examples for which selective inference solutions are required: inference on a "winner", inference on the mean of a region in a regression tree, and inference on the difference in means between a pair of clusters. We argue that conditional guarantees are of scientific interest in such settings. We then review and draw connections between several approaches that provide such guarantees. Finally, we illustrate these approaches in simulation and through an application to single-cell RNA sequencing data.

The summer of 2023 brought record-breaking heat extremes across the globe. In this work, we investigate how much worse these extremes might have become, given identical large-scale conditions. Numerical weather prediction (NWP) ensembles simulate alternate heatwave storylines, but their limited sample size diminishes their statistical coverage of observed heatwaves and ability to represent worst-case events. Using the Spherical Fourier Neural Operator machine learning (ML) weather model, we generated a huge ensemble of 7,424 storyline forecasts of summer temperature extremes. Despite being constrained by the same initial and boundary conditions and using training data from the same modeling system, the ML ensemble produced extreme heatwave scenarios exceeding temperatures from NWP ensembles. For 30% of global land areas, ML-generated extreme events were outside the prediction envelope from extrapolating smaller NWP ensembles with extreme value theory. This research demonstrates that very large ensembles substantially improve upon NWP-based characterization of worst-case heat extremes.

We study harmonic map regression, a nonparametric estimator for manifold-valued responses, that penalizes the empirical Fréchet risk by the Dirichlet energy. By connecting penalized regression to the theory of harmonic maps, the estimator acquires a structural theory that parallels the classical Euclidean smoothing spline. The Euler-Lagrange equation characterizes the solution as a piecewise-geodesic spline, an equivalent kernel controls pointwise risk at the rate $n^{-2/3}$, and the infinite-dimensional variational problem reduces exactly to a finite-dimensional optimization. Such newly established connection reveals a topological phenomenon that has no analogue in Euclidean nonparametric regression and, to our knowledge, has not been studied in the manifold regression literature. On manifolds whose regression curves can wrap around in topologically distinct ways, maps in distinct homotopy classes are separated by energy barriers intrinsic to the geometry of the target, and the Dirichlet penalty makes the estimator sensitive to this structure, recovering the correct topological class with probability tending to one, a phase transition we call topological recovery. A curvature-dependent oracle inequality yields the minimax rate $n^{-2s/(2s+1)}$ for Sobolev order $s$, matching the Euclidean constant on non-positively curved targets, while five geometric obstructions show that the full structural theory is unique to the Dirichlet energy ($s=1$). Simulations on $S^2$, $\mathbb{H}^2$, $SO(3)$, $\mathrm{Sym}^+(2)$, and $T^2$ corroborate the theory, and an application to wind-direction data on $S^1$ illustrates practical advantages.

Multi-arm multi-stage (MAMS) trials have gained popularity, due to their improved efficiency in evaluating multiple treatments. A traditional MAMS trial often decreases the expected sample size of the trial compared to just running a multi-arm approach, but with the drawback of an increase in maximum sample size. For academic led trials this poses a particular challenge, as funding is typically based on the maximum required sample size. To address this, drop-the-loser designs were introduced, where a fixed number of treatments are dropped at each interim stage, thereby reducing the maximum sample size. In this work, we propose an enhanced multi-stage drop-the-loser design that also allows for early stopping of the entire trial for superiority. This approach aims to retain the benefits of a reduced maximum sample size while also lowering the expected sample size. The proposed design is motivated by a trial in atrial fibrillation. We derive analytical expressions for the type I error rate, power, and expected sample size, and compare the proposed design's performance to alternative methods. We outline the key requirements for implementing the proposed design and discuss the contexts in which it should be considered. For the motivating example the results show that the proposed design substantially reduces the expected sample size compared to a standard drop-the-loser design, while lowering the maximum sample size relative to running a traditional MAMS trial or multiple separate trials.

A learning-to-defer (L2D) system decides, for each input, whether to predict on its own or to hand it to one of several available experts. The very well established recipe trains classifier and router jointly by treating the $K$ classes and $J$ experts as competing actions in one shared $(K{+}J)$-action geometry. Subsequent work has proposed a series of incremental fixes within this geometry; we show that each still suffers, to varying severity, from an optimization-level pathology (target distortion, gradient amplification, winner-take-all starvation, set-mass collapse, or class-expert coupling) even under statistical consistency. We step outside the augmented-action family entirely and propose a decoupled surrogate: a softmax classifier head and an independent sigmoid head per expert, mirroring the two natural objects of the problem. We show that per-sample updates are then coordinatewise and the class-expert Hessian block is identically zero, and prove an excess-risk bound with calibration constant $\max\{2\sqrt{2},\sqrt{2J/\lambda}\}$ – to our knowledge the first multi-expert L2D guarantee whose constant does not grow with the expert pool when the per-expert weight is held fixed. On controlled synthetic studies and on CIFAR-10, CIFAR-10H, and Covertype, it is the only method in our comparison that remains stable as the expert pool grows, preserves rare specialists, and improves over a standalone classifier on every real-data benchmark.

We study the population loss landscape of two-layer ReLU networks of the form $\sum_{k=1}^K \mathrm{ReLU}(w_k^\top x)$ in a realisable teacher-student setting with Gaussian covariates. We show that local minima admit an exact low-dimensional representation in terms of summary statistics, yielding a sharp and interpretable characterisation of the landscape. We further establish a direct link with one-pass SGD: local minima correspond to attractive fixed points of the dynamics in summary statistics space. This perspective reveals a hierarchical organisation of minima into discrete families and shows how overparameterisation changes their stability and reachability under gradient-based dynamics. In this overparameterised regime, global minima become increasingly accessible, attracting the dynamics and reducing convergence to spurious solutions. Overall, our results reveal intrinsic limitations of common simplifying assumptions, which may miss essential features of the loss landscape even in minimal neural network models.

Comparing $K$-sample distributions is a fundamental problem in data science that arises in a wide variety of fields and applications. In this article, we introduce a maximum-of-differences approach to make such comparisons. Specifically, we first calculate the pairwise distances from the pooled observations of the $K$ samples. We then define the two observations as connected if their distance is less than a pre-specified threshold value. For each observation, we next calculate the “within" and the “between" probabilities associated with these two types of connections for the given observation, i.e., with other observations within the same sample and between the given observation and the observations in other samples. Subsequently, we propose a maximum-of-differences (MOD) test that finds the maximum value among the standardized squared differences between the “within" and the “between" probabilities of all observations. Accordingly, the proposed test is not only applicable to multivariate data with $K$ samples, but can also be extended to multivariate regression models. Furthermore, we obtain the covariance-adjusted (CA) version of the MOD (CA-MOD) test, which converges to the Type I extreme value distribution under some conditions. Moreover, we demonstrate the asymptotic properties of the two tests under both the null and alternative hypotheses. The performance and usefulness of the tests are illustrated via simulation studies and real examples.

Accurate representation of non-Gaussian distributions of quantities of interest in nonlinear dynamical systems is critical for estimation, control, and decision-making, but can be challenging when forward propagations are expensive to carry out. This paper presents an approach for estimating probability density functions of states evolving under nonlinear dynamics using Seminonparametric (SNP), or Gallant-Nychka, densities. SNP densities employ a probabilists' Hermite polynomial basis to model non-Gaussian behavior and are positive everywhere on the support by construction. We use Monte Carlo to approximate the expectation integrals that arise in the maximum likelihood estimation of SNP coefficients, and introduce a convex relaxation to generate effective initial estimates. The method is demonstrated on density and quantile estimation for the chaotic Lorenz system. The results demonstrate that the proposed method can accurately capture non-Gaussian density structure and compute quantiles using significantly fewer samples than raw Monte Carlo sampling.

Single-cell transcriptomic data approximates the abundance of proteins at a high resolution, but its noisiness necessitates transformation by a pipeline of methods before analysis and inference. In the absence of robust validation of these pipelines and methods, it remains unclear how best to process any particular dataset. To compensate for this, popular visualisation methods, e.g., t-SNE and UMAP, are commonly used to produce descriptions of datasets. Such visualisations are incomplete and provide subjective descriptions of samples rather than statistically meaningful statements about technical noise or biology. In this paper, we introduce the Zero-Inflated Negative-Binomial with Geometric Tail (ZINBGT), a mixture-model-based strategy for producing interpretable visualisations of each gene's expression across cells, along with diagnostic summaries that use Wasserstein distance to highlight outlier genes. These diagnostics are used to reveal an outlier gene within a T. brucei sample. This method is applied to a human immune-cell dataset, highlighting the relationship between sparsity, mean, and spread across genes, as well as revealing an issue with the use of zero-inflated negative-binomial distributions to model single-cell RNA data. An investigation of simulated datasets intended to replicate the immune-cell data revealed discrepancies with the ground truth, establishing purposes for which these simulated datasets are unsuitable. Finally, we list a number of different domains to which this method can be applied.

Testing whether a multivariate distribution is invariant under an orthogonal transformation is a classical problem when the transformation is fixed in advance. We address a harder situation: the transformation is unknown and must be inferred from the data. We frame this as invariance under a group action whose representation is indexed by an unknown parameter. We work in $\mathbb{R}^d$ under a simple-spectrum assumption on the covariance matrix $\Sigma$. Under this assumption, any orthogonal transformation that preserves the distribution must commute with $\Sigma$, and therefore must be a {reflection through a subspace spanned by a subset of the principal directions. This reduces the search over the orthogonal group to a finite family of candidate reflections, one for each subset of principal axes.} For each candidate, we construct a Kolmogorov–Smirnov-type statistic based on projected data and sample splitting. We derive its asymptotic distribution in a triangular-array framework and establish bootstrap validity under suitable regularity conditions. Axial symmetry about an unspecified direction (that is, invariance under reflection across an unknown one-dimensional subspace) and hyperplane (Householder) symmetry about an unspecified normal direction are the two leading particular cases: we treat them in detail, the former driving the simulation study and the latter the real-data application, to show how a concrete problem involving an unknown group representation can be successfully addressed.

The Half-Trek Criterion (HTC) is the primary graphical tool for determining generic identifiability of causal effect coefficients in linear structural equation models (SEMs) with latent confounders. However, HTC is inherently node-wise: it simultaneously resolves all incoming edges of a node, leaving a gap of "inconclusive" causal effects (15-23% in moderate graphs). We introduce Iterative Identification Closure (IIC), a general framework that decouples causal identification into two phases: (1) a seed function S_0 that identifies an initial set of edges from any external source of information (instrumental variables, interventions, non-Gaussianity, prior knowledge, etc.); and (2) Reduced HTC propagation that iteratively substitutes known coefficients to reduce system dimension, enabling identification of edges that standard HTC cannot resolve. The core novelty is iterative identification propagation: newly identified edges feed back to unlock further identification – a mechanism absent from all existing graphical criteria, which treat each edge (or node) in isolation. This propagation is non-trivial: coefficient substitution alters the covariance structure, and soundness requires proving that the modified Jacobian retains generic full rank – a new theoretical result (Reduced HTC Theorem). We prove that IIC is sound, monotone, converges in O(|E|) iterations (empirically <=2), and strictly subsumes both HTC and ancestor decomposition. Exhaustive verification on all graphs with n<=5 (134,144 edges) confirms 100% precision (zero false positives); with combined seeds, IIC reduces the HTC gap by over 80%. The propagation gain is gamma~4x (2 seeds identifying ~3% of edges to 97.5% total identification), far exceeding gamma<=1.2x of prior methods that incorporate side information without iterative feedback.

We propose an adaptive MCMC method that learns a linear preconditioner which is dense in its off-diagonal elements but sparse in its parametrisation. Due to this sparsity, we achieve a per-iteration computational complexity of $O(m^2d)$ for a user-determined parameter $m$, compared with the $O(d^2)$ complexity of existing adaptive strategies that can capture correlation information from the target. Diagonal preconditioning has an $O(d)$ per-iteration complexity, but is known to fail in the case that the target distribution is highly correlated, see \citet[Section 3.5]{hird2025a}. Our preconditioner is constructed using eigeninformation from the target covariance which we infer using online principal components analysis on the MCMC chain. It is composed of a diagonal matrix and a product of carefully chosen reflection matrices. On various numerical tests we show that it outperforms diagonal preconditioning in terms of absolute performance, and that it outperforms traditional dense preconditioning and multiple diagonal plus low-rank alternatives in terms of time-normalised performance.

We propose a Bayesian framework for planning simple step-stress accelerated life tests when items are subject to two independent competing failure modes We assume that the competing risks are independent, with lifetimes following Weibull distributions, and adopt the cumulative exposure model with a log-linear stress-life relationship to connect failure time distributions across stress levels. The optimality criterion is the preposterior variance of the $p$-th quantile of the lifetime distribution at use stress, evaluated without reliance on asymptotic approximations, making the methodology valid regardless of sample size. Building on the idea of quantile-based reparametrisation used in single-mode ALT \citep{zhang2006bayesian}, we extend this approach to the competing risks setting by reparametrising the model parameters for each failure mode to physically interpretable and approximately independent quantities, making it possible to elicit priors directly from engineering knowledge of device behaviour. Posterior inference is carried out using the No-U-Turn Sampler implemented in Stan, and the optimal design is located via Monte Carlo simulation over a grid of candidate designs. The methodology is illustrated on a real step-stress dataset for a solar lighting device subject to capacitor and controller failure modes. A comprehensive sensitivity analysis with respect to the quantile probability, the lower stress level, the prior hyperparameter specification, and the sample size shows that the optimal stress-change time is moderately sensitive to these inputs while the optimal lower stress level consistently favours operation close to use conditions, a finding that holds across all prior specifications considered.

We argue that Bonferroni correction is a better choice for online experimentation than it is commonly given credit for. The case rests on four considerations. First, it is the simplest broadly implementable FWER-controlling method that produces unconditional simultaneous confidence intervals for every metric. Second, in a well-specified decision framework, guardrail and quality metrics use intersection-union logic and cannot inflate the false positive rate, so the Bonferroni denominator is the number of success metrics only, not the total metric count. Third, it is uniquely tractable for pre-experiment sample size calculations. Fourth, we contextualise the power cost empirically. Drawing on a simulation study and an empirical analysis of 1,296 experiments run on Spotify's experimentation platform, Confidence, we show that the power loss relative to more sophisticated FWER methods depends on both how the correction family is specified and how many metrics are truly non-null. When guardrail metrics are incorrectly included in the family, Holm and Hommel are nearly indistinguishable from Bonferroni. When the family is correctly restricted to success metrics only, they gain roughly 4–5 percentage points in ship rate (the fraction of experiments where the treatment is deployed). When few metrics are truly non-null, the gap narrows to near zero regardless of method.

We develop a predictive-first optimisation framework for streaming hidden Markov models. Unlike classical approaches that prioritise full posterior recovery under a fully specified generative model, we assume access to regime-specific predictive models whose parameters are learned online while maintaining a fixed transition prior over regimes. Our objective is to sequentially identify latent regimes while maintaining accurate step-ahead predictive distributions. Because the number of possible regime paths grows exponentially, exact filtering is infeasible. We therefore formulate streaming inference as a constrained projection problem in predictive-distribution space: under a fixed hypothesis budget, we approximate the full posterior predictive by the forward-KL optimal mixture supported on $S$ paths. The solution is the renormalised top-$S$ posterior-weighted mixture, providing a principled derivation of beam search for HMMs. The resulting algorithm is fully recursive and deterministic, performing beam-style truncation with closed-form predictive updates and requiring neither EM nor sampling. Empirical comparisons against Online EM and Sequential Monte Carlo under matched computational budgets demonstrate competitive prequential performance.

We develop a theory of generalization and scaling for Mixture-of-Experts (MoE) Transformers that cleanly separates active per-input capacity from routing combinatorics. By conditioning on fixed routing patterns and union-bounding across them, we derive a sup-norm covering-number bound whose metric entropy scales with the active parameter budget and incurs a MoE-specific routing overhead. Combined with a standard ERM analysis for squared loss, this yields a generalization bound under a $d$-dimensional manifold data model and $C^\beta$ targets, showing that approximation and estimation trade off as in dense networks once active parameters are accounted for appropriately. We further prove a constructive approximation theorem for MoE architectures, showing that, under the approximation construction, error can decrease either by scaling active capacity or by increasing the number of experts, depending on the dominant bottleneck. From these results we derive neural scaling laws for model size, data size, and compute-optimal tradeoffs. Overall, our results provide a transparent statistical reference point for reasoning about MoE scaling, clarifying which behaviors are certified by worst-case theory and which must arise from data-dependent routing structure or optimization dynamics.

This paper introduces a score-driven rating system, a generalization of the classical Elo rating system that employs the score, i.e. the gradient of the log-likelihood, as the updating mechanism for player and team ratings. The proposed framework extends beyond simple win/loss game outcomes and accommodates a wide range of game results, such as point differences, win/draw/loss outcomes, or complete rankings. Theoretical properties of the score are derived, showing that it has zero expected value, sums to zero across all players, and decreases with increasing value of a player's rating, thereby ensuring internal consistency and fairness. Furthermore, the score-driven rating system exhibits a reversion property, meaning that ratings tend to follow the underlying unobserved true skills over time. The proposed framework provides a theoretical rationale for existing dynamic models of sports performance and offers a systematic approach for constructing new ones.

Unobserved confounding is a key challenge when estimating causal effects from a treatment on an outcome. In this work, we assume that we observe a single, potentially multi-dimensional proxy variable of the unobserved confounder and that we know the mechanism that generates the proxy from the confounder. Under an assumption called Single Proxy Identifiability of Causal Effects or simply SPICE, we prove that this error mechanism is complete and causal effects are identifiable. We extend the proxy-based causal identifiability results by Kuroki and Pearl (2014); Pearl (2010) to multi-dimensional continuous settings, more flexible functional relationships and a broader class of distributions. Further, we develop a neural network based estimation framework, SPICE-Net, to estimate causal effects, which is applicable to both discrete and continuous treatments.

The most dangerous error in clinical trial interpretation is equating p > 0.05 with no effect. This review provides a practical, algorithm-based framework for classifying randomized controlled trial (RCT) results into six distinct categories positive, imprecise (+), neutral, inconclusive, negative, and harmful using confidence interval (CI) position relative to the minimal clinically important difference (MCID) as the primary tool, augmented by Bayesian posterior probabilities. We demonstrate that the same p > 0.05 result can represent three fundamentally different conclusions (inconclusive, negative, or neutral), show how Bayesian reanalysis can rescue benefit signals missed by frequentist thresholds, and illustrate the framework with real-world examples from critical care and cardiology trials. The framework synthesizes guidance from Altman, Harrell, Pocock, Zampieri, the ASA, and ICH E9 into a single coherent decision algorithm.

We study community detection in stochastic block models under pure node-level differential privacy, a stringent notion that protects the participation of an individual together with all of their incident edges. This setting is substantially more challenging than edge-private community detection, since modifying a single node can affect linearly many observations. On the algorithmic side, we analyze a node-private estimator based on the exponential mechanism combined with an extension lemma, and show that exact recovery remains achievable. In the standard sparse regime with logarithmic average degree and a fixed number of communities, our results imply that a logarithmic privacy budget suffices to obtain nontrivial recovery guarantees. On the lower bound side, we show that this logarithmic scaling is in fact unavoidable: any pure node-private method must fail to achieve polynomially small exact-recovery error, or polynomially small expected mismatch, unless the privacy budget is at least of this order. Moreover, in the regime of super-logarithmic privacy budgets, our upper and lower bounds yield a matching two-term characterization of the minimax risk, with one term governed by the non-private statistical signal and the other by the privacy budget; these match up to universal constants in the exponents. Taken together, our results identify an inherent logarithmic privacy cost in node-private community detection, absent under edge differential privacy, and provide a precise rate-level characterization of the tradeoff between node privacy and SBM recovery.

Constructing valid inferential methods for constrained parameters in normal and Poisson distributions represents two fundamental and important problems in applied statistics, for which there is currently no unified framework for statistical inference. Most existing studies assume that the nuisance parameters of the model are known, an assumption that is often impractical in real-world applications. However, under the more realistic scenario where nuisance parameters are unknown, the available Bayesian interval estimation methods fail to guarantee nominal coverage and thus cannot provide exact inference. To address these limitations, this paper develops prior-free inferential model (IM) approaches for parameters of interest in constrained normal and Poisson models and demonstrates that the confidence intervals (CIs) obtained from these novel IM methods can achieve exact nominal coverage. Furthermore, considering the discrete nature of the Poisson distribution, we employ random weighting techniques to improve the conservative coverage performance of the IM CIs. Simulation studies show that the coverage probabilities of the improved nonrandomized inferential model (NIM) CIs are closest to the prespecified nominal levels, with corresponding expected lengths shorter than those of Bayesian intervals in weak signal scenarios, whereas the shorter expected lengths of Bayesian intervals in strong signal scenarios come at the cost of sacrificing coverage guarantees. Therefore, the proposed IM and NIM CIs are superior to the Bayesian CIs. Finally, the advantages of the proposed methods are confirmed through an analysis of two experimental datasets on neutrinos in high-energy physics.

AI-based weather forecasting now rivals traditional physics-based ensembles, but state-of-the-art (SOTA) models rely on specialized architectures and massive computational budgets, creating a high barrier to entry. We demonstrate that such complexity is unnecessary for frontier performance. We introduce \ours, a probabilistic forecaster built on a standard U-Net backbone trained with a simple recipe: deterministic pre-training on Mean Absolute Error followed by short probabilistic fine-tuning on the Continuous Ranked Probability Score (CRPS) using Monte Carlo Dropout for stochasticity. As a result, our model matches or exceeds the probabilistic skill of GenCast and IFS ENS at $1.5^\circ$ resolution while reducing training compute by over $10\times$ compared to leading CRPS-based models and inference latency by over $10\times$ compared to diffusion-based models. U-Cast trains in under 12 H200 GPU-days and generates a 15-day ensemble forecast in 3 seconds. These results suggest that scalable, general-purpose architectures paired with efficient training curricula can match complex domain-specific designs at a fraction of the cost, opening the training of frontier probabilistic weather models to the broader community.

Although distributed lag non-linear models (DLNMs) are commonly used to quantify delayed and non-linear exposure-response relationships, most existing applications assume that these relationships are constant across space. However, in many geographical and environmental studies, local characteristics vary substantially across areas, making a spatially varying effect more realistic. Extending DLNMs to allow for spatial heterogeneity remains challenging, and only a limited number of modelling strategies have been proposed in literature. The most popular extension is a two-stage meta-analysis approach, which requires sufficiently large sample sizes at each location. Therefore, its usefulness is limited when working with sparse count data in small area data analyses. Although a number of alternative one-stage approaches have been introduced, their computational burden restricts their applicability in real-life data applications. In this paper, we introduce a computationally efficient Bayesian one-stage spatially-varying DLNM for count data. We define four model variants, differing in the assumed spatial dependence structure and the flexibility of the DLNM spline specification. To address the computational burden typically associated with these flexible models, we use Laplace approximations, offering an efficient alternative to classically used Markov Chain Monte Carlo (MCMC) approaches. Model comparison criteria are provided to facilitate the selection of a suitable model in a real-life data application. The proposed methods are evaluated through simulation studies, and their practical usefulness is illustrated through a real-life data application, investigating the temperature-mortality relationship in every municipality of Sicily, Italy.

Direct effect analyses usually require deciding whether a focal variable is a pre-exposure confounder or a post-exposure mediator. In observational studies, that distinction may be unclear because timing is measured coarsely or the variable reflects an evolving process. Considering the average treatment effect (ATE) and the natural direct effect (NDE) as a common notion of the direct effect when the focal variable is a confounder and a mediator, respectively, we show that, in general, no single observed-data estimand recovers both the ATE when the focal variable is a confounder and the NDE when it is a mediator. Consequently, if a practitioner applies an NDE estimator when the variable is actually pre-exposure, the resulting estimate may have no clear causal interpretation. We identify a no-additive-interaction condition under which these quantities coincide, develop sensitivity bounds for departures from that condition, and propose an alternative model-robust estimand. This estimand equals the ATE when the variable is pre-exposure and an interventional direct effect when it is post-exposure. Moreover, within a natural class of outcome-free stochastic direct effects, it is the unique observed-data functional that remains causally interpretable under both structural roles of the focal variable. We derive an efficient influence function and a doubly robust estimator, yielding robustness at two levels: the estimand is model-robust across the two causal scenarios, and the estimator is doubly robust with respect to nuisance estimation. In simulations and in an NHANES application on elevated PFAS burden, kidney function, and uric acid, mediation-based analyses yielded materially different reported estimates.

This paper introduces a projected functional gradient descent algorithm (P-FGD) for training nonparametric additive quantile regression models in online settings. This algorithm extends the functional stochastic gradient descent framework to the pinball loss. An advantage of P-FGD is that it does not need to store historical data while maintaining $O(J_t\ln J_t)$ computational complexity per step where $J_t$ denotes the number of basis functions. Besides, we only need $O(J_t)$ computational time for quantile function prediction at time $t$. These properties show that P-FGD is much better than the commonly used RKHS in online learning. By leveraging a novel Hilbert space projection identity, we also prove that the proposed online quantile function estimator (P-FGD) achieves the minimax optimal consistency rate $O(t^{-\frac{2s}{2s+1}})$ where $t$ is the current time and $s$ denotes the smoothness degree of the quantile function. Extensions to mini-batch learning are also established.

Causal discovery from observational data remains a fundamental challenge in machine learning and statistics, particularly when variables represent inherently positive quantities such as gene expression levels, asset prices, company revenues, or population counts, which often follow multiplicative rather than additive dynamics. We propose the Hybrid Moment-Ratio Scoring (H-MRS) algorithm, a novel method for learning directed acyclic graphs (DAGs) from positive-valued data by combining moment-based scoring with log-scale regression. The key idea is that for positive-valued variables, the moment ratio $\frac{\mathbb{E}[X_j^2]}{\mathbb{E}[(\mathbb{E}[X_j \mid S])^2]}$ provides an effective criterion for causal ordering, where $S$ denotes candidate parent sets. H-MRS integrates log-scale Ridge regression for moment-ratio estimation with a greedy ordering procedure based on raw-scale moment ratios, followed by Elastic Net-based parent selection to recover the final DAG structure. Experiments on synthetic log-linear data demonstrate competitive precision and recall. The proposed method is computationally efficient and naturally respects positivity constraints, making it suitable for applications in genomics and economics. These results suggest that combining log-scale modeling with raw-scale moment ratios provides a practical framework for causal discovery in positive-valued domains.

Randomized experiments have long been the gold standard for scientists seeking to learn about cause and effect. When randomized experiments are infeasible, scientists often resort to observational studies, which are widely available and often large but rely on untestable assumptions that, when violated, may result in biased estimates. Uncertainty about bias leads to a phenomenon known as the illusion of learning from observational research (Gerber, Green and Kaplan, 2004a): absent prior information about bias, observational results cannot meaningfully contribute to the estimation of a causal parameter. To shatter the illusion, we take an empirical Bayes perspective. We show that the distribution of observational biases can be learned from calibration studies-experiments that target a causal effect that is known a priori to be zero. Calibration identifies the distribution of observational bias and allows observational studies to inform the estimation of causal parameters via empirical Bayes shrinkage. We formalize the illusion phenomenon in an empirical Bayes setting and show that, with an increasing number of calibration and observation studies, both the bias distribution and the causal effect can be consistently recovered. We illustrate our method through a simulation study and a semi-synthetic application based on Ferraro and Miranda (2013)'s water-usage experiment.