This is an expository paper that reviews recent developments on optimal estimation of structured high-dimensional covariance and precision matrices. Minimax rates of convergence for estimating several classes of structured covariance and precision matrices, including bandable, Toeplitz, sparse, and sparse spiked covariance matrices as well as sparse precision matrices, are given under the spectral norm loss. Data-driven adaptive procedures for estimating various classes of matrices are presented. Some key technical tools including large deviation results and minimax lower bound arguments that are used in the theoretical analyses are discussed. In addition, estimation under other losses and a few related problems such as Gaussian graphical models, sparse principal component analysis, factor models, and hypothesis testing on the covariance structure are considered. Some open problems on estimating high-dimensional covariance and precision matrices and their functionals are also discussed.

Professors Cai, Ren and Zhou ought to be congratulated for writing such a wonderful expository paper on optimal estimation of high-dimensional covariance and precision matrices. Nearly all optimality results on large matrix estimation were established by the authors (and their co-authors). Thus, they are the most appropriate team to write this much needed review article. My discussion contains three sections.

In this discussion, we present a brief overview of recent works on the behavior of summary statistics for high-dimensional observations that are time-dependent, and the inference on parameters associated with high-dimensional time series, with emphasis on covariance and auto-covariance matrices.

Transformation models, like the Box-Cox transformation, are widely used in regression to reduce non-additivity, non-normality, and heteroscedasticity. The question of whether one may or may not treat the estimated transformation parameter as fixed in inference about other model parameters has a long and controversial history (Bickel and Doksum, 1981, Hinkley and Runger, 1984). While the frequentist wisdom is that uncertainty regarding the true value of the transformation parameter cannot be ignored, in practice, difficulties in interpretation arise if the transformation is regarded as random and not fixed. In this paper, we suggest a golden mean methodology which attempts to reconcile these philosophies. Penalized estimation yields oracle estimates of transformations that enable treating the transformation parameter as known when the data indicate one of a prespecified set of transformations of scientific interest. When the true transformation is outside this set, rigorous frequentist inference is still achieved. The methodology permits multiple candidate values for the transformation, as is common in applications, as well as simultaneously accommodating variable selection in regression model. Theoretical issues, such as selection consistency and the oracle property, are rigorously established. Numerical studies, including extensive simulation studies and real data examples, illustrate the practical utility of the proposed methods.

We introduce a general method for variable selection in a large-scale regression setting where both the number of parameters and the number of samples are extremely large. The proposed method is based on careful combination of penalized estimators, each applied to a random projection of the sample space into a low-dimensional space. In one special case that we study in detail, the random projections are divided into non-overlapping blocks, each consisting of only a small portion of the original data. Within each block we select the projection yielding the smallest out-of-sample error. Our random ensemble estimator then aggregates the results according to a new maximal-contrast voting scheme to determine the final selected set. Our theoretical results illustrate the effect on performance of increasing the number of non-overlapping blocks. Moreover, we demonstrate that statistical optimality is retained along with the computational speedup. The proposed method achieves minimax rates for approximate recovery over all estimators, using the full set of samples. Furthermore, our theoretical results allow the number of subsamples to grow with the subsample size and do not require irrepresentable condition. The estimator is also compared empirically with several other popular high-dimensional estimators via an extensive simulation study, which reveals its excellent finite-sample performance.

We prove that the convex least squares estimator (LSE) attains a $n^{-1/2}$ pointwise rate of convergence in any region where the truth is linear. In addition, the asymptotic distribution can be characterized by a modified invelope process. Analogous results hold when one uses the derivative of the convex LSE to perform derivative estimation. These asymptotic results facilitate a new consistent testing procedure on the linearity against a convex alternative. Moreover, we show that the convex LSE adapts to the optimal rate at the boundary points of the region where the truth is linear, up to a log-log factor. These conclusions are valid in the context of both density estimation and regression function estimation.

Mode clustering is a nonparametric method for clustering that defines clusters using the basins of attraction of a density estimator’s modes. We provide several enhancements to mode clustering: (i) a soft variant of cluster assignment, (ii) a measure of connectivity between clusters, (iii) a technique for choosing the bandwidth, (iv) a method for denoising small clusters, and (v) an approach to visualizing the clusters. Combining all these enhancements gives us a complete procedure for clustering in multivariate problems. We also compare mode clustering to other clustering methods in several examples.

Motivated by applications in personalized web services and clinical research, we consider a multi-armed bandit problem in a setting where the mean reward of each arm is associated with some covariates. A multi-stage randomized allocation with arm elimination algorithm is proposed to combine the flexibility in reward function modeling and a theoretical guarantee of a cumulative regret minimax rate. When the function smoothness parameter is unknown, the algorithm is equipped with a histogram estimation based smoothness parameter selector using Lepski’s method, and is shown to maintain the regret minimax rate up to a logarithmic factor under a “self-similarity” condition.

This paper studies identifiability and convergence behaviors for parameters of multiple types, including matrix-variate ones, that arise in finite mixtures, and the effects of model fitting with extra mixing components. We consider several notions of strong identifiability in a matrix-variate setting, and use them to establish sharp inequalities relating the distance of mixture densities to the Wasserstein distances of the corresponding mixing measures. Characterization of identifiability is given for a broad range of mixture models commonly employed in practice, including location-covariance mixtures and location-covariance-shape mixtures, for mixtures of symmetric densities, as well as some asymmetric ones. Minimax lower bounds and rates of convergence for the maximum likelihood estimates are established for such classes, which are also confirmed by simulation studies.

