The determination of an optimal design for a given regression problem is an intricate optimization problem, especially for models with multivariate predictors. Design admissibility and invariance are main tools to reduce the complexity of the optimization problem and have been successfully applied for models with univariate predictors. In particular, several authors have developed sufficient conditions for the existence of minimally supported designs in univariate models, where the number of support points of the optimal design equals the number of parameters. These results generalize the celebrated de la Garza phenomenon (Ann. Math. Statistics 25 (1954) 123–130), which states that for a polynomial regression model of degree $\mathit{k}-1$ any optimal design can be based on k points.

This paper provides—for the first time—extensions of these results for models with a multivariate predictor. In particular, we study a geometric characterization of the support points of an optimal design to provide sufficient conditions for the occurrence of the de la Garza phenomenon in models with multivariate predictors and characterize properties of admissible designs in terms of admissibility of designs in conditional univariate regression models.

Valid instrumental variables enable treatment effect inference even when selection into treatment is biased by unobserved confounders. When multiple candidate instruments are available, but some of them are possibly invalid, the previously proposed reinforced design enables one or more nearly independent valid analyses that depend on very different assumptions. That is, we can perform evidence factor analysis. However, the validity of the reinforced design depends crucially on the order in which multiple instrumental variable analyses are conducted. Motivated by the orthogonality of balanced factorial designs, we propose a balanced block design to offset the possible violation of the exclusion restriction by balancing the instruments against each other in the design, and demonstrate its utility for constructing approximate evidence factors under multiple analysis strategies free of the order imposition. We also propose a novel stratification method using multiple, nested candidate instruments, in which case the balanced block design is not applicable. We apply our proposed methods to evaluate (a) the effect of education on future earnings using instrumental variables arising from the disruption of education during World War II via the balanced block design, and (b) the causal effect of malaria on stunting among children in Western Kenya using three nested instruments.

Estimating the orders of the autoregressive fractionally integrated moving average (ARFIMA) model has been a long-standing problem in time series analysis. This paper tackles this challenge by establishing the consistency of the Bayesian information criterion (BIC) for ARFIMA models with independent errors. Since the memory parameter of the model can be any real number, this consistency result is valid for short memory, long memory and nonstationary time series. This paper further extends the consistency of the BIC to ARFIMA models with conditional heteroscedastic errors, thereby extending its applications to encompass many real-life situations. Finite-sample implications of the theoretical results are illustrated via numerical examples.

Inferring causal relationships or related associations from observational data can be invalidated by the existence of hidden confounding. We focus on a high-dimensional linear regression setting, where the measured covariates are affected by hidden confounding and propose the doubly debiased lasso estimator for individual components of the regression coefficient vector. Our advocated method simultaneously corrects both the bias due to estimation of high-dimensional parameters as well as the bias caused by the hidden confounding. We establish its asymptotic normality and also prove that it is efficient in the Gauss–Markov sense. The validity of our methodology relies on a dense confounding assumption, that is, that every confounding variable affects many covariates. The finite sample performance is illustrated with an extensive simulation study and a genomic application. The method is implemented by the $\mathtt{DDL}$ package available from CRAN.

While deep learning approaches to survival data have demonstrated empirical success in applications, most of these methods are difficult to interpret and mathematical understanding of them is lacking. This paper studies the partially linear Cox model, where the nonlinear component of the model is implemented using a deep neural network. The proposed approach is flexible and able to circumvent the curse of dimensionality, yet it facilitates interpretability of the effects of treatment covariates on survival. We establish asymptotic theories of maximum partial likelihood estimators and show that our nonparametric deep neural network estimator achieves the minimax optimal rate of convergence (up to a polylogarithmic factor). Moreover, we prove that the corresponding finite-dimensional estimator for treatment covariate effects is $\sqrt{\mathit{n}}$-consistent, asymptotically normal and attains semiparametric efficiency. Extensive simulation studies and analyses of two real survival data sets show the proposed estimator produces confidence intervals with superior coverage as well as survival time predictions with superior concordance to actual survival times.

