The Nonparanormal: Semiparametric Estimation of High Dimensional Undirected Graphs

Introduction

The linear model is a mainstay of statistical inference that has been extended in several important ways. An extension to high dimensions was achieved by adding a sparsity constraint, leading to the lasso (Tibshirani, 1996). An extension to nonparametric models was achieved by replacing linear functions with smooth functions, leading to additive models (Hastie and Tibshirani, 1999). These two ideas were recently combined, leading to an extension called sparse additive models (SpAM) (Ravikumar et al., 2008b, a). In this paper we consider a similar nonparametric extension of undirected graphical models based on multivariate Gaussian distributions in the high dimensional setting. Specifically, we use a high dimensional Gaussian copula with nonparametric marginals, which we refer to as a nonparanormal distribution.

While Gaussian graphical models can be useful, a reliance on exact normality is limiting. Our goal in this paper is to weaken this assumption. Our approach parallels the ideas behind sparse additive models for regression (Ravikumar et al., 2008b, a). Specifically, we replace the Gaussian with a semiparametric Gaussian copula. This means that we replace the random variable $X=(X_{1},\ldots,X_{p})$ by the transformed random variable $f(X)=\left(f_{1}(X_{1}),\ldots,f_{p}(X_{p})\right)$, and assume that $f(X)$ is multivariate Gaussian. This semiparametric copula results in a nonparametric extension of the normal that we call the nonparanormal distribution. The nonparanormal depends on the functions $\{f_{j}\}$, and a mean $\mu$ and covariance matrix $\Sigma$, all of which are to be estimated from data. While the resulting family of distributions is much richer than the standard parametric normal (the paranormal), the independence relations among the variables are still encoded in the precision matrix $\Omega=\Sigma^{-1}$. We propose a nonparametric estimator for the functions $\{f_{j}\}$, and show how the graphical lasso can be used to estimate the graph in the high dimensional setting. The relationship between linear regression models, Gaussian graphical models, and their extensions to nonparametric and high dimensional models is summarized in Figure 1.

Most theoretical results on semiparametric copulas focus on low or at least finite dimensional models (Tsukahara, 2005). Models with increasing dimension require a more delicate analysis; in particular, simply plugging in the usual empirical distribution of the marginals does not lead to accurate inference. Instead we use a truncated empirical distribution. We give a theoretical analysis of this estimator, proving consistency results with respect to risk, model selection, and estimation of $\Omega$ in the Frobenius norm.

In the following section we review the basic notion of the graph corresponding to a multivariate Gaussian, and formulate different criteria for evaluating estimators of the covariance or inverse covariance. In Section 3 we present the nonparanormal, and in Section 4 we discuss estimation of the model. We present a theoretical analysis of the estimation method in Section 5, with the detailed proofs collected in an appendix. In Section 6 we present experiments with both simulated data and gene microarray data, where the problem is to construct the isoprenoid biosynthetic pathway.

Estimating Undirected Graphs

Let $X=(X_{1},\ldots,X_{p})$ denote a random vector with distribution $P=N(\mu,\Sigma)$. The undirected graph $G=(V,E)$ corresponding to $P$ consists of a vertex set $V$ and an edge set $E$. The set $V$ has $p$ elements, one for each component of $X$. The edge set $E$ consists of ordered pairs $(i,j)$ where $(i,j)\in E$ if there is a edge between $X_{i}$ and $X_{j}$. The edge between $(i,j)$ is excluded from $E$ if and only if $X_{i}$ is independent of $X_{j}$ given the other variables $O_{\backslash\{i,j\}}\equiv(X_{s}:\ 1\leq s\leq p,\ \ s\neq i,j)$, written

It is well-known that, for multivariate Gaussian distributions, (1) holds if and only if $\Omega_{ij}=0$ where $\Omega=\Sigma^{-1}$.

is the sample covariance, with $\overline{X}$ the sample mean. The zeroes of $\Omega$ can then be estimated by applying hypothesis testing to $\widehat{\Omega}$ (Drton and Perlman, 2007, 2008).

When $p>n$, maximum likelihood is no longer useful; in particular, the estimate $\widehat{\Sigma}$ is not positive definite, having rank no greater than $n$. Inspired by the success of the lasso for linear models, several authors have suggested estimating $\Sigma$ by minimizing

