Papers

Stein's method constructs computable discrepancies between a target distribution and a candidate distribution without requiring the target distribution's normalizing constant. These discrepancies support goodness-of-fit tests for model assessment as well as sampling tools for empirical approximation. The R package steinsampling provides the first unified R workflow for applying score-based Stein methods to kernel goodness-of-fit testing of independent or serially dependent observations, point transport, greedy point construction, and sample compression. High-level functions carry out each task in a single call, while the kernel, calibration, optimization, and transition components are provided separately so that users can replace any one of them. A single score and kernel setup can therefore be reused across sampling and testing, making these methods easier to reproduce, compare, and extend.

We propose an online testing framework for detecting distributional changes in serially dependent data with values in a separable metric space. Based on two-sample $U$-statistics, the framework encompasses sequential analogs of energy distance and maximum mean discrepancy (MMD) procedures while accommodating temporal dependence. We establish asymptotic theory for finite and open-ended monitoring horizons that characterizes the full asymptotic run-length distribution under $H_0$ and yields asymptotic false-alarm control. We further establish new spectral approximation results for kernel matrices formed from serially dependent observations, and use them to construct a feasible Monte Carlo calibration procedure. Our flexible window construction encompasses classical, Page-type, and full-scan historical-baseline monitoring and can achieve short detection delays for both early and late changepoints, without requiring sub-Gaussianity or high-order moments of the raw observations. Simulations show reliable false-alarm control across linear, nonlinear, high-dimensional, and functional time-series models and further demonstrate that, over a broad range of alternatives and changepoint locations, the proposed method can achieve substantially shorter delays than recent procedures designed specifically for rapid detection. Applications to foreign exchange rates, electricity-market curves, and daily air transportation networks illustrate the methodology across scalar, functional, and network-valued time series.

Temperature is the most important meteorological driver of electricity demand, but its causal effect has not been estimated under the time-varying confounding that characterizes weather processes. Current temperature is associated with weather conditions, including precipitation, snow cover, and cloud cover, that also influence electricity demand, and past temperature shapes future weather, producing treatment-confounder feedback under which standard regression adjustment is biased, due to conditioning on a time-varying confounder. We estimate the causal effect of temperature on daily Ontario electricity demand from 2018 to 2019 using a marginal structural model with inverse-probability weighting, formulated for a single observed time series rather than a panel of independent subjects, and we extend the longitudinal outcome-adaptive LASSO and adaptive fused LASSO, previously developed for binary treatments, to a continuous treatment using density-ratio weights and a weighted-covariance balance criterion. In a Monte Carlo study, unadjusted regression is biased toward the null with coverage of $0.12$ to $0.15$, whereas the stabilized outcome-adaptive estimators are nearly unbiased with coverage near $0.92$, and the cumulative three-day estimators are approximately unbiased but less efficient. In the Ontario data all estimators identify a positive and highly significant quadratic temperature effect; the unadjusted estimate is $9.34$, the stabilized and outcome-adaptive single-lag estimators give $9.8$ to $10.1$, and the cumulative estimators give smaller values that coincide with a sharp fall in effective sample size and near-unit air-density collinearity. The temperature effect on demand is therefore large and robust to single-lag confounding adjustment, while the cumulative estimates are compromised by positivity limitation that outcome-adaptive selection cannot remove.

We consider the analysis of dyadic network data on ties between individuals in different clusters where the presence of a directed tie is reported by each individual in a dyad and the cluster-level network is of interest. A generalisation of the Social Relations Model (SRM) is proposed which includes actor, partner and dyad effects at the individual and cluster levels. The model additionally uses “double-sampling” of ties to estimate a measurement model which adjusts for and quantifies the extent of reporter effects. The model can be viewed as a type of multilevel structural equation model, with multiple cross-classified random effects, which can be estimated using Markov chain Monte Carlo (MCMC) methods in Bayesian software. Using parameter estimates from this multilevel SRM, we then propose two alternative ways of deriving the between-cluster network that are based on predictions of the strength of between-cluster ties. Our approach is illustrated using data on social support networks in a rural community in Nicaragua where individual reports of bidirectional exchanges of support with individuals from other households are used to derive the between-household network.

