Papers

Statistical validation of synthetic multivariate data requires assessing whether a generator preserves the joint dependence structure of the target population without merely reproducing observed records. We develop a model-agnostic framework based on full conditional distributions. For each coordinate, we normalize the conditional probability assigned to the observed value by the largest conditional probability available in the same record context; averaging this quantity yields a one-sided MAP-alignment statistic that can be estimated using a conditional model fitted on held-out real data. The mathematical contribution is twofold: under strict positivity and compatibility, the complete normalized conditional profile identifies the joint distribution, and its integrated L1 difference defines a metric on finite-state generative processes; we also establish consistency and finite-sample concentration for the corresponding empirical estimators. Because high conditional alignment alone can arise from copying or concentration on conditional modes, we pair it with nearest-real similarity as a separate record-level novelty diagnostic. We evaluate the framework on NSHAP health and aging data, influenza B genomic surveillance, and 34 General Social Survey waves. In GSS, the Large Science Model matched the original-data control in mean conditional alignment while retaining substantial novelty, indicating preservation of conditional structure without row reuse. In influenza B, a Chow-Liu generator matched the control alignment but had almost no novelty, revealing near-reproduction of observed records. The framework therefore distinguishes three statistically different failure modes: loss of dependence, record reuse, and mode concentration, and provides a principled basis for validating synthetic health, surveillance, and population data.

Detection and attribution (D&A) analyses provide a statistical framework for quantifying the contribution of external forcings to observed climate change. Optimal fingerprinting, the primary approach for D&A, is formulated as an errors-in-variables regression with a high-dimensional covariance structure. Reliable inference is challenging because covariance matrices must be estimated from a limited number of control simulations, and climate models may exhibit variability patterns that differ from those of the observed climate system. We develop a spiked optimal fingerprinting framework that exploits the spiked structure of dominant climate variability modes while providing stable covariance estimation in high-dimensional settings. The proposed method develops bias-corrected spiked covariance estimation for constructing adaptive weight matrices and incorporates variability inflation adjustments to account for model–observation differences in internal variability. We further develop a valid uncertainty quantification procedure for the scaling factor estimators and a residual consistency test for assessing model adequacy. Numerical studies demonstrate improved estimation accuracy and uncertainty quantification compared with existing practical approaches. Applied to annual mean near-surface temperature data, the proposed framework produces shorter and better-calibrated confidence intervals and leads to different detection and attribution conclusions in several regions, providing new insights into the interpretation of climate attribution results.

Reconstructing stationary space-time Gaussian processes at unobserved locations is costly when many sites share regular temporal records. We develop a spectral likelihood combining a temporal discrete Fourier transform (DFT) with a Hilbert-space Gaussian process (HSGP) representation of frequency-specific spatial covariance. We derive the exact covariance of the finite-record DFT coefficients and use a Whittle likelihood that approximates distinct frequencies as independent spatial problems. At each temporal frequency, HSGP approximates spatial covariance by evaluating the sampled spectral multiplier at retained Laplacian eigenfrequencies. For models specified by a joint spectral density whose half-spectrum lacks a convenient closed form, this construction avoids repeated Fourier inversion. The fixed spatial basis also permits cached feature projections to be reused in likelihood fitting and held-site reconstruction, with approximation accuracy depending on domain extension and basis size. In a simulation study of such a model, HSGP achieved reconstruction accuracy comparable to a high-accuracy quadrature reference while reducing mean fitting time by 69%. Additional applications to wind-field reconstruction tasks examine the effects of basis rank, temporal record length, and separability, and demonstrate accurate inference on held-out test sites when sufficiently rich bases are used. Taken together, the results indicate that computational savings are attainable when the spectral multiplier can be evaluated directly, numerical spatial inversion is costly, and an adequate basis has rank lower than the number of fitting sites.

