Papers

We develop a new approach to Personalized Federated Learning across heterogeneous clients using Nonparametric Empirical Bayes (NPEB). Leveraging the asymptotic normality of local parameter estimates obtained from Empirical Risk Minimization or M-estimation, our method formulates these estimates as noisy observations to estimate an unknown shared prior via Nonparametric Maximum Likelihood. A key challenge in applying NPEB in this setting is that existing approaches assume known fixed variances, which is not true in practice. To address this, we introduce a Variance-Aware Nonparametric Empirical Bayes (VANEB) framework that leverages the parameter-dependent asymptotic variance of local M-estimators. A key technical contribution is a generalized Tweedie's formula for this heteroskedastic setting. We then establish non-asymptotic error rates for density estimation in the average squared Hellinger distance and derive an oracle denoising inequality that provides error bounds for our estimator. While our theoretical guarantees are rooted in the asymptotic regime of M-estimators, we empirically explore heuristic extensions of VANEB to modern federated learning settings involving Deep Neural Networks (DNNs). For DNNs, we propose VANEB-head and VANEB-FT, which personalize the last fully connected layer via an NPEB step using an approximate diagonal variance estimator. We show that our method has strong performance on popular vision datasets MNIST and CIFAR-10, using a convolutional neural network architecture.

Sequential outcomes in longitudinal studies and multi-wave surveys may be missing not at random at both earlier and later occasions. We study graphical models in which at least one outcome is self-censoring and the response indicator for a later outcome may depend on either the earlier response indicator or the realized earlier outcome. These restrictions define two candidate families under which the relevant full-data distributions are identifiable and whose observed-data models overlap; graphs containing both dependencies form a broader class outside the prespecified comparison. For each candidate family, we establish identification of the full-data distribution under rank or completeness conditions and develop likelihood-based estimation. We then propose a two-stage Vuong-type procedure. The first stage determines whether the candidate models are observationally distinguishable; only after distinguishability is established does the second stage compare their Kullback–Leibler divergences from the true observed-data distribution. We also show that ordinary Wald inference remains asymptotically valid for the selected model-specific functional when the selected model has a fixed positive expected log-likelihood advantage. Simulations evaluate the two-stage procedure across graph classes. We finally apply the procedure to compare candidate missingness models in the Job Corps data and perform downstream functional estimation under the selected model.

Standard win statistics methods determine a win, loss, or tie for a pair of subjects based on their worst outcomes (up to the end of study) that may not fully utilize all patients' conditions or disease experience throughout the follow-up period. While the newly developed win-time statistics fully utilize all patients' longitudinal information, these statistics have been limited to time-to-event endpoints and require monotonic pattern of the events. As such, they are not applicable to any type nor number of hierarchical longitudinal endpoints. We propose the win time in favor of treatment (WINFT), a general measure for any hierarchical longitudinal endpoints, that summarizes the total time a subject in the treatment group spends in a more favorable health state than a subject in the control group. Unlike existing win time methods, the WINFT does not require the component outcomes to be monotone and does not rely on modeling assumptions for estimating state probabilities. This flexibility allows analysis of a complex and diverse set of endpoints, and includes existing win time methods as special cases. Moreover, the WINFT is estimated based on U-statistics, which provide direct framework for variance estimation and confidence interval derivation, under independent censoring and missing at random assumptions, without expensive bootstrapping. We examine the performance of the proposed WINFT estimation method through simulation studies, and illustrate the method using data from the ACTT-1 COVID-19 and HF-ACTION trials. Overall, the WINFT offers a flexible and interpretable estimand for assessing treatment in clinical trial data with complex longitudinal outcomes.

