Missing data is frequently encountered in many areas of statistics. One popular approach to address this issue is through the use of propensity score weighting. However, correctly specifying the statistical model can be a daunting task. Doubly robust estimation is attractive, as the consistency of the estimator is guaranteed when either the outcome regression model or the propensity score model is correctly specified. In this paper, we first employ information projection to develop an efficient and doubly robust estimator via indirect model calibration. The resulting propensity score estimator can be equivalently expressed as a doubly robust regression imputation estimator by imposing the internal bias calibration condition in estimating the regression parameters. In addition, using the γ-divergence measure, we generalize the information projection to allow for outlier-robust propensity score estimation. The study includes the presentation of certain asymptotic properties and findings from a simulation study, which demonstrate that the proposed method enables robust inference, not only in cases of various model assumptions being violated but also in the presence of outliers. A real-life application is also presented using data from the Conservation Effects Assessment Project.

Robustness studies of black-box models is recognized as a necessary task for numerical models based on structural equations and predictive models learned from data. These studies must assess the model’s robustness to possible misspecification of regarding its inputs (e.g., covariate shift). The study of black-box models, through the prism of uncertainty quantification (UQ), is often based on sensitivity analysis involving a probabilistic structure imposed on the inputs, while ML models are solely constructed from observed data. Our work aim at unifying the UQ and ML interpretability approaches, by providing relevant and easy-to-use tools for both paradigms. To provide a generic and understandable framework for robustness studies, we define perturbations of input information relying on quantile constraints and projections with respect to the Wasserstein distance between probability measures, while preserving their dependence structure. We show that this perturbation problem can be analytically solved. Ensuring regularity constraints by means of isotonic polynomial approximations leads to smoother perturbations, which can be more suitable in practice. Numerical experiments on real case studies, from the UQ and ML fields, highlight the computational feasibility of such studies and provide local and global insights on the robustness of black-box models to input perturbations.

The generalized Wendland covariance model is a flexible compactly supported covariance model that allows for a continuous parameterization of smoothness of the underlying Gaussian random field, and includes the celebrated Matérn as a special limit case. However, the generalized Wendland model does not cover the full range of validity of the smoothness parameter of the Matérn model. In this paper, we provide new necessary and sufficient conditions of validity of the generalized Wendland model that allows to fill this gap. The effectiveness of our proposal is illustrated through a simulation study and a re-analysis of a large geo-referenced dataset of yearly total precipitation anomalies.

We investigate the problem of clustering bipartite graphs using a simple spectral method within the framework of the Bipartite Stochastic Block Model (BiSBM), a popular model for bipartite graphs having a community structure. Our focus lies in the high-dimensional setting where the number $n_1$ of rows, and $n_2$ of columns, of the associated adjacency matrix differ significantly. A recent study by [4] has established a sufficient and necessary condition related to the sparsity level $p_{\text{max}}$ of the bipartite graph, enabling the recovery of the latent partition of the rows. In their work, [4] introduces an iterative method that extends the approach proposed by [26] to achieve the stated recovery goal. However, empirical results suggest that the subsequent refinement algorithm does not significantly enhance the performance of the spectral method, indicating that the spectral method achieves exact recovery within the same regime as the refinement method. We establish this claim by deriving new entrywise bounds on the eigenvectors of the similarity matrix utilized by the spectral method. Our analysis extends the framework of [23], which is limited to symmetric matrices with restricted dependencies. As a critical technical step, we also derive an improved concentration inequality tailored for similarity matrices.

We study the well known difficult problem of prediction in measurement error models. By targeting directly at the prediction interval instead of the point prediction, we construct a prediction interval by providing estimators of both the center and the length of the interval which achieves a pre-determined prediction level. The constructing procedure requires a working model for the distribution of the variable prone to error. If the working model is correct, the prediction interval estimator obtains the smallest variability in terms of assessing the true center and length. If the working model is incorrect, the prediction interval estimation is still consistent. We further study how the length of the prediction interval depends on the choice of the true prediction interval center and provide guidance on obtaining minimal prediction interval length. Numerical experiments are conducted to illustrate the performance and we apply our method to predict concentration of $A\mathit{\beta }1-42$ in cerebrospinal fluid in an Alzheimer’s disease data.

The block maxima method is a classical and widely applied statistical method for time series extremes. It has recently been found that respective estimators whose asymptotics are driven by empirical means can be improved by using sliding rather than disjoint block maxima. Similar results are derived for general non-degenerate U-statistics of arbitrary order, in the multivariate time series case. Details are worked out for selected examples: the empirical variance, the probability weighted moment estimator and Kendall’s tau statistic. The results are also extended to the case where the underlying sample is piecewise stationary. The finite-sample properties are illustrated by a Monte Carlo simulation study.

In the context of high-dimensional Gaussian linear regression for ordered variables, we study the variable selection procedure via the minimization of the penalized least-squares criterion. We focus on model selection where the penalty function depends on an unknown multiplicative constant commonly calibrated for prediction. We propose a new proper calibration of this hyperparameter to simultaneously control predictive risk and false discovery rate. We obtain non-asymptotic bounds on the False Discovery Rate with respect to the hyperparameter and we provide an algorithm to calibrate it. This algorithm is based on quantities that can typically be observed in real data applications. The algorithm is validated in an extensive simulation study and is compared with several existing variable selection procedures. Finally, we study an extension of our approach to the case in which an ordering of the variables is not available.

We consider a model where a signal (discrete or continuous) is observed with an additive Gaussian noise process. The signal is issued from a linear combination of a finite but increasing number of translated features. The features are continuously parameterized by their location and depend on some scale parameter. First, we extend previous prediction results for off-the-grid estimators by taking into account here that the scale parameter may vary. The prediction bounds are analogous, but we improve the minimal distance between two consecutive features locations in order to achieve these bounds.

Next, we propose a goodness-of-fit test for the model and give non-asymptotic upper bounds of the testing risk and of the minimax separation rate between two distinguishable signals. In particular, our test encompasses the signal detection framework. We deduce upper bounds on the minimal energy, expressed as the ${\ell }_2$-norm of the linear coefficients, to successfully detect a signal in presence of noise. The general model considered in this paper is a non-linear extension of the classical high-dimensional regression model. It turns out that, in this framework, our upper bound on the minimax separation rate matches (up to a logarithmic factor) the lower bound on the minimax separation rate for signal detection in the high dimensional linear model associated to a fixed dictionary of features. We also propose a procedure to test whether the features of the observed signal belong to a given finite collection under the assumption that the linear coefficients may vary, but have prescribed signs under the null hypothesis. A non-asymptotic upper bound on the testing risk is given.

We illustrate our results on the spikes deconvolution model with Gaussian features on the real line and with the Dirichlet kernel, frequently used in the compressed sensing literature, on the torus.

A fundamental aspect of statistics is the integration of data from different sources. Classically, Fisher and others were focused on how to integrate homogeneous (or only mildly heterogeneous) sets of data. More recently, as data are becoming more accessible, the question of if data sets from different sources should be integrated is becoming more relevant. The current literature treats this as a question with only two answers: integrate or don’t. Here we take a different approach, motivated by information-sharing principles coming from the shrinkage estimation literature. In particular, we deviate from the do/don’t perspective and propose a dial parameter that controls the extent to which two data sources are integrated. How far this dial parameter should be turned is shown to depend, for example, on the informativeness of the different data sources as measured by Fisher information. In the context of generalized linear models, this more nuanced data integration framework leads to relatively simple parameter estimates and valid tests/confidence intervals. Moreover, we demonstrate both theoretically and empirically that setting the dial parameter according to our recommendation leads to more efficient estimation compared to other binary data integration schemes.

