Papers

Japan's 2020 tourism shock coincided with 27% fewer buyers and 35% lower gross merchandise value at one second-hand luxury resale intermediary. We compared 2019 and 2020 co-purchase projections of recorded brand labels with a Holm-corrected battery of five label-invariant metrics and Portrait Divergence against two permutation nulls rebuilding both graphs. Neither test separated the full years (smallest Holm p 0.407; 0.097 under the design-matched paired-swap null). For one random two-block split of brands, the battery detected a planted redirection of 30% of multi-brand baskets (178 of 10,399 second-group baskets, 1.7%) in 83 of 100 simulations (83-92 of 100 across 11 splits, all at least 80); Portrait Divergence reached 80% power only at full redirection. Observed falls in mean degree and clustering resembled those of 20-30% redirection, but the planted rise in degree assortativity was absent. A composition-only model of the observed attrition would be detected in only 15 of 100 simulations. Results moved with construction: without residual and material labels (post hoc), mean degree separated the years under the paired-swap null (Holm p = 0.030), and Portrait Divergence fell below 0.05, unadjusted, at the widest construction and within April-December. Shock-period network comparisons need rebuilt-graph nulls and a per-metric detection floor.

Ignoring unobserved heterogeneity in duration models biases parameter estimates and invalidates inference, but testing for it is non-regular: the null hypothesis lies on the boundary of the parameter space and some parameters are unidentified under the null. These features render standard asymptotic theory inapplicable. This paper develops an EM test for unobserved heterogeneity in censored Weibull duration models, building on the EM approach of Li, Chen, and Marriott (2009). The test statistic has an asymptotic null distribution equal to the square of max{0, N(0,1)}, hence critical values require neither simulation nor bootstrap, and the test accommodates covariate-dependent censoring of arbitrary form. Monte Carlo simulations compare the EM test with the likelihood ratio test (LRT) of Cho and White (2010), information matrix tests, and Lagrange multiplier tests. The EM test has empirical size close to the nominal level for sample sizes of 500 or more, where the LRT remains markedly conservative and the other tests over-reject. Its size-adjusted power is comparable to that of the LRT and higher than that of the other tests. Because size adjustment requires knowledge of the data-generating process and is unavailable in practice, the EM test attains higher power than the LRT in most designs as the tests would actually be applied. In an application to the Stanford Heart Transplant data, the EM test rejects homogeneity in every covariate specification, whereas the LRT's conclusion depends on a user-chosen set of admissible parameter values and on the specification.

Spherical data appear in various applications, including cosmology and the Earth sciences. A standard statistical model for such data employs isotropic spherical random fields. The isotropy assumption may be unrealistic for many real-world datasets. Also, in various applications such as cosmology, only a single realisation of the field is observed, making isotropy impossible to assess through repeated sampling. This paper considers an alternative anisotropic model, referred to as diagonal anisotropy. To test isotropy against diagonal anisotropy, we develop score and cumulative-sum-type tests. Their asymptotic properties and Monte Carlo-based alternatives suitable for moderate sample sizes are established. The performance of the proposed methods is illustrated by numerical studies via applications to simulated data and to actual Planck cosmic microwave background radiation observations. The tests are also applied to three Planck SMICA cosmic microwave background radiation maps and demonstrate how results can vary for different data releases and map construction procedures.

The article identifies the desired properties of the reliability risk metric for the Australian Energy Sector. It proposes a comprehensive set of metrics designed to measure various aspects of the reliability risk when a large portion of power generation is composed of variable renewable energy (VRE). The suggested methods aim to balance the tradeoff between economic benefit and risk management, effectively address the tail risk, and interpret the severity of the outage in a direct way. Properties of the considered reliability metrics are investigated by using simulation studies that use time series of expected unserved energy values. Sample distributions, dependency structures and other statistical properties of the metrics are studied. The results suggest new tools and approaches that can be used to increase the efficiency of the Australian Energy Sector.

