Papers

Minimizing the Mean Squared Error (MSE) is a key objective in machine learning and is commonly used for imputing missing values. While this approach provides accurate point estimates, it introduces systematic biases in downstream analyses. These biases affect key parameters such as variance, prevalence, correlation, slope, and explained variance. The root cause is that imputed values optimized for MSE are averages, which reduce the natural variability in the data. This paper demonstrates that adding noise to imputed values can effectively eliminate these biases. The required noise level is proportional to the MSE. Using a toy example in a multivariate normal setting, we compare two methods: predictive imputation, which minimizes MSE, and stochastic imputation, which incorporates random noise. Simulation results show that predictive methods systematically introduce bias, while stochastic methods preserve the data's natural variability and produce unbiased estimates. We also evaluate three popular imputation tools – missForest, softImpute, and mice – and observe consistent biases in predictive methods. These findings highlight that MSE is an inadequate measure of imputation quality, as it prioritizes accuracy over variability. Incorporating noise into imputation methods is essential to prevent biases and ensure valid downstream analyses, underscoring the importance of stochastic approaches for handling incomplete data.

The rise of large language models (LLMs) has created an urgent need to distinguish between human-written and LLM-generated text to ensure authenticity and societal trust. Existing detectors typically provide a binary classification for an entire passage; however, this is insufficient for human–LLM co-authored text, where the objective is to localize specific segments authored by humans or LLMs. To bridge this gap, we propose algorithms to segment text into human- and LLM-authored pieces. Our key observation is that such a segmentation task is conceptually similar to classical change point detection in time-series analysis. Leveraging this analogy, we adapt change point detection to LLM-generated text detection, develop a weighted algorithm and a generalized algorithm to accommodate heterogeneous detection score variability, and establish the minimax optimality of our procedure. Empirically, we demonstrate the strong performance of our approach against a wide range of existing baselines. The python implementation of our proposal is available at https://github.com/Mamba413/DetectLLMSegmentation.

Training-free conditional diffusion provides a flexible alternative to task-specific conditional model training, but existing samplers often allocate computation inefficiently: independent guided trajectories can vary widely in quality, and additional function evaluations along a single trajectory may not recover from poor early decisions. We propose Tempered Guided Diffusion (TGD), an annealed sequential Monte Carlo framework for training-free conditional sampling with diffusion priors. TGD targets tempered posterior distributions over the clean signal, using noisy diffusion states only as auxiliary variables for proposing reconstructions and propagating particles. Particles are reweighted by incremental likelihood ratios, resampled, and propagated across noise levels, concentrating computation on trajectories plausible under both the prior and observation. Under idealized exact-reconstruction assumptions, full TGD yields a consistent particle approximation to the posterior as the number of particles grows. For expensive reconstruction tasks, Accelerated TGD (A-TGD) retains early particle exploration but prunes to a single high-likelihood trajectory partway through sampling. Experiments on a controlled two-dimensional inverse problem and image inverse problems show improved posterior approximation and favorable wall-clock speed-quality tradeoffs over independent multi-trajectory baselines.

Bayesian predictive inference propagates parameter uncertainty to quantities of interest through the posterior-predictive distribution. In practice, this is typically performed using a two-stage procedure: first approximating the posterior distribution of model parameters, and then propagating posterior samples through the predictive model via Monte Carlo simulation. This sequential workflow can be computationally demanding, particularly for high-fidelity models such as those governed by partial differential equations. We propose a variational Bayesian framework that directly targets the posterior-predictive distribution and jointly learns variational approximations of both the posterior and the corresponding predictive distribution. The formulation introduces a variational upper bound on the Kullback–Leibler divergence together with moment-based regularization terms. The variational distributions are trained in an amortized manner, shifting computational effort to an offline stage and enabling efficient online inference. Numerical experiments ranging from analytical benchmarks to a finite-element solid mechanics problem demonstrate that the proposed method achieves more accurate predictive distributions than conventional two-stage variational inference, while substantially reducing the cost of online predictive inference.