Statistical models of ecological and evolutionary dynamics often include latent variables that are either continuous (e.g., average body size) or discrete (e.g., numerical abundance). Mixed-type hierarchical models containing both are typically fitted using Markov chain Monte Carlo (MCMC), which can be prohibitively slow for large models. Here, we introduce an alternative in the R package RTMB that automates the sequential reduction of small groups of related discrete variables, allowing them to be efficiently marginalized. Sequential reduction is combined with automatic differentiation and the Laplace approximation to estimate parameters and predict both continuous and discrete variables. We demonstrate speed and flexibility using demographic examples (occupancy, dynamic occupancy, N-mixture, and open dynamic N-mixture models), benchmarking RTMB against JAGS and unmarked. We then develop two novel applications. The first is a multi-site open N-mixture model with a spatial latent variable governing site-specific initial abundance and recruitment, which shows that continuous Gaussian Markov random fields can be estimated jointly with discrete abundance dynamics in under a minute. The second is phylogenetic trait imputation for a published data set of female Liolaemus lizards, where we jointly impute a binary trait (viviparity), estimate its state-switching rates, and estimate its effect on a continuous trait (body size) during ancestral state reconstruction. This indicates that phylogenetic comparative methods can estimate linkages among discrete and continuous traits. We envision that intuitive and efficient specification of mixed-type models will allow more expressive representation of ecological and evolutionary dynamics.

Computing matrix inverses (or, e.g., left inverses) using known algorithms is in general computationally prohibitive in high dimension, such as in an operational weather forecasting context. In this article, a swift, scalable algorithm to calculate the left inverse (through singular values decompositions techniques) of a square root of a horizontally localized ensemble covariance matrix is introduced. In this algorithm, calculations involving model variables located in different grid vertical columns can be performed separately. This is the aspect that makes this algorithm scalable by allowing for parallel computing. The algorithm is developed in the case when the number of ensemble members equals the number of model variables in each grid vertical column (this number is of order $10^3$ for current high-resolution models). The algorithm is tested in a synthetic experiment on a three-dimensional grid, and it is found to have very good accuracy. Additionally, a version of the algorithm outputting the left-multiplication of a generic vector by that above-mentioned left inverse is also introduced. For operational applications, if executed on a single processor, that version of the algorithm is already $9$ orders of magnitude faster than a traditional algorithm, and, with multiple, parallel processors, it scales until the number of processors reaches the order of $10^6$. In addition, an application of the algorithm is pointed out: performing the transform from model variables to an alternative set of variables, defined herein as standardized variables, being climatologically uncorrelated and each having unit climatological error variance. Furthermore, it is pointed out that the standardized variables allow constructing a new Hybrid covariance model, which is positive definite and climatologically unbiased, and it features a geographically dependent climatological error covariance matrix.

When providing forecasted probabilities with a predictive model, the ideal model offers perfect calibration: the true probability of the outcome (i.e., the probability that $Y=1$) exactly matches the forecasted probability $f(X)$. In practice, models inevitably exhibit calibration error, and it is therefore important to be able to measure this miscalibration to assess a model's reliability. The Expected Calibration Error (ECE) is the most widely used measure of miscalibration, but is known to be impossible to estimate the ECE with guaranteed accuracy in an assumption-free setting. In this work, we propose an alternative measure, the rankECE, that is based on comparing points with neighboring values of the predicted probability $f(X)$. Our theoretical guarantees and empirical results establish that rankECE provides a better proxy for ECE as compared to binned approximations to ECE, which are the most commonly-used approximations in practice.

We develop a unified framework for $c$-sample homogeneity testing on the circle. The proposed class of $c$-sample Sobolev tests generalizes the two-sample Sobolev tests based on uniform scores and encompasses the existing multisample tests as particular cases. We further embed this class into a broader aggregation framework, where the samplewise Sobolev components are combined through general merging functions, including average-type, maximum-type, and interpolating families. Within this class, we introduce the first Anderson-Darling-type homogeneity test for circular data, and we further propose two new tests designed to detect multimodal departures from homogeneity, constructed from softmax and Poisson kernels. We derive the asymptotic null distribution of the class, prove its consistency against a broad family of fixed alternatives, and obtain the asymptotic distribution under shift-type local alternatives. The tests are distribution-free and therefore do not require resampling. A comprehensive simulation study demonstrates the strong power of the Anderson-Darling-type test compared with Cramér-von Mises-type competitors across several scenarios, and the effectiveness of the softmax and Poisson tests under several alternatives. The testing toolbox is applied to analyze the nursing patterns of polar bears.