We study multiple change-point detection in multivariate time series whose distributions change in a piecewise constant manner. Distributional changes can manifest across different moment orders, from shifts in the mean and covariance to changes in higher-order moments. Higher-order moments capture increasingly rich distributional features but become difficult to estimate in high dimensions. Our tensor representation unifies moments of different orders within a common linear algebraic framework, enabling a new method to detect changes in moments of all orders up to a prescribed fixed order $p$. The resulting procedure accommodates temporal dependence and allows the dimension of the time series to grow with the sample size. Under suitable regularity conditions, the proposed procedure achieves a localization error rate that matches a newly developed minimax lower bound. We further derive limiting distributions under both nonvanishing and vanishing moment jumps and construct asymptotically valid confidence intervals in the vanishing-jump regime. Numerical experiments and real-data analyses demonstrate the method's effectiveness in detecting moment changes and a range of distributional shifts.

As black-box machine learning models become increasingly common, extracting interpretations with uncertainty quantification has become a critical challenge. One popular type of interpretation is leave-one-covariate-out (LOCO) feature importance, while prior LOCO inference methods often require data-splitting or model-refitting. A recent ensemble framework, LOCO-MP, addresses these challenges using minipatches that subsample both observations and features, but massive feature subsampling can hurt prediction in high-dimensional sparse settings. Motivated by this limitation, we consider minipatch ensembles with adaptive feature sampling guided by LOCO importance, and propose LOCO-AdaMP, which enables free LOCO inference for the resulting adaptive minipatch ensemble. We show that LOCO-AdaMP yields substantially improved predictive models while retaining asymptotically valid feature importance inference without data-splitting, despite the complex dependence between the adaptive sampling distribution and the LOCO importance statistics. Our analysis relies on a careful leave-two-out perturbation bound for the iteratively updated sampling probabilities together with the stability of LOCO scores induced by observation subsampling. Empirical results on synthetic and real datasets demonstrate advantages of LOCO-AdaMP over existing methods in predictive performance, inferential power, and stability. Overall, LOCO-AdaMP provides a flexible ensemble framework (agnostic to base models) that delivers both strong predictive performance and asymptotically valid, powerful feature importance inference for regression.

Offline reinforcement learning seeks optimal decision rules from previously collected data. In some applications, a decision can be an entire function, such as a fluence map in radiation therapy or a smooth movement trajectory in robotics. In this paper, we study the finite-sample theory for fitted Q-iteration (FQI) with functional actions in a discounted infinite-horizon setting. Three major difficulties arise in this setting: first, the absence of a Lebesgue probability density for functional actions complicates coverage descriptions; second, conventional coverage requirements can be restrictive; and third, the large functional action space makes greedy optimization in FQI challenging. To address these difficulties, we study smoothness-regularized policy search under a critic-relative coverage condition. This condition measures how well logged data distinguish relevant action-value differences without requiring an action density. Our main theorem gives finite-sample guarantees for learned-policy regret relative to the best value within a fixed smooth class of functional-action policies. The results allow trajectory lengths to be either bounded or growing and Q-functions to be fitted by either functional-input kernel ridge regression or adaptive functional neural networks. For a few examples, we can obtain polynomially decaying regret bounds in the number of logged transitions, up to logarithmic factors, with logarithmically many FQI iterations. Numerical experiments show gains of learned functional-action policies over constant-action policies and support our adoption of a critic-relative coverage condition.

In some randomized trials the primary outcome is costly or slow to measure, while auxiliary variables that predict it are available for everyone. The outcome can then be measured in a probability sample. We study an adaptive design in which the statistician who selects outcomes to measure is blinded to treatment assignment. Outcomes from an initial random sample are used to fit pooled models for the outcome and its residual variance. These set the sampling probabilities for the remaining participants and may be refitted as outcomes accumulate. After unblinding, arm means are estimated by augmented inverse probability weighting. The estimator is unbiased for the complete-data treatment difference for any working models, and a martingale central limit theorem gives Wald intervals under repeated updating. We bound the variance lost by estimating the sampling rule. The bound is linear in the error of the fitted residual variance, quadratic when the optimal probabilities are not truncated, and relates the initial sample size to the learning rate of the models. Relative to designs using treatment assignment, the blinded design loses a term due to unequal residual variances in the two arms and a term due to the conditional treatment effect, which is second order near the null. In simulations, coverage was near nominal. Adaptive sampling was 12% to 34% more efficient than simple random sampling and needed 10% to 26% fewer measured outcomes for the same precision. In a resampling study of an antifungal trial, adaptive sampling reduced the sampling variance by 39%.