The data augmentation (DA) algorithm is considered a useful Markov chain Monte Carlo algorithm that sometimes suffers from slow convergence. It is often possible to convert a DA algorithm into a sandwich algorithm that is computationally equivalent to the DA algorithm, but converges much faster. Theoretically, the reversible Markov chain that drives the sandwich algorithm is at least as good as the corresponding DA chain in terms of performance in the central limit theorem and in the operator norm sense. In this paper, we use the sandwich machinery to compare two DA algorithms. In particular, we provide conditions under which one DA chain can be represented as a sandwich version of the other. Our results are used to extend Hobert and Marchev’s (2008) results on the Haar PX-DA algorithm and to improve the collapsing theorem of Liu et al. (1994) and Liu (1994). We also illustrate our results using Brownlee’s (1965) stack loss data.

Various association measures have been proposed in the literature that equal zero when the associated random variables are independent. However many measures, (e.g., Kendall’s tau), may equal zero even in the presence of an association between the random variables. In order to overcome this drawback, Bergsma and Dassios (2014) proposed a modification of Kendall’s tau, (denoted as $\tau^{*}$), which is non-negative and zero if and only if independence holds. In this article, we investigate the robustness properties and the asymptotic distributions of $\tau^{*}$ and some other well-known measures of association under null and contiguous alternatives. Based on these asymptotic distributions under contiguous alternatives, we study the asymptotic power of the test based on $\tau^{*}$ under contiguous alternatives and compare its performance with the performance of other well-known tests available in the literature.

This paper develops a systematic theory for high-dimensional linear models with dependent errors and/or dependent covariates. To study properties of estimates of the regression parameters, we adopt the framework of functional dependence measures ([43]). For the covariates two schemes are addressed: the random design and the deterministic design. For the former we apply the constrained $\ell_{1}$ minimization approach, while for the latter the Lasso estimation procedure is used. We provide a detailed characterization on how the error rates of the estimates depend on the moment conditions that control the tail behaviors, the dependencies of the underlying processes that generate the errors and the covariates, the dimension and the sample size. Our theory substantially extends earlier ones by allowing dependent and/or heavy-tailed errors and the covariates. As our main tools, we derive exponential tail probability inequalities for dependent sub-Gaussian errors and Nagaev-type inequalities for dependent non-sub-Gaussian errors that arise from linear or non-linear processes.

Scoring functions are used to evaluate and compare partially probabilistic forecasts. We investigate the use of rank-sum functions such as empirical Area Under the Curve (AUC), a widely used measure of classification performance, as a scoring function for the prediction of probabilities of a set of binary outcomes. It is shown that the AUC is not generally a proper scoring function, that is, under certain circumstances it is possible to improve on the expected AUC by modifying the quoted probabilities from their true values. However with some restrictions, or with certain modifications, it can be made proper.

We study parameter identifiability of directed Gaussian graphical models with one latent variable. In the scenario we consider, the latent variable is a confounder that forms a source node of the graph and is a parent to all other nodes, which correspond to the observed variables. We give a graphical condition that is sufficient for the Jacobian matrix of the parametrization map to be full rank, which entails that the parametrization is generically finite-to-one, a fact that is sometimes also referred to as local identifiability. We also derive a graphical condition that is necessary for such identifiability. Finally, we give a condition under which generic parameter identifiability can be determined from identifiability of a model associated with a subgraph. The power of these criteria is assessed via an exhaustive algebraic computational study for small models with 4, 5, and 6 observable variables, and a simulation study for large models with 25 or 35 observable variables.

A key question in modern statistics is how to make fast and reliable inferences for complex, high-dimensional data. While there has been much interest in sparse techniques, current methods do not generalize well to data with nonlinear structure. In this work, we present an orthogonal series estimator for predictors that are complex aggregate objects, such as natural images, galaxy spectra, trajectories, and movies. Our series approach ties together ideas from manifold learning, kernel machine learning, and Fourier methods. We expand the unknown regression on the data in terms of the eigenfunctions of a kernel-based operator, and we take advantage of orthogonality of the basis with respect to the underlying data distribution, $P$, to speed up computations and tuning of parameters. If the kernel is appropriately chosen, then the eigenfunctions adapt to the intrinsic geometry and dimension of the data. We provide theoretical guarantees for a radial kernel with varying bandwidth, and we relate smoothness of the regression function with respect to $P$ to sparsity in the eigenbasis. Finally, using simulated and real-world data, we systematically compare the performance of the spectral series approach with classical kernel smoothing, k-nearest neighbors regression, kernel ridge regression, and state-of-the-art manifold and local regression methods.