An expert who supplies examples of a quantity often also signals how plausible each one is; when is that signal worth using? We study eliciting a Bayesian prior from an expert who provides example points together with their approximate likelihoods. We propose fitting the prior by least squares–minimizing the squared discrepancy between a parametric density and the elicited likelihoods–which defines an M-estimator that remains well posed even for families whose moments do not exist. We establish consistency and asymptotic normality, and prove–under explicit regularity conditions, comprising a well-separation and a uniform-concentration requirement that we verify for the families considered–a non-asymptotic O(1/sqrt(n)) Berry–Esseen bound on its sampling distribution, uniform and nonuniform, by extending a result of Pinelis for maximum-likelihood estimators to the M-estimation setting. We then relax the assumptions that most limit the method in practice. Experts need not report on the density's own scale: an unknown reporting scale can be profiled out in closed form and estimated jointly, at no asymptotic cost for location families. The theory extends to multivariate parameters, where a directional Berry–Esseen bound follows from the multivariate delta method applied to a smooth implicit proxy for the estimator. An additive error floor in the noise model removes a degeneracy in the optimal design, making optimal designs interior. Simulations for normal and beta families confirm the predicted n^-1 error rate and the accuracy of the normal approximation at moderate sample sizes. Finally, we compare the estimator with the sample-only maximum-likelihood baseline, derive an explicit threshold on the expert's reporting noise below which the elicited likelihoods provably reduce estimation error, and calibrate that threshold against eleven datasets of human frequency judgments.

Reinforcement learning (RL) seeks to optimize sequential decisions to maximize population-level benefits over time. However, when deployed in high-stakes settings such as healthcare, RL decisions might systematically restrict some subpopulation's access to valuable services in a manner contrary to the values and goals of stakeholders. Counterfactual fairness (CF) offers a promising framework to address this problem based on causal reasoning. This paper develops a data preprocessing algorithm that, when used in tandem with policy learning, enables CF in RL. Our algorithm relies on a novel quantile distribution mapping method for sequentially estimating the counterfactual states and rewards in the data preprocessing step, subsuming common additivity assumptions used for counterfactual prediction as a special case. We theoretically prove that the per-step level of counterfactual unfairness and infinite-horizon suboptimality gap can be bounded under mild regularity conditions. We also empirically test our algorithm in numerical experiments as well as in application to a real-world interventional digital health dataset.

We study variable screening for high-frequency continuous-time regression, targeting the support of the spot regression coefficient. We establish uniform entrywise bounds for local estimates of the marginal correlation vector and covariate correlation matrix under coordinatewise truncation and finite-variation jumps. Under a local averaged $\alpha$-Hölder condition, the rate is $\{(\log p)/n\}^{\alpha/(2\alpha+1)}$, up to logarithmic factors. This yields sure screening for the marginal support. We also derive a list-size minimax lower bound that matches the localization and multiplicity order when the retained list size $d$ satisfies $\log(p/d)\asymp\log p$. For regression-support recovery under local dependence, we propose partial-residual spot iterative screening with pathwise Lepski adaptation. Under sparse eigenvalues and standardized beta-min, the procedure has the sure screening property and requires only correlation columns associated with the current active set. Simulation studies and a high-frequency factor application illustrate its performance.

