Prevalent cohort studies are commonly conducted in many areas of research when incident cohort studies are deemed infeasible due to logistic or other constraints. While such studies are cost effective, it is known that survival data collected on prevalent cases do not form a representative sample from the target population. When the incidence (e.g. onset of disease) arise from a stationary Poisson process, it allows developing a more efficient methodology. While the stationarity assumption holds in many applications, to the best of our knowledge, the problem of establishing uniform confidence bands using data arisen in such settings has not been addressed in the current literature. We devise a method for obtaining uniform confidence bands for the cumulative hazard and the survival function built on their nonparametric maximum likelihood estimators (NPMLEs). To attain this objective, we first present results on uniform strong consistency, weak convergence and asymptotic efficiency of the NPMLE of the cumulative hazard function. Given the intractable forms of the limiting processes in this case, the idea is to numerically approximate the functionals of the asymptotic processes of the normalized NPMLEs. Our simulation studies reveal the efficiency of the estimators for finite samples. The proposed procedures are illustrated using a set of real data on patients with dementia from the Canadian Study of Health and Aging.

Although there is an extensive literature on the eigenvalues of high-dimensional sample covariance matrices, much of it is specialized to independent components (IC) models—in which observations are represented as linear transformations of random vectors with independent entries. By contrast, less is known in the context of elliptical models, which violate the independence structure of IC models and exhibit quite different statistical phenomena. In particular, very little is known about the scope of bootstrap methods for doing inference with spectral statistics in high-dimensional elliptical models. To fill this gap, we show how a bootstrap approach developed previously for IC models can be extended to handle the different properties of elliptical models. Within this setting, our main theoretical result guarantees that the proposed method consistently approximates the distributions of linear spectral statistics, which play a fundamental role in multivariate analysis. We also provide empirical results showing that the proposed method performs well for a variety of nonlinear spectral statistics.

We study frequentist risk properties of predictive density estimators for mean mixtures of multivariate normal distributions, involving an unknown location parameter $\mathit{\theta }\in {\mathbf{R}}^d$, and which include multivariate skew normal distributions. We provide explicit representations for Bayesian posterior and predictive densities, including the benchmark minimum risk equivariant (MRE) density, which is minimax and generalized Bayes with respect to an improper uniform density for θ. For four dimensions or more, we obtain Bayesian densities that improve uniformly on the MRE density under Kullback-Leibler loss. We also provide plug-in type improvements, investigate implications for certain type of parametric restrictions on θ, and illustrate and comment the findings based on numerical evaluations.

In many modern statistical problems, the limited available data must be used both to develop the hypotheses to test, and to test these hypotheses—that is, both for exploratory and confirmatory data analysis. Reusing the same dataset for both exploration and testing can lead to massive selection bias, leading to many false discoveries. Selective inference is a framework that allows for performing valid inference even when the same data is reused for exploration and testing. In this work, we are interested in the problem of selective inference for data clustering, where a clustering procedure is used to hypothesize a separation of the data points into a collection of subgroups, and we then wish to test whether these data-dependent clusters in fact represent meaningful differences within the data. Recent work by Gao, Bien and Witten (2022) provides a framework for doing selective inference for this setting, where a hierarchical clustering algorithm is used for producing the cluster assignments, which was then extended to k-means clustering by Chen and Witten (2022). Both these works rely on assuming a known covariance structure for the data, but in practice, the noise level needs to be estimated—and this is particularly challenging when the true cluster structure is unknown. In our work, we extend this work to the setting of noise with unknown variance, and provide a selective inference method for this more general setting. Empirical results show that our new method is better able to maintain high power while controlling Type I error when the true noise level is unknown.

In this paper, we develop an epsilon admissible subsets (EAS) model selection approach for performing group variable selection in the high-dimensional multivariate regression setting. This EAS strategy is designed to estimate a posterior-like, generalized fiducial distribution over a parsimonious class of models in the setting of correlated predictors and/or in the absence of a sparsity assumption. The effectiveness of our approach, to this end, is demonstrated empirically in simulation studies, and is compared to other state-of-the-art model/variable selection procedures. Furthermore, assuming a matrix-Normal linear model we show that the EAS strategy achieves strong model selection consistency in the high-dimensional setting if there does exist a sparse, true data generating set of predictors. In contrast to Bayesian approaches for model selection, our generalized fiducial approach completely avoids the problem of simultaneously having to specify arbitrary prior distributions for model parameters and penalize model complexity; our approach allows for inference directly on the model complexity. Implementation of the method is illustrated through yeast data to identify significant cell-cycle regulating transcription factors.

Bayes estimators are well known to provide a means to incorporate prior knowledge that can be expressed in terms of a single prior distribution. However, when this knowledge is too vague to express with a single prior, an alternative approach is needed. Gamma-minimax estimators provide such an approach. These estimators minimize the worst-case Bayes risk over a set Γ of prior distributions that are compatible with the available knowledge. Traditionally, Gamma-minimaxity is defined for parametric models. In this work, we define Gamma-minimax estimators for general models and propose adversarial meta-learning algorithms to compute them when the set of prior distributions is constrained by generalized moments. Accompanying convergence guarantees are also provided. We also introduce a neural network class that provides a rich, but finite-dimensional, class of estimators from which a Gamma-minimax estimator can be selected. We illustrate our method in two settings, namely entropy estimation and a prediction problem that arises in biodiversity studies.