High-dimensional time series has diverse applications in econometrics and finance. Recent models for capturing temporal dependence have employed a bilinear representation for matrix time series, or the Tucker-decomposition based representation in case of tensor time series. A Tucker-decomposition based temporal effect is difficult to interpret on many occasions, along with its computational complexity due to the non-convex nature of the underlying optimization problem. Moreover, the existing tensor models have not sufficiently explored the possibilities of imposing any lower-dimensional pattern on the transition matrices. In this work, we propose a regularized additive tensor autoregressive model with additive interaction of row-wise, column-wise and tube-wise temporal dependence, that offers more interpretability, less computational burden due to its convex nature and estimation of the underlying low rank plus sparse pattern of its transition matrices. We address the issue of identifiability of the various components in our model and subsequently develop a scalable alternating block minimization algorithm for estimating the parameters. We provide a finite sample error bound under high-dimensional scaling for the model parameters. Finally, the efficacy of the proposed model is demonstrated on synthetic and real data.

Many large-scale phase III trials with recurrent event endpoints include a pre-planned interim analysis to evaluate early futility. Conditional power (CP), which quantifies the probability of achieving statistical significance at the final analysis given the interim data, is a commonly used tool to support such decisions. The standard negative binomial model with an offset term, widely adopted for analyzing recurrent events, implicitly assumes that event rates and treatment effects remain constant over the study period. At the interim analysis, however, a substantial proportion of patients have incomplete follow-up, and when the treatment effect is delayed in onset or diminishes over time, the constant-rate assumption introduces systematic bias into the interim estimate and can lead to incorrect futility decisions. In this paper, we propose FLEX-CP-DT, a piecewise negative binomial framework that captures temporal trends in both event rates and treatment effects without imposing the constant-rate assumption. The framework yields a formula-based conditional power calculation that does not require resampling or trial simulation at the interim stage. Through extensive simulations spanning constant-effect and delayed-onset scenarios, we demonstrate that FLEX-CP-DT performs comparably to the standard approach when the constant-rate assumption holds and improves interim futility decision-making when it is violated. A case study calibrated to a published phase 3 bronchiectasis trial further illustrates the practical advantage of the proposed method in reducing the probability of falsely terminating an efficacious drug with delayed treatment onset.

Semicontinuous outcomes characterized by an exact zero mass and heavy-tailed continuous distributions challenge empirical risk modeling and machine learning. Compound Poisson-Gamma Tweedie processes impose a rigid mean-variance coupling that distorts zero-mass probabilities under extreme tail dispersion, flattening continuous density and reducing extreme-risk discrimination. Furthermore, heuristic post-hoc target capping introduces systematic bias while leaving loss gradients vulnerable to tail leverage. We present a unified, statistically coherent likelihood framework for semicontinuous modeling under arbitrary zero-inflation and explicit upper boundary constraints. Optimization and likelihood evaluations operate strictly in a latent space defined by the power-scaling mapping $y_i^* = \min\{(y_i/s)^\lambda, U^*\}$, where $s \in (0, U)$ anchors scale and $\lambda \in (0,1)$ regularizes tail curvature. The architecture partitions the response into three regimes: an unconstrained Bernoulli hurdle, an analytical continuous Gamma body, and an upper boundary point mass governed by incomplete Gamma survival ratios. Rather than an arbitrary truncation, the boundary functions as an intrinsic tail accumulator bounding score and Hessian dynamics. We derive the closed-form deviance objective, score vectors, and block-diagonal Hessian matrix, proving that parameter spaces are information-orthogonal to enable concurrent classification and regression updates. Asymptotic boundary analysis demonstrates that score and curvature scale linearly in the deep tail, establishing a quadratic restoring force that prevents gradient collapse. The framework breaks the Tweedie parameter rigidity paradox while maintaining structural fidelity across both zero occurrence and extreme severity differentiation.