Statistical depth functions provide measures of the outlyingness, or centrality, of the elements of a space with respect to a distribution. It is a nonparametric concept applicable to spaces of any dimension, for instance, multivariate and functional. Liu and Singh (1993) presented a multivariate two-sample test based on depth-ranks. We dedicate this paper to improving the power of the associated test statistic and incorporating its applicability to functional data. In doing so, we obtain a more natural test statistic that is symmetric in both samples. We derive the null asymptotic of the proposed test statistic, also proving the validity of the testing procedure for functional data. Finally, the finite sample performance of the test for functional data is illustrated by means of a simulation study and a real data analysis on annual temperature curves of ocean drifters is executed.

Comparing the survival times among two groups is a common problem in time-to-event analysis, for example if one would like to understand whether one medical treatment is superior to another. In the standard survival analysis setting, there has been a lot of discussion on how to quantify such difference and what can be an intuitive, easily interpretable, summary estimand. In the presence of subjects that are immune to the event of interest (‘cured’), we illustrate that it is not appropriate to just compare the overall survival functions. Instead, it is more informative to compare the cure fractions and the survival of the uncured subpopulations separately from each other. Our research is mainly driven by the question: if the cure fraction is similar for two available treatments, how else can we determine which is preferable? To this end, we estimate the mean survival times in the uncured fractions of both treatment groups and develop both permutation and asymptotic tests for inference. We first propose a nonparametric approach which is then extended to account for covariates by means of the semi-parametric logistic-Cox mixture cure model. The methods are illustrated through practical applications to breast cancer and leukemia data.

We introduce two new tools to assess the validity of statistical distributions. Both the simple and composite null hypothesis contexts are considered. These tools are based on components derived from a new statistical quantity, the comparison curve, which can provide a detailed appraisal of validity. The first tool is a graphical representation of these components on a bar plot (B-plot) accompagnied with related local acceptance regions. These allow getting some ideas and building some confidence about where and to which extent the data contradict the model. The knowledge such gained could also suggest an existing goodness-of-fit test to supplement this assessment with a control of the type I error. Otherwise, a new test may be preferable and the second tool is is the combination of these components to produce a powerful ${\mathit{\chi }}^2$-type goodness-of-fit test. Because the number of these components can be large, we introduce new selection rules to decide on their number. In simulations, our new adaptive goodness-of-fit tests are powerwise competitive with the best solutions recommended. Practical examples show how to use these tools to derive principled information regarding if and possibly where the model departs from the data.

High-dimensional linear models have been widely studied, but the developments in high-dimensional generalized linear models, or GLMs, have been slower. In this paper, we propose an empirical or data-driven prior leading to an empirical Bayes posterior distribution which can be used for estimation of and inference on the coefficient vector in a high-dimensional GLM, as well as for variable selection. We prove that our proposed posterior concentrates around the true/sparse coefficient vector at the optimal rate, provide conditions under which the posterior can achieve variable selection consistency, and prove a Bernstein–von Mises theorem that implies asymptotically valid uncertainty quantification. Computation of the proposed empirical Bayes posterior is simple and efficient, and is shown to perform well in simulations compared to existing Bayesian and non-Bayesian methods in terms of estimation and variable selection.

Taking subset samples from the original data set is an efficient and popular strategy to handle massive data that is too large to be directly modeled. To optimize inference and prediction accuracy, it is crucial to employ a subsampling scheme to collect observations intelligently. In this paper, we propose a space-filling subsampling method that uses distance metric-based strata to select subsamples from high-volume data sets. To minimize the maximal distance from pairs of samples that locate in the same stratum, Voronoi cells of thinnest covering lattices are used to partition the input space. In addition, subsamples that are space-filling according to the response are collected from each stratum. With the help of an algorithm to quickly identify the cell an observation locates in, the computational cost of our subsampling method is proportional to the number of observations and irrelevant to the number of cells, which makes our method applicable to extremely large data sets. Results from simulated studies and real data analysis show that the new method is remarkably better than existing approaches when used in conjunction with Gaussian process models.

Fréchet means on non-Euclidean spaces may exhibit nonstandard asymptotic rates rendering quantile-based asymptotic inference inapplicable. We show here that this affects, among others, all circular distributions whose support is not contained in a closed half circle. We exhaustively describe this phenomenon on the circle and introduce a new concept on metric spaces which we call finite samples smeariness (FSS). In the presence of FSS, for which we provide a test, it turns out that quantile-based tests for equality of Fréchet means systematically feature effective levels higher than their nominal level. This effect can persevere asymptotically for increasing sample size, or may not, depending on the type of FSS. In contrast, suitable bootstrap-based tests correct for FSS and asymptotically attain the correct level. For illustration of relevance we apply our method to directional wind data from two European cities. It turns out that quantile based tests, not correcting for FSS, find a multitude of significant wind changes. This multitude condenses to a few years featuring significant wind changes when our bootstrap tests are applied, correcting for FSS.

We investigate predictive densities for multivariate normal models with unknown mean vectors and known covariance matrices. Bayesian predictive densities based on shrinkage priors often have complex representations, although they are effective in various problems. We consider extended normal models with mean vectors and covariance matrices as parameters, and adopt predictive densities that belong to the extended models including the original normal model. We adopt predictive densities that are optimal with respect to the posterior Bayes risk in the extended models. The proposed predictive density based on a superharmonic shrinkage prior is shown to dominate the Bayesian predictive density based on the uniform prior under a loss function based on the Kullback–Leibler divergence when the variance of future samples is sufficiently large. Our method provides an alternative to the empirical Bayes method, which is widely used to construct tractable predictive densities.

In this manuscript, we study scalar-on-distribution regression; that is, instances where subject-specific distributions or densities are the covariates, related to a scalar outcome via a regression model. In practice, only repeated measures are observed from those covariate distributions and common approaches first use these to estimate subject-specific density functions, which are then used as covariates in standard scalar-on-function regression. We propose a simple and direct method for linear scalar-on-distribution regression that circumvents the intermediate step of estimating subject-specific covariate densities. We show that one can directly use the observed repeated measures as covariates and endow the regression function with a Gaussian process prior to obtain a closed form or conjugate Bayesian inference. Our method subsumes the standard Bayesian non-parametric regression using Gaussian processes as a special case, corresponding to covariates being Dirac-distributions. The model is also invariant to any transformation or ordering of the repeated measures. Theoretically, we show that, despite only using the observed repeated measures from the true density-valued covariates that generated the data, the method can achieve an optimal estimation error bound of the regression function. The theory extends beyond i.i.d. settings to accommodate certain forms of within-subject dependence among the repeated measures. To our knowledge, this is the first theoretical study on Bayesian regression using distribution-valued covariates. We propose numerous extensions including a scalable implementation using low-rank Gaussian processes and a generalization to non-linear scalar-on-distribution regression. Through simulation studies, we demonstrate that our method performs substantially better than approaches that require an intermediate density estimation step especially with a small number of repeated measures per subject. We apply our method to study association of age with activity counts.