is the average log-likelihood and $S$ is the sample covariance matrix. The estimator $\widehat{\Omega}$ can be computed efficiently using the glasso algorithm (Friedman et al., 2007), which is a block coordinate descent algorithm that uses the standard lasso to estimate a single row and column of $\Omega$ in each iteration. Under appropriate sparsity conditions, the resulting estimator $\widehat{\Omega}$ has been shown to have good theoretical properties (Rothman et al., 2008; Ravikumar et al., 2009).

There are several different ways to judge the quality of an estimator $\widehat{\Sigma}$ of the covariance or inverse covariance $\widehat{\Omega}$. We discuss three in this paper, persistency, norm consistency, and sparsistency. Persistency means consistency in risk, when the model is not assumed to be correct. Suppose the true distribution is $P$ has mean $\mu_{0}$, and that we use a multivariate normal $p(x;\mu_{0},\Sigma)$ for prediction. We do not assume that $P$ is normal. We observe a new vector $X\sim P$ and define the prediction risk to be

where $\Sigma_{0}$ is the covariance of $X$ under $P$. If ${\cal S}$ is a set of covariance matrices, the oracle is defined to be the covariance matrix $\Sigma_{*}$ that minimizes $R(\Sigma)$ over ${\cal S}$:

Thus $p(x;\mu_{0},\Sigma_{*})$ is the best predictor of a new observation among all distributions in $\{p(x;\mu_{0},\Sigma):\ \Sigma\in{\cal S}\}$. In particular, if ${\cal S}$ consists of covariance matrices with sparse graphs, then $p(x;\mu_{0},\Sigma_{*})$ is, in some sense, the best sparse predictor. An estimator $\widehat{\Sigma}_{n}$ is persistent if

as the sample size $n$ increases to infinity. Thus, a persistent estimator approximates the best estimator over the class ${\cal S}$, but we do not assume that the true distribution has a covariance matrix in ${\cal S}$, or even that it is Gaussian. Moreover, we allow the dimension $p=p_{n}$ to increase with $n$. On the other hand, norm consistency and sparsistency require that the true distribution is Gaussian. In this case, let $\Sigma_{0}$ denote the true covariance matrix. An estimator is norm consistent if

where $\|\cdot\|$ is a norm. If $E(\Omega)$ denotes the edge set corresponding to $\Omega$. An estimator is sparsistent if

Thus, a sparsistent estimator identifies the correct graph consistently. We summarize known results on these properties for the multivariate normal in Section 5, before presenting our theoretical analysis of the nonparanormal.

The Nonparanormal

We say that a random vector $X=(X_{1},\ldots,X_{p})^{T}$ has a nonparanormal distribution if there exist functions $\{f_{j}\}_{j=1}^{p}$ such that $Z\equiv f(X)\sim N(\mu,\Sigma)$, where $f(X)=(f_{1}(X_{1}),\ldots,f_{p}(X_{p}))$. We then write

When the $f_{j}$’s are monotone and differentiable, the joint probability density function of $X$ is given by

. The nonparanormal distribution $\mathop{\textit{NPN\,}}(\mu,\Sigma,f)$ is a Gaussian copula when the $f_{j}$’s are monotone and differentiable.

Proof. By Sklar’s theorem (Sklar, 1959), any joint distribution can be written as

where the function $C$ is called a copula. For the nonparanormal we have

where $\Phi_{\mu,\Sigma}$ is the multivariate Gaussian cdf and $\Phi$ is the univariate standard Gaussian cdf. Thus, the corresponding copula is

This is exactly a Gaussian copula with parameters $\mu$ and $\Sigma$. If each $f_{j}$ is differentiable then the density of $X$ has the same form as (4). $\;\;\scriptstyle\Box$

Note that the density in (4) is not identifiable; to make the family identifiable we demand that $f_{j}$ preserve means and variances:

Let $F_{j}(x)$ denote the marginal distribution function of $X_{j}$. Then

The following basic fact says that the independence graph of the nonparanormal is encoded in $\Omega=\Sigma^{-1}$, as for the parametric normal.

. If $X\sim\mathop{\textit{NPN\,}}(\mu,\Sigma,f)$ is nonparanormal and each $f_{j}$ is differentiable, then $X_{i}\amalg X_{j}{\,|\,}O_{\backslash\{i,j\}}$ if and only if $\Omega_{ij}=0$, where $\Omega=\Sigma^{-1}$.