The field of distribution-free predictive inference provides tools for provably valid prediction without any assumptions on the distribution of the data, which can be paired with any regression algorithm to provide accurate and reliable predictive intervals. The guarantees provided by these methods are typically marginal, meaning that predictive accuracy holds on average over both the training data set and the test point that is queried. However, it may be preferable to obtain a stronger guarantee of training-conditional coverage, which would ensure that most draws of the training data set result in accurate predictive accuracy on future test points. This property is known to hold for the split conformal prediction method. In this work, we examine the training-conditional coverage properties of several other distribution-free predictive inference methods, and find that training-conditional coverage is achieved by some methods but is impossible to guarantee without further assumptions for others.

The composite quantile regression (CQR) was introduced by Zou and Yuan [Ann. Statist. 36 (2008) 1108–1126] as a robust regression method for linear models with heavy-tailed errors while achieving high efficiency. Its penalized counterpart for high-dimensional sparse models was recently studied in Gu and Zou [IEEE Trans. Inf. Theory 66 (2020) 7132–7154], along with a specialized optimization algorithm based on the alternating direct method of multipliers (ADMM). Compared to the various first-order algorithms for penalized least squares, ADMM-based algorithms are not well-adapted to large-scale problems. To overcome this computational hardness, in this paper we employ a convolution-smoothed technique to CQR, complemented with iteratively reweighted ${\ell }_1$-regularization. The smoothed composite loss function is convex, twice continuously differentiable, and locally strong convex with high probability. We propose a gradient-based algorithm for penalized smoothed CQR via a variant of the majorize-minimization principal, which gains substantial computational efficiency over ADMM. Theoretically, we show that the iteratively reweighted ${\ell }_1$-penalized smoothed CQR estimator achieves near-minimax optimal convergence rate under heavy-tailed errors without any moment constraint, and further achieves near-oracle convergence rate under a weaker minimum signal strength condition than needed in Gu and Zou (2020). Numerical studies demonstrate that the proposed method exhibits significant computational advantages without compromising statistical performance compared to two state-of-the-art methods that achieve robustness and high efficiency simultaneously.

Nonregular designs are attractive, as compared with regular designs, not just because they have flexible run sizes but also because of their performances in terms of generalized resolution, projectivity, and hidden projection property. In this paper, we conduct a comprehensive study on three classes of designs that are obtained from Paley’s two constructions of Hadamard matrices. In terms of generalized resolution, we complete the study of Shi and Tang [15] on strength-two designs by adding results on strength-three designs. In terms of projectivty and hidden projection property, our results substantially expand those of Bulutoglu and Cheng [2]. For the purpose of practical applications, we conduct an extensive search of minimum G-aberration designs from those with maximum generalized resolutions and results are obtained for strength-two designs with 36, 44, 48, 52, 60, 64, 96 and 128 runs and strength-three designs with 72, 88 and 120 runs.

Variable screening for ultrahigh-dimensional data has attracted extensive attention in the past decade. In many applications, researchers learn from previous studies about certain important predictors or control variables related to the response of interest. Such knowledge should be taken into account in the screening procedure. The development of variable screening conditional on prior information, however, has been less fruitful, compared to the vast literature for generic unconditional screening. In this paper, we propose a model-free variable screening paradigm that allows for high-dimensional controls and applies to either continuous or categorical responses. The contribution of each individual predictor is quantified marginally and conditionally in the presence of the control variables as well as the other candidates by reproducing-kernel-based $R^2$ and partial $R^2$ statistics. As a result, the proposed method enjoys the sure screening property and the rank consistency property in the notion of sufficiency, with which its superiority over existing methods is well-established. The advantages of the proposed method are demonstrated by simulation studies encompassing a variety of regression and classification models, and an application to high-throughput gene expression data.

In this paper we present and analyze random number generators for the Poisson Kernel-Based Distribution (PKBD) on the sphere. We show that the only currently available sampling scheme presented in Golzy and Markatou (2020) can be improved by a better selection of hyper-parameters but still yields an unbounded rejection constant as the concentration parameter approaches 1. Furthermore, we introduce two additional and superior sampling methods for which boundedness in the above mentioned case can be obtained. The first method proposes initial draws from angular central Gaussian distribution and offers uniformly bounded rejection constants for a significant part of the PKBD parameter space. The second method uses adaptive rejection sampling and the results of Ulrich (1984) to sample from the projected Saw distribution (Saw, 1978). Finally, both new methods are compared in a simulation study.

Besov priors are nonparametric priors that can model spatially inhomogeneous functions. They are routinely used in inverse problems and imaging, where they exhibit attractive sparsity-promoting and edge-preserving features. A recent line of work has initiated the study of their asymptotic frequentist convergence properties. In the present paper, we consider the theoretical recovery performance of the posterior distributions associated to Besov-Laplace priors in the density estimation model, under the assumption that the observations are generated by a possibly spatially inhomogeneous true density belonging to a Besov space. We improve on existing results and show that carefully tuned Besov-Laplace priors attain optimal posterior contraction rates. Furthermore, we show that hierarchical procedures involving a hyper-prior on the regularity parameter lead to adaptation to any smoothness level.

We investigate the frequentist guarantees of the variational sparse Gaussian process regression model. In the theoretical analysis, we focus on the variational approach with spectral features as inducing variables. We derive guarantees and limitations for the frequentist coverage of the resulting variational credible sets. We also derive sufficient and necessary lower bounds for the number of inducing variables required to achieve minimax posterior contraction rates. The implications of these results are demonstrated for different choices of priors. In a numerical analysis we consider a wider range of inducing variable methods and observe similar phenomena beyond the scope of our theoretical findings.