Varying coefficient (VC) regression models have become indispensable tools in spatial data analysis, providing unmatched flexibility in capturing complex, nonlinear relationships and spatially-varying effects of predictors on responses. Although hierarchical Bayesian approaches offer a rigorous probabilistic framework for uncertainty quantification in VC modeling, their practical application to large-scale spatial datasets remains severely hindered by the scalability limitations of Markov chain Monte Carlo (MCMC) algorithms. In response, the past decade has seen substantial advances in developing more efficient hierarchical Bayesian VC models, primarily through substituting traditional Gaussian processes (GP) with computationally efficient stochastic surrogates for estimating unknown coefficient functions. This article introduces a fundamentally different approach: the Geostatistical Variational Auto-Encoder (GeoVAE), a hierarchical deep generative framework built specifically for joint estimation of multiple spatially varying coefficient functions. GeoVAE departs from both classical GP-based Bayesian models and standard variational auto-encoders (VAEs) in two key ways. First, it constructs a coefficient-specific auto-encoder for each coefficient function, allowing each to capture its own spatial resolution and smoothness. Second, a hierarchical synthesis layer integrates information across these auto-encoders to a shared auto-encoder, explicitly modeling cross-coefficient dependencies arising from the shared spatial domain. This design enables GeoVAE to recover complex spatial patterns at computational cost that scales favorably with sample size, without MCMC or explicit GP covariance representations. We characterize the advantages and limitations of GeoVAE relative to hierarchical Bayesian spatial models, demonstrating its strong potential for large-scale spatial analysis.

Education agencies often summarize test score distributions by counting how many students scored in 3 to 5 different bins. The HETOP model transforms bin counts into estimated means and standard deviations, assuming that scores follow a normal distribution within each school or district. Past HETOP implementations ran slowly, taking 3–60 minutes, if they finished, when given bin counts for all 1,151 districts in Texas. Our new function, \code{fast\_hetop()} in the R package \pkg{binest}, runs all Texas districts in less than a second. Users can choose between maximum likelihood, empirical Bayes, sample-based or population-based estimates. Estimates are strongly correlated with true values, but have bias when the score distribution is skewed and scores are concentrated in the lowest or highest bin.

Crash counts on urban road networks exhibit spatial dependence driven by road connectivity rather than Euclidean proximity alone. We propose a Bayesian negative binomial model that represents spatial variation through a parsimonious combination of spectral components defined directly on road segments. The model combines the normalized random edge neighborhood Gaussian (RENeGe) operator with a continuous spike-and-slab prior, with mixture indicators marginalized out. Eigenvalue-dependent slab variances link spectral regularization to network-supported smoothness, while posterior slab-membership probabilities quantify support for individual components and uncertainty in effective spectral complexity. Segment length is used as an exposure offset, so fitted rates describe crash frequency per unit road length rather than traffic-adjusted risk. In a simulation study with 1,200 datasets and 4,800 Bayesian fits, the proposed model achieved the lowest mean spatial-field recovery error and highest mean test log predictive density under sparse RENeGe truth. Under alternative generating mechanisms, its predictive performance remained close to dense spectral models and exceeded that of a nonspatial model. Applications to Barcelona and central Bogota showed distinct posterior spectral patterns. Exploratory leave-one-segment-out comparisons found similar predictive accuracy for the proposed model and a low-rank edge CAR model, both outperforming a nonspatial negative binomial model, while a spectral BYM2 model achieved the highest exploratory predictive score. These results support probabilistic summaries of spectral complexity without a detectable predictive penalty relative to the corresponding CAR model.