We consider survival data in the presence of a cure fraction, meaning that some subjects will never experience the event of interest. We assume a mixture cure model consisting of two sub-models: one for the probability of being uncured (incidence) and one for the survival of the uncured subjects (latency). Various approaches, ranging from parametric to nonparametric, have been used to model the effect of covariates on the incidence, with the logistic model being the most common one. We propose a monotone single-index model for the incidence and introduce a new estimation method that is based on the profile maximum likelihood approach and techniques from isotonic regression. The monotone single-index structure relaxes the parametric logistic assumption while maintaining interpretability of the regression coefficients. We investigate the consistency of the proposed estimator and show through a simulation study that, when the monotonicity assumption is satisfied, it performs better compared to the non-constrained single-index/Cox mixture cure model. To illustrate its practical use, we use the new method to study melanoma cancer survival data.

The volume function $V(t)$ of a compact set $S\in {\mathbb{R}}^d$ is just the Lebesgue measure of the set of points within a distance to S not larger than t. According to some classical results in geometric measure theory, the volume function turns out to be a polynomial, at least in a finite interval, under a quite intuitive, easy to interpret, sufficient condition (called “positive reach”) which can be seen as an extension of the notion of convexity. However, many other simple sets, not fulfilling the positive reach condition, have also a polynomial volume function. To our knowledge, there is no general, simple geometric description of such sets. Still, the polynomial character of $V(t)$ has some relevant consequences since the polynomial coefficients carry some useful geometric information. In particular, the constant term is the volume of S and the first order coefficient is the boundary measure (in Minkowski’s sense). This paper is focused on sets whose volume function is polynomial on some interval starting at zero, whose length (that we call “polynomial reach”) might be unknown. Our main goal is to approximate such polynomial reach by statistical means, using only a large enough random sample of points inside S. The practical motivation is simple: when the value of the polynomial reach, or rather a lower bound for it, is approximately known, the polynomial coefficients can be estimated from the sample points by using standard methods in polynomial approximation. As a result, we get a quite general method to estimate the volume and boundary measure of the set, relying only on an inner sample of points. This paper explores the theoretical and practical aspects of this idea.

We propose center-outward superquantile and expected shortfall functions, with applications to multivariate risk measurements, extending the standard notion of value at risk and conditional value at risk from the real line to ${\mathbb{R}}^d$. Our new concepts are built upon the recent definition of Monge-Kantorovich quantiles based on the theory of optimal transport, and they provide a natural way to characterize multivariate tail probabilities and central areas of point clouds. They preserve the univariate interpretation of a typical observation that lies beyond or ahead a quantile, but in a meaningful multivariate way. We show that they characterize random vectors and their convergence in distribution, which underlines their importance. Our new concepts are illustrated on both simulated and real datasets.

Volatility modeling is a challenging topic in high-frequency financial data analysis. In this paper, we propose a novel Bayesian framework for modeling and forecasting spot volatility by assuming a latent GARCH structure is embedded into the volatility process at a series of unobserved “anchor” time points, which can well describe the evolving volatility of financial assets in high frequency. We introduce an ideal approximation of latent anchors, which shares similar posterior distribution with true latent anchors. Furthermore, we develop an efficient two-stage inference framework with its corresponding two-stage MCMC sampling algorithm. The simulation study and real data analysis both show our method outperforms the existing alternatives in explanation of latent anchors and the estimation and forecasting of volatility.

Model averaging (MA), a technique for combining estimators from a set of candidate models, has attracted increasing attention in machine learning and statistics. In the existing literature, there is an implicit understanding that MA can be viewed as a form of shrinkage estimation that draws the response vector towards the subspaces spanned by the candidate models. This paper explores this perspective by establishing connections between MA and shrinkage in a linear regression setting with multiple nested models. We first demonstrate that the optimal MA estimator is the best linear estimator with monotonically non-increasing weights in a Gaussian sequence model. The Mallows MA (MMA), which estimates weights by minimizing the Mallows’ $C_p$ over the unit simplex, can be viewed as a variation of the sum of a set of positive-part Stein estimators. Indeed, the latter estimator differs from the MMA only in that its optimization of Mallows’ $C_p$ is within a suitably relaxed weight set. Motivated by these connections, we develop a novel MA procedure based on a blockwise Stein estimation. The resulting Stein-type MA estimator is asymptotically optimal across a broad parameter space when the variance is known. Numerical results support our theoretical findings. The connections established in this paper may open up new avenues for investigating MA from different perspectives. A discussion on some topics for future research concludes the paper.

Space-filling designs are widely used in physical and computer experiments when the model between the response and input factors is uncertain. Recently, Chen and Tang (2022, Ann. Statist. 50, 2925–2949) justified the use of strong orthogonal arrays (SOAs) under a broad class of space-filling criteria. However, when allowable level permutations are applied to an SOA, a class of SOAs can be obtained with different geometrical structures and it is not clear which one should be selected for practical use. In this paper, we address this issue by considering a representative subset of allowable level permutations, called linear allowable level permutations. These special level permutations offer theoretical convenience in classifying various geometrically non-isomorphic SOAs. Based on these results, construction methods are provided to obtain SOAs that are more space-filling than those in the literature.

Inference for functional linear models in the presence of heteroscedastic errors has received insufficient attention given its practical importance; in fact, even a central limit theorem has not been studied in this case. At issue, conditional mean estimates have complicated sampling distributions due to the infinite dimensional regressors, where truncation bias and scaling issues are compounded by non-constant variance under heteroscedasticity. As a foundation for distributional inference, we establish a central limit theorem for the estimated conditional mean under general dependent errors, and subsequently we develop a paired bootstrap method to provide better approximations of sampling distributions. The proposed paired bootstrap does not follow the standard bootstrap algorithm for finite dimensional regressors, as this version fails outside of a narrow window for implementation with functional regressors. The reason owes to a bias with functional regressors in a naive bootstrap construction. Our bootstrap proposal incorporates debiasing and thereby attains much broader validity and flexibility with truncation parameters for inference under heteroscedasticity; even when the naive approach may be valid, the proposed bootstrap method performs better numerically. The bootstrap is applied to construct confidence intervals for centered projections and for conducting hypothesis tests for the multiple conditional means. Our theoretical results on bootstrap consistency are demonstrated through simulation studies and also illustrated with a real data example.

Generalized Bayesian inference replaces the likelihood in the Bayesian posterior with the exponential of a loss function connecting parameter values and observations. As a loss function, it is possible to use Scoring Rules (SRs), which evaluate the match between the observation and the probabilistic model for given parameter values. In this work, we leverage this Scoring Rule posterior for Bayesian Likelihood-Free Inference (LFI). In LFI, we can sample from the model but not evaluate the likelihood; hence, we use the Energy and Kernel SRs in the SR posterior, as they admit unbiased empirical estimates. While traditional Pseudo-Marginal (PM) Markov Chain Monte Carlo (MCMC) can be applied to the SR posterior, it mixes poorly for concentrated targets, such as those obtained with many observations. As such, we propose to use Stochastic Gradient (SG) MCMC, which improves performance over PM-MCMC and scales to higher-dimensional setups as it is rejection-free. SG-MCMC requires differentiating the simulator model; we achieve this effortlessly by implementing the simulator models using automatic differentiation libraries. We compare SG-MCMC sampling for the SR posterior with related LFI approaches and find that the former scales to larger sample sizes and works well on the raw data, while other methods require determining suitable summary statistics. On a chaotic dynamical system from meteorology, our method even allows inferring the parameters of a neural network used to parametrize a part of the update equations.