In survival analysis, the time-to-event variable T is frequently subject to right censoring. Individuals may withdraw from the study for various reasons, or may not experience the event of interest before the end of follow-up. In this paper, we distinguish between two types of censoring: a potentially dependent censoring time C, which may be stochastically related to T, and an independent administrative censoring time A. In addition, the data may exhibit a cure fraction, meaning that some individuals will never experience the event. We build upon a recent work about a fully parametric mixture cure model, which accounts for dependent censoring through copulas. The proposed extension incorporates administrative censoring and allows covariates to affect all model parameters. This framework enables a more accurate modelling of the dependence between survival and censoring times while providing greater flexibility through covariate effects, leading to more individualised estimation of the cure fraction, the dependence structure, and other clinically relevant quantities. Moreover, the presence of covariates allows for weaker identification conditions.

Local differential privacy (LDP) protects individuals in a dataset by perturbing each measurement before release. For real-valued data, a widely used mechanism is additive Laplace noise. We study the problem of estimating the distribution of the latent confidential data from the privatized observations via the nonparametric maximum likelihood estimator (NPMLE) under an i.i.d. sampling model. We first show that under the Laplace convolution model, the NPMLE admits a finite-dimensional reformulation in which the support is restricted to the observation set. This reduces the original infinite-dimensional optimization over all mixing distributions to an $n$-dimensional convex optimization over mixture weights, where $n$ is the sample size. We then study the statistical convergence of the NPMLE under the 1-Wasserstein distance, and explicitly connect its convergence rate with the privacy noise scale. Allowing the privacy noise level to change with the sample size, our analysis shows that the NPMLE remains consistent when the Laplace noise grows at a rate slower than $n^{3/16}$. Conversely, when the Laplace noise is of the order $\sqrt n$ or larger, no estimator can achieve uniformly consistent recovery of the latent distribution.

Self-controlled study designs eliminate confounding by individual-level characteristics that remain constant during the observation time and are thus widely used in pharmacoepidemiology and vaccine safety research. However, existing methods remain vulnerable to time effects, including temporal trends and seasonality in the exposure or outcome, unless these are explicitly modeled or controlled through the study design. We introduce the symmetric pair matching design (SPM), a novel self-controlled method that combines design-based adjustment for time effects with automatic control of time-invariant confounding. Unlike previous design-based approaches that account for time effects, SPM uses observation time both before and after event occurrence, thereby retaining a larger proportion of the available information. We derive the theoretical properties of the method and evaluate its finite-sample performance through simulations. In the simulations, SPM produced unbiased estimates in the presence of temporal trends in both the outcome and exposure, while maintaining greater or comparable statistical efficiency than existing approaches. SPM extends the family of self-controlled methods by providing design-based control of time effects without requiring explicit modeling thereof. By combining robustness to temporal confounding with efficient use of observation time, SPM offers a practical alternative for observational studies with transient exposures. An accompanying R package is provided to facilitate its usage in applied research.

Control variables are widely used in statistical modelling to account for omitted variable bias of known confounders. However, they have largely been underexplored in deep learning. This is surprising, given that deep learning models encode image-inferable covariates, such as demographic variables, into their predictions when these covariates are correlated with the outcome—a form of omitted variable bias referred to as 'shortcut learning'. While many existing confound-control or fairness methods try to restrict the correlation of such covariates with model predictions, we show that this fails to correct for omitted variable bias. We therefore propose a control variable approach for deep learning models, based on generalised additive modelling of the effects of model inputs and covariates. As flexible additive models can suffer from concurvity, we introduce an estimation procedure that refits the final layer of a pre-trained network to include covariate effects, using cross-fitting with ridge penalisation. We show how these effects can be orthogonalised with respect to covariates to exclude their mediated effects and that model predictions can be marginalised over the covariate distribution to control for their effect. This yields unbiased, interpretable predictions and offers flexibility to model the desired effects depending on the scientific or fairness objective. We verify our approach using simulated images, and demonstrate consistent estimation of true effects. Existing methods either require more data or fail to recover the true effects. We apply our method to real neuroimaging data with experimentally induced confounding, where it recovers prediction performance to near the level of a model trained on unconfounded data. Code is available at https://github.com/mpff/cocodeel.

