Papers

Multimodal time-to-event prediction often requires integrating sensitive data distributed across multiple parties, making centralized model training impractical due to privacy constraints. At the same time, most existing multimodal survival models produce single deterministic predictions without indicating how confident the model is in its estimates, which can limit their reliability in real-world decision making. To address these challenges, we propose BVFLMSP, a Bayesian Vertical Federated Learning (VFL) framework for multimodal time-to-event analysis based on a Split Neural Network architecture. In BVFLMSP, each client independently models a specific data modality using a Bayesian neural network, while a central server aggregates intermediate representations to perform survival risk prediction. To enhance privacy, we integrate differential privacy mechanisms by perturbing client side representations before transmission, providing formal privacy guarantees against information leakage during federated training. We first evaluate our Bayesian multimodal survival model against widely used single modality survival baselines and the centralized multimodal baseline MultiSurv. Across multimodal settings, the proposed method shows consistent improvements in discrimination performance, with up to 0.02 higher C-index compared to MultiSurv. We then compare federated and centralized learning under varying privacy budgets across different modality combinations, highlighting the tradeoff between predictive performance and privacy. Experimental results show that BVFLMSP effectively includes multimodal data, improves survival prediction over existing baselines, and remains robust under strict privacy constraints while providing uncertainty estimates.

Generative AI (GAI) reveals an irreducible human core at the center of data science: advances in GAI should sharpen, rather than diminish, the focus on human reasoning in data science education. GAI can now execute many routine data science workflows, including cleaning, summarizing, visualizing, modeling, and drafting reports. Yet the competencies that matter most remain irreducibly human: problem formulation, measurement and design, causal identification, statistical and computational reasoning, ethics and accountability, and sensemaking. Drawing on Donoho's Greater Data Science framework, Nolan and Temple Lang's vision of computational literacy, and the McLuhan-Culkin insight that we shape our tools and thereafter our tools shape us, this paper traces the emergence of data science through three converging lineages: Tukey's intellectual vision of data analysis as a science, the commercial logic of surveillance capitalism that created industrial demand for data scientists, and the academic programs that followed. Mapping GAI's impact onto Donoho's six divisions of Greater Data Science shows that computing with data (GDS3) has been substantially automated, while data gathering, preparation, and exploration (GDS1) and science about data science (GDS6) still require essential human input. The educational implication is that data science curricula should focus on this human core while teaching students how to contribute effectively within iterative prompt-output-prompt cycles using retrieval-augmented generation, and that learning outcomes and assessments should explicitly evaluate reasoning and judgment.

We consider a stopping problem and its application to the decision-making process regarding the optimal timing of organ transplantation for individual patients. At each decision period, the patient state is inspected and a decision is made whether to transplant. If the organ is transplanted, the process terminates; otherwise, the process continues until a transplant happens or the patient dies. Under suitable conditions, we show that there exists a control limit optimal policy. We propose a smoothed perturbation analysis (SPA) estimator for the gradient of the total expected discounted reward with respect to the control limit. Moreover, we show that the SPA estimator is asymptotically unbiased.

In particle physics, as in many areas of science, parameter inference relies on simulations to bridge the gap between theory and experiment. Recent developments in simulation-based inference have boosted the sensitivity of analyses; however, biases induced by simulation-data mismodeling can be difficult to control within standard inference pipelines. In this work, we propose a Template-Adapted Mixture Model to confront this problem in the context of signal fraction estimation: inferring the population proportion of signal in a mixed sample of signal and background, both of which follow arbitrarily complex distributions. We harness many biased simulations to perform data-driven estimates of each process distribution in the signal region, substantially reducing the bias on the signal fraction due to the domain shift between simulation and reality. We explore different methodological choices, including model selection, feature representation, and statistical method, and apply them to a Gaussian toy example and to a semi-realistic di-Higgs measurement. We find that the presented methods successfully leverage the biased simulations to provide estimates with well-calibrated uncertainties.