Variance estimation is important for statistical inference. It becomes nontrivial when observations are masked by serial dependence structures and time-varying mean structures. Existing methods either ignore or sub-optimally handle these nuisance structures. This paper develops a general framework for the estimation of the long-run variance for time series with nonconstant means. The building blocks are difference statistics. The proposed class of estimators is general enough to cover many existing estimators. Necessary and sufficient conditions for consistency are investigated. The first asymptotically optimal estimator is derived. Our proposed estimator is theoretically proven to be invariant to arbitrary mean structures, which may include trends and a possibly divergent number of discontinuities.

The success of the Lasso in the era of high-dimensional data can be attributed to its conducting an implicit model selection, that is, zeroing out regression coefficients that are not significant. By contrast, classical ridge regression cannot reveal a potential sparsity of parameters, and may also introduce a large bias under the high-dimensional setting. Nevertheless, recent work on the Lasso involves debiasing and thresholding, the latter in order to further enhance the model selection. As a consequence, ridge regression may be worth another look since—after debiasing and thresholding—it may offer some advantages over the Lasso, for example, it can be easily computed using a closed-form expression. In this paper, we define a debiased and thresholded ridge regression method, and prove a consistency result and a Gaussian approximation theorem. We further introduce a wild bootstrap algorithm to construct confidence regions and perform hypothesis testing for a linear combination of parameters. In addition to estimation, we consider the problem of prediction, and present a novel, hybrid bootstrap algorithm tailored for prediction intervals. Extensive numerical simulations further show that the debiased and thresholded ridge regression has favorable finite sample performance and may be preferable in some settings.

Given n independent random vectors with common density f on ${\mathbb{R}}^{\mathit{d}}$, we study the weak convergence of three empirical-measure based estimators of the convex λ-level set ${\mathit{L}}_{\mathit{\lambda }}$ of f, namely the excess mass set, the minimum volume set and the maximum probability set, all selected from a class of convex sets $\mathcal{A}$ that contains ${\mathit{L}}_{\mathit{\lambda }}$. Since these set-valued estimators approach ${\mathit{L}}_{\mathit{\lambda }}$, even the formulation of their weak convergence is nonstandard. We identify the joint limiting distribution of the symmetric difference of ${\mathit{L}}_{\mathit{\lambda }}$ and each of the three estimators, at rate ${\mathit{n}}^{-1/3}$. It turns out that the minimum volume set and the maximum probability set estimators are asymptotically indistinguishable, whereas the excess mass set estimator exhibits “richer” limit behavior. Arguments rely on the boundary local empirical process, its cylinder representation, dimension-free concentration around the boundary of ${\mathit{L}}_{\mathit{\lambda }}$, and the set-valued argmax of a drifted Wiener process.

We adapt Higher Criticism (HC) to the comparison of two frequency tables which may—or may not—exhibit moderate differences between the tables in some unknown, relatively small subset out of a large number of categories. Our analysis of the power of the proposed HC test quantifies the rarity and size of assumed differences and applies moderate deviations-analysis to determine the asymptotic powerfulness/powerlessness of our proposed HC procedure.

Our analysis considers the null hypothesis of no difference in underlying generative model against a rare/weak perturbation alternative, in which the frequencies of ${\mathit{N}}^{1-\mathit{\beta }}$ out of the N categories are perturbed by $\mathit{r}(log\mathit{N})/2\mathit{n}$ in the Hellinger distance; here, n is the size of each sample. Our proposed Higher Criticism (HC) test for this setting uses P-values obtained from N exact binomial tests. We characterize the asymptotic performance of the HC-based test in terms of the sparsity parameter β and the perturbation intensity parameter r. Specifically, we derive a region in the $(\mathit{\beta },\mathit{r})$-plane where the test asymptotically has maximal power, while having asymptotically no power outside this region. Our analysis distinguishes between cases in which the counts in both tables are low, versus cases in which counts are high, corresponding to the cases of sparse and dense frequency tables. The phase transition curve of HC in the high-counts regime matches formally the curve delivered by HC in a two-sample normal means model.

