Papers

A/B tests are standard in firm decision making. In the standard pipeline, experimental data is converted to a deployment decision by applying a t-test of the difference in means (the lift) and deploying the treatment if lift is positive and statistically significant. This common workflow answers the wrong question. We argue that firms need a decision rule for economic payoffs in the future deployment environment, not a test of equality in the experimental sample. We develop an ambiguity-averse decision framework in which each arm is evaluated by its ambiguity-penalized value over distributions close to the experimental outcome distribution. The resulting rule has a simple closed form thanks to the Donsker-Varadhan representation and it requires only the outcome data from a standard A/B test plus one interpretable parameter governing trust in the experiment. Our rule is thus no more difficult to implement than a t-test. A mean-variance approximation shows how the rule penalizes variability, while a connection to utility maximization shows it to be a certainty equivalent. We are able to perform a real-world evaluation of our proposed rule in the context of digital marketing using an archive of 552 advertising experiments from an anonymous US-based online platform. The proposed rule substantially reduces regret relative to conventional hypothesis testing. The results show that economically conservative, distribution-aware deployment rules can outperform statistical-significance rules in digital experimentation.

Dynamic and multilayer networks have been widely studied separately, but their joint analysis remains comparatively underdeveloped. Because the relational information of a dynamic multilayer network can be represented as a tensor at each time point $t$, each layer can be summarized through a set of structural statistics, yielding a matrix-valued observation and, consequently, a matrix-valued time series. To exploit this structure, we propose the use of matrix autoregressive (MAR) models, which simultaneously characterize temporal dependence across relational layers and structural statistics. We apply this framework to the ICEWS dataset, which records international interactions among countries under four relational domains and therefore naturally defines a dynamic multilayer network. The results indicate that negative verbal interactions (Verbal-) play a prominent role in the subsequent structural reconfiguration of the material-interaction layers, while mean strength exhibits the strongest temporal persistence and reciprocity the broadest cross-statistic influence. These findings illustrate the usefulness of MAR models for providing a parsimonious and interpretable characterization of temporal and cross-layer dependence in dynamic multilayer networks.

Marginal Mann-Whitney effects are widely used across various fields of research, and extensions of this estimand have been developed in many directions in statistical methodology. In this paper, we focus on an extensions for repeated measurements and factorial designs subject to randomly missing data. In a previous work by Rubarth et al. (2022a), asymptotically correct tests were developed under the assumption of deterministic missing indicators. In contrast, the approach in the present paper accounts for the stochastic nature of missing values under realistic mechanisms. Thus, the involved covariance matrix incorporates the true variability of missing data. The combination with a quasi-randomization procedure using random permutations within each data point yields asymptotically exact tests and a generally improved type-I error control. Additionally, the tests control the type-I error for finite sample sizes in the special case of exchangeable sampling distributions. Simulations across a wide range of settings demonstrate the benefits of the proposed method in small samples, also for different missingness mechanisms. A real data analysis about school children learning math illustrates several practical aspects of the tests' application.

This paper proposes a new theoretically exact Fourier recursion framework for a broad class of latent Markov models (LMMs), encompassing models widely used across a broad range of fields in economics. It can be viewed as a counterpart of the celebrated Kalman filter for non-Gaussian and nonlinear LMMs. Closed-form recursive updates of Fourier coefficients jointly deliver filtering, likelihood evaluation, and simultaneously accumulate the score and Hessian online. We introduce a unified truncated implementation that ensures uniform error control and numerical stability, preventing approximation errors from accumulating through the recursion. We establish asymptotic properties of the feasible maximum likelihood estimator for LMMs, provide a recursion-based consistent estimator for the Fisher information and discuss the possible model misspecification. These results provide the first general asymptotic theory for feasible approximate maximum likelihood estimation in LMMs. Simulations demonstrate its accuracy and stability, while an application to U.S. bankruptcy data recovers a persistent latent bankruptcy-pressure process.