Panel count data arise when recurrent events are observed periodically in a study. The response variable of interest is the number of recurrent events within different time windows instead of the exact onset times of the events. The gamma frailty Poisson process model has been proposed to accommodate the within-subject correlation and overdispersion in panel count data. Although the existing methods based on the gamma frailty Poisson process model have shown some robustness against frailty distribution misspecifications, they are also found to produce biased estimates in some other cases when the gamma frailty assumption is violated. In this paper, we generalize the gamma frailty Poisson process model to allow an unknown frailty distribution for analyzing panel count data. Specifically the frailty distribution is modeled nonparametrically by assigning a Dirichlet Process Gamma Mixture prior. An efficient Gibbs sampler is developed to facilitate the Bayesian computation. Extensive simulation results suggest that the proposed Bayesian approach has an excellent performance in estimating the regression parameters and the baseline mean function and outperforms the corresponding Bayesian method based on the gamma frailty Poisson model when the gamma frailty distribution is misspecified. The proposed method is applied to a skin cancer dataset for an illustration.

The K-function is arguably the most important functional summary statistic for spatial point processes. It is used extensively for goodness-of-fit testing and in connection with minimum contrast estimation for parametric spatial point process models. It is thus pertinent to understand the asymptotic properties of estimates of the K-function. In this paper we derive the functional asymptotic distribution for the K-function estimator. Contrary to previous papers on functional convergence we consider the case of an inhomogeneous intensity function. We moreover handle the fact that practical K-function estimators rely on plugging in an estimate of the intensity function. This removes two serious limitations of the existing literature.

In this paper, we consider functional data with heterogeneity in time and population. We propose a mixture model with segmentation of time to represent this heterogeneity while keeping the functional structure. The maximum likelihood estimator is considered and proved to be identifiable and consistent. In practice, an EM algorithm is used, combined with dynamic programming for the maximization step, to approximate the maximum likelihood estimator. The method is illustrated on a simulated dataset and used on a real dataset of electricity consumption.

We consider dimension reduction of multiview data, which are emerging in scientific studies. Formulating multiview data as multivariate data with block structures corresponding to the different views, or views of data, we estimate top eigenvectors from multiview data that have two-fold sparsity, elementwise sparsity and blockwise sparsity. We propose a Fantope-based optimization criterion with multiple penalties to enforce the desired sparsity patterns and a denoising step is employed to handle potential presence of heteroskedastic noise across different data views. An alternating direction method of multipliers (ADMM) algorithm is used for optimization. We derive the ${\ell }_2$ convergence of the estimated top eigenvectors and establish their sparsity and support recovery properties. Numerical studies are used to illustrate the proposed method. Our code is available in https://github.com/lxiao665/Sparse-and-Integrative-PCA.

Gaussian processes appear as building blocks in various stochastic models and have been found instrumental to account for imprecisely known, latent functions. It is often the case that such functions may be directly or undirectly evaluated, be it in static or in sequential settings. Here we focus on situations where, rather than pointwise evaluations, evaluations of prescribed linear operators at the function of interest are (sequentially) assimilated. While working with operator data is increasingly encountered in the practice of Gaussian process modelling, mathematical details of conditioning and model updating in such settings are typically by-passed. Here we address these questions by highlighting conditions under which Gaussian process modelling coincides with endowing separable Banach spaces of functions with Gaussian measures, and by leveraging existing results on the disintegration of such measures with respect to operator data. Using recent results on path properties of GPs and their connection to RKHS, we extend the Gaussian process – Gaussian measure correspondence beyond the standard setting of Gaussian random elements in the Banach space of continuous functions. Turning then to the sequential settings, we revisit update formulae in the Gaussian measure framework and establish equalities between final and intermediate posterior mean functions and covariance operators. The latter equalities appear as infinite-dimensional and discretization-independent analogues of Gaussian vector update formulae.

We introduce methods and theory for fractionally cointegrated curve time series. We develop a variance-ratio test to determine the dimensions associated with the nonstationary and stationary subspaces. For each subspace, we apply a local Whittle estimator to estimate the long-memory parameter and establish its consistency. A Monte Carlo study of finite-sample performance is included, along with two empirical applications.

We introduce an estimator for distances in a compact Riemannian manifold based on graph Laplacian estimates of the Laplace-Beltrami operator. We upper bound the error in the estimate of manifold distances, or more precisely an estimate of a spectrally truncated variant of manifold distance of interest in non-commutative geometry (cf. [Connes and Suijelekom, 2020]), in terms of spectral errors in the graph Laplacian estimates and, implicitly, several geometric properties of the manifold. A consequence is a proof of consistency for (untruncated) manifold distances. The estimator resembles, and in fact its convergence properties are derived from, a special case of the Kontorovic dual reformulation of Wasserstein distance known as Connes’ Distance Formula.

This paper proposes estimators for the parameters of an explosive fractional Ornstein-Uhlenbeck process. The asymptotic properties for the diffusion estimators are developed under the in-fill asymptotic scheme, while the asymptotic properties for the drift estimators are developed under the double asymptotic scheme for the full range of the Hurst parameter. The double asymptotic distribution of the estimator of the persistency parameter explicitly depends on the initial condition. Simulation results demonstrate the effectiveness of the proposed estimators, and the asymptotic distributions provide a good approximation in finite samples. An empirical application is presented to demonstrate the model’s usefulness and the practical value of the asymptotic theory.

We propose a computationally efficient and sparsity adaptive procedure for estimating changes in unknown subsets of a high-dimensional data sequence. Assuming the data sequence is Gaussian, we prove that the new method successfully estimates the number and locations of changepoints with a given error rate and under minimal conditions for all sparsities of the changing subset. Our method has computational complexity linear up to logarithmic factors in both the length and number of time series, making it applicable to large data sets. Through extensive numerical studies we show that the new methodology is highly competitive in terms of both estimation accuracy and computational cost. The practical usefulness of the method is illustrated by analysing sensor data from a hydro power plant, and an efficient R implementation is available.

This paper proposes procedures for testing the equality hypothesis and the proportionality hypothesis involving a large number of q covariance matrices of dimension $p\times p$. Under a limiting scheme where p, q and the sample sizes from the q populations grow to infinity in a proper manner, the proposed test statistics are shown to be asymptotically normal. Simulation results show that finite sample properties of the test procedures are satisfactory under both the null and alternatives. As an application, we derive a test procedure for the Kronecker product covariance specification for transposable data. Empirical analysis of datasets from the Mouse Aging Project and the 1000 Genomes Project (phase 3) is also conducted.