Randomization testing with permutations is a very common nonparametric approach to hypothesis testing. However, randomization testing can be done with other group transformations including random rotations. In this work, we consider the problem of invariance in quadratic form statistics under a unified framework with closed form p-values. In particular, we propose a nonparametric variant of the classic Durbin-Watson test for testing for autocorrelation in time series data at arbitrary lags. Our test is performed by integrating over a group of invariances of the test statistic, and easy-to-compute analytic formulae for the p-value are derived from concentration inequalities on compact groups on a per-lag basis. Thus, the usual necessity of large-scale Monte Carlo simulations is rendered unnecessary. Our tests outperform the classic Breusch-Godfrey and Ljung-Box tests on simulated data with respect to statistical power to identify significant autocorrelation. They also can be used to identify the presence of autocorrelation at large lags such as in monthly solar intensity data, which follows an approximate 11 year (132 month) cycle. The general formulation of this approach can be easily adapted to other quadratic form statistics.

This paper proposes a nonparametric Bayesian inference framework for partially identified discrete response models. The key observation is that these models map a reduced-form conditional choice probability to an identified set. Consequently, nonparametric Bayesian inference for the conditional probability mass function leads to Bayesian inference for the identified set. The inference framework nests conditional moment inequalities and linear systems with unknown coefficients as special cases. Importantly, our proposal does not require converting conditional moments into unconditional moments or discretizing covariates. We show that the posterior is consistent for the true identified set when the model is correctly specified, show that the posterior can consistently detect model misspecification, and show posterior consistency for a pseudo-identified set that is valid under misspecification. We also verify the assumptions for a class of priors based on Gaussian processes that we use to implement our proposal. These priors offer similar flexibility to frequentist partial identification methods, and are computationally attractive because posterior sampling can be performed in closed-form. We also show that many of the ideas in this paper extend to continuous responses and aggregated discrete responses (e.g., market shares).

Functional data analysis is an important statistical field that treats data as random functions. In practice, the random functions are often not fully observed but instead measured at discrete times. While simpler problems, such as mean and covariance estimation, have been widely studied for discretely observed data, optimal estimation of linear regression for this data type has remained unsolved for over two decades. To tackle this fundamental challenge, we propose a novel approach, referred to as pooling ridge estimation, which combines the advantages of pooling strategy and RKHS-based method by incorporating the unbiased estimation of operators based on discretely observed measurements from all subjects. This unified estimation framework enables us to achieve minimax optimality in prediction risk in arbitrary sampling schemes ranging from sparse to dense designs, for both scalar-on-function and function-on-function regression models. Such methodological and theoretical advances are obtained for the first time and accurately reveal the influence of discrete sampling. For scalar-on-function regression, the phase transition occurs once, separating the convergence behavior into two distinct regimes. Remarkably, for function-on-function regression, up to three phase transitions may occur, determined by the sampling frequencies of the predictor/response functions. Finally, simulation experiments and two real data examples provide empirical support for the proposed methods.

We establish non-asymptotic Berry–Esseen bounds and bootstrap approximations for the Hilbert-space norm of nondegenerate $U$-statistics. Our triangular-array framework allows the kernel to take values in a separable Hilbert space $H_n$ that may vary with the sample size, thereby covering both infinite-dimensional spaces and Euclidean spaces of increasing dimension. The Gaussian approximation combines a Hilbert-space version of Stein's method with refined control of the higher-order terms in the Hoeffding decomposition. Its error depends on fourth moments of the kernel and on the spectral geometry of the covariance operator of the Hájek projection, particularly the proportion of squared spectral mass lying outside the leading eigendirection. Consequently, the theory accommodates singular and approximately low-rank covariance operators, provided that sufficient spectral mass remains beyond the leading direction. For $H_n=\mathbb{R}^{d_n}$, we obtain dimension-explicit bounds under coordinatewise fourth-moment conditions. We also establish approximation bounds for the empirical, Gaussian-weighted, and jackknife multiplier bootstraps, yielding tests with asymptotically correct size and consistency against alternatives at covariance-dependent separation rates. For testing the vector of pairwise Kendall's tau coefficients, a matching minimax lower bound shows that the resulting separation rate is minimax rate-optimal over dense Gaussian correlation alternatives.

