Papers

Cinelli and Hazlett (2020) develops a sensitivity analysis method for the linear regression model that parameterizes omitted variable bias in terms of two partial $R^2$ parameters capturing the residual variation explained by an omitted confounder in the treatment and outcome respectively. This method is often applied to regressions fit to unit-level data when treatment is assigned at a higher level of aggregation – as in clustered observational designs. This paper shows that despite the numerical equivalence of the unit-level regression and an appropriately weighted cluster-aggregated regression for estimating the treatment effect, the sensitivity analysis procedure yields different conclusions depending on the chosen level of analysis. The outcome-confounder partial $R^2$ reflects both between- and within- group variation but the latter is irrelevant to omitted variable bias as it by construction cannot be explained by a group-level confounder. Straightforward corrections to the robustness value and the extreme scenario analysis from the unit-level regression using Pearson's partial-$\eta$ recover equivalence between these two approaches. The paper concludes with a point of caution when benchmarking against unit-level covariates and recommends always including cluster-level averages of these covariates as regressors (Mundlak, 1978).

Competing events are common in clinical and epidemiologic studies, including semi-competing risks in which a terminal event such as death may follow a nonfatal event but also competes with it beforehand. Standard summaries include the cumulative incidence function (CIF) and the incidence rate (IR), defined as the number of observed events divided by observed event-free person-time. With competing events, the naive IR generally depends on the censoring-time distribution unless intensities are constant. We propose the Average Cause-Specific Hazard (ACSH), a survival-weighted rate per event-free person-time that preserves the interpretation of an incidence rate and is defined purely from the event-time distribution, without involving the censoring-time distribution. We develop nonparametric estimation and inference for ACSH and, for two-sample comparisons, introduce ACSH differences and ratios that provide interpretable contrasts without requiring a strong model assumption between two groups. Simulation studies examine the finite-sample performance, and an analysis of the CANVAS trial illustrates the proposed methods.

Circular (directional) data arise whenever observations are measured as angles on the unit circle, such as wind direction, time of day, or calendar phase, and require statistical methods that respect the periodicity of the domain $[0; 2\pi)$. While circular-linear and linear-circular regression models are well established, regression models for a binary or binomial response observed jointly with a circular predictor remain largely undeveloped, with the sole closely related study restricted to the symmetric logit link. This paper develops and evaluates a circular logistic regression framework in which the linear predictor is expressed through the cosine and sine of the circular covariate, and compares the performance of symmetric link functions (logit, probit) against asymmetric alternatives (complementary log-log, Cauchit, and a skew-logit power link) under a generalized linear model formulation. A Monte Carlo simulation generates circular predictors from the von Mises distribution under two concentration regimes and evaluates model fit using the Akaike Information Criterion (AIC) and deviance. The methodology is illustrated with two real data sets: daily rainfall occurrence and wind direction recorded in Macomb, Illinois, and monthly earthquake counts in Western Anatolia, Turkiye, the latter used to connect the binary circular model to the related circular Poisson regression framework for count outcomes. Results indicate that the choice of link function matters most when the circular predictor is broadly dispersed and the response is markedly unbalanced; under high concentration of the predictor, symmetric links are preferred and asymmetric links are prone to instability. Practical guidelines and directions for future software development are discussed.

We consider the problem of parameter inference for ordinary differential equation (ODE) models with structural non-identifiability. Such models arise in a wide range of scientific fields, including control theory, systems biology, and public health. Structural non-identifiability occurs when distinct parameter values provide identical model outputs, resulting in lower-dimensional manifolds of observationally equivalent solutions in the parameter space. This poses challenges for Bayesian inference and Markov chain Monte Carlo (MCMC) methods, often leading to poor mixing and slow convergence. We develop two MCMC methods that use information from structural identifiability analysis. The first, Identifiability-Aware Geometric MCMC, constructs proposals that move within and between non-identifiable manifolds. The second, Identifiability-Aware Pseudo-Marginal MCMC, performs inference on the space of identifiable parameter combinations and reconstructs full parameter values. We show that both methods target the correct posterior distribution and are ergodic under standard conditions. Numerical examples demonstrate improved sampling efficiency and convergence compared with standard MCMC methods.