We consider regression problems in which the marginal distribution of the covariates is informative for estimation and variable selection, rather than merely auxiliary. Motivated by random-design, high-dimensional, and latent-effect settings, we propose a general design-assisted regression framework in which the estimating criterion depends on both the conditional model for $Y \mid \bfX$ and structured features of the covariate distribution. The framework identifies two roles of design information: stabilizing weak design directions through quadratic regularization and correcting latent-effect distortion through nuisance augmentation. We establish oracle properties for the resulting estimator, separate the effects of stochastic error, shrinkage, and approximation, and compare it with a benchmark sparse procedure that ignores design information. These results show that the proposed framework improves estimation while preserving first-order prediction performance. Numerical studies and two real-data applications illustrate the practical impact of incorporating design information.

Model choice in backward Q-learning is recursive because a later-stage choice changes the response supplied to an earlier regression and can alter its model-comparison statistic. Separate stagewise criteria do not directly assess the target-stage prediction risk of a completed Q-learning fit. We address this mismatch by treating the entire backward-fitting rule as the unit of comparison and propose feedback-aware soft tuning for sequential multiple assignment randomised trials. Each backward fit uses its own generated responses and is assessed at a common prediction target. The risk criterion retains downstream effects on upstream comparisons, while a separate correction accounts for estimating the final exponential weights from the same observations. For a fixed finite library of smooth recursive maps, we establish an exact risk identity under a Gaussian shift model and an oracle inequality with an explicit adaptation remainder. Under coordinate-representation and moment conditions, these guarantees transfer to prediction risk at any prespecified stage, with the number of stages fixed as sample size increases. A two-stage construction supplies an explicit observable implementation. Numerical studies examine risk estimation and finite-sample performance, and a simulated attention-deficit/hyperactivity-disorder trial illustrates the relation between comparison feedback and treatment recommendations.

When reconstructing the hidden state of a partial differential equation (PDE) system from partial observations, aggregate prediction error measures performance on a given data distribution but does not reveal how strongly the observations constrain each predicted direction. We introduce UniPDE-Bench, a direction-wise evaluation protocol that incorporates local observation geometry into state-estimation assessment, providing a reference for prediction recovery and confidence-based selection that is independent of the estimator. The protocol whitens the joint observation Jacobian by the noise covariance and normalizes it by a state metric. Its complete right singular basis represents joint variations of the state, which a relative sensitivity threshold partitions into retained and below-threshold directions. In this common basis, the protocol evaluates recovery and the agreement between prediction claims and the geometric partition, while recovery and abstention curves describe confidence-based selection at different claim coverages. In the simulated tasks and observation configurations studied here, confidence rules whose overall ranking exceeds chance can still exhibit below-random abstention on below-threshold directions at some high claim coverages. By separating empirical recovery from direction-selection quality, the protocol relates prediction performance to local observation sensitivity and provides an evaluation of partially observed state estimators beyond aggregate error.

Generalized Poisson matrix factorization (GPMF) is a matrix factorization method for count data with overdispersion based on the generalized Poisson distribution. While GPMF provides point estimates of the model parameters through maximum likelihood estimation, it does not quantify estimation uncertainty. In this paper, we propose a Bayesian extension of GPMF, referred to as Bayesian GPMF, and develop a Gibbs sampler for posterior inference. The proposed method is based on a compound Poisson representation of the generalized Poisson distribution, which introduces latent variables and yields closed-form full conditional posterior distributions for all model parameters. To efficiently sample the overdispersion parameter, we derive an infinite mixture representation of the exponentially tilted beta (EBeta) distribution and develop a finite approximation based on this representation. Simulation studies demonstrate that the proposed Bayesian GPMF achieves smaller mean squared errors than the conventional GPMF while providing credible intervals with empirical coverage probabilities close to the nominal level. Furthermore, the proposed finite approximation attains estimation accuracy comparable to Sampling/Importance Resampling (SIR) while requiring less computation. An application to football event data further illustrates the usefulness of Bayesian GPMF for extracting interpretable latent structures and quantifying estimation uncertainty. These results demonstrate that the proposed framework provides an effective Bayesian approach to generalized Poisson matrix factorization.