Linear thresholding models postulate that the conditional distribution of a response variable in terms of covariates differs on the two sides of a (typically unknown) hyperplane in the covariate space. A key goal in such models is to learn about this separating hyperplane. Exact likelihood or least squares methods to estimate the thresholding parameter involve an indicator function which make them difficult to optimize and are, therefore, often tackled by using a surrogate loss that uses a smooth approximation to the indicator. In this paper, we demonstrate that the resulting estimator is asymptotically normal with a near optimal rate of convergence: $n^{-1}$ up to a log factor, in both classification and regression thresholding models. This is substantially faster than the currently established convergence rates of smoothed estimators for similar models in the statistics and econometrics literatures. We also present a real-data application of our approach to an environmental data set where $C{O}_2$ emission is explained in terms of a separating hyperplane defined through per-capita GDP and urban agglomeration.

The Box–Cox transformation model has been widely applied for many years. The parametric version of this model assumes that the random error follows a parametric distribution, say the normal distribution, and estimates the model parameters using the maximum likelihood method. The semiparametric version assumes that the distribution of the random error is completely unknown; existing methods either need strong assumptions, or are less effective when the distribution of the random error significantly deviates from the normal distribution. We adopt the semiparametric assumption and propose a maximum profile binomial likelihood method. We theoretically establish the joint distribution of the estimators of the model parameters. Through extensive numerical studies, we demonstrate that our method has an advantage over existing methods when the distribution of the random error deviates from the normal distribution. Furthermore, we compare the performance of our method and existing methods on an HIV data set.

This paper addresses the problem of estimating the Hurst exponent of the fractional Brownian motion from continuous time noisy sample. When the Hurst parameter is greater than $3∕4$, consistent estimation is possible only if either the length of the observation interval increases to infinity or intensity of the noise decreases to zero. The main result is a proof of the Local Asymptotic Normality (LAN) of the model in these two regimes which reveals the optimal minimax estimation rates.

In statistical learning, identifying underlying structures of true target functions based on observed data plays a crucial role to facilitate subsequent modeling and analysis. Unlike most of those existing methods that focus on some specific settings under certain model assumptions, a general and novel framework is proposed for recovering the true structures of target functions by using unstructured M-estimation in a reproducing kernel Hilbert space (RKHS) in this paper. This framework is inspired by the fact that gradient functions can be employed as a valid tool to learn underlying structures, including sparse learning, interaction selection and model identification, and it is easy to implement by taking advantage of some nice properties of the RKHS. More importantly, it admits a wide range of loss functions, and thus includes many commonly used methods as special cases, such as mean regression, quantile regression, likelihood-based classification, and margin-based classification, which is also computationally efficient by solving convex optimization tasks. The asymptotic results of the proposed framework are established within a rich family of loss functions without any explicit model specifications. The superior performance of the proposed framework is also demonstrated by a variety of simulated examples and a real case study.

In the famous least sum of trimmed squares (LTS) estimator [21], residuals are first squared and then trimmed. In this article, we first trim residuals – using a depth trimming scheme – and then square the remaining of residuals. The estimator that minimizes the sum of trimmed and squared residuals, is called an LST estimator.

Not only is the LST a robust alternative to the classic least sum of squares (LS) estimator. It also has a high finite sample breakdown point-and can resist, asymptotically, up to $50\mathrm{\%}$ contamination without breakdown – in sharp contrast to the $0\mathrm{\%}$ of the LS estimator.

The population version of the LST is Fisher consistent, and the sample version is strong, root-n consistent, and asymptotically normal. We propose approximate algorithms for computing the LST and test on synthetic and real data sets. Despite being approximate, one of the algorithms compute the LST estimator quickly with relatively small variances in contrast to the famous LTS estimator. Thus, evidence suggests the LST serves as a robust alternative to the LS estimator and is feasible even in high dimension data sets with contamination and outliers.

Envelope methods offer targeted dimension reduction for various statistical models. The goal is to improve efficiency in multivariate parameter estimation by projecting the data onto a lower-dimensional subspace known as the envelope. Envelope approaches have advantages in analyzing data with highly correlated variables, but their iterative Grassmannian optimization algorithms do not scale very well with high-dimensional data. While the connections between envelopes and partial least squares in multivariate linear regression have promoted recent progress in high-dimensional studies of envelopes, we propose a more straightforward way of envelope modeling from a new principal component regression perspective. The proposed procedure, Non-Iterative Envelope Component Estimation (NIECE), has excellent computational advantages over the iterative Grassmannian optimization alternatives in high dimensions. We develop a unified theory that bridges the gap between envelope methods and principal components in regression. The new theoretical insights also shed light on the envelope subspace estimation error as a function of eigenvalue gaps of two symmetric positive definite matrices used in envelope modeling. We apply the new theory and algorithm to several envelope models, including response and predictor reduction in multivariate linear models, logistic regression, and Cox proportional hazard model. Simulations and illustrative data analysis show the potential for NIECE to improve standard methods in linear and generalized linear models significantly.

We develop a family of doubly robust kernel machines for classification in the presence of missing covariates. We assume that the missingness is missing at random and the missing pattern is homogeneous over a subset of covariates. First, we construct a novel convex augmented loss function using inverse probability weighting, multiple imputation, and surrogacy. It features (i) the double robustness against misspecification of the missing mechanism or the imputation model, and (ii) computation feasibility via a constrained quadratic optimization. Second, we obtain theoretical results for the proposed kernel machine, which include Fisher consistency, an upper bound of the excess risk, and the rate of convergence. We demonstrate the finite sample performance of the proposed kernel machine through simulation and real data analysis.