In Bayesian inference problems with non-Gaussian observation noise, the posterior is only as accurate as the noise density, and gradient-based samplers need that density and its gradient evaluable pointwise, whether from an explicit expression or from code, and without an inner solve. We propose Copula Active Subspaces (CAS) to represent this noise density. A componentwise rank transform isolates the noise law's dependence in its copula, and a rank-$r$ reduction keeps only the directions along which that dependence varies. These directions are the leading eigenvectors of the copula score covariance $\boldsymbol{C} := \mathrm{Cov}_{\pi_{\boldsymbol{Z}}}(\nabla\log c^{Z})$, which is what makes the reduction a copula active subspace. Because $\boldsymbol{C}$ vanishes when the coordinates are independent, these are directions of dependence, which the covariance of the data need not identify. From this construction follow a Gaussian-reference KL divergence bound with the explicit constant $\tfrac{1}{2}$, minimized over all rank-$r$ reductions by exactly this eigenspace; a diagnostic for the error the reduction leaves behind, computable from the samples alone; and, from Hermite score matching, a reduced log-density and gradient in closed form, with the truncation orders and the Stage-2 regularization constants chosen on validation samples. The reduction replaces a $d$-dimensional density estimation problem by an $r$-dimensional one. On a $d=20$ noise law and a Bayesian inference problem with that noise, CAS lowers noise KL divergence more than fivefold and posterior KL divergence more than sevenfold against Gaussian-copula, product-of-marginals, and PCA-subspace baselines, and lowers noise KL divergence by factors of about $3.5$ and $2.7$ on two further $d=20$ examples.

Transporting an outcome relationship from a source population to a target population where the outcome is unobserved requires the conditional outcome law to be stable across populations. When it is not (outcome-model drift), the target mean is not point-identified and existing transportability estimators are biased. We study settings, common in cross-cohort educational and biomedical data integration, in which bridge outcomes are recorded in both populations and can be expected to drift alongside the target outcome. The bridge splits the drift into a component it can detect, fixed by a bridge-matching equation, and a residual orthogonal to the bridge that no observed quantity restricts. Because the target is a mean, the residual acts through a single scalar sensitivity parameter $\kappa$. We prove that the identified set has a closed-form irreducible core, whose width is governed by the residual bound and by how much of the outcome the bridge leaves unexplained, and that its center is first-order invariant to the residual. These results hold for any working exponential family, with point identification arising only as the $\kappa=0$ benchmark. We develop a debiased, drift-augmented estimator that is semiparametrically efficient at anchored sensitivity, doubly robust conditional on the bridge-identified drift, and yields rate-robust Imbens-Manski inference for the set. Unlike a bridge-blind sensitivity analysis, which must bound the entire drift channel, the proposed analysis bounds only its bridge-orthogonal part. Simulations confirm that the drift estimator is unbiased under co-drift, that the core width follows its closed form, and that the set fails visibly once the residual bound is exceeded. The motivating setting is the integration of two survey cohorts whose kindergarten mathematics outcome is available in only one.

In spatial observational studies, treatment assignment and outcomes often exhibit spatial dependence patterns, and treatment effects may vary across space and subpopulations due to both measured and unmeasured spatially structured confounders. Accounting for spatial dependence while estimating heterogeneous treatment effects (HTEs) is a central task in spatial causal inference. Causal Bayesian additive regression tree methods are popular nonparametric methods for modeling and estimating HTEs. Despite their flexibility and uncertainty quantification, the axis-aligned split rules often adopted in these models are not suitable for modeling spatial structures. We propose a spatial structure-aware Bayesian nonparametric method, called Graph-Split Bayesian Causal Forest (GSBCF), that integrates graph-split Bayesian additive regression trees (GS-BART) with the Bayesian causal forest propensity-score regression framework for spatial heterogeneous causal inference. Spatial confounding is accommodated through graph-guided split rules in modeling decision trees of the prognostic and HTE functions. We develop an efficient informed proposal sampling algorithm for posterior computation, enabling full Bayesian inference of the spatial conditional average treatment effect function. Simulations and a real data study demonstrate substantially improved estimation accuracy and uncertainty quantification over existing causal BART methods.