History-based risk factors, such as known family history, recorded personal history, and documented regional history, are often treated as covariates. Because these variables depend on testing, diagnosis, and recording, their observed values are influenced by detection. However, monitoring data are often unavailable, making it difficult to distinguish true history effects from detection-driven associations. We study the feedback loop that arises when observed history affects future monitoring and future monitoring affects which histories become observed. Prostate cancer screening serves as a case in point: knowing a family history may raise awareness and increase testing compared with having no known history. We develop Markov models for true and observed history states and show that differential detection can distort both observed risk ratios and the distribution of the observed history variable, and can create an apparent history effect even when true event risk does not depend on history. With history-dependent true risk, the observed history process is generally not Markovian, although the joint true-observed process is. Uniform incomplete detection can attenuate true history effects by shifting individuals with missed events into less recent observed-history states. We also develop an inverse sensitivity-analysis framework that combines a published observed association, a baseline event probability, and plausible detection probabilities to obtain the implied true risk contrast. Numerical analyses illustrate the distortions, and a prostate cancer family-history example demonstrates the sensitivity calculation using external calibration values. The framework is intended for registry, electronic health record, and surveillance studies that use documented history as a risk factor.

We study the problem of estimating unknown parameters in PDE-constrained inverse problems from noisy observations, where the PDE solution is approximated using Physics-Informed Neural Networks (PINNs). While PINNs have demonstrated remarkable empirical success, existing estimators often inherit the slow nonparametric convergence rate of the neural-network solution, leading to biased and statistically inefficient inference for the finite-dimensional parameters of interest. To address this, we propose a two-step debiased estimation procedure that combines neural-network-based nonparametric estimation with an influence-function-based bias correction. By eliminating the first-order sensitivity of the estimator to errors in the nuisance function, our procedure yields a $\sqrt{n}$-consistent and asymptotically normal estimator without requiring undersmoothing of the neural network component. We further extend this framework to Bayesian inference by replacing the original likelihood with a debiased quasi-likelihood and establish a Bernstein-von Mises theorem showing that the resulting posterior contracts at the $\sqrt{n}$-rate with an asymptotic covariance matching that of the frequentist estimator. As a by-product of our analysis, we establish near-minimax optimal convergence rates for estimating a nonparametric regression function and its derivatives in Sobolev spaces using neural networks. Extensive numerical experiments corroborate our theoretical findings and demonstrate the necessity of the proposed debiasing procedure for valid statistical inference in PDE-constrained inverse problems.

Sampling error affects the finite-sample values of residual-based fit indices. We show that the unweighted comparative fit index ($CFI_u$) rests on a quadratic-kernel maximum mean discrepancy and that a Gaussian-kernel comparative fit index (KCFI) converges to $CFI_u$ as bandwidth increases. Because kernel discrepancies lack an analytic sampling correction, we extend expected-discrepancy bootstrap methods to a fitted-correlation reference state: the sample is transformed to reproduce the fitted correlation matrix, then resampled and refitted to estimate the exact-fit sampling component. Across 18 confirmatory factor analysis conditions, this correction reduced the root mean square error (RMSE) of both indices by 62%-87%. Under normality it matched an analytic residual-covariance estimator of $CFI_u$ within .0001, and in small samples it outperformed ordinary bootstrap correction in accuracy and availability. Corrected KCFI had lower RMSE than identically corrected $CFI_u$ in all 29 nonnormal conditions, even against a common target, but in six heavy-tailed small-sample conditions the analytic estimator was more accurate than both. The results validate the reference-state correction where it is needed and delimit where the analytic estimator is preferable.