Bayesian methods are essential in modern data analysis in ecology and evolutionary biology. They provide a flexible framework for modeling complex data-generating processes, while quantifying uncertainty based on the classical subjective interpretation of probability. However, Bayesian inference may provide misleading measures of uncertainty, particularly when the fitted model fails to adequately represent the true data-generating process. Although nonparametric approaches such as leave-k-out diagnostics and bootstrap resampling offer more robust alternatives under model misspecification, their computational cost is often too demanding because they require repeatedly refitting the same Bayesian model. In this paper, we introduce posterior covariance sensitivity analysis (PosCoSeA), a computationally efficient strategy for approximating leave-k-out diagnostics and bootstrap resampling without repeated model refitting. The performance of the methods is evaluated through both simulation studies and an application to real ecological data. We also provide an R package that implements these methods to facilitate their practical application.

We propose an omnibus nonparametric test for structural changes in the multivariate volatility matrix. The test aggregates bounded generalized quantile scores over a range of quantile levels and has a weighted leave-$q$-out $U$-statistic representation. Deleting nearby index pairs renders the centering effect induced by serial dependence asymptotically negligible. All quantities required for implementation, including the variance estimator used for standardization, are constructed under the null, without specifying volatility dynamics under the alternative. The standardized statistic converges to a standard normal distribution. We establish consistency against fixed alternatives that generate a positive integrated quantile-score signal and derive nontrivial local power against smooth departures and increasingly sharp transitions approaching multiple structural breaks. The bounded-score construction avoids the finite fourth- or eighth-moment conditions commonly imposed by least-squares and quasi-likelihood procedures, while aggregation across quantiles uses more distributional information than single-quantile methods. Monte Carlo results show satisfactory size and favorable power under heavy-tailed innovations, with competitive performance under Gaussian innovations. An application to the Fama–French three-factor model provides evidence against stability of the factor covariance matrix over the full sample and several economically relevant subsamples.

Multinomial choice models allow flexible substitution patterns but become computationally demanding with many alternatives or observations. With a fixed per-observation simulation budget, simulated maximum likelihood introduces simulation bias, while each optimization step requires a full-sample likelihood evaluation. We propose Stochastic Approximation with Unbiased Simulated Scores (SAUSS), an averaged stochastic approximation based on conditionally unbiased mini-batch score estimates. Each iteration uses a fixed mini-batch regardless of sample size. For multinomial probit, accept-reject sampling provides exact conditional draws and unbiased score estimates for any fixed number of accepted draws. Under local conditions, asymptotic theory for the averaged estimator and the partial-sum process of the SAUSS iterates incorporates mini-batch and simulation variability and supports random-scaling and plug-in inference. In simulations and an application, SAUSS gives comparable results in less than 1% of the computation time of simulated maximum likelihood. SAUSS extends to limited dependent variable models with conditional-expectation score representations and exact conditional sampling.

Exponential random graph models (ERGMs) describe dependence among network ties, but inference becomes difficult when the likelihood is intractable and candidate network statistics are strongly correlated. We introduce BERGM Elastic Net, an adaptive empirical-Bayes approach that combines lasso shrinkage with ridge stabilization in a Bayesian ERGM. A latent-variable formulation supports approximate exchange sampling, while empirical-Bayes updates adapt the amount of regularization to the observed network. We connect the proposed prior to elastic-net penalized likelihood and clarify the interpretation of thresholded reporting and coefficient grouping. The method is developed for over-specified network models containing many related structural and covariate effects.

Modern medical research, such as dose-finding studies, seamless trials, and shared control designs, often involves comparing multiple treatments simultaneously. Despite its wide applications, most research focuses on continuous endpoints, leaving the inference for general outcome types in high demand. In this article, we propose a new inference method for conducting multiple-treatment comparisons involving endpoints within the generalized linear model (GLM) framework under covariate-adaptive randomization (CAR). First, we investigate the asymptotic properties of the standard Wald z-statistics (z-scores) in multi-arm trials, highlighting issues when the working model is misspecified, particularly through omitted covariates. Our theoretical findings reveal that these \textcolor{black}{z-scores} do not consistently converge to a standard multivariate normal distribution, leading to either conservative or inflated Type I error rates, depending on the specific GLM endpoint. Second, based on these theoretical results, we develop adjusted test statistics to correct the distributional problems. To appropriately control the family-wise Type I error rate inherent in multi-arm comparisons, we incorporate our adjusted statistics with Simes-type multiple-testing procedures. This robust inference method can effectively control Type I error while potentially improving power. Extensive simulation studies and a real-world application to a metastatic breast cancer trial confirm the effectiveness and practicality of our approach.