Estimating an unknown cumulative distribution function (cdf) from data, either as a statistical object of interest or as an input to a downstream optimization problem, is fundamental in operations. In practice, however, distribution estimation is often complicated by incomplete knowledge of the distribution's structure and limited, censored data. To address the first complication, we study distributions satisfying failure-rate shape constraints, especially increasing failure rate (IFR), rather than assuming a fully specified parametric family. To address the second, we consider noisy quantile data: at finitely many prespecified knots, each observation records only whether an independent sample lies below or above the knot. This combination arises naturally in pricing, reliability, and healthcare applications. We formulate the IFR-constrained maximum likelihood estimator and show that the original problem is infinite-dimensional and non-convex. We then develop a tractable two-step approach that solves a finite-dimensional convex optimization problem over transformed knot values and reconstructs a full cdf through shape-preserving interpolation. We establish finite-sample error bounds and convergence rates, yielding practical guidance for offline data collection. We also extend the framework to failure-rate-average, new-better-than-used, and generalized-failure-rate properties. Numerical experiments and case studies in revenue management and reliability demonstrate strong goodness-of-fit and improved downstream decision quality.

We consider statistical testing on the basis of incomplete observations with values in a separable Hilbert space, where the dimension is possibly large or even infinite. The general Hilbert space setting allows various data types as they arise in modern applications, in particular high dimensional and functional data. Possible Hilbert space testing problems are goodness-of-fit, symmetry, homogeneity and independence. We present an approach for modeling incomplete data that covers several problems in practice, e.g., ultra high dimensional random vectors with missing entries or partially observed stochastic processes. We identify a specific structure (independent and identically distributed) in the incomplete data that enables the analysis of statistical procedures with the help of suitable mathematical results (e.g., laws of large numbers and central limit theorems). Additionally, a general and novel concept for testing different hypotheses in this situation is suggested and sketched for the example of testing goodness-of-fit for normality.

In applications such as neuroscience and environmental science, data are naturally modeled as multivariate functional data and often exhibit overlapping cluster structure. Existing clustering methods for functional data typically impose mutually exclusive memberships and therefore fail to capture such structure. We propose a latent factor model based approach with functional factors and a real-valued loading matrix that encodes potentially overlapping cluster memberships. Under mild conditions, we establish identifiability of the loading matrix up to permutation, ensuring that the overlapping cluster structure is recoverable up to label switching. We develop a procedure for estimating both the number of clusters and the associated cluster memberships. This involves solving an infinite-dimensional regression problem in operators, whose solution is characterized using the inner product on the space of Hilbert-Schmidt operators and expressed in terms of real-valued matrices. This formulation enables rigorous asymptotic analysis, and we establish a central limit theorem to facilitate statistical inference on overlapping cluster memberships. We demonstrate the performance of our method using numerical studies and an application to functional magnetic resonance imaging data.

Passive surveillance systems, in which data routinely collected by medical facilities are used to monitor the caseload of infectious diseases, are relatively straightforward to implement but often result in underestimation of the burden of disease due to under-diagnosis and imperfect testing. Targeted active surveillance can be used to correct these case counts to better reflect the true burden of disease. However, when the active surveillance effort is performed at a subset of hospitals and passive surveillance data is reported at an aggregated regional level, the resulting spatial misalignment must be reconciled to estimate the true rate of hospital-presenting disease at the spatial region level. Motivated by a recent active surveillance project for leptospirosis in four Puerto Rican hospitals, we address this challenge and develop a novel Bayesian spatio-temporal framework to better reflect the true number of hospital-presenting individuals with the disease. In particular, our method extends the Poisson-logistic framework to incorporate spatial heterogeneity in the probability of presenting to the hospitals across the study region. Our framework also accounts for imperfect diagnostic testing within the active surveillance data, addressing a common challenge for infectious diseases, particularly for neglected ones like leptospirosis. The model is assessed via simulation under various scenarios and then applied to the motivating leptospirosis data. Our approach offers a comprehensive framework for integrating spatially misaligned passive and active surveillance data, enabling better estimation of true disease burden.