The categorical Gini correlation proposed by Dang et al. [7] is a dependence measure to characterize independence between categorical and numerical variables. The asymptotic distributions of the sample correlation under dependence and independence have been established when the dimension of the numerical variable is fixed. However, its asymptotic behavior for high dimensional data has not been explored. In this paper, we develop the central limit theorem for the Gini correlation in the more realistic setting where the dimensionality of the numerical variable is diverging. We then construct a powerful and consistent test for the K-sample problem based on the asymptotic normality. The proposed test not only avoids computation burden but also gains power over the permutation procedure. Simulation studies and real data illustrations show that the proposed test is more competitive to existing methods across a broad range of realistic situations, especially in unbalanced cases.

Motivated by applications in cybersecurity and epidemiology, we consider the problem of detecting an abrupt change in the intensity of a Poisson process, characterised by a jump (non transitory change) or a bump (transitory change) from constant. We propose a complete study from the nonasymptotic minimax testing point of view, when the constant baseline intensity is known or unknown. Starting from the totally agnostic case where all the parameters (height, location, length) of the change are unknown, the question of minimax adaptation with respect to each parameter is tackled, leading to a comprehensive overview of the various minimax separation rate regimes. We exhibit three such regimes and identify the factors of the two phase transitions by giving the cost of adaptation to each parameter. For any alternative hypothesis, depending on the knowledge or not of the change parameters, we propose minimax or minimax adaptive tests based on simple linear counting statistics derived from Neyman-Pearson tests, or quadratic statistics. When the change location or length is unknown, minimax adaptation is obtained from a scan aggregation principle combined with a Bonferroni or min-p level correction, and a conditioning trick when the baseline intensity is unknown.

We study the nonparametric estimation of the underlying survival function of a survival time in a study with cross-sectional sampling without any follow-up. Under a stationarity assumption on disease incidence rate in the population, the survival function $S_0$ is related to the observed density of the backward recurrence time, $f_0$, via the relationship $S_0(x)=f_0(x)∕f_0(0)$. As $f_0(x)$ is non-decreasing, it is well-known that the nonparametric maximum likelihood estimator of $f_0$ at $x=0$ is inconsistent. In this article, we establish the asymptotic distributions of the estimators of $S_0(x)$ when different consistent estimators of $f_0(0)$ are used. Such results are currently missing in the literature. Another contribution is the establishment of a local Kiefer-Wolfowitz-type result of the form ${sup}_{x\in [0,y]}|{\hat{F}}_n(x)-{\mathbb{F}}_n(x)|=O_p(n^{-2∕3}{(logn)}^{2∕3})$ that makes use of weaker assumptions than existing results, where ${\mathbb{F}}_n$ and ${\hat{F}}_n$ are the empirical distribution function and its least concave majorant, respectively.

This paper is concerned with inference on finite dimensional parameters in semiparametric moment condition models, where the moment functionals are linear with respect to unknown nuisance functions. By exploiting this linearity, we reformulate the inference problem via the Riesz representer, and develop a general inference procedure based on nonparametric likelihood. For treatment effect or missing data analysis, the Riesz representer is typically associated with the inverse propensity score even though the scope of our framework is much wider. In particular, we propose a two-step procedure, where the first step computes the projection weights to approximate the Riesz representer, and the second step reweights the moment conditions so that the likelihood increment admits an asymptotically pivotal chi-square calibration. Our reweighting method is naturally extended to inference on missing data, treatment effects, and data combination models, and other semiparametric problems. Simulation and real data examples illustrate usefulness of the proposed method. We note that our reweighting method and theoretical results are limited to linear functionals.

Recent works have proposed regression models which are invariant across data collection environments [24, 20, 11, 16, 8]. These estimators often have a causal interpretation under conditions on the environments and type of invariance imposed. One recent example, the Causal Dantzig (CD), is consistent under hidden confounding and represents an alternative to classical instrumental variable estimators such as Two Stage Least Squares (TSLS). In this work we derive the CD as a generalized method of moments (GMM) estimator. The GMM representation leads to several practical results, including 1) creation of the Generalized Causal Dantzig (GCD) estimator which can be applied to problems with continuous environments where the CD cannot be fit 2) a Hybrid (GCD-TSLS combination) estimator which has properties superior to GCD or TSLS alone 3) straightforward asymptotic results for all methods using GMM theory. We compare the CD, GCD, TSLS, and Hybrid estimators in simulations and an application to a Flow Cytometry data set. The newly proposed GCD and Hybrid estimators have superior performance to existing methods in many settings.

This paper studies a class of plug-in estimators of the stationary density of an autoregressive model with autoregression parameter $0<\mathit{\varrho }<1$. These use two types of estimator of the innovation density, a standard kernel estimator and a weighted kernel estimator with weights chosen to mimic the condition that the innovation density has mean zero. Bahadur expansions are obtained for this class of estimators in $L_1$, the space of integrable functions. These stochastic expansions establish root-n consistency in the $L_1$-norm. It is shown that the density estimators based on the weighted kernel estimators are asymptotically efficient if an asymptotically efficient estimator of the autoregression parameter is used. Here asymptotic efficiency is understood in the sense of the Hájek–Le Cam convolution theorem.

This paper deals with the grouped variable selection problem. A widely used strategy is to augment the negative log-likelihood function with a sparsity-promoting penalty. Existing methods include the group Lasso, group SCAD, and group MCP. The group Lasso solves a convex optimization problem but suffers from underestimation bias. The group SCAD and group MCP avoid this estimation bias but require solving a nonconvex optimization problem that may be plagued by suboptimal local optima. In this work, we propose an alternative method based on the generalized minimax concave (GMC) penalty, which is a folded concave penalty that maintains the convexity of the objective function. We develop a new method for grouped variable selection in linear regression, the group GMC, that generalizes the strategy of the original GMC estimator. We present a primal-dual algorithm for computing the group GMC estimator and also prove properties of the solution path to guide its numerical computation and tuning parameter selection in practice. We establish error bounds for both the group GMC and original GMC estimators. A rich set of simulation studies and a real data application indicate that the proposed group GMC approach outperforms existing methods in several different aspects under a wide array of scenarios.