We introduce Distributed Lag Neural Additive Models (DLNAMs), neural-additive analogues of Distributed Lag Non-linear Models (DLNMs) for learning nonlinear effects distributed over lags. DLNAMs replace a prespecified spline cross-basis with neural components that learn exposure–lag response surfaces, avoiding choices of basis family, dimension, and knot placement while preserving additive interpretability and familiar distributed-lag summaries. Exp-centered input layers, smooth activations, and learned subnetwork mixtures produce smooth, locally adaptive representations; pointwise uncertainty combines a conditional last-layer Laplace approximation with between-member ensemble variation. In simulations, DLNAMs generally outperformed DLNM comparators, including penalized and treed variants, in recovering known response functions, with lower bias, stronger boundary recovery, and better-calibrated cumulative intervals; gains were largest for more demanding functions. The architecture performed consistently across sample sizes, outcome families, lag horizons, and jointly fitted multi-exposure settings, retaining recovery performance as exposures were added; fit-specific changes were largely confined to optimization, and applications recovered established empirical patterns.

Multi-fidelity Earth system models provide simulations at different levels of complexity and computational cost, but exhaustive exploration of the parameter space at the highest fidelity is often prohibitively expensive. Multi-fidelity emulators can reduce this burden by combining abundant lower-fidelity simulations with limited high-fidelity evaluations. We compare four Gaussian-process-based multi-fidelity approaches: the Kennedy–O'Hagan autoregressive model (K&O), hierarchical kriging (HK), Bayesian hierarchical emulation for multi-level models (BayHEm), and multi-fidelity deep Gaussian processes (MF-DGP). We evaluate the methods using two contrasting applications: a three-fidelity tsunami simulator and a two-fidelity implementation of the Joint UK Land Environment Simulator (JULES). Performance is assessed using leave-one-out predictive accuracy, uncertainty representation, design requirements, and computational characteristics. In the tsunami application, BayHEm gives the lowest normalised root mean square error (NRMSE = 0.031) and highest SCORE (3.025), while MF-DGP performs worse than the single-fidelity baseline when only 10 high-fidelity simulations are available. In the JULES application, MF-DGP gives the lowest NRMSE (0.079) and highest SCORE (3.032), with all 30 held-out high-fidelity observations lying within their nominal 95% predictive intervals. These contrasting results show that no single multi-fidelity emulator is uniformly superior. Instead, method choice should reflect the complexity of the inter-fidelity relationship, the amount of high-fidelity information available, and the importance placed on predictive accuracy and uncertainty quantification.

Model and variable selection are important topics in Bayesian inference. In particular, the selection of different fixed, random effects and hyperparameters in hierarchical models can be difficult because of their complex structure. In this paper, we introduce the use of approximate Bayesian inference for model and variable selection for hierarchical Bayesian models. The approach is based on the application of Markov chain Monte Carlo methods on the model space. In particular, the Metropolis-Hastings algorithm is used to sample from the set of model indices. In this way, the marginal likelihood is used to compute the acceptance probability so that it is not required to estimate the parameter models directly. For this, the integrated nested Laplace approximation (INLA) is used because it provides accurate estimates of the marginal likelihood and it can also provide estimates of the posterior marginal of the model parameters. This method can be applied not only to variable selection but to a wide range of problems subject to model uncertainty. To illustrate the potential of this approach, examples on variable selection, changepoint models and log-Gaussian Cox processes are developed.

We propose a Bayesian generalized network autoregressive (BGNAR) model for multivariate time series whose component series are associated with the nodes of a known network. The proposed framework combines the parsimonious network structure of the generalized network autoregressive (GNAR) model with structured shrinkage and persistence priors adapted from Bayesian vector autoregressive (BVAR) modeling. We adapt Minnesota-type shrinkage to both own-lag and network-lag coefficients, with prior variances decreasing over temporal lags and, for network effects, neighborhood orders. A hierarchical prior on the own-lag coefficients allows information to be shared across nodes while retaining node-specific heterogeneity. We further adapt the sum-of-coefficients and dummy-initial-observation priors to the GNAR parameterization. Posterior inference is performed using a Gibbs sampler. Simulation studies show that BGNAR can use the same deliberately over-specified temporal and neighborhood structure across datasets without dataset-specific BIC order selection, while maintaining forecasting accuracy comparable to BIC-selected GNAR and outperforming the unrestricted BVAR benchmark in the settings considered. The structured prior regularizes weakly supported coefficients toward zero within this fixed model. Posterior distributions for the dynamic coefficients and posterior predictive distributions for future observations provide direct quantification of parameter and predictive uncertainty. An application to a wind-speed network demonstrates that BGNAR achieves point-forecast performance comparable to GNAR while additionally providing posterior inference on own-lag and network-lag effects and posterior predictive uncertainty.

