Accurate imputation of censored data due to the limit of detection (LOD) is essential in many scientific fields. Existing imputation approaches typically rely on strict distributional assumptions or linear regression models, limiting their ability to capture complex non-linear relationships in multidimensional censored data. To address this limitation, we propose a non-parametric imputation method for censored data, termed NPIC, which iteratively imputes each censored variable using random forests and kernel density estimation (KDE) within a Gibbs sampling framework. Specifically, for each censored variable, NPIC treats it as the response and the remaining variables as predictors, trains a random forest model, computes the corresponding residuals, and estimates their density using KDE. This estimated residual density is then used to construct a truncated density, from which the expectation is calculated and used as the imputed value. Simulations on public datasets show that NPIC outperforms state-of-the-art methods. By integrating NPIC with MissForest, we develop NPICM, a unified non-parametric method for imputing both censored and missing values, and demonstrate its effectiveness on a real-world water quality dataset.
This paper develops an optimal model averaging approach for linear measurement error models, when the response variable is subject to random right censoring. Within this context, we propose a bias-corrected weighted least squares estimation for unknown regression parameters. A novel Mallows-type weight choice criterion that skillfully bypasses the unavailable true covariates is developed for allocating model weights. We provide two theoretical justifications for our proposal. First, we establish the asymptotic optimality property of the resulting model averaging estimator when all candidate models are misspecified. Second, we show our model averaging method estimates the model parameters with root-n rate when the true model is among the candidate models. Particularly, it is demonstrated that the optimality is still valid when a model screening strategy is conducted prior to model averaging. Simulation studies and a real data analysis highlight the superiority of the proposed method.
Kernel two-sample tests have been widely used, and the development of efficient methods for high-dimensional, large-scale data is receiving increasing attention in the big data era. However, existing methods, such as the maximum mean discrepancy (MMD) and recently proposed kernel-based tests for large-scale data, are computationally intensive and/or ineffective for some common alternatives in high-dimensional data. In this paper, we propose a new test that exhibits high power across a wide range of alternatives. Furthermore, the new test is more robust to high dimensions than existing methods and does not require optimization procedures for choosing kernel bandwidth and other parameters through data splitting. Numerical studies demonstrate that the new approach performs well on both synthetic and real-world data.
In this paper, we consider the problem of prediction inference for multiple target datasets under varying covariate shifts, a challenging scenario encountered in various economic and biomedical applications. We propose a two-step generative conformal prediction (GCP) method, which is efficient in handling imbalanced datasets. General results on the coverage validity of the proposed method are established, which elucidates the influence of the conditional generative learning procedure on the coverage validity. Furthermore, we prove the robustness of the proposed GCP intervals and provide theoretical evidence to demonstrate the advantages of our method over density ratio-based methods. Numerical studies, including simulation and real data examples, are conducted to compare the proposed method with existing density ratio-based approaches, to empirically illustrate the effectiveness of our method.
Additive models offer a flexible framework for modeling nonlinear relationships between predictors and a continuous response variable, where nonparametric terms are commonly modeled via penalized splines. However, such models can be sensitive to outliers, leading to potentially unreliable estimation and inferential results. Existing methods and software for robust additive models tend to be computationally burdensome, and can come with restrictions such as permitting only one nonparametric term or not offering uncertainty quantification. In this article, we propose a fast, robust approach to fitting additive models based on the gamma-divergence, which we refer to as gamma-divergence additive models or GDAMs. Specifically, we apply gamma-divergence to the restricted maximum likelihood function of the additive model, based on treating the smoothing coefficients as random effects. This leads to an efficient minorization-maximization algorithm that adaptively downweights the impact of outlying observations via a set of normalized power density weights. Simulation studies and an application to data concerning the distribution of federal grants across the United States confirm that GDAMs perform similarly to or better than many existing (non-)robust additive modeling methods under varying degrees of contamination, while also being computationally faster and more scalable.
Canonical polyadic (CP) decomposition is widely used for modeling multiway high-dimensional data, but it does not explicitly incorporate mode-specific structural information such as temporal smoothness, network cohesion, or functional regularity. We propose a structure-regularized CP (SR-CP) framework that uses a unified quadratic penalty to encode diverse structural priors through positive semidefinite matrices. A soft orthogonality-promoting penalty is further introduced to enhance component distinctiveness and numerical stability. For estimation, we develop a cyclic block-coordinate descent algorithm for both rank-one and rank-R decompositions. Each regularized mode-wise update is reformulated as a ridge-type problem, leading to a tensor-specific generalized cross-validation criterion for automatic selection of regularization parameters. We establish convergence rates, consistency, and whole-tensor reconstruction error bounds under general noise conditions. Simulations show improved factor recovery and tensor reconstruction relative to classical CP. A real-image completion study demonstrates robustness under severe missingness, and the FRED-MD application yields temporally coherent and interpretable macroeconomic factors. Overall, SR-CP offers a flexible and principled framework for incorporating structural priors into CP tensor decomposition.
We propose double machine learning (DML) estimators for a partially linear model (PLM) with endogenous treatments and multivariate sample selection. We prove asymptotic normality of the estimators under mild regularity conditions and study finite sample properties on simulated data. The results demonstrate the importance of addressing sample selection in PLMs and the usefulness of the proposed estimators for avoiding selection bias. Moreover, we extend the proposed estimators to the case of endogenous switching.
The R×C ecological inference problem is widely known to be challenging. While most recent approaches rely on linear programming, maximum-entropy formulations, or Bayesian models, direct maximization of the likelihood has received limited attention due to the computational challenges of solving the associated optimization problem. We formulate the problem with the EM algorithm to maximize the likelihood given the observed data. We show that the M-step admits a closed-form solution, and we derive an explicit recursion for the exact E-step that remains computationally feasible only for small-sized instances. To scale beyond these instances, we introduce three polynomial-time approximation methods for the E-step based on multivariate normal approximations and a single multinomial representation. Projection techniques are used to enforce the accounting identity of the local and global probabilities, and a joint EM algorithm is introduced to obtain congruent estimates in both directions of the contingency structure. We evaluate the proposed methods against state-of-the-art alternatives using both real and simulated datasets. Across all settings, our approaches achieve superior accuracy, reducing the error in real election data by 17.9
Environmental spatio-temporal data often exhibit nonlinear dynamics, nonstationarity, and positive skewness, which can limit the adequacy of Gaussian random-field models. We propose a Bayesian neural field framework that combines a coordinate-based neural representation for the mean with flexible closed skew-normal residuals and a separable space–time correlation structure. This formulation enables scalable posterior inference via variational methods and yields probabilistic predictive uncertainty for large datasets. Simulation experiments show that simpler models are preferred when the data-generating mechanism is linear and Gaussian, whereas the neural-field component improves recovery under nonlinear space–time interactions and the flexible closed skew-normal component improves upper-tail calibration under skewness. In a matched-data benchmark, mean-field variational inference gives point predictions close to No-U-Turn sampler but understates part of the posterior dispersion. An application to monthly averages of the daily PM _2.5 -based AQI sub-index in Tehran (2014–2024), evaluated using completely held-out monitoring stations in 2024, illustrates the practical predictive framework. The proposed approach offers a unified Bayesian formulation for nonlinear spatio-temporal structure and asymmetric residual behavior.
Goodness-of-fit testing is a fundamental problem in statistics and has attracted considerable attention from researchers. However, in recent years, most newly proposed tests have been tailored to specific distributions, with relatively few general-purpose goodness-of-fit tests available. This highlights the need for a versatile test applicable to a broad range of distributions. In this paper, we first introduce a series of statistics for testing goodness-of-fit based on empirical distribution function (EDF), and establish the asymptotic distribution of each EDF statistic. The proposed final test statistic is constructed as the weighted sum of transformed p-values derived from each EDF statistic, which approximately follows standard Cauchy distribution. Extensive simulation studies demonstrate that the proposed method is efficiency-robust, in the sense of maintaining high power under small data contamination, and outperforms other existing methods under a wide range of alternative distributions. We further illustrate the performance of our test on several real datasets.
In this paper, we propose a new and broadly applicable root–finding method, called as the upper–crossing/solution (US) algorithm, which belongs to the category of non-bracketing (or open domain) methods. The US algorithm is a general principle for iteratively seeking the unique root θ ^* of a non-linear equation g(θ )=0 and its each iteration consists of two steps: an upper–crossing step (U-step) and a solution step (S-step), where the U-step finds an upper–crossing function or a U-function U(θ |θ ^(t)) [whose form depends on θ ^(t) being the t-th iteration of θ ^* ] based on a new notion of so-called changing direction inequality, and the S-step solves the simple U-equation U(θ |θ ^(t)) =0 to obtain its explicit solution θ ^(t+1) . The US algorithm holds two key advantages: (i) It strongly stably converges to the root θ ^* ; and (ii) it does not depend on any initial values, in contrast to Newton’s method. The key step for applying the US algorithm is to construct one simple U-function U(θ |θ ^(t)) such that an explicit solution to the U-equation U(θ |θ ^(t)) =0 is available. Based on the first–, second–, third– and block–derivative of g(θ ) , four methods are given for establishing such U-functions. An analysis of the convergence rate of the US algorithm is provided. Furthermore, we develop an acceleration technique for the US algorithm, resulting in a weakly stable convergence. Some numerical experiments and comparisons are also presented. Especially, because of the property of strongly stable convergence, the US algorithm could be one of the powerful tools for solving an equation with multiple roots.
Spatial heterogeneity is a defining feature of geospatial data and has attracted sustained attention. Spatial change points, defined as locations where the association patterns among variables undergoes a statistically significant structural shift, are a key manifestation of such heterogeneity. While change point detection methods are highly effective at identifying structural breaks in time series, their extension to spatial settings faces fundamental challenges, primarily because mainstream approaches inherently rely on the natural temporal ordering of observations, a feature that is absent in spatial data. To address this challenge, we propose two nonparametric change point detection methods requiring only observational data and geographic coordinates, where the kernel bandwidth both smooths the regression trend and defines the local scanning scale. The first constructs local discrepancy statistics based on kernel-smoothed residuals by comparing spatial association patterns across the four quadrants within circular neighborhoods defined by a kernel bandwidth. The second introduces a pseudo-temporal ordering by sorting coordinates within local strip-shaped windows and builds a cumulative sum process from the same residuals. Extensive simulations evaluate their performance across diverse spatial heterogeneity configurations, demonstrating strong detection power across scenarios. An application to air quality monitoring data from China validates their practical utility in real-world spatial analysis.
Granger causality analysis (GCA) via vector autoregression with regularization penalty identifies sparse causality in high-dimensional time series. However, the assumptions of linear relationship and uniform lag orders across variables and sensitivity to regularization parameter selection compromise causal discovery accuracy. To address these limitations, a quantile-adaptive distributional Granger causal model (QA-DGCM) is developed to infer non-causal paths via an m-stage hypothesis testings based on distributional Granger causal effect (DGCE) statistic. QA-DGCM treats variables’ quantiles as historical conditionings and quantifies non-linear causality by DGCE. To control the false discovery rate in hypothesis testing, an adaptive threshold is estimated based on the survival functions following asymptotic chi-square distributions to detect different causal relationships with high probability. Theoretical analysis establishes QA-DGCM’s asymptotic properties, including sure screening property and consistency, under mild conditions. As a model-free framework, QA-DGCM requires no specified distributional assumptions and adapts to diverse causality. Simulations demonstrate QA-DGCM’s superior accuracy over existing Granger methods across diverse causal structures. Applied to real-world time series data, it yields unique neurobiological insights advancing whole-brain effective connectivity analysis.
In high-dimensional regression, the choice of regularization penalty typically forces a rigid assumption upon the underlying signal structure, dichotomizing data into strictly sparse (Lasso, q=1 ) or entirely dense (Ridge, q=2 ) regimes. However, real-world data generating mechanisms frequently exist on a continuum between these extremes, requiring flexible geometries to handle varying degrees of sparsity and considerable multicollinearity. In this work, we propose a data-driven framework to learn the optimal regularization norm by elevating the L_q exponent ( q ∈ (0, 2] ) from a discrete choice to a strictly continuous, learnable hyper-parameter. To overcome the computational bottleneck of evaluating non-convex and non-smooth penalty landscapes, we develop a universal proximal coordinate descent solver that utilizes a safeguarded jumping threshold operator and a novel empirical Karush-Kuhn-Tucker (KKT) verification strategy. This solver is coupled with a stochastic Tree-structured Parzen Estimator (TPE) utilizing randomized internal validation splits, enabling the rapid discovery of optimal penalty geometries without over-fitting. We evaluate the framework on simulated architectures, demonstrating its dynamic adaptivity to structural sparsity, collinearity, and varying signal-to-noise ratios. Applied to four high-dimensional genomic datasets (scaling up to P ≈ 50,000 features), our generalized adaptive bridge regression (GABR) framework successfully identifies optimal, off-grid grouping architectures ( q ≈ 1.63 to 1.80), outperforming purely sparse and purely dense alternatives. These results demonstrate that the exact regression geometry can be efficiently learned from the data, enabling a unified approach to high-dimensional inference without the computational restrictions of exhaustive discrete grid searches.
In many research fields, there is an increased availability of network data arising as replicated networks. However, most statistical models for network data in the literature are designed for a single network. Among these, the Stochastic Block Model is arguably the most popular model to perform vertex clustering and community detection. We propose the Hierarchical Stochastic Block Model, a generalization of the SBM to the setting of replicated networks. This model uses a Hierarchical Pitman-Yor prior for the block allocation vector of each graph, and allows different networks to share the same latent blocks. The number of blocks in each graph and the overall number of blocks need not be specify by the practitioner, hence avoiding complicated model selection procedures. A novel MCMC algorithm to perform posterior inference is derived. To illustrate how the model is able to capture different levels of block sharing, the HSBM is fit to a co-authorship and a brain connectomic network.
To address the issues of unstable parameter estimation and variable selection under limited target-domain observations in spatial point processes (SPPs), this paper proposes an algorithmic framework based on transfer learning. Unlike conventional target-only variable selection methods for SPPs, the proposed framework aims to leverage transferable information from source domains to improve target-domain intensity estimation and sparse variable selection. In scenarios where transferable sources are known, we develop a two-stage transfer algorithm by optimizing a Poisson quasi-likelihood objective model combined with an adaptive ℓ _0 -sparse penalty, employing Iterative Hard Thresholding (IHT) and a Warm-Start strategy to improve computational efficiency. Furthermore, when transferable sources cannot be determined, we design a data-driven source detection algorithm based on spatial block cross-validation. This approach screens candidate domains by comparing empirical loss differences, thereby reducing the risk of negative transfer. Numerical simulations conducted under Poisson and Thomas point process settings demonstrate that the proposed method can improve estimation accuracy and variable selection stability, while exhibiting robustness against clustering effects among spatial points. Finally, we apply this algorithm to analyze vehicle crime data in Nottingham, UK, which further illustrates the practical utility of the proposed method.
Copula state space models (SSMs) provide a nonlinear and non-Gaussian framework and have been effectively applied, yet their observability properties remain unexplored. Instead, proposed estimation methods were directly applied to real-world data, without verifying whether they perform reliably. We introduce a novel definition of observability and a numerical approach for assessing observability of general copula SSMs. Given an observation trajectory, we aim to recover the augmented state, which includes both parameters and the state trajectory. Observability depends on the existence of an appropriate estimator for the augmented state. In nonlinear SSMs, observability is not a global property; such an estimator may not exist for all possible observation and state trajectories. Since it is not possible to check all realizations, we consider selected ones - the point masses of a discrete density approximation (quasi-random, deterministic, low-discrepancy sampling), representing the joint distribution of observation and state trajectory. The point masses are called design trajectories. For computation, the Bayesian MCMC framework Stan is chosen. Successful convergence indicates that the augmented state is recoverable; if this is not the case for any design trajectory, the model is considered unobservable. Additionally, we propose quantifying the degree of observability. We show a high degree of observability for copula SSMs where a series of univariate states describes (i) a single and (ii) d univariate time series.
In this paper, we propose three effect-specific optimal subsampling strategies for estimating direct, total and indirect effects in partially linear mediation models with high-dimensional confounders. We construct two inverse-probability weighted subsampling Neyman-orthogonal score functions for the direct and total effects, respectively, to simultaneously eliminate selection bias and regularization bias. The unconditional asymptotic distributions of three subsample effect estimators are established and then we derive effect-specific optimal subsampling probabilities by minimizing the traces of their asymptotic variances. A two-step procedure is proposed for practical implementation by using flexible machine learning and data splitting techniques to estimate the high-dimensional nuisance functions and address potential overfitting. Extensive simulations demonstrate the superior performance of the proposed subsample estimators and an application to a real-world air pollution dataset further confirms the practical utility.
Network data, characterized by interconnected nodes and edges, is pervasive in various domains and has gained significant popularity in recent years. In network data analysis, testing the presence of community structure in a network is one of the most important research tasks. Existing tests are mainly developed for unweighted networks. In practice, many real networks are weighted and our simulation study shows that the existing methods designed for unweighted networks may not be powerful for testing weighted networks. In this paper, we study the problem of testing the existence of a community structure in general networks that are either unweighted or weighted, and either dense or sparse. We propose two new tests, namely, the weighted signed-triangle test and the empirical likelihood test. We find that both methods outperform the existing tests when the network size is small; the empirical likelihood test may further outperform the weighted signed-triangle test in small networks.