Understanding how the brain's cortical regions respond to stimuli like fingertip tapping is a significant challenge in neuroscience, especially with high-resolution functional magnetic resonance imaging data. To address this, we propose a new statistical method for evaluating how a factor's influence locally varies across a complex graph, such as the brain's cortex. Our approach is designed to handle the complexities of large graph sizes and computational demands. Our method uses a novel Bayesian hidden Markov random field to partition the factor's influence into collapsed states with similar effects on the outcome at each node. This unique model promotes sparsity in the partitions and penalizes large variations across adjacent nodes. The result is a highly detailed local influence map that captures subtle changes in the patterns and magnitudes of a factor's influence across the graph. To ensure efficient analysis of large datasets, we developed a specialized Markov chain Monte Carlo algorithm. Our simulation experiments demonstrate significant improvements over existing techniques in both accuracy and scalability. Ultimately, our method provides a powerful new framework for exploring the intricate relationship between stimuli and brain activity, offering a clearer picture of the cortical mapping of fingertips.

Identification in vector autoregressions involves two distinct choices: how extensively to orthogonalize shocks and whether orthogonality is imposed sequentially or simultaneously. Generalized identification imposes no orthogonality; Sims (1980} imposes full orthogonality sequentially; Francis et al. (2026) impose full orthogonality simultaneously; and Buchwalter et al. (2026a) provide clustered partial orthogonalization sequentially. We fill the remaining case by developing simultaneous clustered orthogonalization (SCO). SCO preserves contemporaneous dependence within economically meaningful clusters while imposing orthogonality across clusters jointly, thereby eliminating dependence on cluster ordering. We formulate the associated correlation-maximizing identification problem and show that, in correlation space, it reduces to a quadratic matrix equation. This yields a closed-form solution, and we prove that the associated identification matrix is the unique global maximizer. SCO is order- and scale-invariant and nests generalized identification and full simultaneous orthogonalization as special cases, yielding a flexible family of structural decompositions indexed by the number and composition of clusters.

Prequential adaptation of frozen tabular foundation models under temporal drift, with each label revealed only after prediction, helps some deployments and harms others, yet current practice does not predict which. We study the sources and limits of these gains. A diagnostic anatomy attributes gains to four recurring mechanisms under a streaming protocol that removes three optimistic biases and quantifies a fourth. Within an agnostic total-variation drift class, the target conditional is only partially identified: its identified-set diameter, the wall, is irreducible from unlabeled data uniformly in sample size. A second, orthogonal $L^2$ projection wall quantifies what the frozen representation cannot express. Two canonical mechanism priors collapse the first wall. Under stated nuisance-rate conditions, the wall can be estimated from labeled historical windows at a $\sqrt N$ rate above the margin threshold $\gamma^\star=d_0/(2\alpha_s)$. At $\gamma=0$, the conditional lower-bound program depends on an open affinity estimate; the positive-margin lower branch also remains open. Semi-synthetic data illustrate the finite-sample mechanism with calibrated exponents. Stream-level proxies on eight industrial streams fall on the difficult side under a stated roughness bound, while the equality case $\gamma=\gamma^\star$ remains unresolved.

A standard spatial regression model $Y = X^\top\beta + U$, with covariates $X$ and a Gaussian random field $U$ with Matérn covariance, decouples covariate effects from spatial dependence. In several applications this is mechanistically implausible, as covariates likely drive the dynamics generating the residual covariance, and it is statistically problematic, since the generalised least-squares estimator of $\beta$ may be biased by spatial self-confounding. We propose operator-matched spatial regression, where covariates enter both pointwise and as deterministic forcings in a stochastic partial differential equation (SPDE) that generates the random effect, partitioning the apparent effect into local and operator-mediated components. The formulation arises as the steady state of a transport–relaxation equation forced by $X$, giving the parameters physical interpretations. We show that self-confounding is resolved by construction and that maximum likelihood inference for the regression coefficients, implementable through standard covariance-based methods, is calibrated when the operator parameters are known and, under expanding-domain asymptotics, when they are jointly estimated. A computationally efficient finite element implementation is introduced and justified via bounds on the discretisation error. Applications to temperature–altitude regression and air-quality monitoring show that the method matches or outperforms alternatives and produces physically meaningful results.