In the context of a large system of $N$ neurons interacting through spike events in a mean-field regime as $N\to \infty$, we characterize the estimation of a multidimensional parameter in the spiking rate, when the neural states are observed over a fixed time horizon. We first prove the local asymptotic normality (LAN) property and leverage classical theory to establish the asymptotic efficiency of the maximum likelihood estimator. While the theory of Ibragimov and Hasminskii yields strong results, up to global asymptotic minimax bound, its applicability appears currently limited to models without state resets at spike times. Following then Höpfner's classical approach, we nevertheless derive, in a general setting including neuron reset, the consistency, asymptotic normality and local asymptotic minimax optimality of the estimator. Keywords: Local Asymptotic Normality (LAN); Mean-field regime; Interacting particle system; Multidimensional parameter estimation; Jump rate estimation; Maximum likelihood estimator (MLE); Asymptotic minimax optimality

We investigate the Poisson regression method for Markov and semi-Markov jump processes from a nonparametric angle, allowing the lengths of the time and duration intervals in the partition to vary with the number of observations. Imposing no structural assumptions on the true intensities, we obtain asymptotic normality of the occurence/exposure rates under appropriate shrinking conditions on the partition lengths. We derive asymptotic normality results for both Markov and semi-Markov models using only classical central limit theorems and elementary results for counting processes. All results are illustrated on both simulated and real data.

We propose a new scalable framework for spatio-temporal data fusion with multi-fidelity Gaussian processes (MFGPs) that enables fully likelihood-based inference for both stationary and non-stationary fidelity integration. The framework is designed for environmental applications, where abundant but noisy low-fidelity data (e.g., satellite or reanalysis products) must be fused with sparse yet accurate high-fidelity in-situ observations to obtain high-resolution reconstructions. Our key methodological contribution is a decomposed multi-fidelity covariance formulation that allows the Vecchia approximation to be applied directly to the latent low-fidelity and discrepancy processes. Combined with a Woodbury-based reconstruction, this yields a numerically stable and computationally efficient evaluation of the joint marginal likelihood without ever forming the full multi-fidelity covariance matrix. In addition, we introduce a generalized least squares (GLS) mean-removal strategy with fidelity-specific offsets, preventing systematic biases from being absorbed into cross-fidelity dependence. We validate the proposed approach through extensive experiments on synthetic data and a large-scale real-world application to wind speed reconstruction in the Lombardy region of Italy. The results show that the proposed Vecchia-based MFGP closely matches exact multi-fidelity inference in controlled settings, while substantially outperforming standard single-fidelity spatio-temporal Gaussian processes in terms of predictive accuracy, correlation, and representation of local variability in realistic large-data scenarios.

Given a statistical model, we propose a novel estimation method that yields randomised estimators for the unknown distribution of an observed random variable. We establish non-asymptotic bounds for the performance of these estimators and demonstrate their robustness to potential model misspecification. Notably, these properties are established by circumventing the use of concentration inequalities and empirical process theory. We provide an illustration of this approach to the problem of estimating the intensity of a Poisson process.

Tools from random matrix theory have become central to deep learning theory, using spectral information to provide mechanisms for modeling generalization, robustness, scaling, and failure modes. While often capable of modeling empirical behavior, practical computations are limited by matrix size, often imposing a restriction to models that are too small to be realistic. This motivates the inference of properties of larger models from the behavior of smaller ones. Free decompression (FD) is a recently proposed method for extrapolating spectral information across matrix sizes, but its utility is currently limited by strong assumptions that preclude its implementation on more realistic machine learning (ML) models. We use algebraic spectral curve theory to provide a general FD methodology for spectral densities whose Stieltjes transform satisfies an algebraic relation, a modeling assumption that is more likely to hold in practice. This recasts FD as an evolution along spectral curves which can be readily integrated. Our framework enables the expansion of spectral densities that have multiple or multi-modal bulks, that exist at multiple scales, and that contain atoms, all characteristic of real-world data and popular ML models. We demonstrate the efficacy of our framework on models of interest in modern ML, including Hessian and activation matrices associated with neural networks and large-scale diffusion models.