In some situations, data is collected under systematical and technical constraints due to uncertainty in experimental reports, intermittent measurements, confidentiality, and non-detects. For this reason, it might not be possible to retrieve or receive the data in a conventional format but rather in a grouped form where only the number of occurrences is known within intervals. The challenge is to estimate the density of the underlying ungrouped data based on the observed grouped data with no information regarding the underlying distribution. To overcome this problem, this study introduces a mean-adjusted log-concave (MALC) density estimation method for univariate grouped data, aiming to provide a bandwidth-free non-parametric approach that does not rely on specific distributional assumptions. The performance of the MALC method is evaluated through simulations across various distributions with different sample sizes and grid widths. The results demonstrate the robustness and effectiveness of the MALC approach in grouped data analysis, offering a broader range of applications over traditional methods.

We develop a general framework for constructing and analyzing overlapping sliding-window statistics for dependence and trend detection. For a fixed window size, the overlapping blocks form a Markov chain, and the asymptotic variance of any centered window statistic is determined by the covariance operator of this chain. Using its spectral structure, we obtain an orthogonal decomposition of the space of window functions into components associated with different overlap levels. This leads to a natural notion of incremental dependence information: the part of a statistic that captures exactly the new information introduced by enlarging the window. We give an explicit procedure for extracting these components and apply the method to two classes of examples. For correlation detection, we study symmetric polynomial window functions and identify their informative projected part. For trend detection, we analyze localized rank-based statistics and isolate the contribution of the largest newly introduced lag. The examples also show that different statistics may exhibit different local detection scales, including nonclassical ones. The same viewpoint leads to a natural quantitative measure of incremental information, which can be used to assess how much new dependence structure is captured as the window size increases, and to guide scale selection. Overall, the paper provides a systematic method for designing overlapping-window tests and deriving their asymptotic normalization and local behavior.

We identify the minimax exponent for constant conditional variance estimation under rough random design. The unknown design density is bounded above and away from zero, with no smoothness assumption, and the conditional error laws may depend on the covariates and have uniformly bounded fourth moments. For an $s$-Hölder regression function with $s>1$ in dimension $d>4s$, the minimax root-mean-square risk lies between $cn^{-2(s+1)/(d+4)}e^{-C\sqrt{\log n}}$ and $Cn^{-2(s+1)/(d+4)}$. These bounds show that the rate proposed by Robins is not uniformly attainable over this model class. For $0<s\le1$, we show the exact minimax rate is $n^{-1/2}\vee n^{-4s/(d+4s)}$; for $s>1$ and $d\le4s$, it is $n^{-1/2}$.

We propose a computationally efficient simulation-based estimation method with a two-step procedure for inhomogeneous bivariate Log-Gaussian Cox Processes. It combines classical Poisson estimation for the first-order parameters with simulation-based inference using neural networks for the latent field parameters. By separating the estimations, it reduces the complexity of high dimensional parameter estimation and the need for the simulation-based method to specify broad parameter ranges in the presence of covariates. In addition, we introduce two dimensional image inputs that enable the model to learn spatial information directly. Simulation results demonstrate that the proposed approach provides accurate estimates of the latent field parameters. We further illustrate the method's practical applicability using the gorilla dataset.