Bayesian quantile regression based on the asymmetric Laplace (AL) distribution can be sensitive to extreme observations because of its exponentially decaying tails. We propose a robust error distribution constructed as a finite mixture of the AL distribution and a log-Pareto scale mixture of asymmetric Laplace distributions (LPAL). Unlike a direct log-Pareto extension of the normal location-scale representation of the AL distribution, the proposed AL-LPAL mixture preserves the prescribed quantile and exhibits log-regularly varying behavior in both tails. The LPAL component also has an unbounded density at the target quantile, yielding a distribution that combines sharp central concentration with super-heavy tails. We establish posterior robustness under arbitrarily extreme contamination and provide sufficient conditions for the existence of posterior moments of the regression coefficients and scale parameter. For posterior computation, we develop a Gibbs sampler using latent-variable augmentations and a computationally efficient mean-field variational Bayes approximation. Simulation studies show that the proposed method is competitive under moderate contamination and maintains stable point estimation with comparatively concentrated posterior intervals, particularly when severe contamination affects the quantile of interest. Applications to carbon dioxide and Boston housing data, using the same preprocessing as existing robust Bayesian quantile regression analyses, show favorable predictive performance across nearly all quantile levels and loss criteria considered.

Principal component analyses are often applied to spatial data towards inference on latent modes of spatial variation. These analyses are widespread across domains including spatial transcriptomics and environmental sciences, where the modes of spatial variation are represented by corresponding factors of gene expression or remotely sensed time series measurements. Many methods have been proposed for incorporating spatial information into a probabilistic PCA framework; however, there are three main drawbacks to currently available approaches. First, the loadings matrices are not orthogonal, and subsequent orthogonalization of those loadings corrupts the original prior spatial information. Furthermore, currently proposed methods assume stationarity in their spatial prior. Finally, current methods typically do not achieve linear-time computational complexity with respect to the number of spatial locations. To resolve these problems, we first parameterize the model directly with orthogonal loadings. For the prior distribution, we derive the sampling distribution of an SVD transformation with $k$ unique and $m-k$ repeated singular values. We then show under this model that the maximum a posteriori estimator for the orthogonal loadings is the eigendecomposition of $S + \frac{1}{n}\Sigma$, where $S$ is the empirical covariance matrix and $\Sigma$ is the prior spatial covariance. We develop a minorization-maximization-within-EM algorithm that is linear in computational complexity with respect to the number of spatial locations. We further extend our MM-EM algorithm to handle held-out locations and develop a validation strategy for optimizing the nonstationary prior covariance. Our methodology is used to infer the spatial distribution of direction-specific length scales in a human brain spatial transcriptomics case study, as well as a continental-scale phenology case study in sub-Saharan Africa.

Missing values, atypical observations, and heterogeneity across latent groups are common sources of complexity in regression data. The contaminated Gaussian cluster-weighted model (CG-CWM) provides a natural framework for handling atypical observations, including outliers and leverage points, in model-based clustering. We extend the CG-CWM to data with missing-at-random (MAR) values in both the response and covariate spaces. The proposed model provides clustering in regression analysis while distinguishing typical observations, outliers, and good and bad leverage points. By treating covariates as random, the model preserves assignment dependence, allowing them to contribute directly to cluster formation. Maximum likelihood estimation is performed through an expectation-conditional maximization (ECM) algorithm that accounts for four sources of incomplete information: missing responses and covariates, unknown component memberships, and latent contamination indicators. Conditional on these indicators, the joint distribution of responses and covariates is multivariate Gaussian, yielding closed-form conditional distributions for missing values and incorporating missingness uncertainty directly into parameter updates. Thus, missing values are handled within model fitting rather than by preliminary imputation. The framework provides clustering, clusterwise regression, model-based treatment of MAR values, and detection of atypical observations. Performance is assessed through numerical studies under varying levels of contamination and missingness patterns, and a real data application.