Multivariate functional principal component analysis (MFPCA) is a powerful dimension reduction technique for analyzing multiple functional variables simultaneously. However, existing MFPCA methods assume that all functional observations are recorded over a common, fixed domain. This assumption is often violated in practical applications where the observation period varies across subjects, leading to what is known as variable domain functional data. We propose a novel approach for MFPCA that explicitly accommodates variable domains by extending existing multivariate functional principal component analysis to the variable domain setting. Our methodology involves performing univariate variable domain FPCA for each functional variable separately, stacking the resulting univariate scores, and then smoothing the empirical covariance matrix of these stacked scores over the domain length. This allows us to estimate multivariate eigenfunctions and scores that properly account for varying observation periods. We demonstrate through extensive simulation studies that our proposed method outperforms approaches that ignore the variable domain structure and rely on binning strategies. The practical utility of our method is illustrated through an application analyzing body temperature and capillary oxygen saturation (SpO$_2$) trajectories in COVID-19 hospital admitted patients, where patients experienced varying lengths of stay and monitoring periods.

The study of infectious disease epidemiology for multi-type disease pathogens requires modelling techniques that account for the complex interactions existing between strains across geography and time. In this paper, we propose a novel multi-type spatio-temporal infectious disease model to better support the understanding of these pathogens. We formulate a joint state-space for all epidemics arising for a given multi-type pathogen as well as biologically informed representations of how these epidemic states may interact. We introduce the use of several copula models to uncover the dependence structure of epidemics between strains. We develop a computationally efficient Markov chain Monte Carlo (MCMC) sampling scheme for all proposed models. We also provide robust model comparison techniques using bridge sampling and importance sampling to evaluate model evidence in high-dimensional space. We demonstrate the performance of our proposed models using simulated datasets, where simulated epidemics were successfully identified and associated parameters correctly inferred. The proposed models were also fitted to monthly multi-type incidence data on invasive meningococcal disease from 26 European countries. The accompanying software is freely available as a R package at https://github.com/Matthewadeoye/MultiOutbreaks.

Empirical Bayes methods are widely used for large-scale estimation and inference in the Poisson means problem. Existing results establish theoretical properties of the nonparametric maximum likelihood estimator (NPMLE) for optimal posterior mean estimation, but comparatively less is known about uncertainty quantification (i.e., construction of confidence sets). Two main challenges in constructing confidence sets for the latent parameters based on the NPMLE are its discreteness and its slow rate of prior estimation. We resolve these limitations by introducing a smooth NPMLE that models the prior as a Gamma mixture, which is a flexible class capable of approximating a wide range of continuous priors on $(0,\infty)$. This procedure preserves the convex optimization structure of the classical NPMLE. The smooth NPMLE achieves the optimal nearly parametric rate for posterior mean estimation. Moreover, it achieves a polynomial convergence rate for prior and posterior density estimation under a compact support assumption on the mixing distribution. Based on the smooth NPMLE, we construct plug-in empirical Bayes confidence sets that mimic the oracle optimal (in terms of expected length) marginal coverage sets. We show theoretically and empirically that these sets achieve asymptotically exact marginal coverage and are substantially shorter than existing methods.

We develop high-dimensional goodness-of-fit tests for elliptical models by testing radial–directional independence after affine standardization. The method forms coordinatewise correlations between the log-radius and directional components, using a sum statistic for dense departures, a max statistic for sparse departures, and a Cauchy combination for adaptation. We derive oracle null limits, prove asymptotic independence of the sum and max components under both the null and a balanced local alternative, and establish validity of high-dimensional Hettmansperger–Randles plug-in standardization under explicit perturbation rates. Simulations and data analyses show stable size control, dense–sparse power complementarity, and interpretable coordinate-level diagnostics.