Probabilistic forecasts must sum to unity and cannot express “I don't know.” Possibility theory relaxes this constraint: a subnormal distribution explicitly measures how much of the plausibility budget remains unassigned, ignorance signal that probability cannot represent. This paper develops a verification framework for such forecasts, centred on a five-number scorecard that separately diagnoses whether the forecast pointed at the right outcome (depth-of-truth), how sharply (diffuseness, support margin), how confidently (ignorance), and how dominantly (conditional necessity). A possibility-to-probability conversion preserves ignorance for familiar frequency-based scoring; categorical threshold scores (POD, FAR, CSI, etc.) connect to operational practice. Together, these three complementary facets – possibilistic, probabilistic, and categorical – expose failure modes invisible to any single metric. Storm Prediction Center convective outlook categories serve as the running example throughout; a synthetic reforecast demonstrates diagnostic visualisations and scorecard interpretation. Ignorance is better expressed than repressed.

We study how sampling geometry contributes to uncertainty in modeling spatial geophysical observations as sampled random fields characterized by stationary, isotropic, parametric covariance functions. We incorporate the signature of discrete spatial sampling patterns into an asymptotically unbiased spectral maximum-likelihood estimation method along with analytical uncertainty calculation. We illustrate the broad applicability of our modeling through synthetic and real data examples with sampling patterns that include irregularly bounded contiguous region(s) of interest, structured sweeps of instrumental measurements, and missing observations dispersed across the domain of a field, which spur behaviors from the estimator. We find through asymptotic studies that allocating samples following a growing-domain strategy rather than a densifying, infill scheme best reduces estimator bias and (co)variance, whether the field has been sampled regularly or not. As our modeling assumptions, too, shape how (well) an observed random field can be characterized, we study the effect of covariance parameters assumed a priori. We demonstrate the desirable behavior of the general Matern class and show how to interrogate goodness-of-fit criteria to detect departures from the null hypothesis of Gaussianity, stationarity, and isotropy.

Samplets are data adapted multiresolution analyses of localized discrete signed measures. They can be constructed on scattered data sites in arbitrary dimension such that they exhibit vanishing moments with respect to any prescribed set of primitives. We consider the samplet construction in a probabilistic framework and show that, if choosing polynomials as primitives, the resulting samplet basis converges to signed measures with broken polynomial densities in the infinite data limit. These densities amount to multiwavelets with respect to a hierarchical partition of the region containing the data sites. As a byproduct, we therefore obtain a construction of general multiwavelets that allows for a flexible prescription of vanishing moments going beyond tensor product constructions. For congruent partitions we particularly recover classical multiwavelets with scale- and partition- independent filter coefficients. The theoretical findings are complemented by numerical experiments that illustrate the convergence results in case of random as well as low-discrepancy data sites.

Variational inference (VI) has become a widely used approach for scalable Bayesian inference, but its performance strongly depends on the flexibility of the chosen variational family. In this work, we propose a novel variational family that combines wavelet-based representations for marginal posterior densities with copula functions to model dependence structures. The marginal distributions are constructed using coefficients from the discrete wavelet transform, providing a flexible and adaptive framework capable of capturing complex features such as asymmetry. The joint distribution is then obtained through a copula, allowing for explicit modeling of dependence among parameters, including both independence and Gaussian copula structures. We develop an efficient estimation procedure based on Monte Carlo approximations of the evidence lower bound (ELBO) and automatic differentiation, enabling scalable optimization using gradient-based methods. Through extensive simulation studies, including logistic regression, sparse linear models, and hierarchical models, we demonstrate that the proposed approach achieves posterior mean estimates comparable to Markov chain Monte Carlo (MCMC) methods, while providing improved uncertainty quantification relative to standard variational approaches. Applications to hierarchical logistic regression and Bayesian conditional transformation models further illustrate the practical advantages of the method in complex, high dimensional settings. The proposed wavelet copula variational family offers a flexible and computationally efficient alternative for Bayesian inference.

Many Bayesian inference problems involve high-dimensional models where the performance of standard importance sampling (IS) methods often degrades rapidly as the dimensionality increases. Classical analyses of IS typically rely on the assumption that observations are arbitrary but fixed (i.e., deterministic), thereby neglecting the probabilistic structure that the Bayesian model induces on the data. In this paper, we adopt the perspective that observations are themselves random variables whose distribution is governed by the underlying model. Within this probabilistic framework, we identify a model-dependent function, referred to as the link function, which connects the fixed- and random-observation formulations. We provide a characterization of the $L^2$ Monte Carlo estimation error: specifically, we show that the $L^2$ error bounds are finite and converge at the standard Monte Carlo rate $O(N^{-1/2})$, for arbitrarily large dimension, if and only if the link function is Bochner integrable. This result reveals the fundamental quantity controlling the approximation error and establishes a mechanism to manage the dependence on the model state dimension. Consequently, our approach provides a principled way to alleviate the challenges of high dimensionality, offering insights that transcend worst-case analyses dominant in the existing literature. Finally, we derive explicit analytical examples of the dimensional scaling of the associated errors for several model classes, including linear-Gaussian systems and models with bounded observation functions.