Temporal dependence is frequently encountered in large-scale structured noisy data, arising from scientific studies in neuroscience and meteorology, among others. This challenging characteristic may not align with existing theoretical frameworks or data analysis tools. Motivated by multi-session $\mathrm{fMRI}$ time series data, this paper introduces a novel semi-parametric inference procedure suitable for a broad class of “non-stationary, non-Gaussian, temporally dependent” noise processes in time-course data. It develops a new test statistic based on a tapering-type estimator of the large-dimensional noise auto-covariance matrix and establishes its asymptotic chi-squared distribution. Our method not only relaxes the consistency requirement for the noise covariance matrix estimator but also avoids direct matrix inversion without sacrificing detection power. It adapts well to both stationary and a wider range of temporal noise processes, making it particularly effective for handling challenging scenarios involving very large scales of data and large dimensions of noise covariance matrices. We demonstrate the efficacy of the proposed procedure through simulation evaluations and real $\mathrm{fMRI}$ data analysis.

Heterogeneous effect estimation is crucial in causal inference, with applications across medicine and social science. Many methods for estimating conditional average treatment effects (CATEs) have been proposed, but there are gaps in understanding if and when such methods are optimal. This is especially true when the CATE has nontrivial structure (e.g., smoothness or sparsity). Our work contributes in several ways. First, we study a two-stage doubly robust CATE estimator and give a generic error bound, which yields rates faster than much of the literature. We apply the bound to derive error rates in smooth nonparametric models, and give sufficient conditions for oracle efficiency. Along the way we give a general error bound for regression with estimated outcomes; this is the second main contribution. The third contribution is aimed at understanding the fundamental statistical limits of CATE estimation. To that end, we propose and study a local polynomial adaptation of double-residual regression. We show that this estimator can be oracle efficient under even weaker conditions, and we conjecture that they are minimal in a minimax sense. We go on to give error bounds in the non-trivial regime where oracle rates cannot be achieved. Some finite-sample properties are explored with simulations.

Given a zero-mean Gaussian random field with a covariance function that belongs to a parametric family of covariance functions, we introduce a new notion of likelihood approximations, termed truncated-likelihood functions. Truncated-likelihood functions are based on direct functional approximations of the presumed family of covariance functions. For compactly supported covariance functions, within an increasing-domain asymptotic framework, we provide sufficient conditions under which consistency and asymptotic normality of estimators based on truncated-likelihood functions are preserved. We apply our result to the family of generalized Wendland covariance functions and discuss several examples of Wendland approximations. For families of covariance functions that are not compactly supported, we combine our results with the covariance tapering approach and show that ML estimators, based on truncated-tapered likelihood functions, asymptotically minimize the Kullback-Leibler divergence, when the taper range is fixed.

We consider varying-coefficient models for mixed synchronous and asynchronous longitudinal covariates, where asynchronicity refers to the misalignment of longitudinal measurement times within an individual. We propose three different methods of parameter estimation and inference. The first method is a one-step approach that estimates non-parametric regression functions for synchronous and asynchronous longitudinal covariates simultaneously. The second method is a two-step approach in which synchronous longitudinal covariates are regressed with the longitudinal response by centering the synchronous longitudinal covariates first and, in the second step, the residuals from the first step are regressed with asynchronous longitudinal covariates. The third method is the same as the second method except that in the first step, we omit the asynchronous longitudinal covariate and include a non-parametric intercept in the regression analysis of synchronous longitudinal covariates and the longitudinal response. We further construct simultaneous confidence bands for the non-parametric regression functions to quantify the overall magnitude of variation. Extensive simulation studies provide numerical support for the theoretical findings. The practical utility of the methods is illustrated on a dataset from the ADNI study.

Hypothesis testing procedures are developed to assess linear operator constraints in function-on-scalar regression when incomplete functional responses are observed. The approach enables statistical inferences about the shape and other aspects of the functional regression coefficients within a unified framework encompassing three incomplete sampling scenarios; (i) partially observed response functions as curve segments over random sub-intervals of the domain, (ii) discretely observed functional responses with additive measurement errors, and (iii) the composition of former two scenarios, where partially observed response segments are observed discretely with measurement error. The latter scenario has been little explored to date, although such structured data is increasingly common in applications. For statistical inference, deviations from the constraint space are measured via integrated $L^2$-distance between estimates from the constrained and unconstrained model spaces. Large sample properties of the proposed test procedure are established, including the consistency, asymptotic distribution, and local power of the test statistic. The finite sample power and level of the proposed test are investigated in a simulation study covering a variety of scenarios. The proposed methodologies are illustrated by applications to U.S. obesity prevalence data, analyzing the functional shape of its trends over time, and motion analysis in a study of automotive ergonomics.

The ratio of the hazard functions of two populations or two strata of a single population plays an important role in time-to-event analysis. Cox regression is commonly used to estimate the hazard ratio under the assumption that it is constant in time, which is known as the proportional hazards assumption. However, this assumption is often violated in practice, and when it is violated, the parameter estimated by Cox regression is difficult to interpret. The hazard ratio can be estimated in a nonparametric manner using smoothing, but smoothing-based estimators are sensitive to the selection of tuning parameters, and it is often difficult to perform valid inference with such estimators. In some cases, it is known that the hazard ratio function is monotone. In this article, we demonstrate that monotonicity of the hazard ratio function defines an invariant stochastic order, and we study the properties of this order. Furthermore, we introduce an estimator of the hazard ratio function under a monotonicity constraint. We demonstrate that our estimator converges in distribution to a mean-zero limit, and we use this result to construct asymptotically valid confidence intervals. Finally, we conduct numerical studies to assess the finite-sample behavior of our estimator, and we use our methods to estimate the hazard ratio of progression-free survival in pulmonary adenocarcinoma patients treated with gefitinib or carboplatin-paclitaxel.