Extremal dependence is naturally described by the angular law of large multivariate observations. We introduce anchored geodesic component analysis (AGCA), a dimension-reduction method for extremal angular laws on the positive unit sphere. AGCA approximates angular variation by great subspheres constrained to pass through a chosen reference direction, with balanced complete dependence as the default anchor. Under a bounded sine-squared geodesic loss, the population and empirical problems reduce exactly to eigenanalysis of a second-moment matrix of anchored tangent departures. The resulting scores, loadings, residual risks and explained-variation summaries describe departures from the benchmark and remain well defined for face and near-axis extremes. Low-rank AGCA reconstructions also support tail simulation: bounded Lipschitz functionals and homogeneous tail scores, including portfolio capped excesses and value-at-risk, inherit explicit error bounds from the AGCA residual risk. We establish top-\(k\) consistency for oracle and rank-Pareto AGCA summaries and an oracle central limit theorem whose covariance is that of an independent sample from the limiting angular law. In daily equity-portfolio losses, AGCA finds concentrated benchmark-relative tail directions: ten components explain about \(91\%\) of anchored variation and approximate capped-excess and normalized value-at-risk summaries with about \(1.25\%\) average relative error.

We study online control of a known linear dynamical system with adversarial costs and bounded disturbances, measuring regret against a general class of benchmark policies. We introduce counterfactual tracking, which separates the challenge of learning from the challenge of controlling the system. An online learner builds a reference trajectory by selecting or averaging the trajectories that the benchmark policies would have generated under the realized costs and disturbances, and a corrective law steers the system toward that reference. Charging each change in the reference its recovery cost (the cost of steering the system onto the new reference) reduces the problem to online learning with switching costs. Conversely, under additional natural assumptions, we show that this reduction is tight: the two problems have the same minimax regret up to a system-dependent factor, uniformly over horizons and policy classes. The reduction gives sharp regret guarantees for policy classes beyond standard finite-memory parameterizations. For a class of $N$ possibly nonlinear or history-dependent policies, it achieves $O(\sqrt{T\log N})$ regret over $T$ rounds, provided their trajectories remain within a bounded distance of one another and recovery costs are bounded. For the full $\ell_1$ ball of disturbance-response controllers, it achieves $O(\sqrt{T\log T})$ regret, which is minimax optimal in $T$, without assuming a common decay rate for disturbance effects. The framework also improves the best known regret bounds for linear state-feedback policies.

Factorial trials can be conducted when statistical interactions between two or more treatment factors are not anticipated, but also when they are. A k-in-1 factorial trial answers k single-factor questions with the same number of units as one parallel-group trial if there are no interactions. Literature on k-in-1 factorial trials has dominated trialists understanding of factorial trials in the UK. However, factorial experiments originated from the Design of Experiments field to enable interactions to be robustly estimated. Seemingly conflicting guidance from these literatures poses a source of confusion and misunderstanding for trialists. We bring these literatures together to provide clarity on the arguments that have been used to recommend use of factorial trials in the presence and absence of interactions. We outline motivating examples. We summarise the rationales for using factorial trials, the treatment contrasts of interest, and the properties of their estimators, for the two schools of thought. We describe the debate, going back to 1935, and use an empirical example to illustrate the impact of different analysis approaches. We conclude that it is vital that trialists carefully and clearly specify their objectives. Estimands of interest and treatment contrasts follow, with properties of estimators dictated by this choice.

Data assimilation estimates the state of a dynamical system from model forecasts and incoming observations. Many observation mechanisms, however, are many-to-one, implicit, non-smooth, or accessible only through simulation, and need not provide the residual structures or likelihood guidance required by existing ensemble filters. We introduce implicit data assimilation, in which the analysis law is defined as an energy tilt of the forecast distribution. We then propose the Ensemble Controlled-flow Filter (EnCF), which realizes this update through a stochastic controlled flow and learns the observation-dependent control by adjoint matching from terminal energy gradients. For simulator-defined observations, EnCF-LF learns a surrogate conditional energy from samples and applies the same controlled-flow solver. We prove ideal exactness, derive a one-step error decomposition, and establish non-accumulation of local errors under filter stability. Numerical results show that Kalman-type filters remain preferable for smooth additive-Gaussian observations, while the proposed methods are better suited to non-Gaussian, many-to-one, multimodal, and implicit observation models.