In certain clinical settings, patients diagnosed with a disease of interest are at risk for death as well as a serious nonfatal event, and interest lies in event-free survival (EFS), a composite endpoint defined as the time from disease onset until the earlier of the nonfatal event and death. Component-wise censoring of EFS arises when each component is subject to a different censoring mechanism. For example, the nonfatal event may be interval censored between assessments, and death is right-censored. Further, in studies where individuals enroll after disease onset (including prevalent cohort studies), EFS is left truncated. Methods to estimate EFS probability with left-truncated and right-censored data are available in the literature, but they cannot handle component-wise censoring. We propose a kernel smoothing method to non-parametrically estimate EFS in this setting. Our method can also estimate and test for differences in the restricted mean event-free survival time, and can leverage two types of supplemental data that may be available: data from participants followed for death only (not followed for the nonfatal event), and incident cohort data, which arises when there is no delay between disease onset and study enrollment. We assess the proposed method using simulations and demonstrate the method using data from the Atherosclerosis Risk in Communities (ARIC) Study to estimate dementia-free survival probability following a myocardial infarction.

Generative artificial intelligence changes how firms reach customers, but standard marketing data do not record how often users see and notice a firm's name in generated answers. We develop Generative Marketing Mix Modeling (GMMM) to estimate the causal effects of Generative Engine Optimization (GEO) and Generative Engine Marketing (GEM). For GEO, GMMM combines repeated generated answers with question counts, shares of use across generative systems, and notice probabilities. For GEM, it combines records of sponsored placements with notice probabilities. GMMM compares expected business responses under alternative treatment sequences and establishes sufficient conditions for identifying the resulting effects. We investigate the empirical performance of the proposed method using simulated answers to product recommendation in English and Japanese.

Privacy and data-governance constraints often prevent pooling individual-level data across studies, limiting the use of conventional approaches for causal media- tion analysis in multicenter settings. We propose a federated causal meta-mediation framework for right-censored time-to-event outcomes that enables collaborative es- timation without sharing individual-level data. Our framework targets natural indirect effects in a prespecified population by combining information on mediator and outcome mechanisms across distributed data sources. A site-by-site identifi- cation strategy further allows heterogeneity across data sources to be character- ized, with a variance decomposition separating outcome-related, mediator-related, and interaction components. We develop federated one-step and targeted maxi- mum likelihood estimators that accommodate data-adaptive and machine-learning methods for nuisance-function estimation. The finite-sample performance of the proposed estimators is evaluated through numerical simulations. To illustrate the practical utility of the framework, we apply it on data from the French National Health Data System to evaluate the role of methotrexate coprescription in explain- ing the effect of TNFi versus IL-12/23 inhibitor therapy on treatment persistence among psoriatic patients.

High-dimensional two-sample mean testing is challenging when the dimension exceeds the sample size. The random projection method proposed by Lopes et al. (2011) addresses this difficulty by mapping the data to a lower dimension space where Hotelling's $T^2$ statistic can be applied, while retaining useful covariance information and gaining power when the variables exhibit non-negligible covariance structure. However, single projection tests may be sensitive to the realized projection matrix, whereas existing multiple projection procedures often rely on resampling or simulation for calibration, with limited theoretical understanding. Moreover, projection-based Hotelling tests may be less sensitive to sparse mean differences. To address these limitations, we propose a Cauchy-combined random projection test (CRPT), which applies Hotelling's $T^2$ test after multiple independent random projections and combines the projected $p$-values through the Cauchy transformation. The proposed method retains the ability of random projection methods to incorporate covariance information while reducing reliance on any single projection. Under the Gaussian assumption, we establish the null tail behavior of the proposed statistic and further investigate its asymptotic power under suitable alternatives. To improve sensitivity to sparse alternatives, we further develop a power-enhanced version of CRPT. Simulation studies and real data analysis are conducted to examine the performance and practical applicability of the proposed procedures.