Adaptive Conformal Inference (ACI) extends conformal prediction to non-exchangeable settings by adjusting the nominal miscoverage level online in response to recent coverage errors. When forecasts are issued with horizon $\tau$, however, the outcome needed to evaluate a prediction is observed only $\tau$ steps later, so these adaptive updates must rely on delayed feedback. We study this setting and ask how the effect of delay depends on the persistence of the residual process. First, we show that the $\tau$-delayed ACI recursion can be decomposed into $\tau$ interleaved ACI-like sequences. This representation yields a finite-sample bound on long-run empirical coverage with explicit dependence on $\tau$. We also derive an approximate marginal coverage bound that relates coverage deviation to changes in the underlying environment across the forecast horizon and to the adaptation rate $\gamma$. We then introduce the delay-to-memory ratio $r=\tau/L$, where $L$ is the time scale over which the temporal signal driving non-exchangeability decays. Simulation results show that the usefulness of this ratio depends on the form of temporal dependence: it strongly organizes performance under AR(1) dependence, is less predictive of overall performance under GARCH(1,1) and Markov switching, but more clearly characterizes when scale normalization remains useful in those settings. Abrupt mean- and variance-shift experiments further show that the preferred adaptation rate depends on the residual dynamics. Overall, the results show that the effect of forecast delay is best understood relative to the time scale over which past residual information remains relevant.

Structural MRI provides a noninvasive view of the neuroanatomical differences associated with cognitive impairment and dementia. Using data from the Alzheimer's Disease Neuroimaging Initiative, we investigate which anatomical regions exhibit widespread gray-matter differences and how these regional patterns vary across the clinical spectrum. Addressing this goal requires translating spatially dependent voxel-level evidence into regional conclusions while controlling multiplicity across anatomical regions. We propose the Structure-adaptive E-value FilTer (SEFT), which uses flexible working models to construct spatially adaptive voxel-level scores and aggregates them into regional partial-conjunction e-values. When combined with the e-value Benjamini–Hochberg (e-BH) procedure, these e-values provide finite-sample control of the set-wise false discovery rate under arbitrary interregional dependence. The ADNI analysis reveals a coherent neuroanatomical pattern: the exploratory analysis shows that differences between normal cognition and mild cognitive impairment are concentrated in medial-temporal regions, whereas the differences between mild cognitive impairment and dementia extend more broadly into temporal–limbic and posterior association regions. Both patterns largely overlap the normal-cognition–dementia benchmark, identifying a shared anatomical core across the clinical comparisons.

We develop a filter for time series, defined at each time $t$ as the minimizer of a discounted convex combination of observed and expected losses. The filter can be estimated by simulation to an arbitrary level of accuracy in $O(1)$ flops at each time point $t$ and can be run for all values $t=1,...,T$ in parallel. These methods are applied to robustly compute a preaveraged price process from the more than 1.5 million trades made on a single financial asset in a single day where the noise's variance is infinite. It yields a flat "volatility signature" plot, down to the 1 second level, so the microstructure noise no longer biases the volatility estimate. This is not true when linear methods are employed.

Species and feature sampling problems arise naturally whenever each observed unit is associated with one or more labels from a countable alphabet, and inference focuses on the unobserved portion of the distribution over such an alphabet. Within this framework, recent contributions have shifted attention from inference on the total probability mass over unseen labels to distribution-free confidence intervals for the largest unobserved label probability (i.e., the modal missing probability), thereby providing more refined information on whether the missing mass is concentrated on a few high-prevalence unseen labels or distributed across many negligible ones. Besides lacking a unified framework for species and features, these approaches employ a worst-case perspective that causes substantial information loss and produces overly-conservative intervals. We address these limitations through a unified model-based framework for inference on modal missing probabilities in both species and feature settings, which leverages a flexible Bayesian nonparametric formulation to localize uncertainty around models compatible with the observed data. This leads to sharper closed-form credible intervals that effectively exploit prior information and observed data, while preserving the theoretical frequentist properties and robustness of distribution-free intervals. Simulation studies confirm these improvements, while an organized crime application illustrates how our contribution has the potential to reshape law-enforcement decision-making in investigations.