Proof. From the form of the density (4), it follows that the density factors with respect to the graph of $\Omega$, and therefore obeys the global Markov property of the graph. $\;\;\scriptstyle\Box$

Next we show that the above is true for any choice of identification restrictions.

and let $\Lambda$ be the covariance matrix of $h(X)$. Then $X_{j}\amalg X_{k}{\,|\,}O_{\backslash\{j,k\}}$ if and only if $\Lambda_{jk}^{-1}=0$.

Proof. We can rewrite the covariance matrix as

where $D$ is the diagonal matrix with $\text{diag}(D)=\sigma$. The zero pattern of $\Lambda^{-1}$ is therefore identical to the zero pattern of $\Sigma^{-1}$. $\;\;\scriptstyle\Box$

Thus, it is not necessary to estimate $\mu$ or $\sigma$ to estimate the graph.

Figure 2 shows three examples of 2-dimensional nonparanormal densities. In each case, the component functions $f_{j}(x)$ take the form

where the constants $a_{j}$ and $b_{j}$ are set to enforce the identifiability constraints (5). The covariance in each case is $\Sigma=\binom{1\ .5}{.5\ 1}$ and the mean is $\mu=(0,0)$. The exponent $\alpha_{j}$ determines the nonlinearity. It can be seen how the concavity of the density changes with the exponent $\alpha$, and that $\alpha>1$ can result in multiple modes.

Estimation Method

where $\widetilde{F}_{j}$ is an estimator of $F_{j}$. A natural candidate for $\widetilde{F}_{j}$ is the marginal empirical distribution function

Now, let $\theta$ denote the parameters of the copula. Tsukahara (2005) suggests taking $\widehat{\theta}$ to be the solution of

where $\phi$ is an estimating equation and $\widetilde{F}_{j}(t)=n\widehat{F}_{j}(t)/(n+1)$. In our case, $\theta$ corresponds to the covariance matrix. The resulting estimator $\widehat{\theta}$, called a rank approximate $Z$-estimator, has excellent theoretical properties. However, we are interested in the high dimensional scenario where the dimension $p$ is allowed to increase with $n$; the variance of $\widehat{F}_{j}(t)$ is too large in this case. Instead, we use the following truncated or WinsorizedAfter Charles P. Winsor, whom John Tukey credited with converting him from topology to statistics (Mallows, 1990). estimator:

where $\delta_{n}$ is a truncation parameter. Clearly, there is a bias-variance tradeoff in choosing $\delta_{n}$. In what follows we use

This provides the right balance so that we can achieve the desired rate of convergence in our estimate of $\Omega$ and the associated undirected graph $G$.

Given this estimate of the distribution of variable $X_{j}$, we then estimate the transformation function $f_{j}$ by

and $\widehat{\mu}_{j}$ and $\widehat{\sigma}_{j}$ are the sample mean and the standard deviation:

Now, let $S_{n}(\widetilde{f})$ be the sample covariance matrix of $\widetilde{f}(X^{(1)}),\ldots,\widetilde{f}(X^{(n)})$; that is,

Theoretical Results

In this section we present our theoretical results on risk consistency, model selection consistency, and norm consistency of the covariance $\Sigma$ and inverse covariance $\Omega$. From Lemma 3.3, the estimate of the graph does not depend on $\sigma_{j},\;j\in\{1,\ldots,p\}$ and $\mu$, so we assume that $\sigma_{j}=1$ and $\mu=0$. Our key technical result is an analysis of covariance of the Winsorized estimator defined in (7), (9), and (10). In particular, we show that under appropriate conditions,

where $S_{n}(\widetilde{f})_{jk}$ denotes the $(j,k)$ entry of the matrix. This result allows us to leverage the recent analysis of Rothman et al. (2008) and Ravikumar et al. (2009) in the Gaussian case to obtain consistency results for the nonparanormal. More precisely, our main theorem is the following.

. Suppose that $p=n^{\xi}$ and let $\widetilde{f}$ be the Winsorized estimator defined in (9) with $\delta_{n}=\displaystyle\frac{1}{4n^{1/4}\sqrt{\pi\log n}}$. Define

for $M,\xi>0$. Then for any $\displaystyle\epsilon\geq C(M,\xi)\sqrt{\frac{\log p\log^{2}n}{n^{1/2}}}$ and sufficiently large $n$, we have

where $c_{1},c_{2},c_{3},c_{4}$ are positive constants.