We study the convergence of $Z$-estimators $\widehat{\theta}(\eta)\in \mathbb{R}^{p}$ for which the objective function depends on a parameter $\eta$ that belongs to a Banach space $\mathcal{H}$. Our results include the uniform consistency over $\mathcal{H}$ and the weak convergence in the space of bounded $\mathbb{R}^{p}$-valued functions defined on $\mathcal{H}$. When $\eta$ is a tuning parameter optimally selected at $\eta_{0}$, we provide conditions under which $\eta_{0}$ can be replaced by an estimated $\widehat{\eta}$ without affecting the asymptotic variance. Interestingly, these conditions are free from any rate of convergence of $\widehat{\eta}$ to $\eta_{0}$ but require the space described by $\widehat{\eta}$ to be not too large in terms of bracketing metric entropy. In particular, we show that Nadaraya-Watson estimators satisfy this entropy condition. We highlight several applications of our results and we study the case where $\eta$ is the weight function in weighted regression.

Regression models with functional responses and covariates constitute a powerful and increasingly important model class. However, regression with functional data poses well known and challenging problems of non-identifiability. This non-identifiability can manifest itself in arbitrarily large errors for coefficient surface estimates despite accurate predictions of the responses, thus invalidating substantial interpretations of the fitted models. We offer an accessible rephrasing of these identifiability issues in realistic applications of penalized linear function-on-function-regression and delimit the set of circumstances under which they are likely to occur in practice. Specifically, non-identifiability that persists under smoothness assumptions on the coefficient surface can occur if the functional covariate’s empirical covariance has a kernel which overlaps that of the roughness penalty of the spline estimator. Extensive simulation studies validate the theoretical insights, explore the extent of the problem and allow us to evaluate their practical consequences under varying assumptions about the data generating processes. A case study illustrates the practical significance of the problem. Based on theoretical considerations and our empirical evaluation, we provide immediately applicable diagnostics for lack of identifiability and give recommendations for avoiding estimation artifacts in practice.

This paper considers functional models for longitudinal data with subject and group specific trends modelled using Gaussian processes. Fitting Gaussian process regression models is a computationally challenging task, and various sparse approximations to Gaussian processes have been considered in the literature to ease the computational burden. This manuscript builds on a fast non-standard variational approximation which uses a sparse spectral representation and is able to treat uncertainty in the covariance function hyperparameters. This allows fast variational computational methods to be extended to models where there are many functions to be estimated and where there is a hierarchical model involving the covariance function parameters. The main goal of this paper is to implement this idea in the context of functional models for longitudinal data by allowing individual specific smoothness related to covariates for different subjects. Understanding the relationship of smoothness to individual specific covariates is of great interest in some applications. The methods are illustrated with simulated data and a dataset of streamflow curves generated by a rainfall runoff model, and compared with MCMC. It is also shown how these methods can be used to obtain good proposal distributions for MCMC analyses.

We derive the explicit form of expectation propagation for approximate deterministic Bayesian inference in a simple statistical model. The model corresponds to a random sample from the Normal distribution. The explicit forms, and their derivation, allow a deeper understanding of the issues and challenges involved in practical implementation of expectation propagation for statistical analyses. No auxiliary approximations are used: we follow the expectation propagation prescription exactly. A simulation study shows expectation propagation to be more accurate than mean field variational Bayes for larger sample sizes, but at the cost of considerably more algebraic and computational effort.

Stationary determinantal point processes are proved to be Brillinger mixing. This property is an important step towards asymptotic statistics for these processes. As an important example, a central limit theorem for a wide class of functionals of determinantal point processes is established. This result yields in particular the asymptotic normality of the estimator of the intensity of a stationary determinantal point process and of the kernel estimator of its pair correlation.

We study nonparametric Bayesian inference with location mixtures of the Laplace density and a Dirichlet process prior on the mixing distribution. We derive a contraction rate of the corresponding posterior distribution, both for the mixing distribution relative to the Wasserstein metric and for the mixed density relative to the Hellinger and $L_{q}$ metrics.

We study the performance of nonparametric Bayes procedures for one-dimensional diffusions with periodic drift. We improve existing convergence rate results for Gaussian process (GP) priors with fixed hyper parameters. Moreover, we exhibit several possibilities to achieve adaptation to smoothness. We achieve this by considering hierarchical procedures that involve either a prior on a multiplicative scaling parameter, or a prior on the regularity parameter of the GP.

Generalized linear ($GL$-) statistics are defined as functionals of an $U$-quantile process and unify different classes of statistics such as $U$-statistics and $L$-statistics. We derive a central limit theorem for $GL$-statistics of strongly mixing sequences and arbitrary dimension of the underlying kernel. For this purpose we establish a limit theorem for $U$-statistics and an invariance principle for $U$-processes together with a convergence rate for the remaining term of the Bahadur representation.

An application is given by the generalized median estimator for the tail-parameter of the Pareto distribution, which is commonly used to model exceedances of high thresholds. We use subsampling to calculate confidence intervals and investigate its behaviour under independence and under strong mixing in simulations.

In the present paper, we constructed an estimator of a delta contaminated mixing density function $g(\lambda)$ of an intensity $\lambda$ of the Poisson distribution. The estimator is based on an expansion of the continuous portion $g_{0}(\lambda)$ of the unknown pdf over an overcomplete dictionary with the recovery of the coefficients obtained as the solution of an optimization problem with Lasso penalty. In order to apply Lasso technique in the, so called, prediction setting where it requires virtually no assumptions on the dictionary and, moreover, to ensure fast convergence of Lasso estimator, we use a novel formulation of the optimization problem based on the inversion of the dictionary elements.

