Papers

Interference occurs when one individual's treatment or exposure affects another individual's outcome. In particular, we assume partial interference, where individuals are divided into groups such that there is no interference between individuals in different groups. In observational studies, inverse probability weighting (IPW) based on propensity scores is often used for causal effect estimation. However, under partial interference, the group-level propensity score must be estimated, and it is more likely to take extreme values than the usual individual-level propensity score. As a result, IPW estimators may have large variances. This problem can become more serious when many covariates are available. In this study, we propose an Outcome-Adaptive Lasso based on a mixed-effects logistic regression model to stably estimate causal effects under partial interference. The proposed method performs covariate selection and estimation in the propensity score model simultaneously while accounting for unobserved group-level heterogeneity in treatment assignment. Under regularity conditions, we show that the proposed method has the oracle property and that the IPW estimators based on the proposed method are consistent and asymptotically normal. Through Monte Carlo simulations, we demonstrate that the proposed method tends to select confounders and prognostic factors at high frequencies, while excluding instrumental variables and spurious variables. The results further suggest that the proposed method improves the finite-sample efficiency of IPW estimators. We evaluate the performance of the proposed method using malaria data from the Democratic Republic of the Congo Demographic and Health Survey (DHS).

In multi-institutional studies, different parties hold distinct feature blocks for partially overlapping sets of individuals. Responses may also be missing for some records. In such settings, we propose Assisted Learning with Block-Missing Data (ALB) for sparse high-dimensional linear estimation and coordinatewise inference without pooling records or relying on a coordinating server. ALB minimizes a regularized available-case quadratic loss using cyclic block updates. Each cycle communicates $O(n)$ scalars through sample-level linear summaries, regardless of data dimension $p$, and the iterates converge geometrically to the centralized solution. We derive estimation rates that separate statistical and optimization errors. For inference on a target coefficient, ALB estimates the corresponding precision column and uses a sample-level variance estimator that accounts for dependence among moments computed from overlapping samples. Under sparsity and overlap conditions, the studentized estimator is asymptotically standard normal at the $\sqrt n$ rate, even when there are no complete cases. We also study one-time perturbed covariate and response releases that reduce direct disclosure by replacing unperturbed sample-level quantities with noisy versions. Simulations and an analysis of multimodal Alzheimer's Disease Neuroimaging Initiative data indicate that ALB approximates its centralized benchmark and improves upon complete-case Lasso by incorporating partially observed records.

We propose a regression model for the extreme tail of a response variable, in which covariates rescale the tail without changing its shape. A single covariate-dependent function then characterizes the entire conditional tail, in contrast to extreme quantile regression, which targets a quantile at a pre-specified level. The tail shape itself is left unrestricted: heavy-, light- and short-tailed responses are covered by the same framework. We specify the function through a link function and a linear combination of the covariates, which is in the spirit of a generalized linear model. In estimation, we match the parametric specification to the underlying tail function under a Bregman divergence, over a region localized at the largest observations. The resulting loss is convex, and an $\ell_1$-penalty allows the number of covariates to exceed the effective sample size. The tail localization makes the asymptotic theory deviate from that for classical penalized generalized linear models. Only the tail observations selected by a random threshold are used in the statistical analysis, making them dependent. We derive the convergence rate of the penalized estimator and propose a debiased estimator that is asymptotically normal, yielding confidence intervals for individual coefficients. Its asymptotic variance is determined by the covariance of the score, which under tail localization differs from the Hessian and must be estimated separately. We apply the method to automobile insurance claims data.