Natural and anthropogenic disturbances are impacting the health of forests worldwide. Monitoring forest disturbances at scale is important to inform conservation efforts. Here, we present a scalable approach for country-wide mapping of forest greenness anomalies at the 10 m resolution of Sentinel-2. Using relevant ecological and topographical context and an established representation of the vegetation cycle, we learn a predictive quantile model of the normalised difference vegetation index (NDVI) derived from Sentinel-2 data. The resulting expected seasonal cycles are used to detect NDVI anomalies across Switzerland between April 2017 and August 2025. Goodness-of-fit evaluations show that the conditional model explains 65% of the observed variations in the median seasonal cycle. The model consistently benefits from the local context information, particularly during the green-up period. The approach produces coherent spatial anomaly patterns and enables country-wide quantification of forest browning. Case studies with independent reference data from known events illustrate that the model reliably detects different types of disturbances.

Almost every numerical task can be cast as extrapolation with respect to the fidelity or tolerance parameters of a consistent numerical method. This perspective enables probabilistic uncertainty quantification and optimal experimental design functionality to be deployed, and also unlocks the potential for the convergence of numerical methods to be accelerated. Recent work established Probabilistic Richardson Extrapolation as a proof-of-concept, demonstrating how parallel multi-fidelity simulation can be used to accelerate simulation from a whole-heart model. However, the number of simulations was required to increase super-exponentially in $d$, the number of tolerance parameters appearing in the numerical method. This paper develops a refined notion of 'extrapolation dimension', drastically reducing this simulation requirement when multiple tolerance parameters feature in the numerical method. Sparsity-exploiting methodology is developed that is simultaneously simpler and more powerful compared to earlier work, and this is accompanied by sharp theoretical guarantees and substantial empirical support.

Demographic parity (DP) is a widely studied fairness criterion in regression, enforcing independence between the predictions and sensitive attributes. However, constraining the entire distribution can degrade predictive accuracy and may be unnecessary for many applications, where fairness concerns are localized to specific regions of the distribution. To overcome this issue, we propose a new framework for regression under DP that focuses on the tails of target distribution across sensitive groups. Our methodology builds on optimal transport theory. By enforcing fairness constraints only over targeted regions of the distribution, our approach enables more nuanced and context-sensitive interventions. Leveraging recent advances, we develop an interpretable and flexible algorithm that leverages the geometric structure of optimal transport. We provide theoretical guarantees, including risk bounds and fairness properties, and validate the method through experiments in regression settings.

We study a random model of deep multi-head self-attention in which the weights are resampled independently across layers and heads, as at initialization of training. Viewing depth as a time variable, the residual stream defines a discrete-time interacting particle system on the unit sphere. We prove that, under suitable joint scalings of the depth, the residual step size, and the number of heads, this dynamics admits a nontrivial homogenized limit. Depending on the scaling, the limit is either deterministic or stochastic with common noise; in the mean-field regime, the latter leads to a stochastic nonlinear Fokker–Planck equation for the conditional law of a representative token. In the Gaussian setting, the limiting drift vanishes, making the homogenized dynamics explicit enough to study representation collapse. This yields quantitative trade-offs between dimension, context length, and temperature, and identifies regimes in which clustering can be mitigated.