Representation learning plays a crucial role in automated feature selection, particularly in the context of high-dimensional data, where non-parametric methods often struggle. In this study, we focus on supervised learning scenarios where the pertinent information resides within a lower-dimensional linear subspace of the data, namely the multi-index model. If this subspace were known, it would greatly enhance prediction, computation, and interpretation. To address this challenge, we propose a novel method for joint linear feature learning and non-parametric function estimation, aimed at more effectively leveraging hidden features for learning. Our approach employs empirical risk minimisation, augmented with a penalty on function derivatives, ensuring versatility. Leveraging the orthogonality and rotation invariance properties of Hermite polynomials, we introduce our estimator, named RegFeaL. By using alternative minimisation, we iteratively rotate the data to improve alignment with leading directions. We establish that the expected risk of our method converges in high-probability to the minimal risk under minimal assumptions and with explicit rates. Additionally, we provide empirical results demonstrating the performance of RegFeaL in various experiments.

We investigate the high-dimensional linear regression problem in the presence of noise that is correlated with Gaussian covariates. This type of correlation, known as endogeneity in regression models, often results from unobserved variables and other factors. It poses a significant challenge in causal inference and econometrics. In cases where covariates are high-dimensional, it is common to assume sparsity in the true parameters and to estimate them using regularization techniques, even with endogeneity. However, when sparsity is not applicable, controlling both endogeneity and high dimensionality simultaneously has not been well understood. This study demonstrates that an estimator, even without regularization, can achieve consistency, or benign overfitting, under certain assumptions about the covariance matrix. Specifically, our results indicate that the error of this estimator converges to zero when the covariance matrices of the correlated noise and the instrumental variables meet specific conditions related to their eigenvalues. We explore several extensions that relax these conditions and conduct experiments to validate our theoretical findings. As a technical contribution, we employ the convex Gaussian minimax theorem (CGMT) in our dual problem and expand upon CGMT itself.

Regression analysis is commonly conducted in survey sampling. However, existing methods fail when regression models vary across different clusters of domains. In this paper, we propose a unified framework to study the cluster-wise covariate effect under complex survey sampling based on pairwise penalties, and the associated objective function is solved by the alternating direction method of multipliers. Theoretical properties of the proposed method are investigated under regularity conditions. Numerical experiments demonstrate that the proposed method outperforms its alternatives in terms of identifying the cluster structure and estimation efficiency for both linear regression and logistic regression models. American Community Survey is used as an example to illustrate the advantages of the proposed approach.

Nowadays, several data analysis problems are high-dimensional, requiring a complexity reduction for their modeling. Under the sparsity assumption, variable selection is feasible, removing the non-influential explanatory variables. When factors are present, with their levels being dummy coded, the number of parameters included in the model grows rapidly, leading to high-dimensional problems even in cases with moderate number of factors. This fact emphasizes the need for a drastical parameter reduction, not only through variable selection but also through fusion of levels of factors. The levels fused are those not differentiating significantly in terms of their influence on the response variable. Such fusions, beyond reducing the dimension of the model, propose scale adjustments for categorical predictors. In this work a new regularization technique is introduced, called $L_0$-fused group lasso ($L_0$-FGL) for binary logistic regression. It uses a group lasso penalty for factor selection and for the fusion part it applies a $L_0$ penalty on the differences among the levels’ parameters of a categorical predictor. Using adaptive weights, the adaptive version of the $L_0$-FGL method is derived. Theoretical properties, such as existence, $\sqrt{n}$ consistency and oracle properties under certain conditions, are established. In addition, it is shown that even in the diverging case where the number of parameters $p_n$ grows with the sample size n, $\sqrt{n}$ consistency and a consistency in variable selection result are achieved, as well as a respective result on asymptotic normality for an approximate $L_0$-FGL solution. Two computational methods, the penalized iteratively reweighted least squares (PIRLS) and a block coordinate descent (BCD) approach using quasi Newton, are developed and implemented. A simulation study supports the outstanding performance of $L_0$-FGL, especially in cases with a large number of factors. Finally, we apply our method on a real dataset corresponding to breast cancer recurrence events.

This paper discusses regression analysis of interval-censored failure time data, and a new deep learning approach is proposed under the partially linear Cox model. For the analysis, we need to overcome theoretical and computational challenges arising from complex data structure where the partial likelihood function is no longer available. We propose to use a deep neural network and a B-spline function for approximating the nonlinear component and the baseline cumulative hazard function in the model, respectively. The proposed approach is flexible and able to circumvent the curse of dimensionality. At the same time, it facilitates the interpretability of covariate effects. The asymptotic properties of the resulting estimators are established. In particular, the finite-dimensional estimator of covariate effects is asymptotically normal and attains the semiparametric efficiency, while the deep nonparametric estimator achieves the minimax optimal rate of convergence. A simulation study is conducted to assess the finite-sample performance of the proposed approach and indicates that it works well in practical situations. Finally, the proposed method is applied to a set of real data that motivated this study.

Isolation forest is a popular method for anomaly detection introduced in Liu, Ting and Zhou (2008, 2012). Nonetheless, its statistical properties are little understood. We study the scoring function that is induced by the isolation forest over a finite sample at the limit when the number of trees tends to infinity, based on an analytical expression that we derive. The isolation forest method is proved to be effective at detecting geometrically isolated points within a finite sample. We then study the large sample limit of the scoring function in random designs as well as in sequences of regular designs and we find that the isolation forest method performs as a detector of the support of the underlying distribution. We also find that dense clustered anomalies are not detected asymptotically by the isolation forest method, a phenomenon known as the masking effect, but that isolation forest anomaly detection is robust to training with normal data sparsely contaminated by anomalies. Numerical examples are provided that confirm the theoretical results.

The analysis of paths in undirected graph models can be used to quantify the relevance of the strength of association in multiple paths connecting a pair of vertices of the graph. Some results are available in multivariate Gaussian settings as the covariance of two variables can be decomposed into the sum of measures related to paths joining the variables of the underlying graph. This paper studies the analysis of paths in undirected graph models for binary data, with special focus on Ising models, where the propagation of the variable status through multiple paths joining a pair of vertices is an aspect of interest. A novel logistic regression approach for baseline events in multi-way tables is proposed to show that a parameter of pairwise association can be computed by the sum of components related to paths. These components are based on products of odds ratios which are typically used to measure the dependence represented by the edges in Ising models. Specifically, two parametric decompositions are developed to gain insight on a twofold aspect of interest: the relevance of the multivariate dependence within each path connecting a pair of vertices and the interaction between the multivariate dependence in each path and in the rest of the graph. The results are illustrated through an application to cyber-security risk assessment in industrial networks.

In this paper, we investigate score function-based tests to check the significance of an ultrahigh-dimensional sub-vector of the model coefficients when the nuisance parameter vector is also ultrahigh-dimensional in linear models. We first reanalyze and extend a recently proposed score function-based test to derive, under weaker conditions, its limiting distributions under the null and local alternative hypotheses. As it may fail to work when the correlation between testing covariates and nuisance covariates is high, we propose an orthogonalized score function-based test with two merits: debiasing to make the non-degenerate error term degenerate and reducing the asymptotic variance to enhance power performance. Simulations evaluate the finite-sample performances of the proposed tests, and a real data analysis illustrates its application.