We formulate conditions on the dictionary and the unknown mixing density that yield a sharp oracle inequality for the norm of the difference between $g_{0}(\lambda)$ and its estimator and, thus, obtain a smaller error than in a minimax setting. Numerical simulations and comparisons with the Laguerre functions based estimator recently constructed by [8] also show advantages of our procedure. At last, we apply the technique developed in the paper to estimation of a delta contaminated mixing density of the Poisson intensity of the Saturn’s rings data.

Nonparametric smoothing methods have been widely used in trend analysis. However, the inference procedure usually requires the crucial assumption that the underlying trend function is smooth. This paper considers the situation where the trend function has potential jumps in addition to smooth changes. In order to determine the existence of jumps, we propose a nonparametric test that can survive under dependent and nonstationary errors, where existing tests assuming independence or stationarity can fail. When the existence of jumps is affirmative, we further consider the problem of estimating the number, location and size of jumps. The results are illustrated via both Monte Carlo simulations and a real data example.

In this paper, we analyze methods for estimating the density of a conditional expectation. We compare an estimator based on a straightforward application of kernel density estimation to a bias-corrected estimator that we propose. We prove convergence results for these estimators and show that the bias-corrected estimator has a superior rate of convergence. In a simulated test case, we show that the bias-corrected estimator performs better in a practical example with a realistic sample size.

Inference for the stochastic blockmodel is currently of burgeoning interest in the statistical community, as well as in various application domains as diverse as social networks, citation networks, brain connectivity networks (connectomics), etc. Recent theoretical developments have shown that spectral embedding of graphs yields tractable distributional results; in particular, a random dot product latent position graph formulation of the stochastic blockmodel informs a mixture of normal distributions for the adjacency spectral embedding. We employ this new theory to provide an empirical Bayes methodology for estimation of block memberships of vertices in a random graph drawn from the stochastic blockmodel, and demonstrate its practical utility. The posterior inference is conducted using a Metropolis-within-Gibbs algorithm. The theory and methods are illustrated through Monte Carlo simulation studies, both within the stochastic blockmodel and beyond, and experimental results on a Wikipedia graph are presented.

Principal support vector machine was proposed recently by Li, Artemiou and Li (2011) to combine L1 support vector machine and sufficient dimension reduction. We introduce the principal L$q$ support vector machine as a unified framework for linear and nonlinear sufficient dimension reduction. By noticing that the solution of L1 support vector machine may not be unique, we set $q>1$ to ensure the uniqueness of the solution. The asymptotic distribution of the proposed estimators are derived for $q>1$. We demonstrate through numerical studies that the proposed L2 support vector machine estimators improve existing methods in accuracy, and are less sensitive to the tuning parameter selection.

Graphical models are widely used to model stochastic dependences among large collections of variables. We introduce a new method of estimating undirected conditional independence graphs based on the score matching loss, introduced by Hyvärinen (2005), and subsequently extended in Hyvärinen (2007). The regularized score matching method we propose applies to settings with continuous observations and allows for computationally efficient treatment of possibly non-Gaussian exponential family models. In the well-explored Gaussian setting, regularized score matching avoids issues of asymmetry that arise when applying the technique of neighborhood selection, and compared to existing methods that directly yield symmetric estimates, the score matching approach has the advantage that the considered loss is quadratic and gives piecewise linear solution paths under $\ell_{1}$ regularization. Under suitable irrepresentability conditions, we show that $\ell_{1}$-regularized score matching is consistent for graph estimation in sparse high-dimensional settings. Through numerical experiments and an application to RNAseq data, we confirm that regularized score matching achieves state-of-the-art performance in the Gaussian case and provides a valuable tool for computationally efficient estimation in non-Gaussian graphical models.

Varying coefficient models are useful generalizations of parametric linear models. They allow for parameters that depend on a covariate or that develop in time. They have a wide range of applications in time series analysis and regression. In time series analysis they have turned out to be a powerful approach to infer on behavioral and structural changes over time. In this paper, we are concerned with high dimensional varying coefficient models including the time varying coefficient model. Most studies in high dimensional nonparametric models treat penalization of series estimators. On the other side, kernel smoothing is a well established, well understood and successful approach in nonparametric estimation, in particular in the time varying coefficient model. But not much has been done for kernel smoothing in high-dimensional models. In this paper we will close this gap and we develop a penalized kernel smoothing approach for sparse high-dimensional models. The proposed estimators make use of a novel penalization scheme working with kernel smoothing. We establish a general and systematic theoretical analysis in high dimensions. This complements recent alternative approaches that are based on basis approximations and that allow more direct arguments to carry over insights from high-dimensional linear models. Furthermore, we develop theory not only for regression with independent observations but also for local stationary time series in high-dimensional sparse varying coefficient models. The development of theory for local stationary processes in a high-dimensional setting creates technical challenges. We also address issues of numerical implementation and of data adaptive selection of tuning parameters for penalization.The finite sample performance of the proposed methods is studied by simulations and it is illustrated by an empirical analysis of NASDAQ composite index data.

We consider the problem of recovering a low-rank signal matrix in high-dimensional situations. The main issue is how to estimate the signal matrix in the presence of huge noise. We introduce the power spiked model to describe the structure of singular values of a huge data matrix. We first consider the conventional PCA to recover the signal matrix and show that the estimation of the signal matrix holds consistency properties under severe conditions. The conventional PCA is heavily subjected to the noise. In order to reduce the noise we apply the noise-reduction (NR) methodology and propose a new estimation of the signal matrix. We show that the proposed estimation by the NR method holds the consistency properties under mild conditions and improves the error rate of the conventional PCA effectively. Finally, we demonstrate the reconstruction procedures by using a microarray data set.