A common principle in model diagnostics and forecast evaluation is that fitted or predicted distributions ought to be reliable, ideally in the sense of auto-calibration, where the outcome is a random draw from the posited distribution. For binary responses, auto-calibration is the universal concept of reliability. For real-valued outcomes, a general theory of calibration has been elusive, despite a recent surge of interest in distributional regression and machine learning. We develop a framework rooted in probability theory, which gives rise to hierarchies of calibration, and applies to both predictive distributions and stand-alone point forecasts. In a nutshell, a prediction is conditionally T-calibrated if it can be taken at face value in terms of an identifiable functional T. We introduce population versions of T-reliability diagrams and revisit a score decomposition into measures of miscalibration, discrimination, and uncertainty. In empirical settings, stable and efficient estimators of T-reliability diagrams and score components arise via nonparametric isotonic regression and the pool-adjacent-violators algorithm. For in-sample model diagnostics, we propose a universal coefficient of determination that nests and reinterprets the classical $R^2$ in least squares regression and its natural analog $R^1$ in quantile regression, yet applies to T-regression in general.

Spatial analysis commonly relies on the estimation of a covariance matrix associated with a random field. This estimation strongly impacts the prediction where the process has not been observed, which in turn influences the construction of more sophisticated models. If some of the distances between all the possible pairs of observations in the plane are small, then we may have an ill-conditioned problem that results in a nearly singular covariance matrix. In this paper, we suggest a covariance matrix estimation method that works well even when there are very close pairs of locations on the plane. Our method is an extension to a spatial case of a method that is based on the estimation of eigenvalues of the unitary matrix decomposition of the covariance matrix. Several numerical examples are conducted to provide evidence of good performance in estimating the range parameter of the correlation structure of a spatial regression process. In addition, an application to macroalgae estimation in a restricted area of the Pacific Ocean is developed to determine a suitable estimation of the effective sample size associated with the transect sampling scheme.

In this paper, we consider a simple estimator for tail dependence coefficients of a max-stable time series and show its asymptotic normality under a mild condition. The novelty of our result is that this condition does not involve mixing properties that are common in the literature. More importantly, our condition is linked to the transition between long and short range dependence (LRD/SRD) for max-stable time series. This is based on a recently proposed notion of LRD in the sense of indicators of excursion sets which is meaningfully defined for infinite-variance time series. In particular, we show that asymptotic normality with standard rate of convergence and a function of the sum of tail coefficients as asymptotic variance holds if and only if the max-stable time series is SRD.

The problem of determining whether a given point, or set of points, lies within the convex hull of another set of points in d dimensions arises naturally in the context of certain exponential family models in statistics. This article discusses the general convex hull problem and its application to the particular problem of modelling network data using an exponential-family random graph model (ERGM). While the convex hull question may be solved via a simple linear program, this approach is not well known in the statistical literature. The article also details several substantial improvements to the convex hull-testing algorithm currently implemented in the widely used ergm package for network modeling. It provides direct numerical comparisons of two linear programming packages for R that can be called by ergm and offers several illustrative examples.

It is more and more frequently the case in applications that the data we observe come from one or more random variables taking values in an infinite dimensional space, e.g. curves. The need to have tools adapted to the nature of these data explains the growing interest in the field of functional data analysis. The model we study in this paper assumes a linear dependence between a quantity of interest and several covariates, at least one of which has an infinite dimension. To select the relevant covariates in this context, we investigate adaptations of the Lasso method. Two estimation methods are defined. The first one consists in the minimization of a Group-Lasso criterion on the multivariate functional space H. The second one minimizes the same criterion but on a finite dimensional subspaces of H whose dimension is chosen by a penalized least squares method. We prove oracle inequalities of sparsity in the case where the design is fixed or random. To compute the solutions of both criteria in practice, we propose a coordinate descent algorithm. A numerical study on simulated and real data illustrates the behavior of the estimators.

How can increasingly available observational data be used to improve the design of randomized controlled trials (RCTs)? We seek to design a prospective RCT, with the intent of using an Empirical Bayes estimator to shrink the causal estimates from our trial toward causal estimates obtained from an observational study. We ask: how might we design the experiment to better complement the observational study in this setting?

We show that the risk of such shrinkage estimators can be computed efficiently via numerical integration. We then propose three algorithms for determining the best allocation of units to strata given the estimator’s plannned use: Neyman allocation; a “naïve” design assuming no unmeasured confounding in the observational study; and a robust design accounting for the imperfect parameter estimates we would obtain from the observational study with unmeasured confounding. We propose guardrails on the designs, so that our experiment could be reasonably analyzed without shrinkage if desired.

We demonstrate the viability of these experimental designs through a simulation study involving a rare, binary outcome. Lastly, we deploy our methods on real data from the Women’s Health Initiative, a 1991 study estimating the health effects of hormone therapy on postmenopausal women. In particular, we determine how many units should be allocated to each treatment arm in each stratum of interest in order to maximally reduce estimation risk given the planned use of the shrinkage estimator. We find improved design provides further benefits over and above the benefit of the shrinkage estimator itself.