Estimating optimal individualized treatment rules (ITRs) via outcome weighted learning (OWL) often relies on observed rewards that are noisy or optimistic proxies for the true latent utility. Ignoring this reward uncertainty leads to the selection of policies with inflated apparent performance, yet existing OWL frameworks lack the finite-sample guarantees required to systematically embed such uncertainty into the learning objective. To address this issue, we propose PAC-Bayesian Reward-Certified Outcome Weighted Learning (PROWL). Given a one-sided uncertainty certificate, PROWL constructs a conservative reward and a strictly policy-dependent lower bound on the true expected value. Theoretically, we prove an exact certified reduction that transforms robust policy learning into a unified, split-free cost-sensitive classification task. This formulation enables the derivation of a nonasymptotic PAC-Bayes lower bound for randomized ITRs, where we establish that the optimal posterior maximizing this bound is exactly characterized by a general Bayes update. To overcome the learning-rate selection problem inherent in generalized Bayesian inference, we introduce a fully automated, bounds-based calibration procedure, coupled with a Fisher-consistent certified hinge surrogate for efficient optimization. Our experiments demonstrate that PROWL achieves improvements in estimating robust, high-value treatment regimes under severe reward uncertainty compared to standard methods for ITR estimation.

This paper addresses the problem of clustering measurement vectors that are heteroscedastic in that they can have different covariance matrices. From the assumption that the measurement vectors within a given cluster are Gaussian distributed with possibly different and unknown covariant matrices around the cluster centroid, we introduce a novel cost function to estimate the centroids. The zeros of the gradient of this cost function turn out to be the fixed-points of a certain function. As such, the approach generalizes the methodology employed to derive the existing Mean-Shift algorithm. But as a main and novel theoretical result compared to Mean-Shift, this paper shows that the sole fixed-points of the identified function tend to be the cluster centroids if both the number of measurements per cluster and the distances between centroids are large enough. As a second contribution, this paper introduces the Wald kernel for clustering. This kernel is defined as the p-value of the Wald hypothesis test for testing the mean of a Gaussian. As such, the Wald kernel measures the plausibility that a measurement vector belongs to a given cluster and it scales better with the dimension of the measurement vectors than the usual Gaussian kernel. Finally, the proposed theoretical framework allows us to derive a new clustering algorithm called CENTRE-X that works by estimating the fixed-points of the identified function. As Mean-Shift, CENTRE-X requires no prior knowledge of the number of clusters. It relies on a Wald hypothesis test to significantly reduce the number of fixed points to calculate compared to the Mean-Shift algorithm, thus resulting in a clear gain in complexity. Simulation results on synthetic and real data sets show that CENTRE-X has comparable or better performance than standard clustering algorithms K-means and Mean-Shift, even when the covariance matrices are not perfectly known.

Prognostic covariate adjustment (PROCOVA) is a two-sample two-stage estimation method for covariate adjustment in randomized controlled trials. In the first stage, a prognostic score, defined as the conditional expectation of an outcome given covariates under the control treatment, is estimated using historical data. In the second stage, analysis of covariance with the estimated prognostic score and treatment assignment as explanatory variables is performed, and the average treatment effect is estimated. Although the prognostic score is estimated in this procedure, the variance estimator, which treats the prognostic score as known, has been used. Furthermore, the difference in the asymptotic variance between cases where the prognostic score is known versus where it is estimated has not been previously clarified. In this study, we derived these two asymptotic variances and showed that they are equal. This result also holds when the prognostic score is estimated using machine learning with $L_2$ consistency. We also constructed two variance estimators: one that treats the prognostic score as known, and another that accounts for its estimation, and compared their performance through simulation studies and data applications. For PROCOVA, since both variance estimators are asymptotically valid, it is generally recommended to use a variance estimator that treats the prognostic score as known, as it is simpler to derive and implement. When historical data is small, a variance estimator that explicitly accounts for prognostic score estimation is recommended if conservative inference is preferred.

Modern neural networks of the transformer family require the practitioner to decide, before training begins, how many attention heads to use, how deep the network should be, and how wide each component should be. These decisions are made without knowledge of the task, producing architectures that are systematically larger than necessary: empirical studies find that a substantial fraction of heads and layers can be removed after training without performance loss. This paper introduces DDCL-INCRT, an architecture that determines its own structure during training. Two complementary ideas are combined. The first, DDCL (Deep Dual Competitive Learning), replaces the feedforward block with a dictionary of learned prototype vectors representing the most informative directions in the data. The prototypes spread apart automatically, driven by the training objective, without explicit regularisation. The second, INCRT (Incremental Transformer), controls the number of heads: starting from one, it adds a new head only when the directional information uncaptured by existing heads exceeds a threshold. The main theoretical finding is that these two mechanisms reinforce each other: each new head amplifies prototype separation, which in turn raises the signal triggering the next addition. At convergence, the network self-organises into a hierarchy of heads ordered by representational granularity. This hierarchical structure is proved to be unique and minimal, the smallest architecture sufficient for the task, under the stated conditions. Formal guarantees of stability, convergence, and pruning safety are established throughout. The architecture is not something one designs. It is something one derives.