Score-based diffusion can be defined intrinsically on the manifold of quantum pure states, $\mathbb{CP}^{d-1}$ with the Fubini–Study metric, but no closed-form transition density is available, so the score must be supervised by a local-time teacher taken from the Euclidean limit of the diffusion in normal coordinates. This paper is about what makes that teacher work, and where it stops working. Three training choices turn out not to be incidental: the increment must be divided by the diffusion clock rather than by the elapsed time, since the published expression assumes unit diffusion and a non-unit schedule introduces a time-change mismatch varying by a factor of $400$ across the horizon; the logarithm and exponential maps should be the closed-form Fubini–Study ones, which is the largest single effect we measure; and the global phase must be randomised, because horizontal projection alone does not make a score network descend to the quotient. With these choices the model improves on the published Riemannian local-time baseline in every cell of an eight-benchmark, four-metric comparison over ten seeds, significantly on five of eight after Holm correction, and beats an ambient Euclidean baseline by an order of magnitude everywhere. We then bound what the approximation costs by replacing it with the exact heat kernel of $\mathbb{CP}^{d-1}$, computable up to complex dimension seven, where the local-time teacher loses a factor of $1.2$ to $2.7$.

Over time, clinical trials have increasingly incorporated complex design and analysis elements such as interim analyses, adaptations, multiple endpoints, and sophisticated multiplicity schemes for multiple endpoints and/or treatment arms following the paradigm of frequentist inference. In frequentist clinical trials multiplicity can come from (at least) four sources: multiple looks at the data, multiple endpoints, multiple populations, or multiple treatment comparisons. Normally, Type 1 error control across the multiple hypotheses is implemented to control chance of false positive decisions. To achieve this advanced techniques such as adaptive designs or graphical multiple testing procedures have been developed and are used in the design of clinical trials. However, these methods focus on hypothesis testing while subsequent estimation remains crucial to allow for a benefit-risk assessment and further use of the results by various stakeholders. Through examples, we illustrate challenges in estimation and transparent communication. In general, there are no simple solutions to this conceptual and communicational challenge. The purpose of this paper is to generate awareness of these issues and initiate a discussion about how to address them moving forward.

Self-supervised learning (SSL) excels at finding general-purpose latent representations from complex data, yet lacks a unifying theoretical framework that explains the diverse existing methods and guides the design of new ones. We cast SSL as latent distribution matching (LDM): learning representations that maximize their log-probability under an assumed latent model (alignment), while maximizing latent entropy to prevent collapse (uniformity). This view unifies independent component analysis with contrastive, non-contrastive, and predictive SSL methods, including stop gradient approaches. Leveraging LDM, we derive a nonlinear, sampling-free Bayesian filtering model with a Kalman-based predictor for high-dimensional timeseries. We further prove that predictive LDM yields identifiable latent representations under mild assumptions, even with nonlinear predictors. Overall, LDM clarifies the assumptions behind established SSL methods and provides principled guidance for developing new approaches.

We study expected generalization bounds for the Hierarchical Federated Learning (HFL) setup using Wasserstein distance. We introduce a generalized framework in which data is sampled hierarchically, and we model it with a multi-layered tree structure that induces dependencies among the clients' datasets. We derive generalization bounds in terms of Wasserstein distance under the Lipschitz assumption on the loss function, by applying a supersample construction that allows us to measure the sensitivity of the algorithm to the change of a single node in the sampling tree. By leveraging the FL structure, we recover and strictly imply existing state-of-the-art conditional mutual information (CMI) bounds in the case of bounded losses. We also show that our bound can be applied together with Differential Privacy assumptions, to recover generalization bounds based on algorithmic privacy. To assess the tightness of our bounds, we study the Gaussian Location Model (GLM) and show that we recover the actual asymptotic rate of the generalization error.

We consider function optimization as a sequential decision making problem under budget constraint. This constraint limits the number of objective function evaluations allowed during the optimization. We consider an algorithm inspired by a continuous version of a multi-armed bandit problem which attacks this optimization problem by solving the tradeoff between exploration (initial quasi-uniform search of the domain) and exploitation (local optimization around the potentially global maxima). We introduce the so-called Simultaneous Optimistic Optimization (SOO), a deterministic algorithm that works by domain partitioning. The benefit of such approach are the guarantees on the returned solution and the numerical efficiency of the algorithm. We present this machine learning approach to optimization, and provide the empirical assessment of SOO on the CEC'2014 competition on single objective real-parameter numerical optimization test-suite.