This note studies the geometry of full conformal prediction (FullCP) regions generated by an empirical energy-form pairwise score. Candidate-score convexity alone does not guarantee connected FullCP regions, even for empirical averages of losses convex in the candidate argument. For the energy-form score, each leave-one-out training comparison reduces exactly to a pairwise-dissimilarity sublevel condition. Under symmetry, a constant diagonal, a diagonal lower bound, and attainment of the associated Fréchet-type objective, all comparison regions share a minimizer; if they are convex, every nontrivial exact conformal region is star-shaped about that point. For power distances $\rho_\beta(x,y)=\|x-y\|^\beta$, this geometry holds for $\beta\ge1$, while the conventional energy score is strictly proper for $0<\beta<2$. For $d=1,\beta=1$, every nontrivial empirical-CRPS FullCP region is a nonempty closed interval (possibly $\mathbb R$ when $m=1$). For $1<\beta<2$ and $m\ge2$, explicit data-checkable derivative bounds give Lipschitz control of the radial exits and exact conformal radial function. Combined with directional root search and classical Lipschitz extensions, they yield certified inner and outer radial envelopes of width at most $\delta+2\widetilde Lh_{\mathcal U}$ and same-ray Hausdorff guarantees. An analytic two-dimensional example shows why preserving star-shaped but nonconvex geometry can matter. A staged two-dimensional study finds modest but systematic tightening of the generic certificate and frequent robust nonconvexity witnesses, with detected normalized radial departures typically small. The method is intended for low-dimensional multivariate outputs rather than high-dimensional scaling or runtime improvement.

Probabilistic models of Directed Acyclic Graphs (DAGs) with latent variables impose equality constraints on the observed data distribution beyond ordinary conditional independencies. These so-called Verma constraints arise in nested Markov models associated with Acyclic Directed Mixed Graphs, the latent projection of latent-variable DAGs. While nested Markov models have been extensively studied from the perspectives of graphical representation and causal identification, their implications for semiparametric efficiency theory remain less understood. We develop results toward establishing a semiparametric framework for statistical models defined by Verma constraints. Our key observation is that nested Markov constraints admit weighted conditional-moment representations under post-fixing distributions induced by graphical fixing operations. We show that fixing induces weighted orthogonality relations in L2(P), thereby converting Verma constraints into explicit tangent-space restrictions. Building on this representation, we characterize the tangent-space orthocomplement for models defined by a single nested Markov constraint through residualized weighted moment functions. This geometric formulation yields Hilbert-space characterizations of semiparametric efficient influence functions and efficiency bounds via orthogonal projection and equivalent minimum-variance formulations. We further discuss extensions to models involving multiple nested Markov constraints, for which we characterize a subspace of the orthocomplement as sums of the corresponding weighted orthogonality relations, while leaving the complete tangent-space characterization open. More broadly, our results connect nested graphical structure with semiparametric Hilbert-space geometry and provide a foundation for a general efficiency theory for nested Markov models. We illustrate the framework through several latent-variable DAGs.

The embedding dimension of categorical predictors is usually selected through heuristic tuning, although it directly affects model complexity, approximation quality, and finite-sample generalization. This paper formulates embedding dimension selection as a constrained allocation problem. The main contribution is to show that embedding capacity can be allocated across heterogeneous categorical predictors according to an explicit approximation-estimation tradeoff. We characterize approximation error through the singular-value tail of the latent category representation, while estimation error increases with total embedding complexity. Under a fixed global embedding budget and a tractable approximation model, this leads to a closed-form allocation rule in which the dimension assigned to each predictor is proportional to the square root of its approximation value relative to its parameter cost. Simulation experiments support the proposed approximation-estimation interpretation and show that the allocation rule improves budget efficiency relative to standard uniform and cardinality-based heuristics, particularly when the budget is binding and predictor heterogeneity is substantial. A real-data healthcare application further shows improvements in predictive accuracy and probabilistic calibration. Overall, the results establish embedding-dimension allocation as a principled finite-sample optimization problem rather than a purely heuristic modeling choice.