Fast multiple change-point segmentation methods, which additionally provide faithful statistical statements on the number, locations and sizes of the segments, have recently received great attention. In this paper, we propose a multiscale segmentation method, FDRSeg, which controls the false discovery rate (FDR) in the sense that the number of false jumps is bounded linearly by the number of true jumps. In this way, it adapts the detection power to the number of true jumps. We prove a non-asymptotic upper bound for its FDR in a Gaussian setting, which allows to calibrate the only parameter of FDRSeg properly. Moreover, we show that FDRSeg estimates change-point locations, as well as the signal, in a uniform sense at optimal minimax convergence rates up to a log-factor. The latter is w.r.t. $L^{p}$-risk, $p\ge 1$, over classes of step functions with bounded jump sizes and either bounded, or even increasing, number of change-points. FDRSeg can be efficiently computed by an accelerated dynamic program; its computational complexity is shown to be linear in the number of observations when there are many change-points. The performance of the proposed method is examined by comparisons with some state of the art methods on both simulated and real datasets. An R-package is available online.

We present a novel method for controlling the $k$-familywise error rate ($k$-FWER) in the linear regression setting using the knockoffs framework first introduced by Barber and Candès. Our procedure, which we also refer to as knockoffs, can be applied with any design matrix with at least as many observations as variables, and does not require knowing the noise variance. Unlike other multiple testing procedures which act directly on $p$-values, knockoffs is specifically tailored to linear regression and implicitly accounts for the statistical relationships between hypothesis tests of different coefficients. We prove that knockoffs controls the $k$-FWER exactly in finite samples and show in simulations that it provides superior power to alternative procedures over a range of linear regression problems. We also discuss extensions to controlling other Type I error rates such as the false exceedance rate, and use it to identify candidates for mutations conferring drug-resistance in HIV.

The first Bayesian results for the sparse normal means problem were proven for spike-and-slab priors. However, these priors are less convenient from a computational point of view. In the meanwhile, a large number of continuous shrinkage priors has been proposed. Many of these shrinkage priors can be written as a scale mixture of normals, which makes them particularly easy to implement. We propose general conditions on the prior on the local variance in scale mixtures of normals, such that posterior contraction at the minimax rate is assured. The conditions require tails at least as heavy as Laplace, but not too heavy, and a large amount of mass around zero relative to the tails, more so as the sparsity increases. These conditions give some general guidelines for choosing a shrinkage prior for estimation under a nearly black sparsity assumption. We verify these conditions for the class of priors considered in [12], which includes the horseshoe and the normal-exponential gamma priors, and for the horseshoe+, the inverse-Gaussian prior, the normal-gamma prior, and the spike-and-slab Lasso, and thus extend the number of shrinkage priors which are known to lead to posterior contraction at the minimax estimation rate.

Uniform and nonuniform Berry–Esseen (BE) bounds of optimal orders on the rate of convergence to normality in the delta method for vector statistics are obtained. The results are applicable almost as widely as the delta method itself – except that, quite naturally, the order of the moments needed to be finite is generally $3/2$ times as large as that for the corresponding central limit theorems. Our BE bounds appear new even for the one-dimensional delta method, that is, for smooth functions of the sample mean of univariate random variables. Specific applications to Pearson’s, noncentral Student’s and Hotelling’s statistics, sphericity test statistics, a regularized canonical correlation, and maximum likelihood estimators (MLEs) are given; all these uniform and nonuniform BE bounds appear to be the first known results of these kinds, except for uniform BE bounds for MLEs. The new method allows one to obtain bounds with explicit and rather moderate-size constants. For instance, one has the uniform BE bound $3.61\mathbb{E}(Y_{1}^{6}+Z_{1}^{6})\,(1+\sigma^{-3})/\sqrt{n}$ for the Pearson sample correlation coefficient based on independent identically distributed random pairs $(Y_{1},Z_{1}),\dots,(Y_{n},Z_{n})$ with $\mathbb{E} Y_{1}=\mathbb{E}Z_{1}=\mathbb{E}Y_{1}Z_{1}=0$ and $\mathbb{E}Y_{1}^{2}=\mathbb{E}Z_{1}^{2}=1$, where $\sigma:=\sqrt{\mathbb{E}Y_{1}^{2}Z_{1}^{2}}$.

The dependence structure in a $d$-variate continuous random vector $\mathbf{X}$ is characterized by its unique copula. Starting from the fact that many copulas can be extracted from the global $d$-dimensional copula of $\mathbf{X}$, a very general framework is proposed here for testing that a given collection of induced $p$-dimensional copulas from a multivariate distribution are identical. Many hypotheses of interest in copula modeling fall into this category, including bivariate symmetry (diagonal, radial, joint), exchangeability, as well as various types of equality of copulas. Here, a broad class of test statistics is defined around a matrix representation of the null hypothesis and quadratic functionals including Cramér–von Mises and characteristic function mappings. Since the null hypotheses to be tested are composite by nature, the computation of $\mathrm{P}$-values is achieved using multiplier bootstrap versions of the test statistics. The sample properties of the method are investigated when testing for several types of bivariate symmetry, exchangeability, equality of non-overlapping and overlapping copulas and equality of all bivariate copulas. The general conclusion is that the tests are good at keeping their nominal level and are powerful against a wide variety of alternatives, showing the relevance and reliability of the methodology for the modeling of multivariate datasets with the help of copulas.