We develop graph-based methods for semi-supervised learning based on label propagation on a data similarity graph. When data is abundant or arrive in a stream, the problems of computation and data storage arise for any graph-based method. We propose a fast approximate online algorithm that solves for the harmonic solution on an approximate graph. We show, both empirically and theoretically, that good behavior can be achieved by collapsing nearby points into a set of local representative points that minimize distortion. Moreover, we regularize the harmonic solution to achieve better stability properties. We also present graph-based methods for detecting conditional anomalies and apply them to the identification of unusual clinical actions in hospitals. Our hypothesis is that patient-management actions that are unusual with respect to the past patients may be due to errors and that it is worthwhile to raise an alert if such a condition is encountered. Conditional anomaly detection extends standard unconditional anomaly framework but also faces new problems known as fringe and isolated points. We devise novel nonparametric graph-based methods to tackle these problems. Our methods rely on graph connectivity analysis and soft harmonic solution. Finally, we conduct an extensive human evaluation study of our conditional anomaly methods by 15 experts in critical care.

The goal of this thesis is to investigate the structural properties of certain sequential problems in order to bring the solutions closer to a practical use. In the first part, we put a special emphasis on structures that can be represented as graphs on actions. In the second part, we study the large action spaces that can be of exponential size in the number of base actions or even infinite. For graph bandits, we consider the settings of smoothness of rewards (spectral bandits), side observations, and influence maximization. For large structured domains, we cover kernel bandits, polymatroid bandits, bandits for function optimization (including unknown smoothness), and infinitely many-arms bandits. The thesis aspires to be a survey of the author's contributions on graph and structured bandits.

Sequential hypothesis tests are widely adopted as a principled way to perform multiple tests on data that arrives over time. In particular, researchers frequently utilize group sequential hypothesis tests (GST) to test the same hypotheses at K times or "groups" while data arrives sequentially. In this setting, many methods have been proposed to allow researchers to uniformly control type-1 error across K checks (often known as various alpha-spending budgets). Although these methods are all successfully valid in controlling uniform type-1 error, it is not clear which of these methods are optimal when trying to reject the null as soon as possible. In this paper, we directly optimize the rejection criterion in the GST setting under the same constraints of controlling type-1 and type-2 errors. We use a sample average approximation combined with mixed integer linear programming (S-MILP) approach for this problem and show how our S-MILP approach dominates classical GST procedures such as Lan-DeMets, Pocock, and O'Brien-Fleming methods. We also find that the optimal solution typically aggressively spends the alpha-budget early, shedding insight to the long-standing debate of which alpha-spending budgets are more efficient. We finally apply our optimal S-MILP approach to a recent study on acute kidney injury interventions and find our optimal S-MILP approach can reach the same statistically significant conclusion faster than the original study and other GST methods.

Contextual MDPs are powerful tools with wide applicability in areas from biostatistics to machine learning. However, specializing them to offline datasets has been challenging due to a lack of robust, theoretically backed methods. Our work tackles this problem by introducing a new approach towards adaptive estimation and cost optimization of contextual MDPs. This estimator, to the best of our knowledge, is the first of its kind, and is endowed with strong optimality guarantees. We achieve this by overcoming the key technical challenges evolving from the endogenous properties of contextual MDPs; such as non-stationarity, or model irregularity. Our guarantees are established under complete generality by utilizing the relatively recent and powerful statistical technique of $T$-estimation (Baraud, 2011). We first provide a procedure for selecting an estimator given a sample from a contextual MDP and use it to derive oracle risk bounds under two distinct, but nevertheless meaningful, loss functions. We then consider the problem of determining the optimal control with the aid of the aforementioned density estimate and provide finite sample guarantees for the cost function.