Modeling the joint distribution of extreme values in high-dimensional financial time series is challenging because extremes are sparse and locally extreme observations are not necessarily extreme relative to their full marginal distribution. To address this, we introduce a time-dependent network Hüsler-Reiss model in which market-informed adjacency matrices determine how strongly observations contribute to the estimation. We propose binary and weighted specifications, including the Joint Extremes Adjacency Matrix (JEAM) which combines information about individual extremeness with historical patterns of joint extreme movements. In the forecasting evaluation part, covering one-minute stock returns from three sectors of the S&P 100, JEAM achieves the best out-of-sample log scores for both tail directions; improving scores by 12.5-13.6% in the lower tail and 11.4-14.9% in the upper tail. The results show that incorporating market-informed network structures in the estimation, improves forecast evaluation of extremes across time series.

Finite mixture models characterize heterogeneous populations and are increasingly fitted to distributed data using split-and-conquer procedures that aggregate local mixture estimates at a central server. The aggregation step can be seriously compromised when transmitted local mixture estimates are partially or completely corrupted. To guard against Byzantine failures, existing robust aggregation methods have been developed for settings in which a local mixture estimate is either entirely authentic or entirely corrupted. Such methods can discard useful information when only some component estimates are corrupted. We consider component-wise Byzantine failure, in which some component estimates may be corrupted, but for each mixture component, a majority of the corresponding local estimates remain authentic. We propose component-wise filtered mixture reduction (CFMR), which selects a data-driven anchor, aligns transmitted components, filters each aligned cluster by a majority radius, and aggregates the retained estimates through mixture reduction. By filtering out only unreliable components, CFMR preserves authentic information from partially corrupted machines without requiring knowledge of the failure rates. We establish an adaptive convergence bound for CFMR and show that, under suitable conditions, it attains the oracle rate that would be achieved if the authentic component estimates were known in advance. Simulations and a real-data application show that CFMR remains close to the component-level oracle, whereas whole-machine filtering and unprotected aggregation can deteriorate substantially.

Cluster-based permutation tests are widely used for analyzing MEG data, even though they are limited to cluster-level inference and do not provide spatially resolved effect estimates. We propose a regression-based framework for modeling brain activity directly on the cortical surface. Spatial effects are represented using Wendland radial basis functions, and model-based gradient boosting is employed for data-driven selection and estimation of localized activation patterns. This yields interpretable, spatially resolved effect estimates while mitigating the need for extensive multiple testing correction. In a simulation study with heterogeneous signal structures, the proposed approach recovers localized effects that are difficult to detect using cluster-based methods. An application to experimental MEG data further illustrates its ability to reveal spatially specific activation patterns. Overall, this unified framework provides a flexible and interpretable approach for modeling cortical surface data within a unified statistical setting.

Methods for areal data commonly describe spatial structure through global autocorrelation measures or local indicators based on neighboring values. In many applications, however, it is also relevant to know whether high- or low-valued areas form connected spatial structures. Motivated by topological data analysis, we develop global and local indicators of spatial connectivity for observations defined on an areal adjacency graph. We construct thresholded spatial graphs and use the Betti-0 curve to record their number of connected components. We further derive a local decomposition in which each area's contribution to the Betti-0 curve is expressed as an activation term minus a merging term. Each component fusion contributes to the merging term through a symmetric allocation between the entering area and the neighboring areas through which the fusion occurs. Integrating the merging contributions yields area-level indicators of participation in the connectivity of high- or low-valued regions. Global and local inference is based on random relabeling, with a global test for departures from spatial exchangeability and a conditional test for the local indicators. We illustrate the framework using COVID-19 incidence in Italian provinces during two epidemic waves. Significant departures toward greater connectivity were found for both high- and low-incidence areas, with the largest discrepancy corresponding to high incidence during the first wave. The proposed indicators complement conventional measures of spatial association and are applicable whenever the spatial coherence of threshold exceedances or deficits is of interest.