Stochastic convex optimization is a classical problem with well-understood guarantees under first-order feedback. In contrast, for zero-order optimization with noisy function evaluations, a logarithmic gap has persisted between known upper bounds and the $\Omega(1/\sqrt{T})$ lower bound, even in the one-dimensional case. In this work, we study the problem of minimizing a convex function $f : [0,1] \to [0,1]$ using a zero-order oracle with subGaussian noise. We propose a computationally efficient algorithm that achieves the optimal $O(1/\sqrt{T})$ convergence rate, matching the lower bound. The result closes the existing gap in one dimension, providing the first sharp rate guarantee in this setting.

Stochastic-process models are, as a rule, far easier to simulate than to condition. Non-linear observations, non-Gaussian likelihoods, black-box information, and global constraints all induce intractable conditional laws, requiring bespoke, model-specific constructions. We introduce LatentFlow, a single framework for conditioning stochastic processes, with no learned neural approximations and no training. Our starting point is to write the stochastic process as the deterministic image of a tractable latent innovation, $f_0 = T_{\vartheta}(\xi_0)$, with $\xi_0$ sampled from a simple reference distribution. This reduces process-level conditioning to latent-space inference: pull the likelihood back through $T_{\vartheta}$, sample the resulting latent law with a tractable guided probability flow, and push the samples forward. This construction is provably exact at the level of the target law; in practice, approximation enters only through finite terminal noising, Monte Carlo guidance, and time discretisation of the continuous-time dynamics, each of which is explicit and systematically reducible. As LatentFlow is training-free, conditioning reduces to solving a single reverse-time SDE. This enables conditional sampling in seconds on a single desktop CPU across model classes that have never shared a scalable method: classical spatial priors, nonlinear stochastic dynamics, mechanistic models from the physical and life sciences, stochastic PDEs, heavy-tails and extremes, point and discrete-state processes, and neural or simulator-defined processes.

We show the Randomized Hamiltonian Monte Carlo (RHMC) algorithm has accelerated mixing time guarantees for sampling from log-concave probability distributions. RHMC proceeds by repeatedly simulating the continuous-time Hamiltonian dynamics for some random integration times, and resetting the velocity to be an independent Gaussian random variable between each simulation. We show that when the target distribution is log-concave and satisfies an $\alpha$-Talagrand inequality (for example, if the target distribution is $\alpha$-strongly log-concave), if we use a random integration time from either the triangular or the exponential distribution with mean $\Theta(\alpha^{-1/2})$, then RHMC converges exponentially fast in KL divergence, and the total integration time to reach error $\varepsilon$ in KL divergence scales as $O(\alpha^{-1/2} \log(\varepsilon^{-1}))$. We also show that when the target distribution is log-concave, if we use a sequence of random integration times from the triangular distribution with exponentially increasing means, then the total integration time to reach error $\varepsilon$ in KL divergence scales as $O(\varepsilon^{-1/2})$. Our analysis relies on a bound on the average KL divergence along Hamiltonian dynamics, which is inspired by an analogous result on accelerated optimization methods based on Hamiltonian dynamics.

Interference, under which a unit's outcome is affected by the treatment of other units through network connections, is often present when units interact on a network. When the network of interactions is measured, researchers are often interested in the spillover effect from first-order neighbors. When this is the case, the prevailing approach often involves the neighborhood interference assumption, which is oftentimes overly restrictive. In this paper, we instead rely on a generalized interference assumption, which allows one's potential outcomes to be influenced by the treatment of units from a wider area of the network, referred to as the "interference set". For instance, this can be a community detected through a community detection algorithm, or the set of units that can be reached through a finite network path. Under this assumption, we define new causal estimands to quantify spillover effects from first-order neighbors and, in general, from units at a specific network distance h. We employ two hypothetical Bernoulli distributions with different probabilities for the h-order neighborhood and for the rest of the units in the interference set. We first derive the bias of an approach that relies on a wrong interference set or incorrect exposure mapping function. We then develop new Horvitz-Thompson and Hajek estimators and corresponding weighted regression estimators under the generalized interference assumption. We conduct a series of simulations to assess the bias of OLS estimators – which rely on restrictive interference assumptions and an exposure mapping function – , and the performance of our estimators in different interference scenarios and random graphs. We then apply our estimators to a two-stage randomized trial implemented in Honduras to assess a maternal and child health intervention.