Scoring rules evaluate the performance of probabilistic forecasts. A scoring rule is said to be local if it assigns a score based on the observed outcome and on outcomes that are in some sense “close” to the observed outcome. All scoring rules can be derived from a concave entropy functional and the property of locality follows when the entropy is 1-homogeneous (up to an additive constant). Consequently, except for the log score, a local scoring rule has the remarkable property that it is 0-homogeneous; in other words, it assigns a score that is independent of the normalization of the quoted probability distribution. In many statistical applications, it is not plausible to treat observed outcomes as independent, e.g. time series data or multicomponent measurements. We show that local scoring rules can be easily extended to multidimensional outcome spaces. We also introduce the notion of an extensive scoring rule, i.e. a scoring rule that ensures the score of independent outcomes is a sum of independent scores. We construct local scoring rules that are extensive and show that a scoring rule is a extensive if and only if it is derived from an extensive entropy.

Respondent-Driven Sampling (RDS) is a widely adopted link-tracing sampling design used to draw valid statistical inference from samples of populations for which there is no available sampling frame. RDS estimators rely upon the assumption that each edge (representing a relationship between two individuals) in the underlying network has an equal probability of being sampled. We show that this assumption is violated in even the simplest cases, and that RDS estimators are sensitive to the violation of this assumption.

We aim to estimate multiple networks in the presence of sample heterogeneity, where the independent samples (i.e. observations) may come from different and unknown populations or distributions. Specifically, we consider penalized estimation of multiple precision matrices in the framework of a Gaussian mixture model. A major innovation is to take advantage of the commonalities across the multiple precision matrices through possibly nonconvex fusion regularization, which for example makes it possible to achieve simultaneous discovery of unknown disease subtypes and detection of differential gene (dys)regulations in functional genomics. We embed in the EM algorithm one of two recently proposed methods for estimating multiple precision matrices in Gaussian graphical models. We demonstrate the feasibility and potential usefulness of the proposed methods in an application to glioblastoma subtype discovery and differential gene network analysis with a microarray gene expression data set. We also conduct realistic simulation studies to evaluate and compare the performance of various methods.

In many longitudinal studies, the covariates and response are often intermittently observed at irregular, mismatched and subject-specific times. Last observation carried forward (LOCF) is one of the most commonly used methods to deal with such data when covariates and response are observed asynchronously. However, this can lead to considerable bias. In this paper, we propose a weighted LOCF estimation using asynchronous longitudinal data for the generalized linear model. We further generalize this approach to utilize previously observed covariates in addition to the most recent observation. In comparison to earlier methods, the current methods are valid under weaker assumptions on the covariate process and allow informative observation times which may depend on response even conditional on covariates. Extensive simulation studies provide numerical support for the theoretical findings. Data from an HIV study is used to illustrate our methodology.

The functional delta-method provides a convenient tool for deriving the asymptotic distribution of a plug-in estimator of a statistical functional from the asymptotic distribution of the respective empirical process. Moreover, it provides a tool to derive bootstrap consistency for plug-in estimators from bootstrap consistency of empirical processes. It has recently been shown that the range of applications of the functional delta-method for the asymptotic distribution can be considerably enlarged by employing the notion of quasi-Hadamard differentiability. Here we show in a general setting that this enlargement carries over to the bootstrap. That is, for quasi-Hadamard differentiable functionals bootstrap consistency of the plug-in estimator follows from bootstrap consistency of the respective empirical process. This enlargement often requires convergence in distribution of the bootstrapped empirical process w.r.t. a nonuniform sup-norm. The latter is not problematic as will be illustrated by means of examples.

We consider a population of $N$ individuals, of which we observe the number of actions until time $t$. For each couple of individuals $(i,j)$, $j$ may or not influence $i$, which we model by i.i.d. Bernoulli$(p)$-random variables, for some unknown parameter $p\in(0,1]$. Each individual acts autonomously at some unknown rate $\mu>0$ and acts by mimetism at some rate proportional to the sum of some function $\varphi$ of the ages of the actions of the individuals which influence him. The function $\varphi$ is unknown but assumed, roughly, to be decreasing and with fast decay. The goal of this paper is to estimate $p$, which is the main characteristic of the graph of interactions, in the asymptotic $N\to\infty$, $t\to \infty$. The main issue is that the mean field limit (as $N\to\infty$) of this model is unidentifiable, in that it only depends on the parameters $\mu$ and $p\varphi$. Fortunately, this mean field limit is not valid for large times. We distinguish the subcritical case, where, roughly, the mean number $m_{t}$ of actions per individual increases linearly and the supercritical case, where $m_{t}$ increases exponentially. Although the nuisance parameter $\varphi$ is non-parametric, we are able, in both cases, to estimate $p$ without estimating $\varphi $ in a nonparametric way, with a precision of order $N^{-1/2}+N^{1/2}m_{t}^{-1}$, up to some arbitrarily small loss. We explain, using a Gaussian toy model, the reason why this rate of convergence might be (almost) optimal.