Group Independent Component Analysis (gICA) is widely used to decompose high-dimensional functional MRI data into interpretable brain networks. However, conventional gICA primarily identifies components shared across subjects. This group-level assumption can limit the recovery of networks present only in individuals or subject subsets, reducing sensitivity to intersubject heterogeneity in clinical neuroimaging datasets. We introduce Copula-Linked Group ICA (CoLiG-ICA), an algorithm in the Group ICA 2.0 framework that jointly estimates template-linked, cohort-only, and subject-only brain networks within a unified model. CoLiG-ICA combines ICA-based spatial decomposition, copula-based dependence modeling, and deep learning optimization to preserve the consistency and interpretability of template-constrained ICA while enabling free components beyond the reference networks. By linking subject decompositions to shared templates and jointly estimating cohort-only and subject-only sources, CoLiG-ICA represents individual variability not captured by conventional group priors. We evaluate CoLiG-ICA using resting-state fMRI data from the UCLA-CNP dataset and compare it with conventional constrained ICA in estimating template-linked components, discovering additional free components, improving component independence, and capturing subject-level variability beyond the shared group prior. Compared with MOO-ICAR, CoLiG-ICA showed significantly lower intercomponent spatial dependence, indicating improved subject-level component independence, and significantly reduced motion-related variance in the template-linked components. Additionally, in a schizophrenia-only group analysis, CoLiG-ICA identified three additional resting-state networks beyond the 53 template-linked NeuroMark components: one sensorimotor and two visual networks.

Estimating high dimensional Generalized Structural Equation Models presents severe computational challenges. Traditional simultaneous estimators frequently suffer from numerical instability and prohibitive computational costs. Moreover, there are no tractable algorithms for families such as Poisson, negative binomial, and gamma. To overcome these limitations, this article introduces a Two Stage Quasi-Likelihood Expectation-Maximization framework. The proposed method isolates the structural model from the measurement model. First, it approximates the conditional distribution of the latent variables given the observed indicators. Second, it employs marginal quasi-likelihood estimating equations to evaluate the structural parameters, deriving the necessary conditional moments either exactly or through Monte Carlo integration. This approach completely avoids the need to evaluate the full joint likelihood. Extensive simulations demonstrate that our method drastically reduces computational runtime, providing a numerically stable framework that minimizes the mean squared error and structural bias to yield a scalable and flexible solution for analyzing complex latent variable models.

Rubin's variance ratio (VR) complements the standardized mean difference by detecting covariate imbalance in spread, but no purpose-built omnibus extension exists for more than two groups. We introduce FVR, a size-weighted quadratic combination of pairwise log-variance-ratios that generalizes VR to k groups, together with geometric-mean and maximum pairwise variance-ratio comparators. FVR has an exact relationship to Rubin's VR at k = 2 and a mathematical structure directly parallel to Cohen's f. In Monte Carlo simulations spanning four variance-driven bias mechanisms, k = 3, 4, 6, and sample sizes from 200 to 10,000, FVR and the geometric-mean statistic generally tracked downstream estimation bias at least as well as the maximum statistic. FVR retained substantially more signal when many groups shared a very small sample, while absolute-bias thresholds were unstable at N = 200. Across mechanisms, FVR < 0.10 was reasonably reassuring, values above approximately 0.30 generally indicated concern, and intermediate values were context-dependent. The statistic is implemented for arbitrary analysis weights in the Stata command varatio and illustrated using a multivalued-treatment application.

We present a computational framework for robust topology optimization that integrates polygonal finite-element discretizations, spatially correlated non-Gaussian material modeling, and non-intrusive polynomial-chaos surrogates. Spatial uncertainty in Young's modulus is represented as a homogeneous non-Gaussian random field obtained via a memoryless transformation of a truncated Karhunen-Loève expansion, ensuring physical admissibility through positivity of stiffness while preserving the prescribed autocovariance. Polygonal finite elements provide a stable discretization for density-based optimization on unstructured meshes and mitigate checkerboard artefacts and mesh bias, while the sparse polynomial-chaos expansion enables efficient estimation of low-order statistical moments required by the robust objective at a fraction of the cost of intrusive or Monte Carlo approaches. Numerical studies on a cantilever and a curved beam show that introducing non-Gaussian material variability leads to systematic load-path redistribution and a reallocation of 6-12% of the structural volume, together with a reduction in compliance scatter. The non-intrusive surrogate reproduces intrusive reference results within 3% using an order of magnitude fewer full finite-element analyses. These results demonstrate that the proposed framework offers a physically consistent and computationally efficient route to topology-optimized designs that remain reliable under realistic material uncertainty.