Nonparametric permutation tests are widely used for statistical analysis. However, exact computation of test p-values can be algorithmically challenging, particularly for custom tests with complex test statistics. In contrast, Monte Carlo sampling can be easily applied to any test statistic, but it suffers from poor relative accuracy when estimating small p-values, interfering with multiple hypothesis testing correction and leading to other issues. In this work, we present a hash-augmented adaptive multilevel splitting Monte Carlo algorithm that enables accurate estimation of arbitrarily small p-values in two-sample permutation tests. Using the Kolmogorov-Smirnov and the Mann-Whitney U tests as examples, we highlight potential pitfalls related to the discreteness of the test statistic distribution and show how to address them. By comparing with an exact algorithm, we demonstrate the accuracy of the p-value estimates provided by the proposed algorithm and the validity of the associated confidence intervals. We provide a reference implementation of the proposed algorithm in the Python package hamstest, which allows p-value estimation for a user-defined statistic.

We determine the exact region $\Omega_{\tau,\phi,\beta}:=\{(\tau(C),\phi(C),\beta(C)):C\in\mathcal{C}\}$ of possible joint values of Kendall's tau, Spearman's footrule and Blomqvist's beta over the class $\mathcal{C}$ of all bivariate copulas. The region consists precisely of all triples $(t,p,b)$ satisfying $-1\le b\le 1$, $\frac{3}{16}(1+b)^2-\frac12\le p\le 1-\frac38(1-b)^2$ and $\frac43 p-\frac13\le t\le \frac23 p+\frac13$. In other words, the known exact $(\phi,\beta)$- and $(\tau,\phi)$-regions already characterize the joint region, so that, once the value of Spearman's footrule is fixed, Blomqvist's beta imposes no additional sharp restriction on the possible values of Kendall's tau. The proof is constructive: two one-parameter families of shuffles of $M$ realize the extreme values of Kendall's tau along the lower boundary of the $(\phi,\beta)$-region, ordinal sums spread these families through the whole region, and the vertical fibres are filled using the biaffinity of the concordance function. We further show that $\Omega_{\tau,\phi,\beta}$ is convex with rectangular fixed-footrule sections, identify an affine symmetry of its fibres about $\tau=\phi$, and compute its volume, which equals $\frac{31}{40}$.

Circular data, representing angles or directions, are frequently encountered in computer vision, biology, geology, and meteorology. Traditional regression targets the conditional mean, which is often geometrically misleading for circular responses under multimodal, skewed, or asymmetric data structures. To address these limitations, a lightweight deep generative framework, namely ANGLE, is introduced for non-parametric distributional regression on the circle. The full conditional distribution of an angular response, given Euclidean and circular covariates, is learned through a generative map optimized via a generalized circular energy score (GCES) loss. Desirable theoretical properties, including the strict propriety of the loss and the rotational equivariance of the estimators, are established. Furthermore, both pre- and post-additive noise models are accommodated. A unified toolbox is provided for advancing previously underexplored challenges in circular statistics: extrapolation, sufficient dimension reduction, and conditional distribution equality testing. The framework's efficacy is demonstrated through extensive simulations and real-world applications. Specifically, the proposal is utilized for object pose estimation from imagery and wind direction prediction, which are integral to surveillance, autonomous vehicles, and energy systems, respectively. Superior predictive performance and robust uncertainty quantification of the proposed method in these tasks are revealed.