There are two classical very different extensions of the well-known Gaussian fractional Brownian motion to non-Gaussian frameworks of heavy-tailed stable distributions: the harmonizable fractional stable motion (HFSM) and the linear fractional stable motion (LFSM). As far as we know, while several articles in the literature, some of which appeared a long time ago, have proposed statistical estimators for parameters of LFSM, no estimator has yet been proposed in the framework of HFSM. Among other things, what makes statistical estimation of parameters of HFSM to be a difficult problem is that, in contrast to LFSM, HFSM is not ergodic. The main goal of our work is to propose a new strategy for dealing with this problem and constructing strongly consistent and asymptotically normal statistical estimators for both parameters of HFSM. The keystone of our new strategy consists in the construction of new transforms of HFSM which allow to obtain, at any dyadic level and also at any two consecutive dyadic levels, sequences of independent stable random variables. This new strategy might allow to make statistical inference for more general non-ergodic harmonizable stable processes and fields than HFSM. Moreover, it could turn out to be useful in study of other issues related to them.

We propose a penalized least-squares method to fit the linear regression model with fitted values that are invariant to invertible linear transformations of the design matrix. This invariance is important, for example, when practitioners have categorical predictors and interactions. Our method has the same computational cost as ridge-penalized least squares, which lacks this invariance. We derive the expected squared distance between the vector of population fitted values and its shrinkage estimator as well as the tuning parameter value that minimizes this expectation. In addition to using cross validation, we construct two estimators of this optimal tuning parameter value and study their asymptotic properties. Our numerical experiments and data examples show that our method performs similarly to ridge-penalized least-squares.

The aim of this paper is to give a complete description of the optimal estimation rates for the Hellinger loss when the square root of the density belongs to a Besov ball ${\mathfrak{B}}_{p,\infty }^{\mathit{\alpha }}(R)$. We make them explicit without further conditions when $p<2$, and under a tail dominance condition when p is larger. We also show that these rates can be improved when the density is assumed to be unimodal.

In this erratum, we correct Theorem 5.1 in [1], by mending the limit distribution of the Aalen-Johansen estimator under discontinuous survival distributions.

The key to successful statistical analysis of bivariate extreme events lies in flexible modelling of the tail dependence relationship between the two variables. In the extreme value theory literature, various techniques are available to model separate aspects of tail dependence, based on different asymptotic limits. Results from Balkema and Nolde [1] and Nolde [31] highlight the importance of studying the limiting shape of an appropriately-scaled sample cloud when characterising the whole joint tail. We now develop the first statistical inference for this limit set, which has considerable practical importance for a unified inference framework across different aspects of the joint tail. Moreover, Nolde and Wadsworth [32] link this limit set to various existing extremal dependence frameworks. Hence, a by-product of our new limit set inference is the first set of self-consistent estimators for several extremal dependence measures, avoiding the current possibility of contradictory conclusions. In simulations, our limit set estimator is successful across a range of distributions, and the corresponding extremal dependence estimators provide a major joint improvement and small marginal improvements over existing techniques. We consider an application to sea wave heights, where our estimates successfully capture the expected weakening extremal dependence as the distance between locations increases.

In a tie-breaker design (TBD), subjects with high values of a running variable are given some (usually desirable) treatment, subjects with low values are not, and subjects in the middle are randomized. TBDs are intermediate between regression discontinuity designs (RDDs) and randomized controlled trials (RCTs). TBDs allow a tradeoff between the resource allocation efficiency of an RDD and the statistical efficiency of an RCT. We study a model where the expected response is one multivariate regression for treated subjects and another for control subjects. We propose a prospective D-optimality, analogous to Bayesian optimal design, to understand design tradeoffs without reference to a specific data set. For given covariates, we show how to use convex optimization to choose treatment probabilities that optimize this criterion. We can incorporate a variety of constraints motivated by economic and ethical considerations. In our model, D-optimality for the treatment effect coincides with D-optimality for the whole regression, and, without constraints, an RCT is globally optimal. We show that a monotonicity constraint favoring more deserving subjects induces sparsity in the number of distinct treatment probabilities. We apply the convex optimization solution to a semi-synthetic example involving triage data from the MIMIC-IV-ED database.

The marginal inference of an outcome variable can be improved by closely related covariates with a structured distribution. This differs from standard covariate adjustment in randomized trials, which exploits covariate-treatment independence rather than knowledge on the covariate distribution. Yet it can also be done robustly against misspecification of the outcome-covariate relationship. Starting with a standard estimating function involving only the outcome, we first use a working regression model to compute its conditional expectation given the covariates, and then remove the uninformative part under the covariate distribution model. This effectively projects the initial function onto the joint tangent space of the full data, thereby achieving local efficiency when the regression model is correct. Importantly, even with a faulty working model, the estimator remains unbiased as the subtracted term is always asymptotically centered. Further improvement is possible if the outcome distribution also has its own structure. To demonstrate the process, we consider three examples: one with fully parametric covariates, one with a covariate following a partial parametric model against others, and another with mutually independent covariates.

Robust Bayesian methods for high-dimensional regression problems under diverse sparse regimes are studied. Traditional shrinkage priors are primarily designed to detect a handful of signals from tens of thousands of predictors in the so-called ultra-sparsity domain. However, they may not perform desirably when the degree of sparsity is moderate. In this paper, we propose a robust sparse estimation method under diverse sparsity regimes, which has a tail-adaptive shrinkage property. In this property, the tail-heaviness of the prior adjusts adaptively, becoming larger or smaller as the sparsity level increases or decreases, respectively, to accommodate more or fewer signals, a posteriori. We propose a global-local-tail (GLT) Gaussian mixture distribution that ensures this property. We examine the role of the tail-index of the prior in relation to the underlying sparsity level and demonstrate that the GLT posterior contracts at the minimax optimal rate for sparse normal mean models. We apply both the GLT prior and the Horseshoe prior to a real data problem and simulation examples. Our findings indicate that the varying tail rule based on the GLT prior offers advantages over a fixed tail rule based on the Horseshoe prior in diverse sparsity regimes.

Identifying the number of communities is a fundamental problem in community detection, which has received increasing attention recently. However, rapid advances in technology have led to the emergence of large-scale networks in various disciplines, thereby making existing methods computationally infeasible. To address this challenge, we propose a novel subsampling-based modified Bayesian information criterion (SM-BIC) for identifying the number of communities in a network generated via the stochastic block model and degree-corrected stochastic block model. We first propose a node-pair subsampling method to extract an informative subnetwork from the entire network, and then we derive a purely data-driven criterion to identify the number of communities for the subnetwork. In this way, the SM-BIC can identify the number of communities based on the subsampled network instead of the entire dataset. This leads to important computational advantages over existing methods. We theoretically investigate the computational complexity and identification consistency of the SM-BIC. Furthermore, the advantages of the SM-BIC are demonstrated by extensive numerical studies.

This paper is concerned with the problem of conditional independence testing for discrete data. In recent years, researchers have shed new light on this fundamental problem, emphasizing finite-sample optimality. The non-asymptotic viewpoint adapted in these works has led to novel conditional independence tests that enjoy certain optimality under various regimes. Despite their attractive theoretical properties, the considered tests are not necessarily practical, relying on a Poissonization trick and unspecified constants in their critical values. In this work, we attempt to bridge the gap between theory and practice by reproving optimality without Poissonization and calibrating tests using Monte Carlo permutations. Along the way, we also prove that classical asymptotic ${\mathit{\chi }}^2$- and G-tests are notably sub-optimal in a high-dimensional regime, which justifies the demand for new tools. Our theoretical results are complemented by experiments on both simulated and real-world datasets. Accompanying this paper is an R package UCI that implements the proposed tests.