In clinical diagnosis, sequential screening, and cascaded prediction systems, information is often acquired in stages, with intermediate predictions guiding whether to act on the available evidence or defer action to gather further information. We study sequential predictive inference in these settings, aiming to construct a sequence of prediction sets that supports timely, resource-efficient decisions while remaining reliable throughout the acquisition process. Our objective is to maximize process-level utility subject to any-stage validity, which controls the probability of miscoverage anywhere in the prediction sequence. We develop COINS, a conformal framework for constructing such sequences, and establish a universality theorem identifying the common exclusion structure underlying all nested, calibration-symmetric, any-stage valid prediction sequences. Building on COINS, we develop Vopt-COINS for data-driven optimization of prediction sequences under task-specific utilities based on prediction-set size, stopping time and acquisition costs. Branchwise and localized extensions further support individualized acquisition and calibration. Numerical studies demonstrate the effectiveness and practical benefits of the proposed framework for reliable, resource-efficient sequential prediction.

Pooling latent subpopulations can obscure relationships and yield misleading regression conclusions, including Simpson's paradox (SP). We propose a detector based on the steady-state dynamics of constant-step stochastic gradient descent (SGD). Unlike likelihood-based mixture tests and confounder-search methods, it requires neither a normal-mixture specification nor observed candidate confounders. Two linear pieces compete under a winner-take-all squared-error loss, and their normalized terminal separation forms the test statistic. Using diffusion approximations, we derive its asymptotic null distribution for general centered scalar covariates and for Gaussian multivariate covariates under symmetric noise. The distribution-dependent null center reduces to the dimension-free constant $4/\pi$ when both covariates and noise are Gaussian. We establish asymptotic size control and consistency against fixed mixture alternatives under regularity conditions. We extend the method to intercept and partial-mixture heterogeneity and study endogeneity, heteroskedasticity, and nonlinear misspecification. Simulations examine calibration, power, and robustness. Finally, a three-stage Detect–Screen–Verify toolkit separates evidence of heterogeneity from its substantive explanation. Across eight public datasets, it recovers four established SP benchmarks and identifies four cases that, to our knowledge, have not been documented previously. The detector requires neither latent-group labels nor a prespecified number of mixture components.

The Brier risk summarizes probability-prediction error but does not explain why two prediction rules differ. We develop Brier-UDAM, a framework for paired binary prediction-rule comparisons that separates the Brier-risk difference into mean-bias, dispersion, and Brier-scale alignment contributions. Although the observable contrast identity follows algebraically from Yates-type moment decompositions, the framework gives these terms a latent-risk interpretation and shows that paired contrasts can be identified from observed outcomes and paired predictions without estimating the conditional event-risk function. For two rules evaluated in the same target population, the exact identity $\Delta R = \Delta M + \Delta D + \Delta L$ holds, and the corresponding plug-in estimators reproduce the empirical paired Brier-risk difference exactly in finite samples. We derive influence-function-based inference for first-order regular contrasts and a projection-based confidence procedure for the mean-bias contrast when first-order inference degenerates. Simulations confirmed interpretable component behavior, recovery of the target contrasts, and generally satisfactory inference in regular settings; projection intervals retained coverage near degeneracy but were conservative. In a large public-data analysis of 30-day hospital readmission, models with similar held-out Brier risks showed substantially larger and opposing dispersion and alignment contributions, revealing cancellation that was hidden in the aggregate score. Brier-UDAM therefore provides a diagnostic account of paired Brier-risk differences, while calibration, discrimination, and information-based assessments continue to address their distinct performance questions.