We study the squared $2$-Wasserstein distance to the standard Gaussian as a non-Gaussianity criterion and use it for linear Independent Component Analysis (ICA) and causal discovery in Linear Non-Gaussian Acyclic Models (LiNGAM). Unlike commonly used parametric contrasts and approximations of information-theoretic quantities, this criterion requires no distributional regularity beyond finite second moments, involves neither approximation nor tuning parameters, and can be computed exactly and efficiently from empirical order statistics. Our analysis relies on a strict subadditivity property of the $2$-Wasserstein distance to the Gaussian. At the population level, we prove exact identification of the ICA unmixing matrix, up to signed permutation, and give an analogous characterization of causal orders through sequential least-squares residuals. We then define empirical plug-in estimators and prove distribution-free uniform convergence under finite-moment assumptions, before detailing three practical solvers: a Picard-style orthogonal optimizer for ICA, an exhaustive dynamic program for causal order search, and a greedy order search variant. Empirically, we demonstrate competitive performance for both tasks and provide open-source implementations for source separation and causal discovery.

Conditional independence testing (CIT) is fundamental to modern statistical inference in areas related to causal discovery and variable selection. While marginal independence is relatively well-understood, despite multiple advances, no existing non-parametric CIT provides a unified, efficient, and statistically guaranteed solution across heterogeneous data. We introduce a graph-based test statistic comparing kernel similarities of the response within composite neighborhoods that use exact matching on discrete components and $k_n$-nearest-neighbor matching on continuous ones. The raw statistic, related to prior constructions, suffices under fully discrete conditioning. However, when at least one conditioning variable is continuous, we instead use a local-polynomial debiased variant that cancels the local smoothing bias. We rigorously establish its asymptotic null distribution across all data-type combinations. We further prove a dimension-free $n^{-1/4}$ detection threshold under local alternatives, eliminating the phase transition that affects geometric estimators in high dimensions. Finally, we develop efficient algorithms with near-quadratic complexity and analytic graph-based calibration, bypassing the cubic bottlenecks of global kernel methods.

In geostatistics, it is commonly assumed that sampling locations are selected independently of the underlying spatial process. In practice, however, this assumption is frequently violated. In fisheries, for example, sampling sites are often chosen to maximize expected catches, creating a stochastic dependence between the abundance process and the sampling design. Such preferential sampling can introduce substantial bias and compromise statistical inference. This study investigates the use of constructed covariates, based on average distances from nearest neighbours observations, that are able to mitigate preferential sampling. The inclusion of such covariate in the geostatistical model might be able to account for the stochastic dependence of sampling locations on the spatial variable. If this inclusion sufficiently captures the dependence, conventional methods of inference may be applied without resorting to more complex models. The proposed methodology is evaluated through an extensive simulation study that explores a variety of sampling scenarios and spatial configurations. Additionally, we demonstrate the practical utility of the approach using two real-world datasets: one on fishery landings provided by the Instituto Português do Mar e da Atmosfera, and another concerning lead pollution biomonitoring in Galicia. Results show that incorporating the constructed covariate can substantially reduce the impact of preferential sampling, enabling reliable inference with standard geostatistical tools. We also discuss practical challenges, limitations, and paths for future methodological development.

Measurements of shear wave velocity (Vs) with uncertainty are critical for site-specific probabilistic seismic hazard studies. However, rigorously quantifying the uncertainty in Vs over large enough areas and great enough depths remains challenging. In 2021, Vantassel and Cox (i.e., VC21) proposed a procedure for developing suites of Vs profiles from surface wave testing whose uncertainty were consistent with the experimental dispersion data's uncertainty. The VC21 procedure was a significant step forward, however, it requires a full dispersion data matrix to compute inter-wavelength phase velocity correlations. While applicable to many practical cases, VC21 could not be applied to the case where multiple surface wave arrays of different sizes were deployed at a site as a means of developing broadband dispersion data and deeper Vs profiles. In response, this work extends the VC21 procedure using two possible approaches for estimating a full dispersion data matrix. Approach 1 uses a selection of theoretical dispersion curves from an initial, traditional surface wave inversion. Approach 2 estimates the full data matrix by combining pieces of the data matrix obtained from the experimental dispersion measurements. Both approaches are evaluated using two synthetic datasets; one relatively-simple, three-layered model and one more-complex, five-layered model. Approach 1 and Approach 2 were able to reasonably estimate the true correlation matrix and recover uncertainty-consistent Vs profiles similar to the true distribution of Vs. While the uncertainty of the recovered Vs profiles were higher than is often assumed, the engineering proxies computed from those Vs profiles, namely the time averaged shear wave velocity in upper 30 m and the fundamental site period, showed substantially less uncertainty indicating the Vs profiles, while uncertain, are effective at capturing a site's ...