We consider a multiplicative deconvolution problem, in which the density f or the survival function $SX$ of a strictly positive random variable X is estimated nonparametrically based on an i.i.d. sample from a noisy observation $Y=X\cdot U$ of X. The multiplicative measurement error U is supposed to be independent of X. The objective of this work is to construct a fully data-driven estimation procedure when the error density $f^U$ is unknown. We assume that in addition to the i.i.d. sample from Y, we have at our disposal an additional i.i.d. sample drawn independently from the error distribution. The proposed estimation procedure combines the estimation of the Mellin transformation of the density f and a regularisation of the inverse of the Mellin transform by a spectral cut-off. The derived risk bounds and oracle-type inequalities cover both – the estimation of the density f as well as the survival function $SX$. The main issue addressed in this work is the data-driven choice of the cut-off parameter using a model selection approach. We discuss conditions under which the fully data-driven estimator can attain the oracle-risk up to a constant without any previous knowledge of the error distribution. We compute convergences rates under classical smoothness assumptions. We illustrate the estimation strategy by a simulation study with different choices of distributions.

Two regeneration-based bootstrap methods, namely, the Regeneration based-bootstrap [3, 23] and the Regenerative Block bootstrap [11] are shown to be valid for the problem of estimating the integral of a function with respect to the invariant measure in a β-null recurrent Markov chain with an accessible atom. An extension of the Central Limit Theorem for randomly indexed sequences is also presented.

When studying survival data in the presence of right censoring, it often happens that a certain proportion of the individuals under study do not experience the event of interest and are considered as cured. It is then common to model the data via a mixture cure model. It depends on a model for the conditional probability of being cured (called the incidence) and a model for the conditional survival function of the uncured individuals (called the latency). This work considers a logistic model for the incidence and a semiparametric accelerated failure time model for the latency part. The estimation of this model is obtained via the maximization of the semiparametric likelihood, in which the unknown error density is replaced by a kernel estimator based on the Kaplan-Meier estimator of the error distribution. Asymptotic theory for consistency and asymptotic normality of the parameter estimators is provided. Moreover, the proposed estimation method is compared with several competitors. Finally, the new method is applied to data coming from a cancer clinical trial. An R package, called kmcure, is developed to facilitate the use of the proposed methodology in practice.

It is known that consistent nonparametric regression for a current status censored (CSC) response and a univariate predictor is possible. The paper, for the first time in the literature, presents sharp minimax theory of mean integrated squared error (MISE) convergence and methodology of adaptive estimation. Rate of the MISE convergence is classical, the sharp constant quantifies the effect of CSC, and the results hold under a mild assumption on smoothness of nuisance functions not tied to smoothness of the regression. Then the setting is extended to a multivariate predictor. Real and simulated examples are presented, as well as an illuminating comparison of theoretical results known for CSC and directly observed data.

The continuous ranked probability score (crps) is the most commonly used scoring rule in the evaluation of probabilistic forecasts for real-valued outcomes. To assess and rank forecasting methods, researchers compute the mean crps over given sets of forecast situations, based on the respective predictive distributions and outcomes. We propose a new, isotonicity-based decomposition of the mean crps into interpretable components that quantify miscalibration (MCB), discrimination ability (DSC), and uncertainty (UNC), respectively. In a detailed theoretical analysis, we compare the new approach to empirical decompositions proposed earlier, generalize to population versions, analyse their properties and relationships, and relate to a hierarchy of notions of calibration. The isotonicity-based decomposition guarantees the nonnegativity of the components and quantifies calibration in a sense that is stronger than for other types of decompositions, subject to the nondegeneracy of empirical decompositions. We illustrate the usage of the isotonicity-based decomposition and miscalibration–discrimination (MCB–DSC) plots in case studies from weather prediction and machine learning.

Many areas of science rely on simulators that implicitly encode intractable likelihood functions of complex systems. Classical statistical methods are poorly suited for these so-called likelihood-free inference (LFI) settings, especially outside asymptotic and low-dimensional regimes. At the same time, popular LFI methods — such as Approximate Bayesian Computation or more recent machine learning techniques — do not necessarily lead to valid scientific inference because they do not guarantee confidence sets with nominal coverage in general settings. In addition, LFI currently lacks practical diagnostic tools to check the actual coverage of computed confidence sets across the entire parameter space. In this work, we propose a modular inference framework that bridges classical statistics and modern machine learning to provide (i) a practical approach for constructing confidence sets with near finite-sample validity at any value of the unknown parameters, and (ii) interpretable diagnostics for estimating empirical coverage across the entire parameter space. We refer to this framework as likelihood-free frequentist inference (LF2I). Any method that defines a test statistic can leverage LF2I to create valid confidence sets and diagnostics without costly Monte Carlo or bootstrap samples at fixed parameter settings. We study two likelihood-based test statistics (ACORE and BFF) and demonstrate their performance on high-dimensional complex data. Code is available at https://github.com/lee-group-cmu/lf2i.

We tackle the extension to the vector-valued case of consistency results for Stepwise Uncertainty Reduction sequential experimental design strategies established in [3]. This leads us in the first place to clarify, assuming a compact index set, how the connection between continuous Gaussian processes and Gaussian measures on the Banach space of continuous functions carries over to vector-valued settings. From there, a number of concepts and properties from [3] can be readily extended. However, vector-valued settings do complicate things for some results, mainly due to the lack of continuity for the pseudo-inverse mapping that affects the conditional mean and covariance function given finitely many pointwise observations. We apply obtained results to the Integrated Bernoulli Variance and the Expected Measure Variance uncertainty functionals employed in [9] for the estimation for excursion sets of vector-valued functions.

In the standard difference-in-differences research design, the parallel trend assumption can be violated when the effect of some unmeasured confounders on the outcome trend is different between the treated and untreated populations. Progress can be made if there is an exogenous variable that (i) does not directly influence the change in outcome (i.e. the outcome trend) except through influencing the change in exposure (i.e. the exposure trend), and (ii) is not related to the unmeasured exposure - outcome confounders on the trend scale. Such exogenous variable is called an instrument for difference-in-differences. For continuous outcomes that lend themselves to linear modelling, so-called instrumented difference-in-differences methods have been proposed. In this paper, we will suggest novel multiplicative structural mean models for instrumented difference-in-differences, which allow one to identify and estimate the average treatment effect that is stable over time on the multiplicative scale, in the whole population or among the treated, when (i) a valid instrument for difference-in-differences is available and (ii) there is no carry-over effect across periods. We discuss the identifiability of these models, then develop efficient semi-parametric estimation approaches that allow the use of flexible, data-adaptive or machine learning methods to estimate the nuisance parameters. We apply our proposal on health care data to investigate the risk of moderate to severe weight gain under sulfonylurea treatment compared to metformin treatment, among new users of antihyperglycemic drugs.