Randomized dose-optimization trials may screen several candidate doses using binary efficacy and toxicity outcomes. A dose is inadmissible if efficacy is insufficient or toxicity is excessive, so each dose-specific null hypothesis is a union region and strong familywise error control must hold across arbitrary mixtures of inadmissible and promising doses. We study exact verification of multistage monitoring rules in which armwise decisions may be coupled through the number of active arms remaining at each interim analysis. For uncoupled rules, we derive an exact product representation and show that a complete-null boundary configuration is least favourable. Under active-count coupling, this factorization no longer holds because a promising dose can remain active and alter future boundaries for null doses. Assuming independent sampling across arms, prespecified per-dose analysis schedules, analysis timing not driven by the observed efficacy or toxicity outcomes, arm-specific monitoring statistics, and stagewise monotonicity, we establish an exact finite-boundary characterization without parametric restrictions on the within-arm joint efficacy-toxicity distribution. Each dose need only be evaluated at an efficacy-null boundary, a toxicity-null boundary, or a maximally favourable alternative. A deterministic finite-state recursion then verifies a fixed decision table without Monte Carlo error. Numerical studies confirmed the prespecified strong familywise error control and showed that active-count coupling can shift the least-favourable configuration away from the complete null while improving joint retention of multiple promising doses. A published randomized dose-ranging trial was used as a clinical illustration of how the framework could be prospectively implemented. The framework separates monitoring-rule construction from rigorous error verification.

The strictly standardized mean difference (SSMD) is an interpretable, unit-free effect size originally developed for quality control and hit selection in high-throughput screening and subsequently extended to effect-size assessment and error-rate control. Existing SSMD theory considers one variable at a time. Here, we extend SSMD to a linear combination of multiple variables and determine the coefficients that maximize the resulting SSMD. We show that this optimization can be formulated as a generalized Rayleigh quotient with a closed-form solution and identify conditions under which the optimal direction coincides with Fisher's linear discriminant. We derive method-of-moments, bias-corrected, and maximum-likelihood estimators of the maximized SSMD and construct confidence intervals based on Hotelling's T-squared statistic and the noncentral F distribution under different correlation and variance structures. We also characterize classification cutoffs for the resulting linear score and establish its relationship with the area under the receiver operating characteristic curve (AUROC). Simulations show that SSMD maximization coincides with Fisher's linear discriminant and logistic regression under equal covariance but diverges under unequal covariance and class imbalance. In paired designs, SSMD maximization uniquely incorporates cross-covariance information. Linear support vector machines are also evaluated as a modern comparator. This framework provides an interpretable approach to constructing multivariable biomarkers and signatures for high-throughput screening, cytokine profiling, metabolomics, salivary diagnostics, and continuous-monitoring applications.

We develop an expected shortfall factor model (ESFM) to estimate and price common variation in the severity of lower-tail losses in large panels of asset returns. Mean factor models describe common variation in average returns, while quantile factor models describe common movements in tail thresholds. ESFM instead captures common variation in the average severity of losses below those thresholds. The model combines observed risk exposures with latent common factors. We estimate ESFM using an orthogonalized two-step procedure under which first-stage quantile estimation error has no first-order effect on the ES coefficient estimates. We establish nonasymptotic error bounds for the ES coefficients, a finite-sample Gaussian approximation, and consistent selection of the number of latent factors. Applied to a large panel of equities, ESFM uncovers common factors that react sharply to market stress and contain information not captured by mean and quantile factors. Average returns increase across portfolios sorted on ESFM exposure; high-minus-low portfolios earn annualized returns of 8.0%–11.7% and Fama–French five-factor alphas of 10.3%–15.0%. These spreads remain positive and statistically significant after conditioning separately and jointly on mean- and quantile-factor exposures. Tail-by-tail spanning tests show that ESFM factors retain significant alphas after controlling for standard traded factors and the corresponding mean and quantile factors. Adding ESFM to these benchmark factor sets increases the maximum attainable Sharpe ratio. These findings identify common loss severity as a distinct and priced dimension of downside risk.