Bell and Cover showed that an investor who multiplies the initial unit of capital by an independent uniform random variable on $(0,2)$, and then uses the log-optimal portfolio, wins a head-to-head wealth comparison with probability at least one half against every independently randomized competitor. We explain very simply how this result transfers to testing by betting: for any composite null $\mathcal P$ and simple alternative $Q$, denoting $E^*$ as the corresponding numeraire e-variable, we show that $UE^*$ exceeds any other e-variable $E$ with probability at least half. Interestingly, we show that this competitive optimality result is actually equivalent to the numeraire inequality $\mathbb E_Q[E/E^*]\leq1$, and in general randomization only helps the numeraire and fails to improve the competitive advantage of an arbitrary e-variable. Under optional stopping with or without knowledge of $U$, we emphasize a key distinction between e-process validity and competitive optimality. We also show that competitive optimality comes at the price of expected log wealth and power: thresholding $UE^*$ at $1/\alpha$ has sharp size at most $\alpha/2$, but the factor of two actually disappears under optional stopping. Even after correcting for this factor of two, the test is dominated in conditional rejection probability by randomizing the testing threshold (randomized Markov's inequality). Thus, Bell-Cover randomization is optimal for a specific competitive objective, at the cost of others.

Objectives: There has been increased appreciation of the differences between marginal and conditional estimates and different types of effect measures regarding their applicability to different target populations. This issue of transportability is of concern in model-based cost-effectiveness analysis (CEA) when treatment effects from (international) trials are applied to (country-specific) baseline risk estimates. The objective of this paper is to create awareness regarding the issues that arise when using different types of treatment effect and baseline risk estimates in a model-based CEA to inform health technology assessment (HTA). Methods: We clarify collapsibility, marginal versus conditional estimation, and transportability; derive the ideal modeling approach implied by a marginal cost-effectiveness estimand; and characterize the issues of common modeling approaches, illustrated with a fictitious state-transition model. Results: An individual-level simulation that predicts outcomes from conditional inputs and averages them over the target population targets the marginal cost-effectiveness estimand. Cohort-model approaches that marginalize inputs early, use conditional inputs at mean covariate values, or combine a conditional effect with a marginal baseline (or vice versa) can misstate cost-effectiveness results even with correct-population inputs; inputs from the wrong population add further error. Conclusions: The most rigorous approach is an individual-level simulation that carries conditional inputs (baseline risk, treatment effect, prognostic effects and effect modifiers) and marginalizes late. Cohort approaches instead marginalize early, relying on aggregated inputs, and do not necessarily target the marginal cost-effectiveness estimand of interest for HTA. Model developers should document, for each input, whether it is marginal or conditional and its population.

Artificial intelligence models deployed as clinical decision support tools often suffer substantial performance degradation over time and across institutions due to covariate shift, concept drift, and cross-system heterogeneity. Fully retraining complex models is frequently infeasible, particularly in EHR settings where labeled outcome data are scarce and regulatory constraints limit model modification. We propose GUARD (Guided and Uncertainty-Aware Robustness to Domain shift), a unified statistical framework for updating an existing deployed model using limited labeled target data and multiple heterogeneous source populations while guarding against future distributional shifts. GUARD formulates robust multi-source transfer learning as a semi-supervised, adversarial optimization problem anchored at the current model, and employs double machine learning with cross-fitting to obtain debiased, efficient estimates that leverage abundant unlabeled target covariates. It then constructs an uncertainty-aware, guided group DRO estimator that combines source and target information while accounting for sampling variability in the source models. Our framework jointly provides principled robustness to domain shift, efficient semi-supervised estimation, and valid statistical inference for multi-source transfer of clinical prediction models. We demonstrate the utility of GUARD via extensive simulation experiments and a real world application to predicting future disease activity in rheumatoid arthritis from EHR data spanning more than a decade. GUARD recalibrates a pretrained model with only a small number of labeled records per year and remains accurate over multi-year horizons where target-only and standard transfer estimators degrade sharply.

Virtual simulators are widely used for studying complex physical phenomena, from particle collisions to rocket propulsion. Such "computer experiments" can be highly time-intensive, and a Bayesian surrogate model can be used for efficient emulation with reliable uncertainty quantification. To train accurate surrogates with a limited sample size $n$, recent work has explored the incorporation of monotonicity (or isotonicity) information, which can often be elicited from physical systems. In practical applications with many input variables, however, existing Bayesian isotonic models can face statistical and computational limitations, which may result in worse performance compared to models that do not incorporate isotonicity. We propose a new transformed additive isotonic model (TAIM), which aims to tame this "curse-of-dimensionality". TAIM makes use of a flexible transformed additive isotonic modeling framework, which leverages a data-estimated link transformation and a monotone basis model with spike-and-slab priors on basis weights. Prediction-wise, TAIM achieves (up to log factors) a posterior contraction rate of $O(n^{-1/3})$ when the true black-box function is in a transformed additive isotonic form with mild smoothness conditions. Such a rate does not depend on the input dimension $d$ for terms involving $n$, which softens the effect of dimensionality on posterior predictions. Computation-wise, TAIM allows for efficient posterior inference via a carefully designed Gibbs sampler, where each sampling iteration requires only linear work in $d$. We further present an extension of TAIM that can model potential deviations from transformed additivity. Numerical experiments and two applications show the effectiveness of TAIM for isotonic surrogate modeling with many input variables.

Price stability remains a pillar in monetary policy practices and carries a special importance within monetary unions. Mainstream economics tried to leverage price stability using price indices and several metrics to shed light on specific dynamics and optimal macroeconomic levels. The wide availability of data led researchers to consider the study of systems using Random Matrix Theory, based on inner correlation patterns. This aims to enhance the multivariate analysis by removing noisy patterns from the signal and improve data quality for further inferences. This work considers the collection of monthly inflation indices in the Eurozone as a system of prices to analyze its eigenvalues' statistical and asymptotic properties and uncover inner country-level insights. Results confirm the system cannot assumed to be randomly generated, and the data exhibit noise-dominated patterns, due to small and persistent variations at the country-level. The latter make the inter-country correlations more dynamic and the separation of the signal from the noise quiet difficult. Findings identified two countries as distorting inflation dynamics besides three other distinct, regional-based groups of countries. Variability sources might stem from economic episodes fueling inflation spikes in some countries, as well as methodological aspects used to ensure data quality and representativeness in the European Union. Despite being complex, the system demonstrates a certain stability, in terms of self-organization; while large monthly fluctuations cannot be considered as rare events, but part of the data-generating process.

We propose a novel joint sampling design for surveys involving two interrelated populations of parents and children in the presence of complex family structures, including shared parental responsibility and blended families. The design combines a graph-theoretic representation of co-parenting relationships with balanced sampling. This framework ensures that, in the first sampling phase, exactly one parent is selected from each co-responsibility pair while maintaining accurate inference for both the parent and child populations. A second sampling phase selects one child per sampled parent using unequal probabilities that produce consistent weighting and unbiased estimation. The proposed framework accommodates overlapping family relationships while balancing auxiliary information for both populations. Results based on the Luxembourg population registers show that the estimated totals are virtually identical to the known population totals for all balancing variables, and other auxiliary variables, demonstrating the excellent precision and efficiency of the proposed sampling design.

We present spsurv, an R package for semi-parametric time-to-event regression based on Bernstein-polynomial estimation of unknown baseline functions. The package provides a unified modelling interface for proportional hazards (PH), proportional odds (PO), and accelerated failure time (AFT) models for right-censored data, with either maximum likelihood or Bayesian estimation via Stan. Smooth baseline hazard, odds-function, or log-time structures are estimated without assuming a parametric baseline family, while retaining familiar hazard-ratio, odds-ratio, and time-ratio interpretations. We describe methodology, implementation, and syntax; evaluate finite-sample behaviour in a Monte Carlo study; and illustrate usage with oncology trials.

In this paper, we formulate a general differencing framework for variance estimation across a range of settings. We demonstrate that conventional difference-based noise variance estimators cannot achieve the desired bias-correcting power in nonparametric regression with repeated measurements. A new high-order bias-corrected differencing scheme, adapted to repeated measurements, is proposed by interlacing inter-group and intra-group differencing. The theoretical properties of the new sequences and estimators are studied. Our proposals are particularly efficient in finite samples and under high signal-to-noise ratio scenarios, where asymptotic convergence has not yet fully taken effect, due to their strong bias-correcting power.

Clustering problems increasingly involve complex objects observed over space, such as distributions, matrices, functions, images, or multivariate data, for which a scientifically meaningful dissimilarity between objects is often easier to specify and computationally more tractable than an object-response likelihood model. We propose DISCCO, a Bayesian framework for clustering spatially indexed complex objects using only a pairwise distance matrix and a spatial adjacency graph, with broad applicability and minimal user modeling requirements. The model combines a hierarchical distance-based likelihood accounting for within-cluster compactness and between-cluster separation with a random spatial graph partition prior, ensuring that posterior clusters are spatially contiguous. A key model feature is a set of node-specific frailty parameters that induce dependence among overlapping within-cluster distances and provide posterior summaries of object-level centrality or peripherality within each inferred cluster. We develop a partially collapsed Markov chain Monte Carlo algorithm for posterior inference. Simulations with distribution- and matrix-valued responses show that the proposed spatial distance-clustering framework improves region recovery relative to existing distance-clustering methods, while the frailty layer provides interpretable centrality summaries. Real applications to Houston Census Block Group racial-composition distributions and Western US county-level cancer mortality matrices illustrate how the method recovers interpretable contiguous clusters and frailty-based centrality maps.

In target trial emulation, the treatment strategies that patients follow are inferred from the treatments they actually receive in routine care. However, in most settings, the outcome may preclude the observation of planned treatment, giving rise to immortal time bias through misclassification of treatment strategy, while markers of treatment response may influence subsequent treatment decisions, giving rise to time-varying confounding. A key step toward unbiased treatment effect estimation is to partition follow-up into sufficiently short time intervals to unfold the feedback relationships involving treatment and represent the resulting causal relations with a directed acyclic graph. In this study, we present the possible within-interval causal orderings induced by this partitioning, discuss their causal implications, and assess their plausibility across clinical settings. For each causal ordering, we derive the corresponding g-formula. Using ancestral multi-world networks and simulations, we show that the standard cloning-censoring-weighting estimator is invalid when treatment affects the outcome within a time interval, and we propose a modified version of the method that restores its validity in this setting. Finally, we analyze the consequences of choosing time intervals that are either excessively wide or excessively narrow, thereby formally establishing the need for time partitioning and providing practical guidance for selecting an appropriate partition based on the clinical setting.

We propose a Bayesian nonparametric factor model that infers the number of factors, induces row-wise sparsity, and merges redundant dictionary elements via exact clustering. A Dirichlet process prior is placed on the columns of an overcomplete loading matrix and fully marginalized to an exact Pólya urn, avoiding stick-breaking and auxiliary variables. A spike-and-slab base measure allows entire factors to be exactly zero. The model uniquely combines exact zeros, exchangeability over columns, and exact merging within a single marginalized Dirichlet process, unlike CUSP, MGP, or the beta process. An exact Gibbs sampler with canonical relabeling and parallel C/MPI implementation is developed. We prove posterior contraction at rate $\sqrt{M s_0 \log n / n}$ for the covariance matrix, and in the fixed-dictionary setting obtain the minimax optimal rate $\sqrt{s_0 \log n / n}$ plus underfitting consistency; the overfitting direction is an open conjecture. The spike-and-slab is essential: without it the effective dimension scales as $pM$, yielding a slower rate. Simulations show the method is the only fully adaptive approach to recover the true rank, achieving the smallest covariance, loading, and signal-reconstruction errors, beating an oracle baseline. On van 't Veer breast cancer data ($n=97$, $p=1213$), the posterior concentrates on eight interpretable programmes; seven pass coherence and two pass Bonferroni-corrected Hallmark enrichment.

Background Decision curve analysis (DCA) evaluates prediction model usefulness using net benefit (NB). Population-level NB is often interpreted as a proxy for population-level expected utility, but this assumes comparability of utility across subgroups. We aimed to characterize when subgroup utility heterogeneity may invalidate population-level DCA conclusions. Methods We compared prediction-driven treatment with default strategies in populations containing subgroups with different utility values. We derived an inconsistency region: combinations of subgroup-specific {\Delta}NB values for which population-level NB and utility favor different strategies. We also developed a practical robustness framework. Results Opposite signs of subgroup-specific {\Delta}NB provide a warning signal for possible inconsistency. The inconsistency region is larger when subgroup sizes are more similar and subgroup-specific \(a-c\) values are more different, where \(a-c\) is the incremental utility of a true positive relative to a false negative. When \(a-c\) differs across subgroups, population-level NB combines quantities on different implicit utility scales and may conflict with population utility. A real-world case study illustrates the problem. Conclusions Using population-level NB as a proxy for population utility implicitly assumes homogeneous \(a-c\) across subgroups. Subgroup DCA and our framework can identify and assess when utility heterogeneity may invalidate population-level conclusions.

In observational causal inference, covariate balancing is widely used to reduce source-target covariate shift, but under weak overlap in high dimensions, stronger balance can induce concentrated weights and increase variance. Balance measures how well the target covariate distribution is represented, but does not by itself determine how reliably the counterfactual mean can be estimated. We develop risk-calibrated balancing for the average treatment effect on the treated, which applies ridge augmentation to any normalised base weights and selects its penalty using conditional prediction risk of the counterfactual mean. Under a random-effects predictive model, we derive an exact finite-sample decomposition of this risk into residual covariate imbalance and weight-induced variance. For design-independent base weights under proportional asymptotics, we characterise how limiting risk depends on source and target covariance geometry, population mean shift, and weight concentration. For covariate-adaptive base weights, we develop a uniformly consistent target-aware risk estimator whose minimiser attains vanishing scaled oracle excess risk. Simulations show that the high-dimensional risk predictions remain informative for adaptive balancing and that target-aware tuning generally reduces excess target risk. Empirical analyses of job-training and single-cell perturbation data show that risk-calibrated balancing generally improves on the corresponding base estimators, with larger gains under weaker overlap.

Building on the James–Stein approach to leading eigenvector estimation (Goldberg and Kercheval, Proc. Natl. Acad. Sci. USA 120, e2207046120, 2023), we develop a data-adaptive augmented James–Stein shrinkage framework for estimating leading eigenvectors and eigenspaces under a generalized spiked population model, in the high-dimensional regime where the dimension $p$ and sample size $n$ grow proportionally. For each spiked eigenvector, we construct an augmented target subspace that combines auxiliary information, either from domain knowledge or prior information, with the remaining sample spiked eigenvectors. This augmentation allows information shared across the sample spiked components to be exploited while retaining a fully data-driven shrinkage rule. We show that the resulting eigenvector estimator strictly improves upon standard PCA whenever the target subspace contains nonvanishing information about the population eigenvector, while asymptotically reverting to PCA when the target is uninformative. The individual estimators further yield a nested sequence of estimators for all leading spiked eigenspaces, with analogous dominance properties. The proposed estimator also strictly improves upon the existing HDLSS-motivated shrinkage estimator in the proportional high-dimensional regime. Simulation studies demonstrate substantial finite-sample gains and robustness to misspecification of the number of spikes.

Linear mixed-effect panel data models are widely used in longitudinal biomedical, environmental, and social studies, but are often sensitive to data contamination and numerous covariates. To address these challenges, we propose a robust variable selection approach based on a least squares approximation (LSA) of the density power divergence (DPD) objective function combined with the Adaptive LASSO penalty. The LSA converts the nonlinear DPD objective into a computationally efficient quadratic approximation while preserving the robustness of DPD estimation. Under suitable regularity conditions, the proposed DPD Adaptive LASSO-LSA estimator is shown to possess oracle properties, including selection consistency and asymptotic normality. Simulation studies demonstrate improved robustness, better sparsity recovery, and significantly enhanced computational efficiency than penalized likelihood methods, while maintaining performance comparable to existing DPD-based approaches. An application to the SWAN bone mineral density dataset illustrates the practical relevance of the proposed methodology for robust estimation and reliable identification of important covariates.

A non-significant coefficient in a dynamic panel model need not imply the absence of a relationship. It may instead reflect heterogeneous effects averaged toward zero, reciprocal dynamics overlooked by a recursive specification, or relationships masked by the omission of correlated covariates. Standard linear estimators cannot distinguish among these possibilities. We develop an inferential workflow for political panel data that resolves this ambiguity by combining flexible autoregressive estimation, forecast-necessity testing, functional characterization, and same-data linear benchmarking. The workflow first identifies relationships required for out-of-sample prediction, then characterizes their functional form across political contexts, and finally, distinguishes differences arising from estimator flexibility from those due to model specification. Applied to the causal sequence model of democratization on the V-Dem panel of 113 countries, the workflow reproduces the model's central finding - the protective belt of civil society, the rule of law, and institutionalized parties - while recovering reciprocal relationships from democracy to its institutional supports that a linear model cannot detect. Most importantly, three weak published direct effects, of which two are null, and one is marginally significant, receive three different diagnoses: one dissolves under the full specification, one reflects heterogeneous effects averaged toward zero, and one was masked by the reduced variable set. The workflow corrects the published record in both directions, removing one relationship and recovering two. More broadly, the workflow provides a framework for evaluating dynamic political theories under a model class capable of representing nonlinear and reciprocal mechanisms while preserving relationship-level interpretation and explicit inferential standards.