We introduce a novel latent space network model for analyzing multivariate time series of neural spike-train data. The methodology is motivated by an experimental study in mice, where neuronal responses were collected under a sequence of visual discrimination tasks. We adopt a latent variable framework to model the firing rates of aggregated brain areas, while simultaneously inferring the interactions between regions via a hidden network structure. This interaction network is embedded in a geometric latent space, enabling interpretable visualizations and novel model-based summaries. The proposed framework provides an intuitive interpretation of the latent variables, which bear a conceptual connection to node eigen-centrality measures. To capture temporal dependence, we incorporate a nested hidden Markov structure that can flexibly represent non-linear shifts that are induced by the changing of experimental conditions. We further establish theoretical properties of the model by deriving sufficient conditions that prevent degeneracy, thereby guiding our model assumptions. Overall, the proposed methodology provides a unified framework to characterize brain activity, its temporal dynamics, and spillover effects through a hidden latent space network model.

Envelope methods improve estimation efficiency in multivariate analysis by isolating low-dimensional structures that contain all the information material to the parameter of interest. In multivariate linear regression with random predictors, predictor envelope models achieve this goal by removing variation in the predictors that is immaterial to the regression. In many applications, observations are naturally divided into several groups, such as treatment groups, regions, or demographic strata, and the regression relationship may differ between groups. Motivated by this setting, we propose a groupwise predictor envelope model for multivariate linear regression. The proposed model assumes that the group-specific regression coefficient matrices are represented through a common predictor envelope subspace while allowing group-specific regression effects and group-specific error covariance matrices. We derive an objective function for estimating the common predictor envelope, obtain the corresponding regression estimators, and establish asymptotic normality together with an explicit asymptotic variance formula. Moreover, we show that the proposed estimator is asymptotically more efficient than the estimator obtained by fitting predictor envelope models separately to each group. This theoretical advantage over the existing work is also demonstrated through simulation studies.

Cox regression with Lasso penalization is widely used for variable selection in time-to-event data, but reliable coefficient inference after selection remains difficult. We investigate bootstrap inference for the debiased Cox estimator after Cox Lasso selection. Two score-based procedures are considered: a wild bootstrap using independent multipliers and an Efron bootstrap using centered multinomial resampling weights. Both procedures keep the original Cox Lasso fit and debiasing matrix fixed and apply bootstrap weights to the score contributions that determine the first-order distribution of the debiased estimator. We establish asymptotic validity for both bootstrap schemes and study their finite-sample behavior in extensive simulations. The bootstrap intervals improve coverage accuracy over first-order debiased Wald intervals in many small- and moderate-sample settings, with the size of the improvement depending on the simulation setting and tuning rule. A SUPPORT2 application illustrates the procedures in practice.

Many time series are not stable over their full observation horizon, but may contain scientifically meaningful periods during which a signal remains stable up to a prescribed tolerance. We formulate this as an existence test for stable intervals in a non-stationary regression model with dependent, locally stationary errors. For a signal $d$ derived from the mean function, including deviations from reference levels and derivatives, stability over duration $\delta$ is encoded by $d_\infty=\inf_{t\in[0,1-\delta]}\sup_{s\in[t,t+\delta]}|d(s)|$. The hypothesis $d_\infty\ge\Delta$ states that no interval of length $\delta$ remains within the tolerance $\Delta$, while rejection provides evidence for the existence of a relevant stable period. We estimate $d$ by local linear regression and construct plug-in tests for $d_\infty$. The asymptotic distribution is determined only by near-extremal windows at which the minimax functional is attained. We formalize this localization through extremal sets and derive Gaussian and extreme value approximations for kernel estimators over possibly shrinking index sets with time-varying long-run variance. The resulting tests are consistent and have asymptotic level control. The theory also extends supremum-based relevant-change tests to time-varying long-run variance and derivative-based hypotheses. Simulations and applications to physiological and engineering time series illustrate the method.