In this paper we develop an online statistical inference approach for high-dimensional generalized linear models with streaming data for real-time estimation and inference. We propose an online debiased lasso method that aligns with the data collection scheme of streaming data. Online debiased lasso differs from offline debiased lasso in two important aspects. First, it updates component-wise confidence intervals of regression coefficients with only summary statistics of the historical data. Second, online debiased lasso adds an additional term to correct approximation errors accumulated throughout the online updating procedure. We show that our proposed online debiased estimators in generalized linear models are asymptotically normal. This result provides a theoretical basis for carrying out real-time interim statistical inference with streaming data. Extensive numerical experiments are conducted to evaluate the performance of our proposed online debiased lasso method. These experiments demonstrate the effectiveness of our algorithm and support the theoretical results. Furthermore, we illustrate the application of our method with a high-dimensional text dataset.

This paper introduces a new class of space-time integer-valued ARMA models referred to as STINARMA. This class arises as the natural space-time extension of the INARMA models and, simultaneously, as the integer-valued counterpart of the conventional STARMA models. In this work, the moving average subclass STINMA$(q_{m_1,\dots ,m_q})$ is studied in detail. Particular attention is given to the derivation of first- and second-order moments, including space-time autocorrelations. Due to its large potential use in real-data applications, the Poisson STINMA($1_1$) process is analyzed in further detail. Estimation methods are also addressed and their performance is demonstrated through a simulation study and by analysing the daily number of hospital admissions observed over time in three Portuguese locations.

We study the problem of detecting latent geometric structure in random graphs. To this end, we consider the soft high-dimensional random geometric graph $G(n,p,d,q)$, where each of the n vertices corresponds to an independent random point distributed uniformly on the sphere ${\mathbb{S}}^{d-1}$, and the probability that two vertices are connected by an edge is a decreasing function of the Euclidean distance between the points. The probability of connection is parametrized by $q\in [0,1]$, with smaller q corresponding to weaker dependence on the geometry; this can also be interpreted as the level of noise in the geometric graph. In particular, the model smoothly interpolates between the hard spherical random geometric graph $G(n,p,d)$ (corresponding to $q=1$) and the Erdős–Rényi model $G(n,p)$ (corresponding to $q=0$). We focus on the dense regime (i.e., p is a constant).

We show that if $nq\to 0$ or $d»n^3{q}^2$, then geometry is lost: $G(n,p,d,q)$ is asymptotically indistinguishable from $G(n,p)$. On the other hand, if $d«n^3{q}^6$, then the signed triangle statistic provides an asymptotically powerful test for detecting geometry. These results generalize those of Bubeck, Ding, Eldan, and Rácz (2016) for $G(n,p,d)$, and give quantitative bounds on how the noise level affects the dimension threshold for losing geometry. We also prove analogous results under a related but different distributional assumption which corresponds to the random dot product graph.

We develop sampling algorithms to fit Bayesian hierarchical models, the computational complexity of which scales linearly with the number of observations and the number of parameters in the model. We focus on crossed random effect and nested multilevel models, which are used ubiquitously in applied sciences. The posterior dependence in both classes is sparse: in crossed random effects models it resembles a random graph, whereas in nested multilevel models it is tree-structured. For each class we identify a framework for scalable computation, building on previous work. Methods for crossed models are based on extensions of appropriately designed collapsed Gibbs samplers, where we introduce the idea of local centering; while methods for nested models are based on sparse linear algebra (SLA) and data augmentation. We provide a theoretical analysis of the proposed algorithms in some simplified settings, including a comparison with previously proposed methodologies and an average-case analysis based on random graph theory. Numerical experiments, including two challenging real data analyses on predicting electoral results and real estate prices, compare with off-the-shelf Hamiltonian Monte Carlo, displaying drastic improvement in performance. The code for replicating the experiments in the article and for implementing our methods can be found at https://github.com/timsf/crossed-effects and https://github.com/timsf/nested-effects.

We establish optimal convergence rates up to a $\log$ factor for a class of deep neural networks in a classification setting under a restraint sometimes referred to as the Tsybakov noise condition. We construct classifiers based on empirical risk minimization in a general setting where the boundary of the Bayes rule can be approximated well by neural networks. Corresponding rates of convergence are proven with respect to the misclassification error using an additional condition that acts as a requirement for the “correct noise exponent”. It is then shown that these rates are optimal in the minimax sense. For other estimation procedures, similar convergence rates have been established. Our first main contribution is to prove that the rates are optimal under the additional condition. Secondly, our main theorem establishes almost optimal rates in a generalized setting. We use this to show optimal rates which circumvent the curse of dimensionality.

Estimation of a conditional mean (linking a set of features to an outcome of interest) is a fundamental statistical task. While there is an appeal to flexible nonparametric procedures, effective estimation in many classical nonparametric function spaces, e.g., multivariate Sobolev spaces, can be prohibitively difficult – both statistically and computationally – especially when the number of features is large. In this paper, we present some sieve estimators for regression in multivariate product spaces. We take Sobolev-type smoothness spaces as an example, though our general framework can be applied to many reproducing kernel Hilbert spaces. These spaces are more amenable to multivariate regression, and allow us to, in-part, avoid the curse of dimensionality. Our estimator can be easily applied to multivariate nonparametric problems and has appealing statistical and computational properties. Moreover, it can effectively leverage additional structure such as feature sparsity.

The setting where an unknown number m of the largest data is missing from an underlying Pareto-type distribution is considered. Solutions are provided for estimating the extreme value index, the number of missing data and extreme quantiles. Asymptotic results of the parameter estimators and an adaptive selection method for the number of top data used in the estimation are proposed for the case where all missing data are beyond the observed data. An estimator of the number of missing extremes spread over the largest observed data is also proposed. To this purpose, a key component is a likelihood solution based on exponential representations of spacings between the largest observations. An effective and fast optimization procedure is established using regularization, and simulation experiments are provided. The methodology is illustrated with a dataset from the diamond mining industry, where large-carat diamonds are expected to be missing.