We study the prophet inequality, a fundamental problem in online decision-making and optimal stopping, in a practical setting where rewards are observed only through noisy realizations and reward distributions are unknown. At each stage, the decision-maker receives a noisy reward whose true value follows a linear model with an unknown latent parameter, and observes a feature vector drawn from a distribution. To address this challenge, we propose algorithms that integrate learning and decision-making via lower-confidence-bound (LCB) thresholding. In the i.i.d.\ setting, we establish that both an Explore-then-Decide strategy and an $\varepsilon$-Greedy variant achieve the sharp competitive ratio of $1 - 1/e$, under a mild condition on the optimal value. For non-identical distributions, we show that a competitive ratio of $1/2$ can be guaranteed against a relaxed benchmark. Moreover, with limited window access to past rewards, the tight ratio of $1/2$ against the optimal benchmark is achieved.

In this paper, we propose a method to analyze correlations in pandemic-related data across different geographical regions, relying on the analysis of correlations for non-stationary time series, which are typical of pandemic data. Unlike traditional epidemiological approaches focused on medical and modeling perspectives during a pandemic, our method emphasizes post-pandemic analysis to assess how societal responses; such as lockdowns, travel restrictions, mobility patterns, and vaccination campaigns, manifest in the collective behavior of regions. These insights can inform future public health strategies and enhance understanding of the complex dynamics underlying pandemic spread and control.

Effectively modeling irregularly sampled longitudinal data is essential for understanding disease progression and improving risk prediction. We propose a two-view mixture model that integrates static baseline covariates and longitudinal biomarker trajectories within a unified probabilistic clustering framework. Temporal patterns are modeled using Neural Ordinary Differential Equations. Model training uses an EM algorithm with a sparsity-inducing log-penalty for interpretable subgroup discovery. Application of the model to an Irish cohort of ANCA-associated vasculitis patients reveals subgroups with heterogeneous serum creatinine trajectories and variation in end-stage kidney disease outcomes.

The increasing availability of large-scale global datasets has generated a demand for scalable spatial prediction methods defined on spherical domains. Classical spatial models that rely on Euclidean distance representations are inappropriate for spherical data because planar projections distort geodesic distances and spatial neighborhood structures, while traditional kriging-based prediction methods are often computationally prohibitive for massive datasets. To address these challenges, we propose a Spherical DeepKriging framework for spatial prediction on $\mathbb{S}^2$. The proposed approach constructs a flexible prediction model by integrating thin-plate spline (TPS) basis functions defined intrinsically on the sphere. Simulation studies and real data analyses are presented to demonstrate the superior predictive performance of the proposed method.

We study the normal mean inference problem, which involves simultaneous testing of the means of many normal distributions. This problem has been extensively studied within the empirical Bayes (EB) framework. However, the reliability of most EB methods heavily depends on two key conditions: (i) the prior distribution is correctly specified, and (ii) it can be accurately estimated. In practice, both conditions are difficult to satisfy, and it is often unclear whether they hold in a given application. To overcome these limitations, we propose a new algorithm, called COIN (COnformal Inference for Normal mean inference problem). Unlike traditional empirical Bayes approaches, COIN produces decision rules whose validity does not depend on the correct specification or accurate estimation of the prior. We theoretically prove that COIN asymptotically controls the false discovery rate at the nominal level, even in the presence of prior misspecification or estimation errors. Since the COIN algorithm requires an external training dataset to estimate the prior distribution and conformity score function, we introduce two data-splitting strategies – sample-splitting and feature-splitting – for the case where such external data are unavailable. We provide theoretical guarantees for the data-splitting strategies and demonstrate their effectiveness through extensive numerical studies and three real data examples.