The proof of the above theorem is given in Section 7. The following corollary is immediate, which specifies the scaling of the dimension in terms of sample size.

. Let $M>\displaystyle 1+\frac{1}{\xi}$. Then

. Suppose that the data are generated as $X^{(i)}\sim\mathop{\textit{NPN\,}}(\mu_{0},\Sigma_{0},f_{0})$, and let $\Omega_{0}=\Sigma_{0}^{-1}$. If the regularization parameter $\lambda_{n}$ is chosen as

then the nonparanormal estimator $\widehat{\Omega}_{n}$ of (11) satisfies

is the number of nonzero off-diagonal elements of the true precision matrix.

To prove the model selection consistency result, we need further assumptions. We follow Ravikumar (2009) and let the $p^{2}\times p^{2}$ Fisher information matrix of $\Sigma_{0}$ be $\Gamma\equiv\Sigma_{0}\otimes\Sigma_{0}$ where $\otimes$ is the Kronecker matrix product, and define the support set $S$ of $\Omega_{0}=\Sigma_{0}^{-1}$ as

We use $S^{c}$ to denote the complement of $S$ in the set $\{1,\ldots,p\}\times\{1,\ldots,p\}$, and for any two subsets $T$ and $T^{\prime}$ of $\{1,\ldots,p\}\times\{1,\ldots,p\}$, we use $\Gamma_{TT^{\prime}}$ to denote the sub-matrix with rows and columns of $\Gamma$ indexed by $T$ and $T^{\prime}$ respectively.

. There exists some $\alpha\in(0,1]$, such that $\left\|\Gamma_{S^{c}S}(\Gamma_{SS})^{-1}\right\|_{\infty}\leq 1-\alpha.$

As in Ravikumar et al. (2009), we define two quantities $K_{\Sigma_{0}}\equiv\|\Sigma_{0}\|_{\infty}$ and $K_{\Gamma}\equiv\|(\Gamma_{SS})^{-1}\|_{\infty}$. Further, we define the maximum row degree as

. The quantities $K_{\Sigma^{0}}$ and $K_{\Gamma}$ are bounded, and there are positive constants $C_{1}$ and $C_{2}$ such that

The proof of the following uses our Theorem 5.1 in place of equation (12) in the analysis of Ravikumar et al. (2009),

. Suppose the regularization parameter is chosen as

Then the nonparanormal estimator $\widehat{\Omega}_{n}$ satisfies

where ${\mathcal{G}}(\widehat{\Omega}_{n},\Omega_{0})$ is the event

Our persistency (risk consistency) result parallels the persistency result for additive models given in Ravikumar et al. (2008a), and allows model dimension that grows exponentially with sample size. The definition in this theorem uses the fact (from Lemma 7.1) that $\sup_{x}\Phi^{-1}\left(\widetilde{F}_{j}(x)\right)\leq\sqrt{2\log n}$ when $\delta_{n}=1/(4n^{1/4}\sqrt{\pi\log n})$.

In the next theorem, we do not assume the true model is nonparanormal and define the population and sample risks as

. Suppose that $p\leq e^{n^{\xi}}$ for some $\xi<1$, and define the classes

Hence the Winsorized estimator of $(f,\Omega)$ with $\delta_{n}=1/(4n^{1/4}\sqrt{\pi\log n})$ is persistent over ${\mathcal{C}}_{n}$ when $L_{n}=o\left(n^{(1-\xi)/2}/\sqrt{\log n}\right)$.

The proofs of Theorems 5.1 and 5.7 are given in Section 7.

Experimental Results

Note that we can reuse the glasso implementation to fit a sparse nonparanormal. In particular, after computing the Winsorized sample covariance $S_{n}(\widetilde{f})$, we pass this matrix to the glasso routine to carry out the optimization

We begin by describing a procedure to generate graphs as in (Meinshausen and Bühlmann, 2006), with respect to which several distributions can then be defined. We generate a $p$-dimensional sparse graph $G\equiv(V,E)$ as follows: Let $V=\{1,\ldots,p\}$ corresponding to variables $X=(X_{1},\ldots,X_{p})$. We associate each index $j$ with a point $(Y_{j}^{(1)},Y_{j}^{(2)})\in^{2}$ where

for $k=1,2$. Each pair of nodes $(i,j)$ is included in the edge set $E$ with probability