The Davis-Kahan theorem is a fundamental tool in spectral analysis, providing quantitative control over the distance between the eigenspaces of a symmetric matrix and its perturbation. However, when the matrix dimension is large, computing leading eigenvectors is computationally expensive, limiting the practical use of spectral methods in modern large-scale applications. This paper addresses this problem by proposing an independent Bernoulli sampling scheme and proves that the leading left singular vectors of the subsampled matrix faithfully approximate the target subspace of a low-rank symmetric matrix. Our main result is a subsampled Davis-Kahan bound that gives an explicit error bound depending directly on the sampling probability. The bound reveals the trade-off: the computational cost scales linearly with the sampling probability, while the statistical error scales as the inverse square root of the sampling probability. Our result thus extends the Davis-Kahan theorem to the subsampled setting, enabling scalable spectral analysis of large-scale symmetric matrices.

We present a test statistic for the comparison of DNA sequences under some of the most popular evolutionary processes available in the literature. Theoretical properties for the test statistic as well as its empirical performance by stochastic simulations are presented. The proposed test statistic is a generalized $U$-statistics built for tests of distributional homogeneity under the null hypothesis. We show that a dicothomous situation exists here. Under the null hypothesis, the $U$-statistics kernel is first-order degenerated, this test statistic falls in the quasi $U$-statistics class and follows an asymptotic normal law, albeit of higher order than the standard case. Under heterogeneity, the asymptotic normality is attained on the more usual first-order asymptotics. Asymptotic normality is proven for the cases: high-dimension/large sample size, high-dimension/small sample size, low-dimension/large sample size. Moreover, the case of local alternatives is discussed, and the contiguity of the test statistic for them is established. Simulation studies are performed to assess some finite-dimensional properties of the test statistic, regarding issues such as balanced/unbalanced samples, dimension and sample size.

Cluster randomized crossover (CRXO) trials are often used when individual randomization is impractical and the number of available clusters is limited. However, statistical analysis of CRXO trials is complex because of the need to account for complex correlation structures over time. It becomes especially challenging when very few clusters are used because standard modeling assumptions may lead to unstable variance estimates, poor confidence interval coverage, and inflated type I error. This study evaluates individual-level mixed-effects and fixed-effects models with and without a cluster-period random effect, cluster-period summary analysis using normal- or \(t\)-based inference, and two-period crossover-difference estimators. Using extensive simulation studies under both nested exchangeable and discrete time decay correlation structures, we compare model performance in terms of bias, root mean squared error, coverage probability, type I error, and convergence. Across scenarios, all models produced approximately unbiased treatment effect estimates, but their inferential performance differed substantially. Models that explicitly accounted for cluster-period heterogeneity generally provided the most reliable control of coverage and type I error, whereas simpler exchangeable models performed adequately only when the true correlation structure closely matched their assumptions. Cluster-period level analysis performance improved with increasing numbers of periods but was unreliable in the sparsest designs. Overall, the findings suggest that in CRXO trials with very few clusters, accurate modeling of cluster-period correlation is more important than the choice between fixed and random cluster intercepts, and that results from extremely sparse designs should be interpreted with caution.

We propose the Hypercomplex Irregular Autoregressive (\HIAR) model, a quaternion extension of the \IAR/\CIAR/\BIAR{} family for four-component observations at irregular times. Temporal dependence is represented by a real power of a quaternion parameter and estimated in state-space form through the Gaussian innovation likelihood of the Kalman filter, conditional on the adopted covariance specifications. Across 12,000 Monte Carlo fits, mean absolute bias decreased with sample size and ranged from 0.0009 to 0.0023 at $N=300$. We applied the model to 30,824 Sentinel-2 pixel series covering the Mata de Santa Genebra ARIE from 2020 to 2023. The optimizer reported successful numerical termination in 98.73% of the fits; the median $\|{\Phihat}\|$ was 0.6248, 3,872 pixels met the operational high-persistence threshold ($\|{\Phihat}\}\geq0.95$), and the median residual RMSE for band B8 was 5.1551 percentage points. The results demonstrate the computational feasibility of \HIAR{} and its ability to generate descriptors of multispectral persistence, vector dynamics, predictive error, and spatial discontinuities potentially associated with forest edges.