Orthogonal array, a classical and effective tool for collecting data, has been flourished with its applications in modern computer experiments and engineering statistics. Driven by the wide use of computer experiments with both qualitative and quantitative factors, multiple computer experiments, multifidelity computer experiments, cross-validation and stochastic optimization, orthogonal arrays with certain structures have been introduced. Sliced orthogonal arrays and nested orthogonal arrays are examples of such arrays. This article introduces a flexible, fresh construction method, which uses smaller arrays and a special structure. The method uncovers the hidden structure of many existing fixed-level orthogonal arrays of given run sizes, possibly with more columns. It also allows fixed-level orthogonal arrays of nearly strength three to be constructed, which are useful as there are not many construction methods for fixed-level orthogonal arrays of strength three, and also helpful for generating Latin hypercube designs with desirable low-dimensional projections. Theoretical properties of the proposed method are explored. As by-products, several theoretical results on orthogonal arrays are obtained.

Predicting extreme events is important in many applications in risk analysis. Extreme-value theory suggests modelling extremes by max-stable distributions. The Bayesian approach provides a natural framework for statistical prediction. Although various Bayesian inferential procedures have been proposed in the literature of univariate extremes and some for multivariate extremes, the study of their asymptotic properties has been left largely untouched. In this paper we focus on a semiparametric Bayesian method for estimating max-stable distributions in arbitrary dimension. We establish consistency of the pertaining posterior distributions for fairly general, well-specified max-stable models, whose margins can be short-, light- or heavy-tailed. We then extend our consistency results to the case where data are samples of block maxima whose distribution is only approximately a max-stable one, which represents the most realistic inferential setting.

We study nonparametric density estimation problems where error is measured in the Wasserstein distance, a metric on probability distributions popular in many areas of statistics and machine learning. We give the first minimax-optimal rates for this problem for general Wasserstein distances for two classes of densities: smooth probability densities on ${[0,1]}^{\mathit{d}}$ bounded away from 0, and sub-Gaussian densities lying in the Hölder class ${\mathcal{C}}^{\mathit{s}}$, $\mathit{s}\in (0,1)$. Unlike classical nonparametric density estimation, these rates depend on whether the densities in question are bounded below, even in the compactly supported case. Motivated by variational problems involving the Wasserstein distance, we also show how to construct discretely supported measures, suitable for computational purposes, which achieve the minimax rates. Our main technical tool is an inequality giving a nearly tight dual characterization of the Wasserstein distances in terms of Besov norms.

The detection of weak and rare effects in large amounts of data arises in a number of modern data analysis problems. Known results show that in this situation the potential of statistical inference is severely limited by the large-scale multiple testing that is inherent in these problems. Here, we show that fundamentally more powerful statistical inference is possible when there is some structure in the signal that can be exploited, for example, if the signal is clustered in many small blocks, as is the case in some relevant applications. We derive the detection boundary in such a situation where we allow both the number of blocks and the block length to grow polynomially with sample size. We derive these results both for the univariate and the multivariate settings as well as for the problem of detecting clusters in a network. These results recover as special cases the sparse signal detection problem (Ann. Statist. 32 (2004) 962–994) where there is no structure in the signal, as well as the scan problem (Statist. Sinica 23 (2013) 409–428) where the signal comprises a single interval. We develop methodology that allows optimal adaptive detection in the general setting, thus exploiting the structure if it is present without incurring a relevant penalty in the case where there is no structure. The advantage of this methodology can be considerable, as in the case of no structure the means need to increase at the rate $\sqrt{log\mathit{n}}$ to ensure detection, while the presence of structure allows detection even if the means decrease at a polynomial rate.

Local Fréchet regression is a nonparametric regression method for metric space valued responses and Euclidean predictors, which can be utilized to obtain estimates of smooth trajectories taking values in general metric spaces from noisy metric space valued random objects. We derive uniform rates of convergence, which so far have eluded theoretical analysis of this method, for both fixed and random target trajectories, where we utilize tools from empirical processes. These results are shown to be widely applicable in metric space valued data analysis. In addition to simulations, we provide two pertinent examples where these results are important: The consistent estimation of the location of properly defined extrema in metric space valued trajectories, which we illustrate with the problem of locating the age of minimum brain connectivity as obtained from fMRI data; and time warping for metric space valued trajectories, illustrated with yearly age-at-death distributions for different countries.