Personalized predictive modelling has been growing rapidly with the increasing availability of Electronic Health Records. This approach aims to improve a model's predictive performance by fitting a unique model to each individual. We train the model on a subset of the training data consisting of individuals similar to the individual being predicted, identified through some similarity metric. Earlier studies show that using a personalized model trained on a customized subset of the data leads to better prediction than using a global model trained on the full dataset. In this work, we develop a new patient similarity metric to improve the prediction of a personalized model for binary response data. Specifically, we introduce a weighted cosine similarity metric that extends the standard cosine similarity metric by assigning predictor-specific weights when computing similarity between participants. These weights are estimated using a supervised approach with the relaxed adaptive group lasso. Results from simulation studies and an analysis of intensive care unit data show that although our proposed similarity metric leads to a slight deterioration in calibration, it produces substantial gains in discrimination. Overall predictive performance measured by the Brier Score improves because the increase in discrimination outweighs the loss in calibration; therefore, our proposed similarity metric more effectively identifies similar participants, resulting in improved predictive accuracy.

Randomized trials in chronic disease settings often measure treatment benefit through recurrent non-fatal events that are truncated by death, where conventional summaries either discard recurrences, conflate the treatment effect with survival, or treat death as censoring and forfeit a causal interpretation. While-alive estimands measure event burden per unit time alive, but doubly robust estimation for the exposure-weighted while-alive rate remains undeveloped, particularly in cluster-randomized trials (CRTs). We develop a doubly robust estimator based on a local Nelson-Aalen representation, with augmented estimating equations targeting the marginal hazard of the terminal event and the weighted recurrent event rate among those alive; the estimator accommodates multiple event types through prespecified clinical weights and remains consistent if either the censoring model or the outcome working models are correctly specified. For CRTs, we define a new pair of individual-average and cluster-average estimands under informative cluster size, with inference based on cluster-level influence functions. We establish component-wise double robustness and asymptotic normality, corroborate the theory in simulations, and illustrate the methods with reanalyses of data from two completed randomized trials.

In transfer-learning settings, a model derived from abundant surrogate labels may be deployed in a target population where gold-standard outcomes are unobserved. Evaluating its target performance is essential for determining whether decisions based on the model remain reliable, yet it is difficult when gold labels are scarce, and covariate distributions differ across data sources. We study a three-sample setting with a small gold-labeled source, a larger surrogate-labeled source, and an unlabeled target. Under conditional transportability, we evaluate the surrogate-derived model against the latent gold-standard outcome in the target population. We propose cross-fitted estimators that transport information from the two labeled sources through source-specific density ratios. We also combine outcome-regression augmentation with a kernel correction for estimating the model near a threshold, accounting for uncertainty from all three samples. We establish asymptotically linear inference for TPR and FPR, consistency and pointwise inference for the ROC curve, and asymptotically normal inference for AUC. Simulations assess bias, coverage, and sensitivity to bandwidth and relative sample sizes. A retrospective temporal validation on Chatbot Arena and a semi-synthetic ACS-Income study provide validation in real-world AI applications.

Understanding whether a digital marketing campaign is effective is central to designing effective customer engagement strategies. We analyze a large-scale, longitudinal promotional email campaign conducted by a U.S.\ retailer to evaluate how value-added incentives, such as free shipping, compare with traditional price discounts in influencing customer purchasing behavior. The analysis is complicated by non-compliance, due to not opening emails, multiple longitudinal mediators, and zero-inflated mediators and purchase outcomes. To address these challenges, we develop a Bayesian causal mediation framework based on enriched Dirichlet process mixture models and estimate the causal estimands using a scalable G-computation algorithm. We show that analyses ignoring email-opening behavior substantially attenuate estimated effects. Value-added incentives consistently outperform price discounts, yielding higher estimated potential purchase amounts, with benefits accumulating over time. We design an individualized sequential emailing strategy that optimizes expected purchase count in the observed data.

Clinical studies generally assume that group-level averages are useful quantities for guiding individual-level decisions in the clinical care of each individual patient. Precision medicine has notably closed the gap towards truly individualized care through highly refined subgrouping. Today, digital health technologies and other modern sources of dense, personal "small data" enable a different approach to treatment individualization—one that seeks to characterize a single person's own recurring health patterns first and foremost, rather than identifying the best subgroup to which they might belong. In this chapter, we review the key concepts underlying n-of-1 studies, single-case designs, and other "multitudinal" approaches for digital health applications, and explore their relationships to other digital health methods. We also share some promising future directions for "esametry", the statistics of the digitized multitudes within each of us.