We consider dimensionality reduction for high-dimensional observations accompanied by a supplied partition into two or more clusters. The objective is not to construct a low-rank projection, but to retain an interpretable subset of the original coordinates that preserves the distributional information distinguishing the clusters. For each coordinate, the proposed procedure compares the cluster-specific empirical distribution functions through a several-sample Kolmogorov-Smirnov separation statistic. We formalize the resulting marginal cluster support and establish simultaneous finite-sample concentration over all coordinates, explicit bounds for false inclusions and omissions, and exact support recovery when the minimum distributional separation dominates the high-dimensional stochastic error. We also quantify the dimension inflation induced by using an unadjusted testing level and give a familywise-error-controlled version. Under a conditional sufficiency condition, sure screening preserves the full-data posterior cluster probabilities, mutual information, and Bayes risk; an additional result characterizes robustness to imperfectly estimated cluster labels. The procedure is invariant to strictly increasing coordinate transformations and can retain low-variance cluster signals that principal components may discard. We further develop average dual information, a criterion combining partition agreement after transformation with structural coverage of cluster-relevant coordinates, and derive its basic properties and consistency. Simulations illustrate the theory, the interpretability of the selected coordinates, and the distinction between cluster-directed screening and variance-directed projection.

Bayes estimation has been extensively studied and widely used in statistics, decision theory, signal processing, machine learning, and system identification. Among its variants, empirical Bayes (EB) estimation has attracted considerable attention due to its favorable estimation performance and computational tractability. However, the direct plug-in dependence of an EB estimator on hyperparameters can make it locally sensitive to hyper-parameter perturbations. This paper considers the linear regression model and focuses on the EB estimator by employing the marginal maximum likelihood hyper-parameter estimator. For conciseness, this estimator is simply referred to as the EB estimator. Given a family of EB weighting functions, a generalized Bayes estimator is constructed with the same excess mean squared error (XMSE) as the corresponding EB estimator. Here, the XMSE is a second-order asymptotic measure of the mean squared error difference between the estimator of interest and the maximum likelihood estimator. Furthermore, the EB estimator is shown to be at most firstorder sensitive to hyper-parameter perturbations, whereas the constructed Bayes estimator is at most second-order sensitive, making it locally more robust. The computational complexities of these two estimators are also analyzed. In some cases, the constructed Bayes estimator can be computationally comparable to, or more efficient than, the EB estimator. These theoretical results are further supported by numerical simulations.

This paper proposes and studies a likelihood-ratio test of the null hypothesis that weak stochastic transitivity (WST) does not hold. This is the reverse of the formulation commonly used in the literature, where WST is set as the null hypothesis. We show that, even when the number of items grows with the number of comparisons per pair, the uniform size converges to the nominal level, with the critical value determined by a chi-bar-square distribution. We further establish that the Type II error converges uniformly to zero under a sufficient signal-strength condition. Simulation studies demonstrate good finite-sample Type I error control and show that power increases with the sample size and signal strength.

Rank-based tests need no distributional assumptions to be valid, but the sample size they require does depend on the shapes of the two outcome distributions, which are unknown at the design stage. Standard practice substitutes a variance calibrated under the null. We ask how wrong that substitution can be. Writing $\Pi$ for the ratio of the true asymptotic variance of the estimated relative effect to the substituted one, we prove that under balanced allocation $\Pi \le 16\theta(1-\theta)/\{2+\theta(1-\theta)\} \le 16/9$ for all ordinal distributions in any number of categories: within the first-order asymptotic calculation the required sample size can exceed the calculated one by at most $77.8\%$, and an explicit two-point family attains the bound at every effect size. The effect-specific envelope, not the constant, is the operative quantity for a design: the ceiling falls to $1.63$ at $\theta=0.65$ and $1.37$ at $\theta=0.75$, so a trial powered for a larger effect is correspondingly less exposed. The proof is elementary. In five extracted trial comparisons and 4,711 generated alternatives the realized cost stays within a few percent, because the extremal configuration is one that shift-type treatment mechanisms do not produce. These facts are complementary, and fix how the bound should be used: it is a design sensitivity certificate, not an inflation factor. A protocol can report the conventional sample size alongside the largest requirement consistent with any ordinal configuration.