Statistical methods for testing aggregate rare-variant genetic associations are typically based on either burden or dispersion tests (or a combination of the two). These methods lack statistical power in the presence of diverse genetic architectures. Moreover, few aggregate rare-variant association methods have been developed specifically for survival data. To address these issues, we propose data-adaptive gene- and pathway-based association tests based on Schoenfeld residuals in Cox proportional hazards models for association studies between an aggregate of rare-variants and survival outcomes. Our methods improve statistical power while maintaining flexibility across various genetic effect sizes and directions. We develop an efficient R package that enables fast computation and supports data simulation as well as gene- and pathway-level testing. Applying our approach to late bladder toxicity following radiotherapy for non-metastatic prostate cancer, we identify biologically relevant genes and pathways, replicate known signals, and capture additional associations. Our method provides a powerful, adaptive framework for survival-based genetic association studies of rare-variants. Keywords: aSPU, time-to-event outcomes, rare-variant associations, Cox regression, Schoenfeld residuals

Optimization over the space of probability measures endowed with the Wasserstein-2 geometry is central to modern machine learning and mean-field modeling. However, traditional methods relying on full Wasserstein gradients often suffer from high computational overhead in high-dimensional or ill-conditioned settings. We propose a randomized coordinate descent framework specifically designed for the Wasserstein manifold, introducing both Random Wasserstein Coordinate Descent (RWCD) and Random Wasserstein Coordinate Proximal{-Gradient} (RWCP) for composite objectives. By exploiting coordinate-wise structures, our methods adapt to anisotropic objective landscapes where full-gradient approaches typically struggle. We provide a rigorous convergence analysis across various landscape geometries, establishing guarantees under non-convex, Polyak-Łojasiewicz, and geodesically convex conditions. Our theoretical results mirror the classic convergence properties found in Euclidean space, revealing a compelling symmetry between coordinate descent on vectors and on probability measures. The developed techniques are inherently adaptive to the Wasserstein geometry and offer a robust analytical template that can be extended to other optimization solvers within the space of measures. Numerical experiments on ill-conditioned energies demonstrate that our framework offers significant speedups over conventional full-gradient methods.

In spatio-temporal analysis, we often record data at specific time intervals but with varying spatial locations between these timepoints. We propose a conditional model to analyze such spatio-temporal data that accommodates the dependencies alongside second-order stationary explanatory variables, which may be infinite-dimensional and accommodate spatio-temporal covariates. Because of the absence of a mixing-type dependence condition in this case, which is typically required by the existing studies, we consider a weaker polynomially decaying moment contraction (PMC) condition on the covariates. In this paper, we obtain nonparametric point estimates of the mean and covariate functions of such a regression model, which we then show to be statistically consistent. We also obtain a simultaneous confidence interval of the mean function using the central limit theorem for the proposed estimator. Such simultaneous inference tools can be used to test for certain specifications of the mean function. Some simulation studies and two real-data analyses have been illustrated to corroborate the findings.

Brownian motion and fractional Brownian motion have been widely applied in statistical modeling in finance, telecommunication, network traffic, neuroscience, physics, and other fields. More realistic models for real time series data, such as multifractional processes, generalize these classical models by allowing their regularity to vary over time. A new class of Gaussian Haar-based multifractional processes, which utilizes the Haar wavelet series representation, was recently introduced. It significantly extends the range of available models by incorporating more general classes of Hurst functions. The Rmfrac package was developed to simulate multifractional time series. The package also comprises several functions for the analysis and visualization of time series. It includes the estimation of the Hurst function and local fractal dimension, clustering realizations and computing various geometric statistics of these time series. The package also offers a Shiny application to visualize simulation and estimation results. The article presents an overview of the Rmfrac package and exemplifies its main functionalities.

The Wallace–Freeman estimator is a classical minimum message length estimator whose relationship with likelihood-based asymptotic theory has not been fully developed. We show that, in regular parametric models, the Wallace–Freeman criterion is equivalent, up to constants, to a penalised likelihood criterion with penalty weight \(n^{-1}\). This representation places the estimator within the standard theory of penalised M-estimation and yields existence, consistency, an asymptotic linear expansion, and asymptotic normality under regularity conditions. We further derive the first-order difference between the Wallace–Freeman estimator and the maximum likelihood estimator, showing that it is an explicit \(O(n^{-1})\) shift determined by the gradient of the Wallace–Freeman penalty. Combining this expansion with the Cox–Snell formula gives a first-order bias expansion for the Wallace–Freeman estimator. The result clarifies its relationship with maximum likelihood, Jeffreys-prior penalisation, and Firth-type bias reduction. We illustrate the theory for the Weibull model, where the penalty modifies the leading bias of the maximum likelihood estimator of the shape parameter.