where $y_{i}\equiv(y^{(1)}_{i},y^{(2)}_{i})$ is the observation of $(Y^{(1)}_{i},Y^{(2)}_{i})$ and $\|\cdot\|_{n}$ represents the Euclidean distance. Here, $s=0.125$ is a parameter that controls the sparsity level of the generated graph. We restrict the maximum degree of the graph to be four and build the inverse covariance matrix $\Omega_{0}$ according to

where the value $0.245$ guarantees positive definiteness of the inverse covariance matrix.

Given $\Omega_{0}$, $n$ data points are sampled from

where $\mu_{0}=(1.5,\ldots,1.5)$, $\Sigma_{0}=\Omega^{-1}_{0}$. For simplicity, the transformations functions for all dimensions are the same $f_{1}=\ldots=f_{p}=f_{0}$. To sample data from the nonparanormal distribution, we also need $g_{0}\equiv f^{-1}_{0}$, two different transformations $g_{0}$ are employed:

Next, we define the following two transformation families:

. (Gaussian CDF Transformation) Let $g_{0}$ be a one-dimensional Gaussian cumulative distribution function with mean $\mu_{g_{0}}$ and the standard deviation $\sigma_{g_{0}}$, i.e.,

We define the transformation function $g_{j}=f_{j}^{-1}$ for the $j$-th dimension as

. (Symmetric Power Transformation) Let $g_{0}$ be the symmetric and odd transformation given by

where $\alpha>0$ is a parameter. We define the power transformation for the $j$-th dimension as

These transformation are constructed to preserve the marginal mean and standard deviation. In the following experiments, we refer to them as the cdf transformation and the power transformation, respectively. For the cdf transformation, we set $\mu_{g_{0}}=0.05$ and $\sigma_{g_{0}}=0.4$. For the power transformation, we set $\alpha=3$.

To visualize these two transformations, we sample $5000$ data points from a one-dimensional normal distribution ${N}(0.5,1.0)$ and then apply the above two transformations; the results are shown in Figure 3. It can be seen how the cdf and power transformations map a univariate normal distribution into a highly skewed and a bi-modal distribution, respectively.

To generate synthetic data, we set $p=40$, resulting in $\binom{40}{2}+40=820$ parameters to be estimated, and vary the sample sizes from $n=200$ to $n=1000$. Three conditions are considered, corresponding to using the cdf transform, the power transform, or no transformation. In each case, both the glasso and the nonparanormal are applied to estimate the graph.

We choose a set of regularization parameters $\Lambda$; for each $\lambda\in\Lambda$, we obtain an estimate $\widehat{\Omega}_{n}$ which is a $40\times 40$ matrix. The upper triangular matrix has 780 parameters; we can vectorize it to get a 780-dimensional parameter vector. A regularization path is trace of these parameters over all the regularization parameters within $\Lambda$. The regularization paths for both methods are plotted in Figure 4. For the cdf transformation and the power transformation, the nonparanormal separates the relevant and the irrelevant dimensions very well. For the glasso, relevant variables are mixed with irrelevant variables. If no transformation is applied, the paths for both methods are almost the same.

A.2 Estimated transformations

For sample size $n=1000$, we plot the estimated transformations for three of the variables in Figure 5. It is clear that Winsorization plays a significant role for the power transformation. This is intuitive due to the high skewness of the nonparanormal distribution resulting from the power transformations.

A.3 Quantitative comparison

To evaluate the performance for structure estimation quantitatively, we use false positive and false negative rates. Let $G=(V,E)$ be a $p$-dimensional graph (which has at most $\binom{p}{2}$ edges) in which there are $|E|=r$ edges, and let $\widehat{G}^{\lambda}=(V,\widehat{E}^{\lambda})$ be an estimated graph using the regularization parameter $\lambda$. The number of false positives at $\lambda$ is

The number of false negatives at $\lambda$ is defined as

The oracle regularization level $\lambda^{*}$ is then

To illustrate the overall performance of these two methods over the full paths, ROC curves are shown in Figure 7, using

The curves clearly show how the performance of both methods improves with sample size, and that the nonparanormal is superior to the Gaussian model in most cases.

A.4 Visualization of typical runs

Figure 8 shows typical runs for the cdf and power transformations. It’s clear that when the glasso estimates the graph incorrectly, the mistakes include both false positives and negatives.

B Gene microarray data