Modelling and forecasting the occurrence of extreme events is especially difficult when the event process is nonstationary, with changes in both the rate at which extremes occur and the magnitude of the extremes when they occur. We approach this task by developing a Bayesian point process model for extreme events, which uses a self-exciting Hawkes process to model the rate at which extremes occur. The Hawkes process has a structure which allows events to occur in clusters, making it realistic for many types of data. We use a flexible Bayesian nonparametric approach based on the Dirichlet process to learn the temporal excitation pattern from the data. Further, we build on Extreme Value Theory by using a Generalised Pareto Distribution (GPD) to model the magnitudes of the extremes, with a hierarchical mark model allowing these magnitudes to vary across Hawkes-induced clusters. A hierarchical specification of the model results in partial pooling, allowing for more accurate GPD estimation even in clusters with only a small number of observations. We develop an MCMC algorithm to sample from the resulting hierarchical model. A simulation study confirms that the two flexible components improve prediction when the corresponding features are present in the data-generating mechanism, and across four real data sets the nonparametric Hawkes model with hierarchical GPD marks gives the best held-out predictive performance among the model variants considered.

In many Phase II statistical process control (SPC) problems, the main concern is not whether a monitored process has ever changed, but whether it is currently operating at an acceptable level. This distinction is especially important when monitoring continues after a signal, or when corrective action may restore the process. We develop Bayesian monitoring procedures for this formulation of the Phase II task. For recoverable processes that may alternate between in-control and out-of-control states, we derive recursions for the posterior probability that the process is presently in control. For sequential tracking problems in which a latent parameter evolves over time, we monitor the posterior probability that the parameter lies inside an acceptable region of behavior. The methods are studied through calibrated time-between-failure experiments, Gaussian and Binomial tracking examples, and a held-out multivariate data illustration using white wine quality measurements.

Sequential multiple assignment randomized trials (SMARTs) provide a systematic framework for constructing and evaluating dynamic treatment regimens (DTRs). In clinical studies, longitudinal biomarkers are routinely collected to monitor disease progression and define treatment response. However, the integration of longitudinal biomarker data into survival analysis for DTR evaluation within a SMART remains unexplored. We propose a joint longitudinal-survival framework to estimate DTR-specific survival outcomes within a two-stage SMART. A linear mixed model specifies the biomarker trajectory, and a relative risk model links the survival process to the current latent biomarker value. To accommodate the time-varying treatment assignment, treatment effects are parameterized through piecewise cumulative exposure terms with a structural change at the decision point. Joint-model parameters are estimated by maximum likelihood using pseudo-adaptive Gauss-Hermite quadrature for random-effect integration. Under standard causal identification assumptions, DTR-specific marginal survival outcomes are obtained through a plug-in parametric g-formula. We compare the joint-model framework to the inverse probability of treatment weighted Kaplan-Meier estimator and implement the multiple comparisons with the best method to identify the optimal DTR with multiplicity control. Through a simulation study and an application to a SMART for androgen-independent prostate cancer, we demonstrate that the joint-model framework produces unbiased estimates under correct model specification, with substantial efficiency gains and higher accuracy in optimal DTR identification. This framework exploits prognostic information from longitudinal biomarkers, enables valid causal inference for DTR-specific survival outcomes, and provides a generalizable analytic tool for developing evidence-based adaptive treatment strategies.

This paper studies the statistical estimation of exact Wasserstein barycenters. Existing non-asymptotic results for empirical barycenters exhibit a severe curse of dimensionality. Motivated by the semi-dual formulation of the barycenter problem and its associated Sobolev optimization geometry, we develop a smoothness-aware approach that combines density estimation with Sobolev geometric structure to estimate the population barycenter. We establish nonparametric convergence rates for estimating both the barycenter functional and its minimizer, demonstrating how smoothness can substantially improve statistical performance.

Detecting observations from a minority class under severe class imbalance is a central challenge in applications such as fraud detection, medical screening, and industrial quality control. In these settings, each positive prediction triggers a costly follow-up action, an MRI scan, a transaction audit, whose execution is subject to real operational constraints. This paper proposes a formal classification framework under capacity constraints: given a user-defined bound limit $b$ on the proportion of observations that can be labeled as belonging to the minority class, the goal is to find the classifier that maximizes sensitivity on that class. We characterize the optimal classifier under this constraint and establish its equivalence with the classical Bayes classifier under a reweighting of the prior probabilities. We also introduce a capacity-adjusted performance metric $M$ that accounts for the effective detection rate when the capacity constraint is binding. The framework is implemented on top of standard learning methods, k-NN, SVM, random forests, and neural networks, and statistical consistency is established for each. We further show that these methods reduce to post-hoc thresholding when no hyperparameters are oriented toward the capacity-constrained objective, and introduce a capacity-aware support vector machine that exploits the constraint during training and achieves the strongest empirical performance. Experiments on the Taiwanese credit card default dataset confirm that capacity-constrained classifiers substantially outperform both classical approaches and SMOTE under high imbalance regimes. The framework extends naturally to multiclass settings and online environments.