Many prediction problems arise with data collected in groups. In this setting, hierarchical conformal prediction (HCP) (Lee et al., 2026) provides distribution-free prediction sets for a new observation from a previously unseen group under hierarchical exchangeability. In many applications, however, prediction is conducted only after a few observations from the group of interest have already been collected. Standard HCP cannot leverage these observations, as its required symmetry conditions do not hold in this setting. At the same time, the initial sample may still be too small for standard conformal prediction applied within the test group to be informative. We develop predictive inference methods for this setting. Our proposed method, Generalized HCP (GHCP), restores the relevant symmetry needed for conformal inference by assigning the test group a randomly "donated" reference group size. GHCP further leverages the initial test group observations to improve the quality of the nonconformity scores for prediction within that group. To improve efficiency, we introduce a variant that restricts the set of eligible donors. We demonstrate the performance of the proposed method through simulations and an illustration on the American Community Survey dataset.

Causal mediation analysis requires reconstruction of counterfactual distributions to estimate natural direct and indirect effects. Inverse probability weighting estimators rely on models for treatment assignment and mediator density ratios, whereas moment balancing approaches require researchers to specify in advance which functions of the covariates and mediators should be balanced. We propose Energy Balancing Weights for Mediation Analysis (EBWMA), which targets the joint mediator-covariate distribution used to identify counterfactual means such as E[Y(1, M(0))]. Under standard identification conditions for natural effects, EBWMA constructs weights whose weighted empirical distribution approximates this target, without modeling treatment assignment, mediator density ratios, or the outcome regression. The weights minimize energy distance through two quadratic programming problems solved sequentially. In simulations with nonlinearly transformed, skewed, or binary covariates and nonlinear mediator and outcome models, EBWMA generally achieved favorable bias and root mean squared error, with uniformly lower Monte Carlo variability than gradient boosting-based inverse probability weighting and moment balancing weights. In an illustrative analysis of the National Health and Nutrition Examination Survey I Epidemiologic Follow-up Study, EBWMA gave the smallest standardized mean differences for most covariates and for the mediator.

Competing-risk cure models describe time-to-event populations with individuals immune to all event types or an event of interest, yet literature is fragmented across model families. We review 26 papers across five axes: cure definition/scope; decomposition/cure mechanism; latency; dependence, censoring, and masked causes; and estimation. We distinguish global from cause-specific cure and incidence–latency mixtures from vertical susceptibility factorizations, latent competing-causes/zero-count constructions, defective-survival models, and zero-inflated mixture or cumulative incidence function (CIF) formulations. We compare parametric, piecewise-constant, PH, AFT, transformation, CIF-based, nonparametric, and partially specified latency models for right/interval censoring, clustering, and masked causes. Mixture formulations dominate, but similar names can mask different estimands, cure mechanisms, latent-risk/censoring assumptions, and regression interpretations. Latent-failure dependence is modeled less often than cure or latency; failure–censoring dependence, within-cluster association, and masked causes occur in smaller subsets. Estimation spans likelihood and expectation-maximization (EM), including neural-network M-steps, estimating equations, inverse-probability-of-censoring weighting, Bayesian computation, and copula-graphic estimation. A reproducible defective-Gompertz analysis of public bone-marrow-transplant data shows that fitted tail probabilities require model- and endpoint-specific interpretation, not all-method comparison. Reproducibility remains limited: most implementations use custom code, few offer repository access, and no widely adopted, clearly licensed R/Python framework unifies the constructions. This taxonomy supports transparent model selection/reporting, estimation-method comparison, and needs for theory, software, benchmarking, and reproducible applications.