Every chi-square-family statistic for a two-way contingency table (Pearson's X^2, the power-divergence members, the variance-stabilized T_root) is referred to an approximate null distribution (chi-square, moment-matched chi-square, saddlepoint, or bootstrap) that miscalibrates on sparse or heterogeneous tables, where the exact distribution is a lattice of atoms rather than a continuous curve. The exact reference conditions on both margins (the multivariate Fisher noncentral hypergeometric law) but is usually treated as uncomputable, because the number of tables sharing a margin is astronomical. We show that the exact conditional moments and the exact conditional distribution of any chi-square-family statistic are computable without enumerating tables, because the statistic is additive over cells and the margin-conditional law factorizes cell by cell: a dynamic program walks the table one cell at a time, carrying per row-capacity state either a few moment accumulators, a map from statistic value to probability, or one complex number (the characteristic function), at a cost set by the number of margin states rather than the number of tables. The moment engine computes the exact conditional moments of a 5x5 table at three per cell (about 79 billion tables on the margin fibre) in three seconds; the distribution engine returns the exact tail, verified against enumeration to machine precision, where a three-moment chi-square fit misses it by up to 0.28. The engine furnishes Monte-Carlo-free ground truth against which any approximate reference can be measured, extends to the non-null alternative, and is the computational substrate beneath the exact conditional interval and reporting tools of the companion papers.

Neural network classifiers trained by cross-entropy minimization are highly sensitive to label noise and adversarial contamination. While robust alternatives offer bounded influence and resistance to corruption, their statistical foundations in the deep learning setting are insufficient due to a fundamental difficulty: neural parameterizations are non-identifiable, so the population loss minimizer is an equivalence class of parameters, not a unique point. We develop a consistency theory for robust neural classifiers based on the S-divergence family that requires no identifiability assumption. Casting training as stochastic optimization over a non-identifiable parameter space, we prove that empirical S-divergence minimizers converge to the population-optimal equivalence class under mild regularity conditions, and verify these conditions for three architecture choices. We further establish that limit points of the robust training algorithm are stationary points of the empirical objective. Experiments on vision and language benchmark datasets confirm that S-divergence training maintains clean-data accuracy while exhibiting performance competitive with existing robust methods.

Computer models of complex engineering systems rely on proper tuning of their model parameters to ensure accurate predictions of the system behavior. The challenge of effectively calibrating many-parameter models is the difficulty of sampling in high-dimensional spaces and the computational expense of generating a large number of samples to characterize the calibrated parameter distributions. The method of active subspaces has been shown to be effective at constructing low-dimensional latent spaces for Bayesian inverse problems when the misfit function (i.e., negative log-likelihood) is treated as the function of interest. On the other hand, works that implement surrogate modeling for inference often focus on approximating the predictive model itself. In this work, an integrated dimension reduction and surrogate modeling framework for efficient and robust model calibration based on the Kennedy O'Hagan framework is proposed, with the following key components. First, an active subspace of the misfit function is identified. Then, a surrogate model for the misfit is constructed in this low-dimensional latent space. Care is taken to ensure that the assumed probabilistic structure of the misfit surrogate is compatible with the structure imposed on the misfit by the observation noise and computer model discrepancy. Further, a generalized likelihood function is defined that can account for the misfit surrogate uncertainty along with the other usual sources of uncertainty, e.g., experimental noise, model inadequacy, etc. This general formulation is shown to be valid for surrogates of any deterministic bijective function of the original likelihood, not just the misfit. Finally, a strategy for incorporating the uncertainty in identifying the active subspace is included.

We study high-dimensional one-sample mean testing for time series with strong common serial dependence driven by latent dynamic factors. After estimating the dynamic factor loading space from lagged autocovariance, we project the data onto its orthogonal complement and construct three factor-adjusted tests: a max test for sparse alternatives, a quadratic test for dense alternatives, and a Cauchy combination test for unknown sparsity. The idiosyncratic component is allowed to be non-Gaussian sub-Gaussian vector white noise. We establish the Gumbel limit of the max statistic, the normal limit and local power function of the quadratic statistic, their asymptotic independence, and the validity of the Cauchy combination. In the strong-factor case, the refined projection expansion shows that the quadratic statistic remains valid for dimensions as large as $p=o(n^2)$. A random-loading residual bootstrap is developed for finite-sample calibration. Simulation studies and a real data application demonstrate reliable size control and competitive power for high-dimensional observations with strong dependence.

Ranking data arise in scientific and machine learning applications, including recommendation systems, information retrieval, voting, marketing, and AI preference ranking from human feedback. Existing statistical work has primarily focused on inference tasks such as preference estimation, rank aggregation, and ranking prediction. However, generating realistic synthetic rankings from an observed population is important for privacy-preserving data sharing, benchmark construction, simulation, and uncertainty quantification. This task is challenging because rankings are high-dimensional combinatorial objects with non-Euclidean dependence structures, while ranking populations often exhibit substantial preference heterogeneity. We propose a framework for population-level generative modeling through a latent preference simplex embedding. It estimates a low-dimensional latent preference simplex through a likelihood-based ranking model, leverages flow matching to learn the population distribution of latent preferences, and generates new rankings through the fitted probabilistic ranking model. We show that ranking generation admits an oracle reduction to latent distribution learning and derive finite-sample generative guarantees that clarify how the number of items, ranking length, and latent dimension affect accuracy. Experiments on synthetic and real datasets demonstrate improved population-level fidelity and provide a statistically interpretable representation of preference heterogeneity.

This article proposes a comprehensive framework for the identification and statistical quantification of impulsive behavior in signals, with a primary focus on condition monitoring. We concentrate on evaluating impulsivity, where such behavior results from normal operation or additional disturbances. Such an evaluation is crucial in the context of local damage detection, as the presence of impulsive disturbances significantly complicates the machine condition monitoring process. To address problem of impulsivity assessment we introduce a two-stage methodology to make processing workflow effective. First, we propose an objective selection criterion for "best-performing" impulsivity measures based on the Mann-Whitney statistic, allowing for the systematic comparison of various classical and advanced metrics across diverse signal scenarios. Second, we establish a formal procedure for assessing statistical significance using bootstrap-driven resampling and define a magnitude index to quantify the intensity of detected impulsivity. The framework is validated through extensive Monte Carlo simulations for three reference signal scenarios and applied to real-world vibration data from an industrial compressor. By systematizing existing measures and providing a statistically grounded pipeline, this research extends prior works, offering a scalable tool for distinguishing between diagnostically useful signals and those corrupted by anomalous interference.

Augmented inverse-probability weighting (AIPW), targeted maximum likelihood estimation (TMLE), and double/debiased machine learning (DML) are three routes to the same efficient influence function for the average treatment effect — settled theory we treat as background. This tutorial's contribution is its worked, shared-nuisance reconciliation on real data: what a practitioner must actually match for the routes, and the software packages, to agree. Working the effect of smoking cessation on weight change in the open NHEFS data (n=1566) with one shared Super Learner library and identical cross-fitting folds, we build all three estimators by hand from one influence function, in full-sample, cross-fit, and double-cross-fit variants; the six resulting doubly-robust estimates span only 3.32–3.42 kg, consistent with the established benchmark. We then reconcile the same estimand across our engine and the tmle, AIPW, DoubleML, and tmle3 packages. At their defaults the estimates span 3.32–3.49 kg. Once the library and folds are matched and single-split noise is averaged out, the three library-sharing implementations agree to within 0.01 kg — so the residual spread traces to the nuisance library, folds, and repetitions, not to the estimator label. Because the estimators share one influence function, they agree under good overlap; under a positivity violation the pooled ATE is not identified without additional extrapolation assumptions, and their finite-sample estimates can then diverge sharply. We therefore place a positivity diagnosis ahead of estimator choice, illustrate the failure on a no-overlap example, and close with a reporting checklist. Open-source R code reproduces every number.

Many functional data analyses reduce random functions to scalar summaries or conditional mean curves. This is limiting when we wish to understand how covariates affect the distribution of entire functional responses, including their shape, timing, or variability. We study the problem of estimating conditional laws of functional outcomes and show that these objects can be estimated and evaluated in a practical nonparametric framework. To do this, we introduce functional distributional random forests, which estimate each conditional law as a covariate-dependent distribution over sampled functions by training a random forest to minimize a kernel-based maximum mean discrepancy within the leaf nodes of the decision tree. This supports inference on arbitrary functionals of the conditional distribution while keeping predictive samples tied to realistic curves. We consider a variety of kernels defined on function spaces, including Sobolev and operator-induced kernels. We also provide conditions for consistency of our estimator and develop scoring rules for comparing it to baseline estimators. In simulations, our method recovers distributional changes that are missed by baseline methods. In an application to NHANES accelerometer data, it identifies interesting covariate-associated changes in both median activity profiles and predictive dispersion.

A predictive correction is a prespecified modification of an existing predictive distribution intended to reflect an anticipated change in future outcomes given their inputs, motivated, for example, by instrument recalibration, assay drift, or a known intervention. We study how to accumulate anytime-valid evidence that such a correction predicts incoming target outcomes better than the uncorrected source predictive distribution. A fixed nonnegative tilt transforms the source predictive into a corrected predictive, and the corrected-to-source predictive likelihood ratio is a conditional e-value whose running product forms an e-process. This process remains valid under optional stopping and arbitrary input sequences, including adaptively selected ones, while its logarithm equals the cumulative predictive log-score advantage of the correction. A conditional drift decomposition characterizes evidence growth under an arbitrary target predictive distribution, and a correction-dependent half-space identifies misspecified target distributions for which the same false-confirmation bound continues to hold. When the predictive likelihood ratio is strictly positive, its reciprocal yields an anytime-valid refutation boundary, while an overshoot identity explains why the realized null crossing probability may fall below the nominal level. Label-shift, conditional mean and variance, subgroup-specific, and exponential-family corrections arise as special cases. Prespecified mixtures accommodate uncertainty over corrections, predictable tilts permit adaptive betting, and beyond-tolerance comparisons target changes large enough to justify action. Cross-family calculations and synthetic experiments show that a boundary crossing supports the proposed correction relative to its reference but does not uniquely identify the mechanism responsible for the shift.

In this paper, we propose a distribution-free test for detecting changepoint in the mean direction of angular data. The uncertainty in angular measurements is quantified through the square of an angle, derived from the intrinsic geometry of the torus. It is established that, under the null hypothesis, the test statistic distributionally converges to the Kolmogorov distribution, while under the alternative hypothesis, both the consistency of the test and the asymptotic properties of the changepoint estimator are established. Through extensive simulations, we compare the empirical performance of the proposed method with two existing approaches for angular data and further benchmark it against a test based on the circular arc length distance. Finally, we demonstrate the practical utility of our approach by analyzing the timestamps of extreme events in Bitcoin, Ethereum, and Gold price datasets, where the continuous, high-frequency nature of the data is modeled in the circular framework.

In the classical normal means problem, independent train–test folds can be constructed by perturbing the data with normal randomization. Averaging over $K$ such folds yields a cross-validation estimator whose bias depends on the marginal distribution of the randomization variables, while its variance depends on their joint distribution. This raises the questions: which joint law is optimal, and how to construct the corresponding randomization scheme? We show that: (i) for smooth estimators, antithetic randomization with pairwise correlation $\rho=-1/(K-1)$ is necessary and sufficient for the reducible variance due to randomization to remain bounded as the bias vanishes; (ii) a general construction yields a class of antithetic schemes, within which the jointly normal scheme is minimax optimal; and (iii) for non-smooth estimators with finitely many jump discontinuities, antithetic randomization improves the asymptotic rate of the reducible variance, while a simple control variate restores bounded variance when the discontinuities are known.

Two-way contingency tables are tested for independence throughout applied statistics (genomics, network and text co-occurrence, pharmacovigilance, ecology, survey cross-tabulation), and the routine test reads asymptotic Gaussian tail probabilities off the table cell by cell. On the tables people actually analyze this is badly miscalibrated: across 5,543 real public two-way tables, two-thirds have counts small enough or margins heterogeneous enough that the asymptotic maximum-cell independence scan false-positives at a mean 33% against a 0.05 target, while an exact margin-conditional anchor holds at about 1%. The failure is sharpest when the departure is sparse, the association concentrated in a few cells, where the omnibus chi-square is inefficient and the sharp object is the detection boundary. In three parts we answer where a sparse signal can be seen, which combiner attains that limit, and how to calibrate at small counts. Part I specializes the sparse-detection theory of Ingster (1997), Donoho and Jin (2004), and Chhor, Mukherjee and Sen (2024) to independence: the table is a heteroscedastic Gaussian sequence whose variance profile is the expected-count table, fixed by the margins via CVe^2 = (1 + CVr^2)(1 + CVc^2) - 1, giving an exact signal map and a closed-form separation radius valid under a moderate-deviation growing-count condition. Part II shows higher criticism, with its closed-form Jaeschke-Eicker null, attains that boundary adaptively over unknown sparsity. Part III measures the small-count calibration failure (size 0.3 to 0.6 even under uniform margins, so the driver is the per-cell tail), removes it with exact margin-conditional laws (size 0.003 to 0.009), and gives a routing rule sending large-count tables to the asymptotic gate and small-count tables to the exact anchor. Every claim is reproduced from openly deposited, deterministically seeded code.

The Cox model remains the default for survival analysis, but the proportional hazards assumption is often violated and hazard ratios can be difficult to interpret. Accelerated failure time (AFT) models provide an intuitive time-scale alternative, yet flexible covariate adjustment while preserving valid inference on a target exposure remains challenging. For the partially linear AFT model under right censoring, a rank-based debiased machine learning (DML) framework remains undeveloped: the rank-based pairwise moment is not Neyman orthogonal and standard cross-fitting does not directly apply to U-statistics. We develop the first such framework by combining an orthogonalized rank-based U-statistic, a censoring-corrected influence function, and block-pairwise cross-fitting, yielding valid inference under flexible nuisance estimation. Simulations and an application to All of Us electronic health record data demonstrate finite-sample performance and practical utility.

Estimating individualized treatment regimes (ITRs) is fundamental in data-driven personalized decision-making problems, such as precision medicine. Most of the ITR literature either focuses on categorical/continuous treatments or assumes no unmeasured confounding. In this paper, we make the first attempt to estimate the optimal individualized functional treatment regime (IFTR) for observational data where the treatment is a function and unmeasured confounding is present. We establish an identification result for a class of IFTRs under the proximal causal inference framework. Based on the identification result, we develop an algorithm of finding the optimal IFTR. The appealing practical performance of the proposed method is demonstrated by a simulation study. The proposed method is applied to an accelerometry dataset collected by the US National Health and Nutrition Examination Survey to find the optimal physical activity distribution for the best of the Triglyceride-Glucose index.

Prediction intervals for multi-modal regression with tabular variables, text, images, or other input sources are difficult to calibrate when those sources disagree or one is missing. A single global quantile averages these regimes together instead of calibrating to the modality pattern observed at test time. We address this through a modality-aware conformal calibration layer. The layer trains or reuses one predictor per modality, computes a disagreement score from their predictions, and uses that score in split conformal calibration under a strict split protocol. We use the score in two complementary ways. First, a continuous disagreement-scaled method reallocates interval width across examples while preserving the usual marginal split-conformal guarantee. Second, a Mondrian (stratified) method calibrates within groups defined by disagreement or modality availability fixed before calibration, giving group guarantees under joint exchangeability of the calibration and test examples. Across four multi-modal datasets, the disagreement-scaled layer matches or improves the marginal conformal baseline in 59 of 60 paired runs for interval continuous ranked probability score (CRPS) and in 52 of 60 for interval width, while keeping empirical coverage near the 95% target. In stress tests with missing modalities, mask-matched recalibration recovers up to 19.5 percentage points of coverage in the hardest fixed-mask regime. The result is a simple, model-agnostic reliability layer for multi-modal regression systems. A project page is available at https://unco3892.github.io/modality-aware-conformal.

Large-scale assessments suppress subgroup achievement estimates below minimum sample size thresholds, such as the National Assessment of Educational Progress (NAEP) rule of 62, disproportionately affecting historically underrepresented groups. This study introduces Small Area Bayesian Dynamic Borrowing (SABDB), a unit-level small area estimation method assigning each regression coefficient its own between-area variance, so cross-area borrowing adapts to each coefficient's heterogeneity. We compare SABDB against the unit-level Hierarchical Bayesian Small Area Estimation (HBSAE) model in a simulation study and an empirical case study. The simulation mimics the NAEP eighth-grade mathematics assessment, and the case study uses the PISA 2018 dataset. Across both studies, SABDB achieved near-nominal coverage, whereas HBSAE's intervals were narrow but severely miscalibrated.

This study investigated the behavior of Principal Component Analysis (PCA) when applied to datasets with extremely large numbers of observations. Although statistical theory suggests that sampling error diminishes and sample estimates converge toward their population values as sample size increases, relatively little empirical evidence exists regarding the behavior of PCA at scales measured in billions or trillions of observations. Three datasets were analyzed: a 10-billion observation random dataset, a 1-trillion observation random dataset, and a 10-billion observation engineered dataset designed to contain three latent factors. Results showed that the PCA solutions obtained from the 10BillionRandom and 1TrillionRandom datasets were nearly identical, indicating substantial stability of PCA at extremely large sample sizes. In contrast, the engineered dataset produced three dominant principal components that accounted for 99.996% of the total standardized variance and successfully recovered the intended latent-factor structure. These findings suggest that PCA solutions converge rapidly at very large sample sizes and suggest that PCA solutions may reach practical convergence well before sample sizes reach the trillions. These findings have implications for large-scale applications in fields such as remote sensing, digital mapping, environmental modeling, and other domains where datasets routinely contain millions or billions of observations.

Hierarchical composite endpoints are commonly used in clinical trials when component out- comes differ in clinical importance. The win ratio compares patients across treatment groups according to a pre-specified order of clinical priority and has an intuitive interpretation. A strict hierarchical rule, however, uses a lower-priority endpoint only when all higher-priority endpoints are tied or non-informative. This may reduce power when the treatment effect is mainly ex- pressed through lower-priority outcomes. Recent threshold-based extensions relax this hierar- chy by allowing lower-priority outcomes to contribute when higher-priority outcomes are within pre-specified margins. We propose a continuous-weighted win ratio for hierarchical composite endpoints. The method replaces hard transitions between endpoint levels with smooth weights, so that a lower-priority endpoint can contribute gradually when the higher-priority endpoint is close to tied. The resulting estimator is a two-sample U-statistic. We derive its large-sample distribution using the Hoeffding projection and provide a consistent plug-in variance estimator. We also describe the role of the tuning parameter through local alternatives and local efficiency. The proposed class includes the strict win ratio and hard-threshold rules as special or limit- ing cases. Simulation studies show that relaxing the strict hierarchy can improve power when treatment effects are concentrated on lower-priority endpoints, but may reduce power when the highest-priority endpoint carries the main treatment signal.

The goal of this paper is to model node heterogeneity and link homophily for dynamic networks. The proposed framework brings new insights on how networks evolve over time. It also provides more sophisticated tools for the prediction of future networks with statistical guarantees. The new model accounts for the link homophily associated with both observed traits and latent traits. The joint modeling of node heterogeneity and both observed and latent homophily effects also poses the significant challenge in statistical inference, resulted from the large number of confounding parameters in the model. To overcome this, we propose a novel normalized squared loss, paving the way for efficient and stable estimation of parameters in a high-dimensional setting. We provide a rigorous theoretical analysis of the estimation method, and demonstrate its effectiveness through extensive simulations and the illustration with some real-world network data.

Understanding the relationship between home and away goal counts in football provides valuable insights into match-level dynamics. While the influence of home advantage is well-established, with historical records indicating roughly 50% of matches are won by home teams (versus about 30% by the away team), properly determining the joint goal distribution while accounting for key match factors remains under-explored. In this paper, we develop a Bayesian bivariate Conditional Poisson (BCP) regression model to explicitly capture the dependence between home and away goal counts, addressing a core limitation of traditional football scoring models that assume independence or only allow for positive correlation. The BCP regression is applied to the English Premier League (EPL) data spanning three seasons, incorporating stadium attendance and committed fouls by both teams as regressors. Inference is conducted within a Bayesian framework, enabling interpretable uncertainty quantification and model validation through posterior predictive checks. Our results reveal a negative correlation between home and away goal counts and highlight an asymmetry in how match attendance influences home versus away scoring.

Many economic and financial relationships may change gradually rather than abruptly. We study panel data models in which the coefficient vector is continuous and piecewise linear in calendar time, with a finite number of unknown kink dates at which its slope changes. We propose a penalised least squares estimator that applies adaptive weighted group penalties to the second differences of the coefficient path, and develop asymptotic theory showing that it recovers both the number and the locations of the kinks with probability approaching one. To our knowledge, this is the first panel framework to estimate an unknown number of common kink dates in a time-varying coefficient path under fixed effects. We establish that endpoint slopes converge at the usual cubic regime-length rate and interior slopes at rates determined by their own and adjacent regime lengths. We also develop a coefficient-by-coefficient extension allowing individual regressors to kink at different dates. Monte Carlo evidence supports the good finite sample properties, and we illustrate the method through an application in macro-finance, specifically the relationship between debt and growth.

This paper introduces Mixtures of Geodesic Factor Analyzers (MGFA) on Riemannian homogeneous spaces. MGFA uses a geodesic factor model within each mixture component, providing greater expressiveness than mixtures of Riemannian radial distributions and enabling clustering of manifold-valued data with anisotropic subpopulations. We establish root-$n$ consistency for the MGFA maximum likelihood estimator (MLE), thereby filling a theoretical gap for mixtures of Riemannian radial distributions as a special case. We also propose an iterative estimation algorithm and implement it on spheres, shape spaces, and hyperbolic spaces. Numerical experiments show that MGFA substantially outperforms competing methods in well-specified regimes while remaining robust under model misspecification. Finally, case studies on corpus callosum and left hippocampus shape datasets demonstrate MGFA's effectiveness for both 2D contour and 3D shape analysis.

Covariate adjustment improves the efficiency of treatment-effect analyses in randomized clinical trials, provided the adjustment targets the correct quantity. For time-to-event endpoints, two marginal targets are of primary interest: the log-rank test for the presence of a treatment effect and the marginal hazard ratio for its magnitude. Existing covariate adjustment approaches reach these targets by different ways. Augmentation adjusts the log-rank score by regressing derived outcomes on the baseline covariates within each arm. Weighting instead reweights the two arms to balance the covariates before the survival comparison is formed: inverse probability weighting does so through a fitted propensity model, while calibration weighting solves directly for weights that match covariate means. In this manuscript, we first develop balancing weighting for time-to-event endpoints, covering both calibration weights (stable balancing weights and entropy balancing) and propensity score weights, and prove that any balancing-regular weighting is first-order equivalent to the augmented log-rank score and to the root of the marginal Cox score. All three routes therefore deliver the same estimator to first order, and calibration reaches it without fitting any model. The weighted procedures thereby inherit the validity and guaranteed efficiency gain of the augmentation approach. In addition, we show that the efficiency gain grows with the prognostic strength of the adjustment covariates, while the practical caveat lies in variance estimation, for which we give recommendations to guard against finite-sample Type I error inflation. We further confirm our results through simulation studies and an analysis of the REWIND cardiovascular trial.

Polypharmacy, commonly defined as the concurrent use of multiple medications, is increasingly prevalent in aging populations and is associated with adverse health outcomes. Motivated by longitudinal studies of aging people with HIV (PWH), we study how medication use evolves over time and characterizes multimorbidity patterns. Longitudinal medication data present substantial methodological challenges, including high dimensionality, sparsity, irregular observation times, and structured pharmacologic relationships among medications. Existing approaches are typically task-specific and lack a unified framework for learning general-purpose representations of medication trajectories. We propose LOPEL (LOngitudinal Polypharmacy Embedding Learning), a two-stage self-supervised framework for learning low-dimensional representations of longitudinal medication data. In the first stage, visit-level embeddings are learned using a Gaussian process model that incorporates pharmacologic similarity through the Anatomical Therapeutic Chemical hierarchy. In the second stage, subject-level embeddings are constructed by modeling trajectories over time and defining similarity through a Wasserstein-based representation capturing temporal dynamics and uncertainty. Simulation studies demonstrate that LOPEL accurately recovers latent structure under realistic conditions with high-dimensional sparsity and irregular sampling. In an application to aging cohorts of PWH, LOPEL identifies clinically meaningful subgroups that differ in the timing, composition, and progression of medication use, highlighting heterogeneity relevant for risk stratification and clinical management.