We propose a minimum distance estimation method for robust regression in sparse high-dimensional settings. Likelihood-based estimators lack resilience against outliers and model misspecification, a critical issue when dealing with high-dimensional noisy data. Our method, Minimum Distance Lasso (MD-Lasso), combines minimum distance functionals customarily used in nonparametric estimation for robustness, with $\ell_{1}$-regularization. MD-Lasso is governed by a scaling parameter capping the influence of outliers: the loss is locally convex and close to quadratic for small squared residuals, and flattens for squared residuals larger than the scaling parameter. As the parameter approaches infinity the estimator becomes equivalent to least-squares Lasso. MD-Lasso is able to maintain the robustness of minimum distance functionals in sparse high-dimensional regression. The estimator achieves maximum breakdown point and enjoys consistency with fast convergence rates under mild conditions on the model error distribution. These hold for any solution in a convexity region around the true parameter and in certain cases for every solution. We provide an alternative set of results that do not require the solutions to lie within the convexity region but where the $\ell_{2}$-norm of the feasible solutions is constrained within a safety radius. Thanks to this constraint, a first-order optimization method is able to produce local optima that are consistent. A connection is established with re-weighted least-squares that intuitively explains MD-Lasso robustness. The merits of our method are demonstrated through simulation and eQTL analysis.

We introduce a general framework for estimation of inverse covariance, or precision, matrices from heterogeneous populations. The proposed framework uses a Laplacian shrinkage penalty to encourage similarity among estimates from disparate, but related, subpopulations, while allowing for differences among matrices. We propose an efficient alternating direction method of multipliers (ADMM) algorithm for parameter estimation, as well as its extension for faster computation in high dimensions by thresholding the empirical covariance matrix to identify the joint block diagonal structure in the estimated precision matrices. We establish both variable selection and norm consistency of the proposed estimator for distributions with exponential or polynomial tails. Further, to extend the applicability of the method to the settings with unknown populations structure, we propose a Laplacian penalty based on hierarchical clustering, and discuss conditions under which this data-driven choice results in consistent estimation of precision matrices in heterogenous populations. Extensive numerical studies and applications to gene expression data from subtypes of cancer with distinct clinical outcomes indicate the potential advantages of the proposed method over existing approaches.

This article proposes a first analysis of kernel spectral clustering methods in the regime where the dimension $p$ of the data vectors to be clustered and their number $n$ grow large at the same rate. We demonstrate, under a $k$-class Gaussian mixture model, that the normalized Laplacian matrix associated with the kernel matrix asymptotically behaves similar to a so-called spiked random matrix. Some of the isolated eigenvalue-eigenvector pairs in this model are shown to carry the clustering information upon a separability condition classical in spiked matrix models. We evaluate precisely the position of these eigenvalues and the content of the eigenvectors, which unveil important (sometimes quite disruptive) aspects of kernel spectral clustering both from a theoretical and practical standpoints. Our results are then compared to the actual clustering performance of images from the MNIST database, thereby revealing an important match between theory and practice.

We propose a comprehensive framework for additive regression models for non-Gaussian functional responses, allowing for multiple (partially) nested or crossed functional random effects with flexible correlation structures for, e.g., spatial, temporal, or longitudinal functional data as well as linear and nonlinear effects of functional and scalar covariates that may vary smoothly over the index of the functional response. Our implementation handles functional responses from any exponential family distribution as well as many others like Beta- or scaled and shifted $t$-distributions. Development is motivated by and evaluated on an application to large-scale longitudinal feeding records of pigs. Results in extensive simulation studies as well as replications of two previously published simulation studies for generalized functional mixed models demonstrate the good performance of our proposal. The approach is implemented in well-documented open source software in the pffr function in R-package refund.

Matrix completion has been well studied under the uniform sampling model and the trace-norm regularized methods perform well both theoretically and numerically in such a setting. However, the uniform sampling model is unrealistic for a range of applications and the standard trace-norm relaxation can behave very poorly when the underlying sampling scheme is non-uniform.

In this paper we propose and analyze a max-norm constrained empirical risk minimization method for noisy matrix completion under a general sampling model. The optimal rate of convergence is established under the Frobenius norm loss in the context of approximately low-rank matrix reconstruction. It is shown that the max-norm constrained method is minimax rate-optimal and yields a unified and robust approximate recovery guarantee, with respect to the sampling distributions. The computational effectiveness of this method is also discussed, based on first-order algorithms for solving convex optimizations involving max-norm regularization.

A new method of hypothesis testing is proposed ensuring that as the sample size grows, the probability of a type I error will become arbitrarily small by allowing the significance level to decrease with the number of observations in the study. Furthermore, the corresponding sequence of hypothesis tests will make only a finite number of errors with probability one, under mild regulatory conditions, including both i.i.d. and strongly mixing processes. It can be used as an alternative to arbitrary fixed significance levels such as $0.05$ or $0.01$.

This approach resolves the Jeffreys-Lindley paradox. It is also robust to multiple comparisons, and to a lesser extent, optional stopping making it somewhat robust to data fishing or p-hacking.