Synthetic data generation enables data curators to more easily share datasets that limits the potential for disclosive inferences about data subjects in confidential datasets. Generating synthetic data requires navigating numerous design choices; however, most existing open source software fails to provide common software infrastructure for making such design choices efficiently. In this paper, we introduce tidysynthesis, a meta-package for synthetic data generation that enables better interoperability between existing modeling frameworks and statistical data privacy methods. tidysynthesis allows users more flexibility to specify and iterate on synthetic data algorithms by providing a common syntax to easily create and modify synthetic data generation pipelines. We demonstrate the features and extensibility of tidysynthesis, as well as provide end-to-end examples for synthetic data generation using data from the American Community Survey

We study the long-time behavior of the Wasserstein gradient flow of the squared Maximum Mean Discrepancy (MMD) between a probability measure $\rho$ and a target measure $\mu$, where the underlying kernel is given by a Coulomb potential. For $L^\infty$ target densities $\mu$, we establish the existence of global weak solutions starting from arbitrary Borel probability measures and prove that the density $\rho_t$ belongs to $L^\infty$ for any $t>0$. We also show that the Hölder norm can grow exponentially in time. On the flat torus ${\mathbb{T}}^\mathsf{d}$, we prove a global metric PL inequality for every finite-Coulomb-energy source and nearly uniform target. For general bounded, uniformly positive targets, we prove exponential decay of the squared MMD without requiring a lower bound on the initial data, using a defective PL inequality. We also prove that the usual PL inequality may fail when the target vanishes only at one point and that, when $\mathsf{d}\ge2$, no PL constant can hold uniformly over all targets satisfying a prescribed lower bound. On ${\mathbb{R}}^\mathsf{d}$, for $\mathsf{d}\ge2$, under radial symmetry, source-support inclusion, and target-positivity assumptions, we establish a PL inequality and exponential convergence. On the unrestricted whole-space class, neither a multiplicative squared-MMD decay modulus uniform over the initial datum nor a global PL inequality can hold. Finally, in every dimension and in both spatial settings, we prove that every Lagrangian critical point coincides with the target when $(\rho-\mu)^+$ is absolutely continuous. In dimension two, the energy supplies uniform tightness. This implies that if our constructed solutions have finite energy at some positive time, then they converge to the target narrowly and strongly in negative-order Sobolev spaces.

Limited patient availability complicates sample size determination in pediatric clinical trials. Although Bayesian methods incorporating external data offer a solution, rigorously controlling the type I error rate remains difficult. Psioda and Ibrahim (2019) proposed a simulation-based framework as a practical solution. However, although their framework was designed to relax the type I error control, this relaxation fails when the external data exhibit a large treatment effect, making it difficult to design clinical trials that incorporate external data. Furthermore, restricting the support of sampling priors can cause trial outcomes to fall outside of this support, leading to lower power. Additionally, their analytic prior formulation may induce bias, and their method is not applicable to hybrid-controlled trials involving two-group comparisons. Thus, we propose modifications to both the sampling and analytic prior specifications and extend the framework to hybrid-controlled trials. We redefine the null sampling prior as a normal distribution centered at the null boundary, ensuring a Bayesian type I error evaluation. For the analytic prior, we employ a weakly informative prior for the second component of a robust mixture prior to mitigate bias under prior-data conflict. Furthermore, we extend this methodology to hybrid-controlled trials. Simulation studies and a pediatric case study of cutaneous lupus erythematosus demonstrate that our method substantially reduces the required sample size compared with both frequentist and original Bayesian methods, while maintaining the target operating characteristics and controlling estimation bias under prior-data conflict. This framework provides a reliable and efficient approach for designing clinical trials that incorporate external information.