The paper analyzes cointegration in vector autoregressive processes (VARs) for the cases when both the number of coordinates, N, and the number of time periods, T, are large and of the same order. We propose a way to examine a VAR of order 1 for the presence of cointegration based on a modification of the Johansen likelihood ratio test. The advantage of our procedure over the original Johansen test and its finite sample corrections is that our test does not suffer from overrejection. This is achieved through novel asymptotic theorems for eigenvalues of matrices in the test statistic in the regime of proportionally growing N and T. Our theoretical findings are supported by Monte Carlo simulations and an empirical illustration. Moreover, we find a surprising connection with multivariate analysis of variance (MANOVA) and explain why it emerges.

Given partially observed pairwise comparison data generated by the Bradley–Terry–Luce (BTL) model, we study the problem of top-k ranking. That is, to optimally identify the set of top-k players. We derive the minimax rate with respect to a normalized Hamming loss. This provides the first result in the literature that characterizes the partial recovery error in terms of the proportion of mistakes for top-k ranking. We also derive the optimal signal to noise ratio condition for the exact recovery of the top-k set. The maximum likelihood estimator (MLE) is shown to achieve both optimal partial recovery and optimal exact recovery. On the other hand, we show another popular algorithm, the spectral method, is in general suboptimal. Our results complement the recent work (Ann. Statist. 47 (2019) 2204–2235) that shows both the MLE and the spectral method achieve the optimal sample complexity for exact recovery. It turns out the leading constants of the sample complexity are different for the two algorithms. Another contribution that may be of independent interest is the analysis of the MLE without any penalty or regularization for the BTL model. This closes an important gap between theory and practice in the literature of ranking.

Given a random sample of points from some unknown density, we propose a method for estimating density level sets, for a positive threshold t, under the r-convexity assumption. This shape condition generalizes the convexity property and allows to consider level sets with more than one connected component. The main problem in practice is that r is an unknown geometric characteristic of the set related to its curvature, which may depend on t. A stochastic algorithm is proposed for selecting its value from data. The resulting reconstruction of the level set is able to achieve minimax rates for Hausdorff metric and distance in measure uniformly on the level t.

This paper establishes a precise high-dimensional asymptotic theory for boosting on separable data, taking statistical and computational perspectives. We consider a high-dimensional setting where the number of features (weak learners) p scales with the sample size n, in an overparametrized regime. Under a class of statistical models, we provide an exact analysis of the generalization error of boosting when the algorithm interpolates the training data and maximizes the empirical ${\ell }_1$-margin. Further, we explicitly pin down the relation between the boosting test error and the optimal Bayes error, as well as the proportion of active features at interpolation (with zero initialization). In turn, these precise characterizations answer certain questions raised in (Neural Comput. 11 (1999) 1493–1517; Ann. Statist. 26 (1998) 1651–1686) surrounding boosting, under assumed data generating processes. At the heart of our theory lies an in-depth study of the maximum-${\ell }_1$-margin, which can be accurately described by a new system of nonlinear equations; to analyze this margin, we rely on Gaussian comparison techniques and develop a novel uniform deviation argument. Our statistical and computational arguments can handle (1) any finite-rank spiked covariance model for the feature distribution and (2) variants of boosting corresponding to general ${\ell }_{\mathit{q}}$-geometry, $\mathit{q}\in [1,2]$. As a final component, via the Lindeberg principle, we establish a universality result showcasing that the scaled ${\ell }_1$-margin (asymptotically) remains the same, whether the covariates used for boosting arise from a nonlinear random feature model or an appropriately linearized model with matching moments.

A new framework is developed to intrinsically analyze sparsely observed Riemannian functional data. It features four innovative components: a frame-independent covariance function, a smooth vector bundle termed covariance vector bundle, a parallel transport and a smooth bundle metric on the covariance vector bundle. The introduced intrinsic covariance function links estimation of covariance structure to smoothing problems that involve raw covariance observations derived from sparsely observed Riemannian functional data, while the covariance vector bundle provides a rigorous mathematical foundation for formulating such smoothing problems. The parallel transport and the bundle metric together make it possible to measure fidelity of fit to the covariance function. They also play a critical role in quantifying the quality of estimators for the covariance function. As an illustration, based on the proposed framework, we develop a local linear smoothing estimator for the covariance function, analyze its theoretical properties and provide numerical demonstration via simulated and real data sets. The intrinsic feature of the framework makes it applicable to not only Euclidean submanifolds but also manifolds without a canonical ambient space.