We present a general method for obtaining posterior samples from a model for which the likelihood function is intractable or otherwise unavailable, but from which simulation is possible. While the inability to evaluate the likelihood impedes the implementation of traditional sampling algorithms, like the Metropolis–Hastings algorithm, our approach, based on the Fourier integral theorem, allows for an exact simulation-based estimate of such likelihoods. Moreover, given its foundations in the Fourier integral theorem, our approach comes with clear guidelines for setting the two tuning parameters on which it relies. Through comparison with alternative methods, specifically synthetic likelihoods and ABC, including both dimension-reduction and full-data approaches, across a variety of examples, culminating in a highly nonlinear state space model based on the Ricker map, we will demonstrate how the Fourier approach is able to provide statistically sensible, full-data likelihood estimation, requiring only that data can be simulated from the model.

Posterior contraction rates with regard to non-intrinsic metrics have been a long-standing challenge in the Bayesian analysis of high-dimensional models. This paper establishes the minimax-optimal posterior contraction rates of the Bayesian sparse spiked covariance matrix model under the spectral norm, a non-intrinsic metric for the Gaussian covariance matrix model. Our proof technique relies on the recent advance in the geometric properties of Euclidean representation for subspaces and low-rank matrices, a local asymptotic normality argument, and the distributional approximation to the asymptotic posterior distribution of the sparse spiked covariance matrix.

While analysis of variance (ANOVA) tests for differences in the means of independent samples, it is unsuitable for evaluating differences in tail behavior, especially when means do not exist or empirical estimation of means or higher moments is not consistent due to heavy-tailed distributions. Here, we propose an ANOVA-like decomposition to analyze tail variability, allowing for flexible representation of heavy tails through a set of user-defined extreme quantiles, possibly located outside the range of observations. Assuming regular variation, we introduce a test for significant tail differences across multiple independent samples and derive its asymptotic distribution. We investigate the theoretical behavior of the test statistics for the case of two samples, each following a Pareto distribution, and explore strategies for setting test hyperparameters. We conduct simulations that highlight generally reliable test behavior for a wide range of finite-sample situations. The test is applied to identify clusters of financial stock indices with similar extreme log-returns and to detect temporal changes in daily precipitation extremes in Germany.

We develop a doubly penalized constrained maximum likelihood (dPCML) method for using summary information from external studies to improve estimation efficiency for an internal study that has individual-level data, in the presence of study population heterogeneity and external information uncertainty. The dPCML method can simultaneously select and incorporate the external information that agrees with the internal study while properly accounting for the uncertainty of the external information. It allows partial information where only some but not all parameter estimates from external models are reported and/or certain parameters are known to be unequal between the internal and external studies. It can still effectively account for the external information uncertainty with only external sample sizes available instead of standard errors of parameter estimates. It covers some existing data integration methods as special cases. A detailed theoretical investigation is carried out to establish asymptotic properties of the dPCML estimator, including estimation consistency, external information selection consistency, and asymptotic normality. We also provide an algorithm for implementation and conduct comprehensive simulation studies. As an application, we build an updated model to study the risk of having high-grade prostate cancer by integrating information from two widely used risk calculators.

We consider distinguishing between two distribution tail models when tails of one model are lighter (or heavier) than those of the other. Two procedures are proposed: one scale-free and one location- and scale-free, and their asymptotic properties are established. We show the advantage of using these procedures for distinguishing between certain tail models in comparison with the tests proposed in the literature by simulation and apply them to data on daily precipitation in Green Bay, US and Jena, Germany.

Cloud-to-ground lightning strikes observed in a specific geographical domain over time can be naturally modeled by a spatio-temporal point process. Our focus lies in the parametric estimation of its intensity function, incorporating both spatial factors (such as altitude) and spatio-temporal covariates (such as field temperature, precipitation, etc.). The events are observed in France over a span of three years. Spatio-temporal covariates are observed with resolution $0.1^{\circ }\times 0.1^{\circ }$ ($\approx 100$km${}^2$) and six-hour periods. This results in an extensive dataset, further characterized by a significant excess of zeroes (i.e., spatio-temporal cells with no observed events). We reexamine composite likelihood methods commonly employed for spatial point processes, especially in situations where covariates are piecewise constant. Additionally, we extend these methods to account for zero-deflated subsampling, a strategy involving dependent subsampling, with a focus on selecting more cells in regions where events are observed. A simulation study is conducted to illustrate these novel methodologies, followed by their application to the dataset of lightning strikes.

We develop a new approach for estimating the risk of an arbitrary estimator of the mean vector in the classical normal means problem. The key idea is to generate two auxiliary data vectors, by adding carefully constructed normal noise vectors to the original data. We then train the estimator of interest on the first auxiliary vector and test it on the second. In order to stabilize the risk estimate, we average this procedure over multiple draws of the synthetic noise vector. A key aspect of this coupled bootstrap (CB) approach is that it delivers an unbiased estimate of risk under no assumptions on the estimator of the mean vector, albeit for a modified and slightly “harder” version of the original problem, where the noise variance is elevated. We prove that, under the assumptions required for the validity of Stein’s unbiased risk estimator (SURE), a limiting version of the CB estimator recovers SURE exactly. We then analyze a bias-variance decomposition of the error of the CB estimator, which elucidates the effects of the variance of the auxiliary noise and the number of bootstrap samples on the estimator’s accuracy. Lastly, we demonstrate that the CB estimator performs favorably in various simulated experiments.

In causal inference there are often multiple reasonable estimators for a given target quantity. For example, one may reasonably use inverse probability weighting, an instrumental variables approach, or construct an estimate based on proxy outcomes if the actual outcome is difficult to measure. Ideally, the practitioner decides on an estimator before looking at the data. However, this might be challenging in practice since a priori it might not be clear to a practitioner how to choose the method. If the final model is chosen after peeking at the data, naive inferential procedures may fail. This raises the need for a model selection tool, with rigorous asymptotic guarantees. Since there is usually no loss function available in causal inference, standard model selection techniques do not apply.

We propose a model selection procedure that estimates the squared ${\ell }_2$-deviation of a finite-dimensional estimator from its target. The procedure relies on knowing an asymptotically unbiased (potentially highly variable) estimate of the parameter of interest. The resulting estimator is discontinuous and does not have a Gaussian limit distribution. Thus, standard asymptotic expansions do not apply. We derive asymptotically valid confidence intervals for low-dimensional settings that take into account the model selection step.

The performance of the approach for estimation and inference for average treatment effects is evaluated on simulated data sets in low-dimensional settings, including experimental data, instrumental variables settings and observational data with selection on observables.

This work considers stationary vector count time series models defined via deterministic functions of a latent stationary vector Gaussian series. The construction is very general and ensures a pre-specified marginal distribution for the counts in each dimension, depending on unknown parameters that can be marginally estimated. The vector Gaussian series injects flexibility into the model’s temporal and cross-dimensional dependencies, perhaps through a parametric model akin to a vector autoregression. We show that the latent Gaussian model can be estimated by relating the covariances of the counts and the latent Gaussian series. In a possibly high-dimensional setting, concentration bounds are established for the differences between the estimated and true latent Gaussian autocovariances, in terms of those for the observed count series and the estimated marginal parameters. The results are applied to the case where the latent Gaussian series is a vector autoregression, and its parameters are estimated sparsely through a LASSO-type procedure.