A multi-regional clinical trial (MRCT) is a single clinical trial conducted in multiple regions simultaneously under a common protocol, which may be used to support parallel submissions to multiple regulatory authorities. For a regional regulatory authority, treatment effects defined specifically for its own region are more relevant to consider than overall treatment effects based on all regions included in an MRCT. A regional treatment effect can be estimated consistently using local data from the region of interest; however, this approach is generally inefficient as it excludes data from other regions and ignores possible similarities between regions. On the other hand, simply pooling data across regions requires strong assumptions and may introduce bias when the required assumptions are not met. Here, we propose a simple and robust approach to estimating a regional treatment effect in a two-arm randomized MRCT. The proposed approach uses a working regression model to incorporate information from baseline covariates as well as data from other regions for improved efficiency. The model accounts for residual regional differences (after adjusting for measured covariates) using interaction terms that describe how the dependence of outcome on treatment and covariates may vary across regions. The adaptive lasso is used to identify null interactions and thus achieve selective borrowing of information from other regions. The resulting regional treatment effect estimator is consistent and asymptotically normal even when the working model is misspecified, and able to improve efficiency over local estimation when there are similarities between regions in the form of null interactions.

In the face of spatial interference, researchers are often interested in estimating treatment effects at specific points located in space. Wang et al. (2025) and Pollmann (2023) provide design-based frameworks for estimating spillover effects on points located across a range of distances from interventions. Although these frameworks are design-based, we show their proposed estimands rely on outcomes that are directly unobservable and therefore, require outcome modeling. When using modeled outcomes in practice, even the typically design-unbiased Horvitz-Thompson estimator can accrue bias as a result of the modeling error. The performance of spatial outcome models depends on the density or resolution of observed outcomes. Through simulation, we find that the bias of the estimators decays with increasing outcome density, but not with increasing numbers of intervention units, and standard errors using modeled outcomes converge to the oracle standard errors. To demonstrate the role of outcome modeling with spatial interference, we reanalyze an experiment from Collazos et al. (2021) on the effect of hot spots policing on crime and provide several suggestions for practice.

Large household surveys support precise national estimates but can become statistically fragile after repeated disaggregation by geography, sector, sex, age, and outcome category. This paper develops an auditable design-based framework for deciding how far such disaggregation can be taken in a stratified multistage survey. The framework is built around a nested contribution ledger that reconstructs each domain total through the first stage probability proportional to size expansion, the certainty-plus-random hamlet-group selection, and the second-stage household expansion. Nonlinear domain parameters are expressed as ratios of these totals and analyzed by first-order linearization. The two independent National Sample Survey subsamples then provide a natural replication variance estimator. A granularity-stability profile combines the resulting relative standard error with replicate support and concentration diagnostics, so that a detailed estimate is accompanied by evidence about whether the design can sustain it. Finite population unbiasedness of the total estimator, asymptotic validity of the ratio linearization, and unbiasedness of the two-subsample variance estimator for linearized totals are established. The method is illustrated with the 71st-round Social Consumption: Education survey, focusing on home language versus medium of instruction and reported reasons for preferring private educational institutions in India and Himachal Pradesh. The application preserves the substantive analysis in the original project while replacing ad hoc calculation with a reproducible inferential workflow.

Nominal LoRA rank is a design parameter; calibrated spectral evidence is a separate inferential quantity. This article develops a finite-sample framework for inferring effective rank structure in public foundation-model adapters. The theoretical core is an exact chi-square divergence for the fixed-dimensional Gaussian rank-one reference experiment, with an unknown signal direction integrated under a rotation-invariant reference prior. The resulting series yields a computable finite-sample Le Cam bound at concrete layer sizes, an explicit remainder bound for numerical truncation, and the rectangular Baik-Ben Arous-Peche (BBP) limit. A compact-manifold Laplace expansion shows that finite-sample likelihood evidence also depends on leading spectral gaps through the factor $s_1^{|m-n|}\prod_{i\ge2}(s_1^2-s_i^2)$, motivating joint calibration of clustered singular values. Building on these results, we introduce an empirical-null workflow for PEFT LoRA adapters: factor reconstruction, Monte Carlo $p$-values, stagewise and block testing, and module-wise and corpus-level BH reporting. In an audit of 26 public adapters, 684 modules, six architecture families, and 31,770 public-checkpoint spectra rows, calibrated effective rank is typically much smaller than nominal rank and differs systematically from 95% energy retention. A measured RoBERTa-RTE slice on $n=24$ examples illustrates the measurement path from calibrated ranks to task evaluation, without treating the slice as a utility study. The main empirical finding is that calibrated effective rank is usually far below nominal rank, and that energy retention and statistical surprise answer different questions.