As an example of the practical applications of this technique, it is used as a lag order selection mechanism in simulations. It performs well relative to other model selection criteria. In another simulation, hypothesis tests and confidence intervals for the mean are investigated, demonstrating the improved performance even in small samples. We also show that under mild regularity conditions, any sequence of two-sided hypothesis tests with fixed significance level will make an infinite number of mistakes with positive probability.

In this paper, we present several estimators of the diagonal elements of the inverse of the covariance matrix, called precision matrix, of a sample of independent and identically distributed random vectors. The main focus is on the case of high dimensional vectors having a sparse precision matrix. It is now well understood that when the underlying distribution is Gaussian, the columns of the precision matrix can be estimated independently form one another by solving linear regression problems under sparsity constraints. This approach leads to a computationally efficient strategy for estimating the precision matrix that starts by estimating the regression vectors, then estimates the diagonal entries of the precision matrix and, in a final step, combines these estimators for getting estimators of the off-diagonal entries. While the step of estimating the regression vector has been intensively studied over the past decade, the problem of deriving statistically accurate estimators of the diagonal entries has received much less attention. The goal of the present paper is to fill this gap by presenting four estimators—that seem the most natural ones—of the diagonal entries of the precision matrix and then performing a comprehensive empirical evaluation of these estimators. The estimators under consideration are the residual variance, the relaxed maximum likelihood, the symmetry-enforced maximum likelihood and the penalized maximum likelihood. We show, both theoretically and empirically, that when the aforementioned regression vectors are estimated without error, the symmetry-enforced maximum likelihood estimator has the smallest estimation error. However, in a more realistic setting when the regression vector is estimated by a sparsity-favoring computationally efficient method, the qualities of the estimators become relatively comparable with a slight advantage for the residual variance estimator.

In this paper, we extend the correspondence between Bayes’ estimation and optimal interpolation in a Reproducing Kernel Hilbert Space (RKHS) to the case of convex constraints such as boundedness, monotonicity or convexity. In the unconstrained interpolation case, the mean of the posterior distribution of a Gaussian Process (GP) given data interpolation is known to be the optimal interpolation function minimizing the norm in the RKHS associated to the GP. In the constrained case, we prove that the Maximum A Posteriori (MAP) or Mode of the posterior distribution is the optimal constrained interpolation function in the RKHS. So, the general correspondence is achieved with the MAP estimator and not the mean of the posterior distribution. A numerical example is given to illustrate this last result.

Probabilistic binary classification typically calls for a vector of marginal probabilities where each element gives the probability of assigning the corresponding case to class 1. Scoring rules are principled ways to assess probabilistic forecasts about any outcome that is subsequently observed. We develop a class of proper scoring rules called linear scoring rules that are specifically adapted to probabilistic binary classification. When applied in competition situations, we show that all linear scoring rules essentially balance the needs of organizers and competitors. Linear scoring rules can also be used to train classifiers. Finally, since scoring rules have a statistical decision theoretic foundation, a linear scoring rule can be constructed for any user-defined misclassification loss function.

A new risk bound is presented for the problem of convex/concave function estimation, using the least squares estimator. The best known risk bound, as had appeared in Guntuboyina and Sen, scaled like $\log(en)\:n^{-4/5}$ under the mean squared error loss, up to a constant factor. The authors in [8] had conjectured that the logarithmic term may be an artifact of their proof. We show that indeed the logarithmic term is unnecessary and prove a risk bound which scales like $n^{-4/5}$ up to constant factors. Our proof technique has one extra peeling step than in a usual chaining type argument. Our risk bound holds in expectation as well as with high probability and also extends to the case of model misspecification, where the true function may not be concave.

When describing adequately complex data structures one is often confronted with the fact that mean as well as variance (or more generally dispersion) is highly influenced by some covariates. Drawbacks of the available methods is that they are often based on approximations and hence a theoretical study should deal with also studying these approximations. This however is often ignored, making the statistical inference incomplete. In the proposed framework of double generalized modelling based on proper dispersion models we avoid this drawback and as such are in a good position to use recent results on Bregman divergence for establishing theoretical results for the proposed estimators in fairly general settings. We also study variable selection when there is a large number of covariates, with this number possibly tending to infinity with the sample size. The proposed estimation and selection procedure is investigated via a simulation study, that includes also a comparative study with competitors. The use of the methods is illustrated via some real data applications.

Bayesian analysis is increasingly popular for use in social science and other application areas where the data are observations from an informative sample. An informative sampling design leads to inclusion probabilities that are correlated with the response variable of interest. Model inference performed on the observed sample taken from the population will be biased for the population generative model under informative sampling since the balance of information in the sample data is different from that for the population. Typical approaches to account for an informative sampling design under Bayesian estimation are often difficult to implement because they require re-parameterization of the hypothesized generating model, or focus on design, rather than model-based, inference. We propose to construct a pseudo-posterior distribution that utilizes sampling weights based on the marginal inclusion probabilities to exponentiate the likelihood contribution of each sampled unit, which weights the information in the sample back to the population. Our approach provides a nearly automated estimation procedure applicable to any model specified by the data analyst for the population and retains the population model parameterization and posterior sampling geometry. We construct conditions on known marginal and pairwise inclusion probabilities that define a class of sampling designs where $L_{1}$ consistency of the pseudo posterior is guaranteed. We demonstrate our method on an application concerning the Bureau of Labor Statistics Job Openings and Labor Turnover Survey.