In this study, we consider a dataset based on Affymetrix GeneChip microarrays for the plant Arabidopsis thaliana, (Wille, 2004). The sample size is $n=118$. The expression levels for each chip are pre-processed by log-transformation and standardization. A subset of 40 genes from the isoprenoid pathway are chosen, and we study the associations among them using both the paranormal and nonparanormal models. Even though these data are generally treated as multivariate Gaussian in the previous analysis (Wille, 2004), our study shows that the results of the nonparanormal and the glasso are very different over a wide range of regularization parameters. This suggests the nonparanormal could support different scientific conclusions.

We first compare the regularization paths of the two methods, in Figure 9. To generate the paths, we select 50 regularization parameters on an evenly spaced grid in the interval $[0.16,1.2]$. Although the paths for the two methods look similar, there are some subtle differences. In particular, variables become nonzero in a different order, especially when the regularization parameter is in the range $\lambda\in[0.2,0.3]$. As shown below, these subtle differences in the paths lead to different model selection behaviors.

B.2 Comparison of the selected graphs

Figure 11 compares the estimated graphs for the two methods at several values of the regularization parameter $\lambda$ in the range $[0.16,0.37]$. For each $\lambda$, we show the estimated graph from the nonparanormal in the first column. In the second column we show the graph obtained by scanning the full regularization path of the glasso fit and finding the graph having the smallest symmetric difference with the nonparanormal graph. The symmetric difference graph is shown in in the third column. The closest glasso fit is different, with edges selected by the glasso not selected by the nonparanormal, and vice-versa. Several estimated transformations are plotted in Figure 11, which are are nonlinear. Interestingly, several of the differences between the fitted graphs are related to these variables.

Proofs

We assume, without loss of generality from Lemma 3.3, that $\mu_{j}=0$ and $\sigma_{j}=1$ for all $j=1,\ldots,p$. Thus, define $\widetilde{f}_{j}(x)\equiv\Phi^{-1}(\widetilde{F}_{j}(x))$ and ${f}_{j}(x)\equiv\Phi^{-1}({F}_{j}(x))$, and let $g_{j}\equiv f_{j}^{-1}$.

We start with some useful lemmas; the first is from Abramovich et al. (2006).

. ​(Gaussian Distribution function vs. Quantile function) Let $\Phi$ and $\phi$ denote the distribution and density functions of a standard Gaussian random variable. Then

. ​(Distribution function of the transformed random variable) For any $\alpha\in(-\infty,\infty)$

which holds for any $t$. $\;\;\scriptstyle\Box$

. ​(Gaussian maximal inequality) Let $W_{1},\ldots,W_{n}$ be independently and identically distributed standard Gaussian random variables. Then for any $\alpha>0$

from which the result follows. $\;\;\scriptstyle\Box$

. For any $\alpha>0$ that satisfies $1-\delta_{n}-\Phi\left(\sqrt{\alpha\log n}\right)>0$ for all $n$, we have

Equation (22) then follows from equation (21). The proof of equation (23) uses the same argument. $\;\;\scriptstyle\Box$

Now let $M>2$ be some constant and set $\displaystyle\beta=\frac{1}{2}$. We split the interval

The behaviors of the function estimates in these two regions is different, and so we first establish bounds on the probability that a sample can fall in the end region ${\mathcal{E}}_{n}$.

. Let $\displaystyle A\equiv\sqrt{\frac{2}{\pi}}(\sqrt{M}-\sqrt{\beta})$. Then

Proof. Using equation (21) and the mean value theorem, we have

The result of the lemma follows directly. $\;\;\scriptstyle\Box$

We next bound the error of the Winsorized estimate of a component function over the end region.

Proof. From Lemma 7.2 and the definition of ${\mathcal{E}}_{n}$, we have

Given the fact that $\delta_{n}=\displaystyle\frac{1}{4n^{1/4}\sqrt{\pi\log n}}$, we have $\widetilde{F}_{j}(t)\in\displaystyle\left(\frac{1}{n},1-\frac{1}{n}\right)$. Therefore, from equation (20),

The result follows from the triangle inequality and $\sqrt{M}+\sqrt{2}\leq\sqrt{2(M+2)}$. $\;\;\scriptstyle\Box$

We only need to analyze the rate for the first term above, since the second one is of higher order (Cai et al., 2008). Let

with $c_{1}$ as a generic positive constant. Therefore