Graphical causal inference supplies a complete theory of efficient covariate adjustment for the average treatment effect: one adjustment set, computable from the graph, is optimal under every compatible distribution. We show that this is a property of the average treatment effect's inverse-prevalence weights, not of causal estimands in general. For the average treatment effect on the treated we index the efficiency bound by the adjustment set and derive exact identities for its change under treatment-side and outcome-side extensions of a valid set. Covariates that predict only the treated-arm outcome are exactly efficiency-neutral, and covariates that predict the control-arm outcome can strictly increase the bound when the propensity is below one half – a reversal of the supplementation lemma whose source is an arithmetic-geometric-mean inequality that holds for the average treatment effect and fails for the treated-population estimand. A construction with two faithful distributions on one graph proves that no graphical optimality criterion exists for the treated-population estimand; under no effect modification the ATE-optimal set is nonetheless optimal among the graphically valid sets, with an exact expression for its advantage. The results extend to weighted average treatment effects with propensity-dependent weights, yielding symmetric thresholds for overlap weights, an estimand-drift phenomenon under instrument adjustment, and a characterization of constant weights as the only smooth positive weights for which outcome-side supplementation never increases the bound. Simulations and the LaLonde data provide illustrations.

Despite its ubiquity, clustering lacks a universally accepted definition of what is a cluster. Kleinberg's Impossibility Theorem formalizes this difficulty by showing that no flat clustering method can simultaneously satisfy three natural axioms: scale invariance, richness, and consistency. In this paper, we ask whether this impossibility persists when the output is a hierarchy rather than a single partition. We show that, in contrast to the flat clustering setting, the hierarchical analog of these axioms are jointly satisfiable. In fact, there exist uncountably many hierarchical clustering methods satisfying these axioms, which we call admissible. We explicitly construct several admissible methods, including methods based on well-separated clusters and a non-binary version of single linkage. For certain pairs of admissible methods, the hierarchy produced by one always refines that produced by the other. This refinement relation defines a partial order on the class of admissible methods. This partially ordered set has no greatest element and contains uncountably many pairwise incompatible maximal elements, revealing substantial diversity among admissible methods. Nevertheless, this diversity is constrained: every admissible method contains a hierarchy of sufficiently well-separated clusters, and every finite collection of admissible methods shares such a nontrivial common backbone.

We propose conditional f-modeling (Cf-modeling), a framework for empirical Bayes inference with covariates. A central identity shows that the conditional marginal score function determines not only the posterior mean through Tweedie's formula, but also the posterior moment-generating function, providing a basis for recovering posterior quantities without explicit prior modeling. Motivated by this observation, we treat the conditional marginal score as the primary object of inference and estimate it directly using an energy-based representation and Hyvärinen score matching, thereby avoiding potentially intractable covariate-dependent normalizing constants. The resulting framework flexibly accommodates covariate effects and heteroscedasticity and provides a practical approach to posterior moment estimation and uncertainty quantification. We demonstrate the effectiveness of the proposed method through simulations and an RNA-seq application.

Applied work often locates a threshold by linearizing a smooth function about a reference point and solving for the crossing. When the function is curved, the linear crossing can be substantively displaced even when standard errors are valid. We introduce the curvature-overstatement parameter $\Theta_{COT}=\log(|h_2^*|/|h_1^*|)$, the log ratio of second- to first-order threshold displacement. On the quadratic branch continuous with the linear solution, $|h_2^*/h_1^*|=2/(1+\sqrt{1-u})$, where $u=2qa/b^2$ is a dimensionless index formed from the local gap, slope, and curvature. Thus $\Theta_{COT}$ is scale-free, depends on the local parameters through one scalar, and has regular-domain range $(-\infty,\log 2)$, with boundary limit $\log 2$ at tangency. We derive regular asymptotic inference, characterize local-to-tangency and weak-slope failures, and give a remainder bound linking the second-order crossing to the true threshold. The main practical contribution is an adequacy test that can affirm that the linear threshold is accurate within a prespecified proportional tolerance, rather than treating failure to detect curvature as evidence of adequacy. Monte Carlo results confirm regular-case calibration, the predicted nonstandard behavior near tangency, and failure under weak slope, while bootstrap diagnostics identify regimes in which regular inference should not be used. COT therefore provides an effect-size scale, adequacy test, and diagnostics for deciding whether a first-order threshold is accurate enough to report.