This paper presents a parametric estimation method for ill-observed linear stationary Hawkes processes. When the exact locations of points are not observed, but only counts over time intervals of fixed size, methods based on the likelihood are not feasible. We show that spectral estimation based on Whittle’s method is adapted to this case and provides consistent and asymptotically normal estimators, provided a mild moment condition on the reproduction function. Simulated data sets and a case-study illustrate the performances of the estimation, notably of the reproduction function even when time intervals are relatively large.

We consider multivariate centered Gaussian models for the random variable $\mathit{Z}=({\mathit{Z}}_1,\dots ,{\mathit{Z}}_{\mathit{p}})$, invariant under the action of a subgroup of the group of permutations on $\{1,\dots ,\mathit{p}\}$. Using the representation theory of the symmetric group on the field of reals, we derive the distribution of the maximum likelihood estimate of the covariance parameter Σ and also the analytic expression of the normalizing constant of the Diaconis–Ylvisaker conjugate prior for the precision parameter $\mathit{K}={\mathrm{\Sigma }}^{-1}$. We can thus perform Bayesian model selection in the class of complete Gaussian models invariant by the action of a subgroup of the symmetric group, which we could also call complete RCOP models. We illustrate our results with a toy example of dimension 4 and several examples for selection within cyclic groups, including a high- dimensional example with $\mathit{p}=100$.

We consider the problem of ranking n players from partial pairwise comparison data under the Bradley–Terry–Luce model. For the first time in the literature, the minimax rate of this ranking problem is derived with respect to the Kendall’s tau distance that measures the difference between two rank vectors by counting the number of inversions. The minimax rate of ranking exhibits a transition between an exponential rate and a polynomial rate depending on the magnitude of the signal-to-noise ratio of the problem. To the best of our knowledge, this phenomenon is unique to full ranking and has not been seen in any other statistical estimation problem. To achieve the minimax rate, we propose a divide-and-conquer ranking algorithm that first divides the n players into groups of similar skills and then computes local MLE within each group. The optimality of the proposed algorithm is established by a careful approximate independence argument between the two steps.

We propose confidence regions with asymptotically correct uniform coverage probability of parameters whose Fisher information matrix can be singular at important points of the parameter set. Our work is motivated by the need for reliable inference on scale parameters close or equal to zero in mixed models, which is obtained as a special case. The confidence regions are constructed by inverting a continuous extension of the score test statistic standardized by expected information, which we show exists at points of singular information under regularity conditions. Similar results have previously only been obtained for scalar parameters, under conditions stronger than ours, and applications to mixed models have not been considered. In simulations our confidence regions have near-nominal coverage with as few as $\mathit{n}=20$ observations, regardless of how close to the boundary the true parameter is. It is a corollary of our main results that the proposed test statistic has an asymptotic chi-square distribution with degrees of freedom equal to the number of tested parameters, even if they are on the boundary of the parameter set.

One fundamental goal of high-dimensional statistics is to detect or recover planted structure (such as a low-rank matrix) hidden in noisy data. A growing body of work studies low-degree polynomials as a restricted model of computation for such problems: it has been demonstrated in various settings that low-degree polynomials of the data can match the statistical performance of the best known polynomial-time algorithms. Prior work has studied the power of low-degree polynomials for the task of detecting the presence of hidden structures. In this work, we extend these methods to address problems of estimation and recovery (instead of detection). For a large class of “signal plus noise” problems, we give a user-friendly lower bound for the best possible mean squared error achievable by any degree-D polynomial. To our knowledge, these are the first results to establish low-degree hardness of recovery problems for which the associated detection problem is easy. As applications, we give a tight characterization of the low-degree minimum mean squared error for the planted submatrix and planted dense subgraph problems, resolving (in the low-degree framework) open problems about the computational complexity of recovery in both cases.