We propose a family of control variate estimators for variance reduction in design-based survey sampling and causal inference, with and without interference. In these settings, inverse probability weighting (IPW) estimators are widely used, but may have large variance when sampling, treatment, or exposure probabilities are small. Building on the observation that several common estimators, including the Hajek estimator, the normalized estimator, the augmented inverse probability weighting (AIPW) estimator, and the targeted maximum likelihood estimator (TMLE), all correct the Horvitz-Thompson estimator by canceling part of its randomness, we provide a unified interpretation of these estimators as special cases of a general control variate estimator. We then construct optimal control variates that can reduce the finite sample variance compared to these common estimators. We parameterize the proposed control variates by their bases and characterize the optimal bases through a stochastic optimization formulation. In survey sampling and causal inference without interference, the optimal bases are characterized by leading eigenvectors of matrices that depend on both the design-based sampling structure and the model-based outcome uncertainty. In causal inference under network interference, the optimal bases solve a nonconvex quadratic optimization problem; we provide a $\frac{1}{2}$-approximate solution and an alternating local search heuristic. We apply the control variate estimators to the Swiss Environmental Panel survey data and the Chinese social network data, and conduct extensive simulations to show that the proposed control variate estimators can achieve substantial variance reduction.

Constraint-based causal discovery like PC and FCI depends on its conditional independence test. Partial correlation and the Generalised Covariance Measure (GCM) detect only the conditional covariance of residuals, so they miss dependence in the mean's nonlinear part, the scale, and the tails. Tests that detect more are biased inside PC, not scalable, only continuous, or not aimed at the tails. Our Generalised Feature Covariance Measure (GFCM) is valid, sensitive beyond covariance, robust inside PC, and applicable to mixed-type data. It runs the GCM template on a configurable set of residual features with conditional mean zero (centered moments and conditional quantile indicators), pooled in blocks and combined by the Cauchy rule, with a growing-knot spline nuisance at regression cost. We contribute (i) a centering result making the scale feature Neyman orthogonal, where the uncentered version is biased; (ii) the orientation asymmetry the mean-quantile construction creates inside PC, and its fix; (iii) a Phi-faithfulness theory under which PC with GFCM recovers the CPDAG of the set's detection class; and (iv) a benchmark of CI tests sensitive beyond covariance on synthetic data, semi-synthetic tail injections, and PC discovery on random DAGs. Under size-corrected power, GFCM recovers the scale and tail edges the covariance family misses and alone keeps power at the deep conditioning sets PC issues. It stays calibrated as n grows, whereas FFCI, the boosted GCM, and the partial copula test do not, and it handles mixed-type data directly. Inside PC at scale it attains the lowest skeleton SHD among tests that stay calibrated, while the others inflate false edges. Validity rests on an additive nuisance, and the tail advantage is shown on simulated and semi-synthetic data, as no fully real benchmark with both heavy tails and known structure exists.

Decorrelating a multivariate system need not destroy source-variable identity. We introduce Constrained Oblique Rotation with Anchored Loadings (CORAL), which minimizes residual cross-correlation while guaranteeing a declared minimum correlation between each transformed variable and its designated source. For a p-variable correlation matrix $R$, we show that every exact decorrelator can be written as $R^{-1/2}Q$ for some orthogonal matrix $Q$, and define $\rho_\star(R)$ as the maximum common source fidelity compatible with exact decorrelation. Constructive lower bounds and rigorous analytical upper bounds tightly bracket $\rho_\star$ at [0.972,0.976], [0.959,0.962], and [0.949,0.950] in simulations with p={6,18,50}, respectively, compared with PCA's largest achievable minimum correlation between distinct principal components and matched source variables of 0.358, 0.329, and 0.172. Corresponding intervals are [0.751,0.758] for World Development Indicators and [0.826,0.835] for wine chemistry data sets. Thus, loss of source-variable identity is not inherent to exact decorrelation but depends on the decorrelator selected. CORAL uses constrained Riemannian optimization and extends to exact support restrictions.