We investigate a general matrix factorization for deviance-based data losses, extending the ubiquitous singular value decomposition beyond squared error loss. While similar approaches have been explored before, our method leverages classical statistical methodology from generalized linear models (GLMs) and provides an efficient algorithm that is flexible enough to allow for structural zeros via entry weights. Moreover, by adapting results from GLM theory, we provide support for these decompositions by (i) showing strong consistency under the GLM setup, (ii) checking the adequacy of a chosen exponential family via a generalized Hosmer-Lemeshow test, and (iii) determining the rank of the decomposition via a maximum eigenvalue gap method. To further support our findings, we conduct simulation studies to assess robustness to decomposition assumptions and extensive case studies using benchmark datasets from image face recognition, natural language processing, network analysis, and biomedical studies. Our theoretical and empirical results indicate that the proposed decomposition is more flexible, general, and robust, and can thus provide improved performance when compared to similar methods. To facilitate applications, an R package dmf with efficient model fitting and family and rank determination is also provided.

The variance of noise plays an important role in many change-point detection procedures and the associated inferences. Most commonly used variance estimators require strong assumptions on the true mean structure or normality of the error distribution, which may not hold in applications. More importantly, the qualities of these estimators have not been discussed systematically in the literature. In this paper, we introduce a framework of equivariant variance estimation for multiple change-point models. In particular, we characterize the set of all equivariant unbiased quadratic variance estimators for a family of change-point model classes, and develop a minimax theory for such estimators.

Because it determines a center-outward ordering of observations in ${\mathbb{R}}^d$ with $d\ge 2$, the concept of statistical depth permits to define quantiles and ranks for multivariate data and use them for various statistical tasks (e.g. inference, hypothesis testing). Whereas many depth functions have been proposed ad-hoc in the literature since the seminal contribution of [50], not all of them possess the properties desirable to emulate the notion of quantile function for univariate probability distributions. In this paper, we propose an extension of the integrated rank-weighted statistical depth (IRW depth in abbreviated form) originally introduced in [40], modified in order to satisfy the property of affine invariance, fulfilling thus all the four key axioms listed in the nomenclature elaborated by [59]. The variant we propose, referred to as the affine invariant IRW depth (AI-IRW in short), involves the precision matrix of the (supposedly square integrable) d-dimensional random vector X under study, in order to take into account the directions along which X is most variable to assign a depth value to any point $x\in {\mathbb{R}}^d$. The accuracy of the sampling version of the AI-IRW depth is investigated from a non-asymptotic perspective. Namely, a concentration result for the statistical counterpart of the AI-IRW depth is proved. Beyond the theoretical analysis carried out, applications to anomaly detection are considered and numerical results are displayed, providing strong empirical evidence of the relevance of the depth function we propose here.

We propose a bootstrap-based test to detect a mean shift in a sequence of high-dimensional observations with unknown time-varying heteroscedasticity. The proposed test builds on the U-statistic based approach in Wang et al. (2022), targets a dense alternative, and adopts a wild bootstrap procedure to generate critical values. The bootstrap-based test is free of tuning parameters and is capable of accommodating unconditional time varying heteroscedasticity in the high-dimensional observations, as demonstrated in our theory and simulations. Theoretically, we justify the bootstrap consistency by using the recently proposed unconditional approach in Bücher and Kojadinovic (2019). Extensions to testing for multiple change-points and estimation using wild binary segmentation are also presented. Numerical simulations demonstrate the robustness of the proposed testing and estimation procedures with respect to different kinds of time-varying heteroscedasticity.

Computer experiments are increasingly being used to build high-quality surrogate models for complex emulation systems. Maximin distance Latin hypercube design is an efficient approach for designing computer experiments. Algorithmic search is commonly used for finding such designs but becomes ineffective when searching for large designs. Theoretical construction of such designs is fast but limited and challenging. In this paper, we propose a series of construction methods for maximin distance Latin hypercube designs. We use a piece-wise linear transformation to obtain balanced designs which are then rotated to generate Latin hypercube designs. Theoretical results guarantee that the generated designs are asymptotically optimal under the maximin distance criterion. The generated designs also exhibit low column correlations and mirror symmetry, which significantly benefits the identification of the main and interaction effects. Moreover, we present numerical comparisons with existing methods to demonstrate the superiority of the proposed methods.

We define Gaussian graphical models on directed acyclic graphs with coloured vertices and edges, calling them RDAG (restricted directed acyclic graph) models. If two vertices or edges have the same colour, their parameters in the model must be the same. We present an algorithm to find the maximum likelihood estimate (MLE) in an RDAG model, and characterise when the MLE exists, via linear independence conditions. We relate properties of a graph, and its colouring, to the number of samples needed for the MLE to exist and to be unique. We also characterise when an RDAG model is equal to an associated undirected graphical model and study connections to groups and invariant theory. We provide examples and simulations to study the benefits of RDAGs over uncoloured DAGs.

In this paper, we deal with nonparametric regression for circular data, meaning that observations are represented by points lying on the unit circle. We propose a kernel estimation procedure with data-driven selection of the bandwidth parameter. For this purpose, we use a warping strategy combined with a Goldenshluger-Lepski type estimator. To study optimality of our methodology, we consider the minimax setting and prove, by establishing upper and lower bounds, that our procedure is nearly optimal on anisotropic Hölder classes of functions for pointwise estimation. The obtained rates also reveal the specific nature of regression for circular responses. Finally, a numerical study is conducted, illustrating the good performances of our approach.