Thus, we only need to carry out our analysis on the event $\mathcal{A}_{n}$. On this event, we have the following decomposition:

We now analyze each of these terms separately.

. On the event $\mathcal{A}_{n}$, let $\beta=1/2$ and $\displaystyle\epsilon\geq C(M,\xi)\sqrt{\frac{\log p\log^{2}n}{n^{1/2}}}$, then

with the same parameter $A$ as in Lemma 7.5. Such a $\theta_{1}$ guarantees that

Using the Bernstein’s inequality, for $\beta=\displaystyle\frac{1}{2}$,

where $c_{1},c_{2},c_{3}>0$ are generic constants.

By adding and subtracting terms $f_{j}(t)$ and $f_{s}(t)$, we have

The first term can further be decomposed to be

Also, from the definition of $\mathcal{E}_{n}$, we have

Since $\displaystyle\epsilon\geq C(M,\xi)\sqrt{\frac{\log p\log^{2}n}{n^{1/2}}}$, we have

The claim of the lemma then follows directly. $\;\;\scriptstyle\Box$

. From the above analysis, we see that the data in the tails doesn’t affect the rate. Using exactly the same argument, we can also show that

. On the event $\mathcal{A}_{n}$, let $\beta=1/2$ and $\displaystyle\epsilon\geq C(M,\xi)\sqrt{\frac{\log p\log^{2}n}{n^{1/2}}}$. There exist generic constants $c_{1},c_{2},c_{3},c_{4}$, such that

Since $\delta_{n}=\displaystyle\frac{1}{4n^{\beta/2}\sqrt{2\pi\beta\log n}}$, using Mill’s inequality we have

From (24) and (25), it is easy to see that

where $c_{3}$ and $c_{4}$ are generic positive constants.

From the definition of $\widetilde{F}_{j}$, we have

From equation (21) and the fact that $1-\delta_{n}\geq\Phi\left(\sqrt{\beta\log n}\right)$, we have that

Finally, using the Dvoretzky-Kiefer-Wolfowitz inequality,

where $c_{1},c_{2},c_{3},c_{4}$ are generic constants. $\;\;\scriptstyle\Box$

The conclusion of Theorem 5.1 follows from Lemma 7.7 and Lemma 7.9.

B Proof of Theorem 5.7

Proof. First note that the population and sample risks are

Therefore, for all $(f,\Omega)\in{\mathcal{M}}_{n}^{p}\oplus{\mathcal{C}}_{n}$, we have

Now, if ${\mathcal{F}}$ is a class of functions, we have

where $\log N_{[\,]}(\epsilon,{\mathcal{F}})$ is the bracketing entropy. For the class of one dimensional, bounded and monotone functions, the bracketing entropy satisfies

for some $K>0$ (van der Vaart and Wellner, 1996).

Now, let ${\mathcal{P}}_{n,p}$ be the class of all functions of the form $m(x)=f_{j}(x_{j})f_{k}(x_{k})$ for $j,k\in\{1,\ldots,p\}$, where $f_{j}\in{\mathcal{M}}_{n}$ for each $j$. Then the bracketing entropy satisfies

and the bracketing integral satisfies $J_{[\,]}(C\sqrt{\log n},{\mathcal{P}}_{n,p})=O(\sqrt{\log n\log p})$. It follows from (26) and Markov’s inequality that

and the conclusion follows. $\;\;\scriptstyle\Box$

Concluding Remarks

In this paper we have introduced the nonparanormal, a type of Gaussian copula with nonparametric marginals that is suitable for estimating high dimensional undirected graphs. The nonparanormal can be viewed as an extension of sparse additive models to the setting of graphical models. We proposed an estimator for the component functions that is based on thresholding the tails of the empirical distribution function at appropriate levels. A theoretical analysis was given to bound the difference between the sample covariance with respect to these estimated functions and the true sample covariance. This analysis was leveraged with the recent work of Ravikumar et al. (2009) and Rothman et al. (2008) to obtain consistency results for the nonparanormal. Computationally, fitting a high dimensional nonparanormal is no more difficult than estimating a multivariate Gaussian, and indeed one can exploit existing software for the graphical lasso. Our experimental results indicate that the sparse nonparanormal can give very different results than a sparse Gaussian graphical model, suggesting that it may be a useful tool for relaxing the normality assumption, which is often made only for convenience.

References