We propose a novel data harmonization approach known as Tensor-ComBat (TC) for structural neuroimaging data. Tensor-Combat is a novel spatially aware harmonization method that aims to estimate and remove unwanted technical variation between voxel-level images from different study sites or scanners. Tensor-ComBat uses a Bayesian tensor response regression (BTRR) model to estimate spatially distributed scanner effects via a low-rank PARAFAC decomposition on the model coefficients, and subsequently removes these scanner effects via a post-hoc ComBat harmonization step. Unlike the classical ComBat method that treats the ROIs or voxels in the image as interchangeable, the Tensor-ComBat approach incorporates the information about the spatial configurations of imaging voxels when estimating the model parameters, resulting an improved harmonization pipeline. The proposed Tensor-ComBat method is fit using a Markov chain Monte Carlo (MCMC) that is not only able to produce point estimates, but also able to quantify uncertainty associated with these estimates. Due to the fully Bayesian implementation of Tensor-ComBat, it is able to perform inference related to significant spatially distributed scanner effects and/or biological effects that is not straightforward under existing methods. The methods is applied to over 2100 T1w-MRI scans in ADNI-1 that illustrates greater scanner effect removal, improved biological prediction, and superior reproducibility compared to state-of-the-art approaches such as ComBat.

Conformal risk control (CRC) provides distribution-free guarantees for controlling the expected loss at a user-specified level. Existing theory typically assumes that the loss decreases monotonically with a tuning parameter that governs the size of the prediction set. However, this assumption is often violated in practice, where losses may behave non-monotonically due to competing objectives such as coverage and efficiency. In this paper, we study CRC under non-monotone loss functions when the tuning parameter is selected from a finite grid, a setting commonly arising in thresholding and discretized decision rules. Revisiting a known counterexample, we show that the validity of CRC without monotonicity depends critically on the relationship between the calibration sample size and the grid resolution. In particular, reliable risk control can still be achieved when the calibration sample is sufficiently large relative to the grid size. We establish a finite-sample guarantee for bounded losses over a grid of size $m$, showing that the excess risk above the target level $\alpha$ scales on the order of $\sqrt{\log(m)/n}$, where $n$ is the calibration sample size. A matching lower bound demonstrates that this rate is minimax optimal. We also derive refined guarantees under additional structural conditions, including Lipschitz continuity and monotonicity, and extend the analysis to settings with distribution shift via importance weighting. Numerical experiments on synthetic multilabel classification and real object detection data illustrate the practical implications of non-monotonicity. Methods that explicitly account for finite-sample uncertainty achieve more stable risk control than approaches based on monotonicity transformations, while maintaining competitive prediction set sizes.

Causal mediation analysis provides techniques for defining and estimating effects that may be endowed with mechanistic interpretations. With many scientific investigations seeking to address mechanistic questions, causal direct and indirect effects have garnered much attention. The natural direct and indirect effects, the most widely used among such causal mediation estimands, are limited in their practical utility due to stringent identification requirements. Accordingly, considerable effort has been invested in developing alternative direct and indirect effect decompositions with relaxed identification requirements. Such efforts often yield effect definitions with nuanced and challenging interpretations. By contrast, relatively limited attention has been paid to relaxing the identification assumptions of the natural direct and indirect effects. Motivated by a secondary aim of a recent non-randomized vaccine prospective cohort study (NCT05168813), we present a set of relaxed conditions under which the natural direct effect is identifiable in spite of unobserved baseline confounding of the exposure-mediator pathway; we use this result to investigate the effect mediated by putative immune correlates of protection. Relaxing the commonly used but restrictive cross-world counterfactual independence assumption, we discuss strategies for evaluating the natural direct effect in non-randomized settings that arise in the analysis of vaccine studies. We revisit prior studies of semi-parametric efficiency theory to demonstrate the construction of flexible, multiply robust estimators of the natural direct effect and discuss efficient estimation strategies that do not place restrictive modeling assumptions on nuisance functions.