In this paper we show how to extend quantile-stratified sampling to produce simulations from various multivariate distributions. These simulations have desirable space-filling and coverage properties relative to simulation using IID sampling. We examine the coverage performance of these simulations against IID sampling by looking at plots of ordered log-density values from the simulations.

Sufficient dimension reduction (SDR) seeks the minimal subspace of the predictors that captures the full conditional distribution of the response, which is known as the central subspace (CS). When the response is multivariate, the problem becomes considerably more challenging, particularly when the sample size is limited. Existing methods face different limitations:inverse regression approaches rely on strong distributional assumptions and matrix inversion, and their multi-response extensions suffer from severe slice sparsity; forward regression methods depend on computationally intensive iterative smoothing whose cost grows with the response dimension; and deep learning-based approaches demand large amounts of labeled data. To circumvent these shortcomings, we propose an SDR framework based on the generalized Stein's lemma. Our method constructs a cross-moment matrix between the multivariate response and the marginal score function of the predictors, and recovers the CS via its singular value decomposition. The proposed method does not rely on the linearity condition, avoids matrix inversion as well as iterative smoothing, and can leverage unlabeled data when available. We establish convergence guarantees for the proposed estimator under standard regularity conditions. Moreover, we propose a practical rank-selection algorithm to estimate the dimension of the CS. Extensive simulation studies and a real data application demonstrate that the proposed methods consistently outperform existing approaches across a variety of settings, particularly in moderate-dimensional, label-scarce scenarios with high noise levels.

In this paper, we study the effects of employing multiplier bootstrap to analyze the asymptotic distributions of the largest eigenvalues of high-dimensional sample covariance matrices in both spiked and non-spiked models. Our findings demonstrate that the multiplier bootstrap establishes several phase transitions in the limiting edge distributions of both unconditional and conditional bootstrapped covariance matrices, provided the different classes of multipliers. In the nonspiked setting, unbounded multipliers lead to Frechet or Gumbel limits for the largest eigenvalue of the bootstrapped covariance matrix, both conditionally on the observed data and unconditionally. For bounded multipliers, the unconditional model exhibits transitions among Tracy-Widom, Gaussian, or Weibull limits, determined jointly by the aspect ratio p/n, the upper-endpoint behavior of the multipliers, and the population covariance matrix. The conditional model displays analogous Gaussian and Weibull regimes; in contrast, the conditional counterpart of the unconditional Tracy-Widom regime collapses to a point mass. In the spiked setting, under suitable signal-strength conditions, the leading eigenvalues of both the unconditional and conditional bootstrapped sample covariance matrices are asymptotically Gaussian for bounded as well as unbounded multipliers, under some mild assumptions. Our theoretical results also clarify the feasibility and adaptability of the multiplier bootstrap for spectral inference in high-dimensional sample covariance models. Numerical simulations confirm the accuracy of our results and the effectiveness of the proposed spectral inference procedures, which may be of independent interest.

We study anomaly detection in temporally dependent network sequences. Methods based only on adjacency matrices, which are widely used for static networks, can miss changes in the way edges evolve. We instead represent each pair of consecutive networks by separate formation and dissolution event matrices, and embed the resulting matrix sequences using unfolded adjacency spectral embedding. Comparing these embeddings against a stationary baseline yields vertex- and network-level statistics that detect an anomalous transition and identify whether it involves formation, dissolution, or both. For autoregressive random dot product graphs, we establish uniform rowwise consistency of the transition-event embeddings and derive high-probability guarantees for vertex- and network-level detection and exact recovery of the anomalous vertex set. The theory separates a vertex's own displacement from interference caused by other vertices and from autoregressive memory. We quantify the memory left by an earlier anomaly, show that it decays geometrically after the transition mechanism returns to baseline, and give conditions under which it is negligible. A degree-corrected stochastic block model extension gives exact community recovery and provides community-reassignment, centre-shift, split, and merge anomaly statistics with high-probability detection and recovery guarantees. Simulations and applications to an international trade dataset and a primary-school contact network illustrate the performance of the proposed method, revealing a dissolution-driven trade decline followed by formation-driven recovery during COVID-19 and temporary merge-type mixing between school classes.