Partial missingness is common in high-dimensional data, but most existing change-point procedures are developed for fully observed sequences. We introduce gMiss, a graph-based framework for testing and localizing a change in the observed-data distribution of a partially observed high-dimensional sequence. The method treats the observed values together with the missingness indicators as the object of inference, so the target alternative is a change in the induced observed data law. It is designed for general distributional changes and requires neither sparsity nor Gaussianity. When the augmented observations are independent, the full permutation test controls type I error in finite samples. The procedure combines graph scans based on elementwise imputation and distance imputation. The two scans capture complementary graph patterns. Simulation results indicate that gMiss maintains accurate null calibration across the MCAR and MAR designs considered, remains competitive under Gaussian location alternatives, and exhibits strong power and localization performance in many non-Gaussian location and scale settings. We further illustrate the practical utility of the method through an application to genomic copy-number data, where gMiss identifies additional candidate boundaries that are visually plausible in the raw heatmap.

Large-scale multiple hypothesis testing is central to the analysis of high-throughput data, where controlling false discoveries is critical. Classical procedures typically rely on theoretical null distributions and often adjust for dependence among test statistics, but these approaches may be misleading when the empirical distribution of null statistics deviates from theoretical assumptions. Motivated by this observation, we investigate the role of the empirical cumulative distribution function (c.d.f.) of null test statistics in controlling the false discovery proportion (FDP). We first show that, under an oracle scenario where the empirical c.d.f. of the test statistics for all null hypotheses is known, FDP control can be achieved optimally regardless of the dependence structure, highlighting that explicit modeling of dependence may be unnecessary. Building on this insight, we propose an empirical c.d.f.-based FDP control (eFDP) method, implemented via a multivariate mixture model framework and a nonparametric estimation procedure for the empirical c.d.f.s, establish its asymptotic convergence, and construct an FDP control procedure that achieves asymptotic FDP control. Extensive simulations demonstrate that eFDP attains more accurate FDP control and higher power than existing approaches, particularly under strong dependence, and analysis of a high-dimensional breast cancer gene expression dataset confirms its practical utility.

Clinical prediction models are increasingly fitted by penalized logistic regression, because collinearity or many candidate predictors makes maximum likelihood unstable or impossible. Calibration is then almost always assessed by a grouped goodness-of-fit test such as the Hosmer-Lemeshow test. We show that this combination is invalid. Under ridge regression the grouped standardized residuals acquire a non-centrality induced by shrinkage, so the reference distribution used in practice is wrong, and at the penalty that most improves the fitted probabilities the test rejects correctly specified models between 92 and 100 per cent of the time. We derive the corrected law and define the shrinkage-corrected Hosmer-Lemeshow test, which subtracts an estimate of that non-centrality, restoring the maximum likelihood reference exactly to first order, and is made valid by prepivoting at a power cost we measure. We also give the attenuation law governing what any such test can detect once the linear predictor must be estimated. In glaucoma diagnosis by confocal laser tomography, where the maximum likelihood estimate does not exist, the corrected test finds the evidence for misfit weaker by more than three orders of magnitude: the fitted risks are too flat rather than mis-ordered, so the model needs recalibration rather than rebuilding.

This work presents DPQBootstrap, a probabilistic resampling method for time series based on a Dirichlet quantile bootstrap and rank-based temporal reconstruction. The method aims to generate pseudo-series that preserve the main characteristics of the original series, including variability, quantiles, extreme values, and part of its temporal dependence. A simulation study compares DPQBootstrap with the Maximum Entropy Bootstrap implemented in the meboot package across eight simulation scenarios. The results show that DPQBootstrap achieves 106 wins over 86 for meboot according to the Interval Score, indicating a strong ability to preserve the statistical properties of the original time series. In addition, an empirical application was conducted within a Bayesian ARMA forecasting framework, with and without DPQBootstrap integration. The results indicate that incorporating DPQBootstrap can improve point forecast accuracy and the quality of predictive intervals, particularly when combined with appropriately specified ARMA models.