Area-specific causal inference is important in many policy and survey applications, where the goal is to evaluate treatment effects for small geographic or demographic domains. Existing causal small area estimation methods, however, typically rely on a strong data requirement that treatment status is observed for all units in the population. This assumption is often unrealistic in practical survey settings, where both treatment and outcome variables are observed only for sampled units, while auxiliary covariates are available for the full population. To address this limitation, we develop a new identification strategy for area-specific treatment effects under this more realistic data structure by combining survey-only covariates with population-level auxiliary information. Based on this result, we propose a doubly robust estimator that remains consistent when either the outcome regression model or the treatment and area assignment models are correctly specified. We further derive the semiparametric efficiency bound for the target parameter and show that the proposed estimator attains this bound under regularity conditions. Simulation studies demonstrate favorable finite-sample performance, particularly in settings with small sample sizes within areas, and an empirical application illustrates the practical relevance of the proposed framework.

We provide a unified theoretical analysis of Linear Discriminant Analysis with simultaneous multilabel scatter matrix formulations and Stiefel orthogonality constraints. Our contributions span both algebraic structure and statistical guarantees. On the algebraic side, we characterize the rank of the multilabel between-class scatter matrix, showing that the effective discriminant dimensionality can strictly exceed the classical single-label bound of $C-1$; we establish a multilabel partition of variance and prove that all four Fisher objectives are equivalent under the $W^\top S_t^{ML} W = I_r$ constraint while characterizing their divergence under the Stiefel constraint; and we prove a two-sided label-distance preservation bound relating projected distances to Hamming distances in label space. On the statistical side, we establish a finite-sample $O(k_{\max}\sqrt{d\log d/n}/gap_r)$ bound on the subspace estimation error under sub-Gaussian noise with a matching $\Omega(\sigma^2 d/(n\,gap_r))$ minimax lower bound, establishing a near-minimax-optimal rate (matching up to logarithmic and $k_{\max}$ factors) for multilabel discriminant subspace estimation. We further provide high-probability distance concentration, robustness guarantees under label interactions, and a regularization analysis preserving the spectral structure when $d \gg n$. All results are verified numerically on synthetic data generated from the linear label-effect model, covering both the algebraic identities and the multilabel-specific quantities ($k_{\max}$, $\kappa(S_t^{ML})$, $\|\Gamma/n\|_2$, $\Delta_r$) that govern the statistical bounds. The numerical experiments are designed as a sanity check for the theorems rather than as an empirical benchmark; evaluation on real multilabel datasets is left to future work targeting application-oriented venues.

Externally controlled trials (ECTs) are increasingly used when randomized controls are infeasible, unethical, or insufficient, including applications in rare diseases, oncology, pediatrics, and post-approval effectiveness research. Although methodological work has expanded rapidly across causal inference, Bayesian dynamic borrowing, and hybrid trial designs, the literature remains fragmented. We adopt a six-step scientific roadmap to organize modern ECT methodology in two primary settings: (i) single-arm trials that evaluate efficacy through comparison with external controls, and (ii) hybrid controlled trials that augment the internal control arm with external controls drawn from real-world data or historical studies. The roadmap clarifies causal estimands, identifiability assumptions, and how statistical parameters arise from identification, and shows how modeling and borrowing strategies trade off efficiency and robustness, especially under covariate shift and outcome drift. Within this framework, we synthesize and evaluate recent Bayesian and frequentist developments, compare their strengths, limitations, operating characteristics, and available software, and emphasize the role of sensitivity analysis. By re-framing ECT methodology through a causal lens, this work establishes a coherent foundation for integrating external data into regulatory and clinical decision-making and highlights core challenges and opportunities for future research.