This paper presents a framework for constructing structure-preserving representations of complex systems through discrete approximation, and demonstrates its use in studying the vertical temperature and salinity structures in the mesopelagic zone across the global ocean using the ARGO dataset. Clustering serves as a means of organizing complexity into a finite set of structures that approximate the overall oceanic conditions, and a color encoding design then integrates these structures into a coherent map, with the three color components derived from interpretable geometric features of a profile: its initial level, its magnitude of variation, and its shape. Instead of focusing on specific depth levels or computing zonal averages within selected regions, our approach preserves the full vertical structure of individual profiles and incorporates each profile in the global ocean, capturing both fine-scale profile detail and large-scale spatial variability. By clustering over one million profiles collected over a decade, we identify and characterize representative profile shapes, which form the basis for a visualization strategy that provides an integrated, comprehensive, and interpretable presentation of the large-scale spatial distributions of these oceanic vertical patterns.

Watermarking provides a principled way to authenticate text generated by large language models (LLMs). In practice, however, the final text may be mixed-source, with watermark evidence surviving at only a subset of token positions after rewriting, insertion, deletion, or paraphrasing. Although prior work has studied global detection of watermark signals, when such signals can be localized remains unclear. We formulate watermark localization as a token-level multiple-testing problem based on pivotal statistics, with a latent indicator recording whether watermark dependence survives at each position. Under an asymptotic regime indexed by exponents for signal sparsity, next-token concentration, and effective-vocabulary growth, we derive a sharp boundary for global detection and phase transitions for discovery and classification within the class of coordinatewise pivot-based localization rules. We show that discovery is strictly harder than detection and that consistent classification is impossible across the parameter regime within this class. We then develop an adaptive thresholding method that does not require knowledge of the exponents or time-varying next-token distributions, but uses a data-driven estimate of the surviving watermark fraction. The method attains the optimal discovery boundary and near-optimal discovery power relative to homogeneous pivot-based rules. Simulations support the theoretical phase transitions, while experiments on model-generated texts demonstrate practical localization performance under common edit mechanisms.

Thinned count autoregressions are popular for modelling infectious disease surveillance data due to their flexibility and interpretability. However, such a model is challenging to fit since it involves high-dimensional and serially correlated integer-valued unknowns. One solution is to consider an analogous continuous-valued surrogate model whose values are post-hoc mapped to integers. This procedure produces biased estimates as inference is performed using samples from such a surrogate model and not the thinned count autoregression itself. In this work, we propose a sequential importance sampling procedure to correct this misspecified model. We demonstrate its validity in a simulation study and its applicability for epidemic curve reconstruction using rotavirus data from Germany and meningococcus data from France.

The Widely Applicable Information Criterion (WAIC) is a cornerstone of Bayesian model selection, but its conditional independence assumption renders it inappropriate for sequential and spatially correlated data - a limitation that leads to systematically underestimated model complexity and over-optimistic predictive assessments. We introduce CC-WAIC, a principled generalization of WAIC that explicitly incorporates the full posterior covariance structure of log-likelihood contributions. CC-WAIC reduces exactly to WAIC when independence holds, making it a natural extension rather than an ad-hoc modification. To overcome the prohibitive computational cost of the full covariance matrix, we develop a linear-time implementation using a banded covariance approximation that reduces complexity, with theoretical guarantees on the approximation error. We further introduce an effective sample size correction to mitigate finite-sample bias in MCMC estimation. Through extensive simulations on Hidden Markov Models and real-world applications to Old Faithful geyser data and S&P 500 volatility modelling, we demonstrate that CC-WAIC substantially outperforms WAIC, LOO-CV, iWAIC, and WAICNF, particularly under strong temporal dependence and limited sample sizes. The proposed criterion offers a computationally scalable and theoretically grounded tool for Bayesian model selection in the dependent data settings that are ubiquitous across modern science. Limitations include reliance on exponential mixing and exact conditional likelihoods; extensions to long-memory processes and approximate inference are discussed.