Regularized rank-based estimation of high-dimensional nonparanormal graphical models

Introduction.

Estimating covariance or precision matrices is of fundamental importance in multivariate statistical methodologies and applications. In particular, when data follow a joint normal distribution, $\mathbf{X}=(X_{1},\ldots,X_{p})\sim N_{p}(\bolds{\mu},\bolds{\Sigma}),$ the precision matrix $\bolds{\Theta}=\bolds{\Sigma}^{-1}$ can be directly translated into a Gaussian graphical model. The Gaussian graphical model serves as a noncausal structured approach to explore the complex systems consisting of Gaussian random variables, and it finds many interesting applications in areas such as gene expression genomics and macroeconomics determinants study [Friedman (2004), Wille et al. (2004), Dobra, Eicher and Lenkoski (2010)]. The precision matrix plays a critical role in the Gaussian graphical models because the zero entries in $\bolds{\Theta}=(\theta_{ij})_{p\times p}$ precisely capture the desired conditional independencies, that is, $\theta_{ij}=0$ if and only if $X_{i}\perp\!\!\!\!\perp X_{j}|\mathbf{X}\setminus\{X_{i},X_{j}\}$ [Lauritzen (1996), Edwards (2000)].

Although the normality assumption can be relaxed if we only focus on estimating a precision matrix, it plays an essential role in making the neat connection between a sparse precision matrix and a sparse Gaussian graphical model. Without normality, we ought to be very cautious when translating a good sparse precision matrix estimator into an interpretable sparse Gaussian graphical model. However, the normality assumption often fails in reality. For example, the observed data are often skewed or have heavy tails. To illustrate the issue of nonnormality in real applications, let us consider the gene expression data to construct isoprenoid genetic regulatory network in Arabidposis thaliana [Wille et al. (2004)], including 16 genes from the mevalonate (MVA) pathway in the cytosolic, 18 genes from the plastidial (MEP) pathway in the chloroplast and 5 encode proteins in the mitochondrial. This dataset contains gene expression measurements of 39 genes assayed on $n=118$ Affymetrix GeneChip microarrays. This dataset was analyzed by Wille et al. (2004), Li and Gui (2006) and Drton and Perlman (2007) in the context of Gaussian graphical modeling after taking the log-transformation of the data. However, the normality assumption is still inappropriate even after the log-transformation. To show this, we conduct the normality test at the significance level of $0.05$ as in Table 1. It is clear that at most 9 out of 39 genes would pass any of three normality tests. Even after log-transformation, at least $60\%$ genes reject the null hypothesis of normality. With Bonferroni correction there are still over $30\%$ genes that fail to pass any normality test. Figure 1 plots the histograms of two key isoprenoid genes MECPS in the MEP pathway and MK in the MVA pathway after the log-transformation, clearly showing the nonnormality of the data after the log-transformation.

Using transformation to achieve normality is a classical idea in statistical modeling. The celebrated Box–Cox transformation is widely used in regression analysis. However, any parametric modeling of the transformation suffers from model mis-specification which could lead to misleading inference. In this paper we take a nonparametric transformation strategy to handle the nonnormality issue. Let $F(\cdot)$ be the CDF of a continuous random variable $X$ and $\Phi^{-1}(\cdot)$ be the inverse of the CDF of $N(0,1)$. Consider the transformation from $X$ to $Z$ by $Z=\Phi^{-1}(F(X))$. Then it is easy to see that $Z$ is standard normal regardless of $F$. Motivated by this simple fact, we consider modeling the data by the following nonparanormal model:

The nonparanormal model: $\mathbf{X}=(X_{1},\ldots,X_{p})$ follows a $p$-dimensional nonparanormal distribution if there exists a vector of unknown univariate monotone increasing transformations, denoted by $\mathbf{f}=(f_{1},\ldots,f_{p})$, such that the transformed random vector follows a multivariate normal distribution with mean 0 and covariance $\bolds{\Sigma}$,

where without loss of generality the diagonals of $\bolds{\Sigma}$ are equal to 1.

Note that model (1) implies that $f_{j}(X_{j})$ is a standard normal random variable. Thus, $f_{j}$ must be $\Phi^{-1}\circ F_{j}$ where $F_{j}$ is the CDF of $X_{j}$. The marginal normality is always achieved by transformations, so model (1) basically assumes that those marginally normal-transformed variables are jointly normal as well. We follow Liu, Lafferty and Wasserman (2009) to call model (1) the nonparanormal model, but model (1) is in fact a semiparametric Gaussian copula model. The semiparametric Gaussian copula model is a nice combination of flexibility and interpretability, and it has generated a lot of interests in statistics, econometrics and finance; see Song (2000), Chen and Fan (2006), Chen, Fan and Tsyrennikov (2006) and references therein. Let $\mathbf{Z}=(Z_{1},\ldots,Z_{p})=(f_{1}(X_{1}),\ldots,f_{p}(X_{p}))$. By the joint normality assumption of $\mathbf{Z}$, we know that $\theta_{ij}=0$ if and only if $Z_{i}\perp\!\!\!\!\perp Z_{j}|\mathbf{Z}\setminus\{Z_{i},Z_{j}\}.$ Interestingly, we have that

Therefore, a sparse $\bolds{\Theta}$ can be directly translated into a sparse graphical model for presenting the original variables.

In this work we primarily focus on estimating $\bolds{\Theta}$ which is then used to construct a nonparanormal graphical model. As for the nonparametric transformation function, by the expression $f_{j}=\Phi^{-1}\circ F_{j}$, we have a natural estimator for the transformation function of the $j$th variable as $\hat{f}_{j}=\Phi^{-1}\circ\hat{F}^{+}_{j}$ where $\hat{F}^{+}_{j}$ is a Winsorized empirical CDF of the $j$th variables. Note that the Winsorization is used to avoid infinity value and to achieve better bias-variance tradeoff; see Liu, Lafferty and Wasserman (2009) for detailed discussion. In this paper we show that we can directly estimate $\bolds{\Theta}$ without estimating these nonparametric transformation functions at all. This statement seems to be a bit surprising because a natural estimation scheme is a two-stage procedure: first estimate $f_{j}$ and then apply a well-developed sparse Gaussian graphical model estimation method to the transformed data $\hat{\mathbf{z}}_{i}=\mathbf{\hat{f}}(\mathbf{x}_{i}),1\leq i\leq n$. Liu, Lafferty and Wasserman (2009) have actually studied this “plug-in” estimation approach. They proposed a Winsorized estimator of the nonparametric transformation function and used the graphical lasso in the second stage. They established convergence rate of the “plug-in” estimator when $p$ is restricted to a polynomial order of $n$. However, Liu, Lafferty and Wasserman (2009) did not get a satisfactory rate of convergence for the “plug-in” approach, because the rate of convergence can be established for the Gaussian graphical model even when $p$ grows with $n$ almost exponentially fast [Ravikumar et al. (2011)]. As noted in Liu, Lafferty and Wasserman (2009), it is very challenging, if not impossible, to push the theory of the “plug-in” approach to handle exponentially large dimensions. One might ask if using a better estimator for the transformation functions could improve the rate of convergence such that $p$ could be allowed to be nearly exponentially large relative to $n$. This is a legitimate direction for research. We do not pursue this direction in this work. Instead, we show that we could use a rank-based estimation approach to achieve the exact same goal without estimating these transformation functions at all.

Our estimator is constructed in two steps. First, we propose using the adjusted Spearman’s rank correlation to get a nonparametric sample estimate of $\bolds{\Sigma}$. As the second step, we compute a sparse estimator $\bolds{\Theta}$ from the rank-based sample estimate of $\bolds{\Sigma}$. For that purpose, we consider several regularized rank estimators, including the rank-based graphical lasso, the rank-based neighborhood Dantzig selector and the rank-based CLIME. The complete methodological details are presented in Section 2. In Section 3 we establish theoretical properties of the proposed rank-based estimators, regarding both precision matrix estimation and graphical model selection. In particular, we are motivated by the theory to consider the adaptive version of the rank-based neighborhood Dantzig selector and the rank-based CLIME, which can select the true support set with an overwhelming probability without assuming a stringent irrepresentable condition required for the oracle and rank-based graphical lasso. Section 4 contains numerical results and Section 5 has some concluding remarks. Technical proofs are presented in an Appendix.

A referee informed us in his/her review report that Liu et al. (2012) also independently used the rank-based correlation in the context of nonparametric Gaussian graphical model estimation. A major focus in Liu et al. (2012) is the numerical demonstration of the robustness property of the rank-based methods using both Spearman’s rho and Kendall’s tau when data are contaminated. In the present paper we provide a systematic analysis of the rank-based estimators, and our theoretical analysis further leads to the rank-based adaptive Dantizg selector and the rank-based adaptive CLIME in order to achieve improved sparsity recovery properties. Our theoretical analysis of the rank-based adaptive Dantizg selector is of independent interest. Although the theory is established for the rank-based estimators using Spearman’s rho, the same analysis can be easily adopted to prove the theoretical properties of the rank-based estimators using Kendall’s tau rank correlation.

Methodology.

Suppose an oracle knows the underlying transformation vector; then the oracle could easily recover “oracle data” by applying these true transformations, that is, $\mathbf{z}_{i}=\mathbf{f}(\mathbf{x}_{i}),1\leq i\leq n$. Before presenting our rank-based estimators, it is helpful to revisit the “oracle” procedures that are defined based on the “oracle data.”

The oracle neighborhood lasso selection. Under the nonparanormal model, for each $k=1,\ldots,p$, the “oracle” variable $Z_{k}$ given $\mathbf{Z}_{(k)}$ is normally distributed as $N(\mathbf{Z}_{(k)}^{T}\bolds{\Sigma}_{(k)}^{-1}\bolds{\sigma}_{(k)},1-\bolds{\sigma}_{(k)}^{T}\bolds{\Sigma}_{(k)}^{-1}\sigma_{(k)}),$ which can be written as $Z_{k}=\mathbf{Z}_{(k)}^{T}\bolds{\beta}_{k}+\varepsilon_{k}$ with $\bolds{\beta}_{k}=\bolds{\Sigma}_{(k)}^{-1}\bolds{\sigma}_{(k)}$ and $\varepsilon_{k}\sim N(0,1-\bolds{\sigma}_{(k)}^{T}\bolds{\Sigma}_{(k)}^{-1}\bolds{\sigma}_{(k)}).$ Notice that $\bolds{\beta}_{k}$ and $\varepsilon_{k}$ are closely related to the precision matrix $\bolds{\Theta}$, that is, $\theta_{kk}=1/\operatorname{Var}(\varepsilon_{k})$ and $\bolds{\theta}_{(k)}=-\bolds{\beta}_{k}/\operatorname{Var}(\varepsilon_{k}).$ Thus for the $k$th variable, $\bolds{\theta}_{(k)}$ and $\bolds{\beta}_{k}$ share the same sparsity pattern. Following Meinshausen and Bühlmann (2006), the oracle neighborhood lasso selection obtains the solution $\hat{\bolds{\beta}}^{o}_{k}$ from the following lasso penalized least squares problem:

and then the sparsity pattern of $\bolds{\Theta}$ can be estimated by aggregating the neighborhood support set of $\hat{\bolds{\beta}}^{o}_{k}=(\hat{\beta}^{o}_{jk})_{j\neq k}$ ($\widehat{ne}_{k}=\{j\dvtx\hat{\beta}^{o}_{jk}\neq 0\}$) via intersection or union.

Then (3) can be written in the following equivalent form:

The oracle neighborhood Dantzig selector. Following Yuan (2010) the lasso least squares in (3) can be replaced with the Dantzig selector

Then the sparsity pattern of $\bolds{\Theta}$ can be similarly estimated by aggregating via intersection or union. Furthermore, we notice that

Then (5) can be written in the following equivalent form:

Cai, Liu and Luo (2011) compared the CLIME with the graphical lasso, and showed that the CLIME enjoys nice theoretical properties without assuming the irrepresentable condition of Ravikumar et al. (2011) for the graphical lasso.

2 The proposed rank-based estimators.

The existing theoretical results in the literature can be directly applied to these oracle estimators. However, the “oracle data” $\mathbf{z}_{1},\mathbf{z}_{2},\ldots,\mathbf{z}_{n}$ are unavailable and thus the above-mentioned “oracle” procedures are not genuine estimators. Naturally we wish to construct a genuine estimator that can mimic the oracle estimator. To this end, we can derive an alternative estimator of $\bolds{\Sigma}$ based on the actual data $\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n}$ and then feed this genuine covariance estimator to the graphical lasso, the neighborhood selection or CLIME. To implement this natural idea, we propose a rank-based estimation scheme. Note that $\bolds{\Sigma}$ can be viewed as the correlation matrix as well, that is, $\sigma_{ij}=\operatorname{corr}(\mathbf{z}_{i},\mathbf{z}_{j})$. Let ($x_{1i},x_{2i},\ldots,x_{ni}$) be the observed values of variable $X_{i}$. We convert them to ranks denoted by $\mathbf{r}_{i}=(r_{1i},r_{2i},\ldots,r_{ni})$. Spearman’s rank correlation $\hat{r}_{ij}$ is defined as Pearson’s correlation between $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$. Spearman’s rank correlation is a nonparametric measure of dependence between two variables. It is important to note that $\mathbf{r}_{i}$ are the ranks of the “oracle” data. Therefore, $\hat{r}_{ij}$ is also identical to the Spearman’s rank correlation between the “oracle” variables $Z_{i},Z_{j}$. In other words, in the framework of rank-based estimation, we can treat the observed data as the “oracle” data and avoid estimating $p$ nonparametric transformation functions. We make a note here that one may consider other rank correlation measures such as Kendall’s tau correlation. To fix the idea we use Spearman’s rank correlation throughout this paper.

The nonparanormal model implies that $(Z_{i},Z_{j})$ follows a bivariate normal distribution with correlation parameter $\sigma_{ij}$. Then a classical result due to Kendall (1948) [see also Kruskal (1958)] shows that

which indicates that $\hat{r}_{ij}$ is a biased estimator of $\sigma_{ij}$. To correct the bias, Kendall (1948) suggested using the adjusted Spearman’s rank correlation

Combining (8) and (9) we see that $\hat{r}^{s}_{ij}$ is an asymptotically unbiased estimator of $\sigma_{ij}$. Naturally we define the rank-based sample estimate of $\bolds{\Sigma}$ as follows:

In Section 3 we show $\hat{\mathbf{R}}^{s}$ is a good estimator of $\bolds{\Sigma}$. Then we naturally come up with the following rank-based estimators of $\bolds{\Theta}$ by using the graphical lasso, the neighborhood Dantzig selector and CLIME:

The rank-based neighborhood Dantzig selector: A rank-based estimate of $\bolds{\beta}_{k}$ can be solved by

The support of $\bolds{\Theta}$ can be estimated from the support of $\hat{\bolds{\beta}}_{1}^{s.nd},\ldots,\hat{\bolds{\beta}}_{p}^{s.nd}$ via aggregation by union or intersection. We can also construct the rank-based precision matrix estimator $\hat{\bolds{\Theta}}_{nd}^{s}=(\hat{\theta}^{s.nd}_{ij})_{1\leq i,j\leq p}$ with

($k=1,\ldots,p$). We can symmetrize $\hat{\bolds{\Theta}}_{nd}^{s}$ by solving the following optimization problem [Yuan (2010)]:

Theoretical analysis of the rank-based neighborhood Dantzig selector in Section 3.2 motivated us to consider using the adaptive Dantzig selector in the rank-based neighborhood estimation in order to achieve better graphical model selection performance. See Section 3.2 for more details.

for $k=1,\ldots,p$. Then $\hat{\bolds{\Theta}}_{c}^{s}$ is exactly equivalent to $(\hat{\bolds{\theta}}^{s.c}_{1},\ldots,\hat{\bolds{\theta}}^{s.c}_{p})$. Note that $\hat{\bolds{\Theta}}_{c}^{s}$ could be asymmetric. Following Cai, Liu and Luo (2011) we consider

with $\breve{\theta}^{s.c}_{ij}=\hat{\theta}^{s.c}_{ij}I_{\{|\hat{\theta}^{s.c}_{ij}|\leq|\hat{\theta}^{s.c}_{ji}|\}}+\hat{\theta}^{s.c}_{ji}I_{\{|\hat{\theta}^{s.c}_{ij}|>|\hat{\theta}^{s.c}_{ji}|\}}.$ In the original paper Cai, Liu and Luo (2011) proposed to use hard thresholding for graphical model selection. Borrowing the basic idea from Zou (2006), we propose an adaptive version of the rank-based CLIME in order to achieve better graphical model selection. See Section 3.3 for more details.

3 Rank-based neighborhood lasso?

One might consider the rank-based neighborhood lasso defined as follows:

However, there is a technical problem for the above definition. The Spearman’s rank correlation matrix $\hat{\mathbf{R}}$ is always positive semidefinite, but the adjusted correlation matrix $\hat{\mathbf{R}}^{s}$ could become indefinite. To our best knowledge, Devlin, Gnanadesikan and Kettenring (1975) were the first to point out the indefinite issue of the estimated rank correlation matrix. Here we also use a toy example to illustrate this point. Consider the $3\times 3$ correlation matrix

Note that $\mathbf{A}$ is positive-definite with eigenvalues $1.99$, $1.00$ and $0.01$, but $2\sin(\frac{\pi}{6}\mathbf{A})$ becomes indefinite with eigenvalues $2.01$, $1.00$ and $-0.01$. The negative eigenvalues will make (14) an ill-defined optimization problem. Fortunately, the positive definite issue does not cause any problem for the graphical lasso, Dantzig selector and CLIME. Notice that the diagonal elements of $\hat{\mathbf{R}}^{s}$ are obviously strictly positive, and thus Lemma 3 in Ravikumar et al. (2011) suggests that the rank-based graphical lasso always has a unique positive definite solution for any regularization parameter $\lambda>0$. The rank-based neighborhood Dantzig selector and the rank-based CLIME are still well defined, even when $\hat{\mathbf{R}}^{s}_{(k)}$ becomes indefinite, and the according optimization algorithms also tolerate the indefiniteness of $\hat{\mathbf{R}}^{s}_{(k)}$. One might consider a positive definite correction of $\hat{\mathbf{R}}^{s}$ for implementing the neighborhood lasso estimator. However, the resulting estimator shall behave similarly to the rank-based neighborhood Dantzig selector because the lasso penalized least squares and Dantzig selector, in general, work very similarly [Bickel, Ritov and Tsybakov (2009), James, Radchenko and Lv (2009)].

Theoretical properties.

In this section we establish theoretical properties for the proposed rank-based estimators. The main conclusion drawn from these theoretical results is that the rank-based graphical lasso, neighborhood Dantzig selector and CLIME work as well as their oracle counterparts in terms of the rates of convergence. We first provide useful concentration bounds concerning the accuracy of the rank-based sample correlation matrix.

Fix any $0<\varepsilon<1$, and let $n\geq\frac{12\pi}{\varepsilon}$. Then there exists some absolute constant $c_{0}>0$, and we have the following concentration bounds:

Lemma 1 is a key ingredient of our theoretical analysis. It basically shows that the rank-based sample estimator of $\bolds{\Sigma}$ works as well as the usual sample covariance estimator of $\bolds{\Sigma}$ based on the “oracle data.”

Element-wise maximal bound: if $\lambda$ is chosen such that

with probability at least $1-p^{2}\exp(-\frac{\kappa^{2}}{16}c_{0}n\lambda^{2})$, the rank-based graphical lasso estimator $\hat{\bolds{\Theta}}_{g}^{s}$ satisfies that $\hat{\theta}^{s.g}_{ij}=0$ for any $(i,j)\in\mathcal{A}^{c}$ and

Graphical model selection consistency: picking a regularization parameter $\lambda$ to satisfy that

then with probability at least $1-p^{2}\exp(-\frac{\kappa^{2}}{16}c_{0}n\lambda^{2})$, $\hat{\bolds{\Theta}}_{g}^{s}$ is sign consistent satisfying that $\operatorname{sign}(\hat{\theta}^{s.g}_{ij})=\operatorname{sign}(\theta_{ij}^{*})$ for any $(i,j)\in\mathcal{A}$ and $\hat{\theta}^{s.g}_{ij}=0$ for any $(i,j)\in\mathcal{A}^{c}$.

Rates of convergence: assume $n\gg d^{2}\log p$, and pick a regularization parameter $\lambda$ such that $d^{-1}\gg\lambda=O(({\log p}/n)^{1/2})$. Then we have

Graphical model selection consistency: assume $\psi_{\min}$ is also fixed and $n\gg d^{2}\log p$. Pick a $\lambda$ such that $d^{-1}\gg\lambda=O(({\log p}/n)^{1/2}).$ Then we have $\operatorname{sign}(\hat{\theta}^{s.g}_{ij})=\operatorname{sign}(\theta_{ij}^{*})$, $\forall(i,j)\in\mathcal{A}$ and $\operatorname{sign}(\hat{\theta}^{s.g}_{ij})=0$, $\forall(i,j)\in\mathcal{A}^{c}$.

Under the same conditions of Theorem 1 and Corollary 1, by the results in Ravikumar et al. (2011), we know that the conclusions of Theorem 1 and Corollary 1 hold for the oracle graphical lasso. In other words, the rank-based graphical lasso estimator is comparable to its oracle counterpart in terms of rates of convergence.

2 Rank-based neighborhood Dantzig selector.

Likewise we can partition $\bolds{\sigma}^{*}_{(k)}$ and $\bolds{\Sigma}_{(k)}^{*}$ with respect to $\mathcal{A}_{k}$.

Pick the $\lambda$ such that $d\lambda=o(1)$ and $bn\lambda\geq 12\pi M$. With probability at least $1-p^{2}\exp(-c_{0}\frac{b^{2}}{M^{2}}n\lambda^{2})$, there exists $C_{b,B,M}>0$ depending on $b$, $B$ and $M$ only such that

Suppose that $b$, $B$ and $M$ are all fixed. Let $n\gg d^{2}\log p$, and pick $\lambda$ such that $d^{-1}\gg\lambda=O(({\log p}/n)^{1/2})$. Then we have

Dantzig selector and the lasso are closely related [Bickel, Ritov and Tsybakov (2009), James, Radchenko and Lv (2009)]. Similarly to the lasso, the Dantzig selector tends to over-select. Zou (2006) proposed the adaptive weighting idea to develop the adaptive lasso which improves the selection performance of the lasso and corrects its bias too. The very same idea can be used to improve the selection performance of Dantzig selector which leads to the adaptive Dantzig selector [Dicker and Lin (2009)]. We can extend the rank-based Dantzig selector to the rank-based adaptive Dantzig selector. Given adaptive weights $\mathbf{w}_{k}$, consider

where $\circ$ denotes the Hadamard product, and $\mathbf{a}_{d\times 1}\leq\mathbf{b}_{d\times 1}$ denotes the set of entrywise inequalities $a_{i}\leq b_{i}$ for ease of notation. In both our theoretical analysis and numerical implementation, we utilize the optimal solution $\hat{\bolds{\beta}}_{k}^{s.nd}$ of the rank-based Dantzig selector to construct the adaptive weights $\mathbf{w}_{k}$ by

For each $k$, we pick $\lambda=\lambda_{d}$ as in (11) satisfying that $\lambda_{d}\geq\frac{12\pi M}{bn}$ and $o(1)=d\lambda_{d}\leq\min\{\frac{\psi_{k}}{2C_{0}},\frac{1}{4C_{0}d(\psi_{k}+2G_{k})}-\frac{1}{C_{0}n}\}$, and pick $\lambda=\lambda_{ad}$ as in (15) such that $\frac{\psi_{k}^{2}}{8G_{k}}\geq\lambda_{ad}\geq\max\{\frac{12\pi}{n},(C_{0}d\lambda_{d}+\frac{1}{n})\frac{H_{k}\psi_{k}}{G_{k}}\},$ and $o(1)=d\lambda_{ad}\leq\min\{\lambda_{\min}(\bolds{\Sigma}^{*}_{\mathcal{A}_{k}\mathcal{A}_{k}}),\frac{1}{2G_{k}},\frac{\psi_{k}}{8G_{k}(\psi_{k}+G_{k})}\}.$ In addition, we also choose $\mathbf{w}_{k}=\mathbf{w}_{k}^{d}$ as in (16) for each $k$. Then with a probability at least $1-p^{2}\exp(-c_{0}n\cdot\min\{\lambda_{ad}^{2},\frac{b^{2}}{M^{2}}\lambda_{d}^{2}\})$, for each $k$, the rank-based adaptive Dantzig selector finds the unique solution $\hat{\bolds{\beta}}^{s.nad}_{k}=(\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}},\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}^{c}})$ with $\operatorname{sign}(\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}})=\operatorname{sign}(\bolds{\beta}^{*}_{\mathcal{A}_{k}})$ and $\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}^{c}}=\mathbf{0}$, and thus the rank-based neighborhood adaptive Dantzig selector is consistent for the graphical model selection.

Suppose $b$, $B$, $M$, $\psi_{k}$, $G_{k}$ and $H_{k}$ ($1\leq k\leq p$) are all constants. Assume that $n\gg d^{4}\log p$ and $\lambda_{\min}(\bolds{\Sigma}^{*}_{\mathcal{A}_{k}\mathcal{A}_{k}})\gg d^{2}({\log p}/n)^{1/2}$. Pick the tuning parameters $\lambda_{d}$ and $\lambda_{ad}$ such that $\frac{1}{d}\gg\lambda_{d}=O(({\log p}/n)^{1/2})$ and $\min\{\frac{1}{d}\cdot\lambda_{\min}(\bolds{\Sigma}^{*}_{\mathcal{A}_{k}\mathcal{A}_{k}}),\frac{1}{d}\}\gg\lambda_{ad}\gg d\lambda_{d}$. Then with probability tending to $1$, for each $k$, the rank-based adaptive Dantzig selector with $\mathbf{w}_{k}=\mathbf{w}_{k}^{d}$ as in (16) finds the unique optimal solution $\hat{\bolds{\beta}}^{s.nad}_{k}=(\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}},\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}^{c}})$ with $\operatorname{sign}(\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}})=\operatorname{sign}(\bolds{\beta}^{*}_{\mathcal{A}_{k}})$ and $\hat{\bolds{\beta}}^{s.nad}_{\mathcal{A}_{k}^{c}}=\mathbf{0}$, and thus the rank-based neighborhood adaptive Dantzig selector is consistent for the graphical model selection.

The sign-consistency of the adaptive Dantzig selector is similar to that of the adaptive lasso [van de Geer, Bühlmann and Zhou (2011)]. Based on Theorem 2 we construct the adaptive weights in (16) which is critical for the success of the rank-based adaptive Dantzig selector in the high-dimensional setting. It is important to mention that the rank-based adaptive Dantzig selector does not require the strong irrepresentable condition for the rank-based graphical lasso to have the sparsity recovery property. Our treatment of the adaptive Dantzig selector is fundamentally different from Dicker and Lin (2009). Dicker and Lin (2009) focused on the canonical linear regression model and constructed the adaptive weights as the inverse of the absolute values of ordinary least square estimator. Their theoretical results only hold in the classical fixed $p$ setting. In our problem $p$ can be much bigger than $n$. The choice of adaptive weights in (16) plays a critical role in establishing the graphical model selection consistency for the adaptive Dantzig selector under the high-dimensional setting where $p$ is at a nearly exponential rate to $n$. Our technical analysis uses some key ideas such as the strong duality and the complementary slackness from the linear optimization theory [Bertsimas and Tsitsiklis (1997), Boyd and Vandenberghe (2004)].

3 Rank-based CLIME.

Compared to the graphical lasso, the CLIME can enjoy nice theoretical properties without assuming the irrepresentable condition [Cai, Liu and Luo (2011)]. This continues to hold when comparing the rank-based graphical lasso and the rank-based CLIME.

Moreover, assume that $n\gg d^{2}\log p$, and suppose $M$ is a fixed constant. Pick a regularization parameter $\lambda$ satisfying $\lambda=O(({\log p}/n)^{1/2})$. Then we have

Theorem 4 is parallel to Theorem 6 in Cai, Liu and Luo (2011) which can be used to establish the rate of convergence of the oracle CLIME.

To improve graphical model selection performance, Cai, Liu and Luo (2011) suggested an additional thresholding step by applying the element-wise hard-thresholding rule to $\hat{\bolds{\Theta}}_{c}^{s}$,

where $\tau_{n}\geq 2M\lambda$ is the threshold, and $\lambda$ is given in Theorem 4. Here we show that consistent graphical model selection can be achieved by an adaptive version of the rank-based CLIME. Given an adaptive weight matrix $\mathbf{W}$ we define the rank-based adaptive CLIME as follows:

where $\mathbf{A}_{p\times p}\leq\mathbf{B}_{p\times p}$ is a simplified expression for the set of inequalities $a_{ij}\leq b_{ij}$ (for all $1\leq i,j\leq p$). Write $\mathbf{W}=(\mathbf{w}_{1},\ldots,\mathbf{w}_{p})$. By Lemma 1 in Cai, Liu and Luo (2011) the above linear programming problem in (18) is exactly equivalent to $p$ vector minimization subproblems,

In both our theory and implementation, we utilize the rank-based CLIME’s optimal solution $\hat{\bolds{\Theta}}_{c}^{s}$ to construct an adaptive weight matrix $\mathbf{W}$ by

The nice theoretical property of the rank-based CLIME allows us to construct the adaptive weights in (19), which is critical for establishing the graphical model selection consistency for the rank-based adaptive CLIME estimator in the high-dimensional setting without the strong ir-representable condition.

Numerical properties.

In this section we present both simulation studies and real examples to demonstrate the finite sample performance of the proposed rank-based estimators.

In the simulation study, we consider both Gaussian data and nonparanormal data. In models 1–4 we draw $n$ independent samples from $N_{p}(0,\bolds{\Sigma})$ with four different $\bolds{\Theta}$: {longlist}[Model 1:]

$\theta_{ii}=1$ and $\theta_{i,i+1}=0.5$;

$\theta_{ii}=1$, $\theta_{i,i+1}=0.4$ and $\theta_{i,i+2}=\theta_{i,i+3}=0.2$;

Randomly choose $16$ nodes to be the hub nodes in $\bolds{\Theta}$, and each of them connects with $5$ distinct nodes with $\bolds{\Theta}_{ij}=0.2$. Elements, not associated with hub nodes, are set as in $\bolds{\Theta}$. The diagonal element $\sigma$ is chosen similarly as that in the previous model.

$\bolds{\Theta}=\bolds{\Theta}_{0}+\sigma I$, where $\bolds{\Theta}_{0}$ is a zero-diagonal symmetric matrix. Each off-diagonal element ${\bolds{\Theta}_{0}}_{ij}$ independently follows a point mass $0.99\delta_{0}+0.01\delta_{0.2}$, and the diagonal element $\sigma$ is set to be the absolute value of the minimal negative eigenvalue of $\bolds{\Theta}_{0}$ to ensure the semi-positive-definiteness of $\bolds{\Theta}$.

In models 1b–4b we first generate $n$ independent data from $N_{p}(0,\bolds{\Sigma})$ and then transfer the normal data using transformation functions

where $f_{1}(x)=x$, $f_{2}(x)=\log(x)$, $f_{3}(x)=x^{\frac{1}{3}},$ $f_{4}(x)=\log(\frac{x}{1-x})$ and $f_{5}(x)=f_{2}(x)I_{\{x<-1\}}+f_{1}(x)I_{\{-1\leq x\leq 1\}}+(f_{4}(x-1)+1)I_{\{x>1\}}$. In all cases we let $n=300$ and $p=100$.

2 Applications to gene expression genomics.

We illustrate our proposed rank-based estimators on a real data set to recover the isoprenoid genetic regulatory network in Arabidposis thaliana [Wille et al. (2004)]. This dataset contains the gene expression measurements of 39 genes (excluding protein GGPPS7 in the MEP pathway) assayed on $n=118$ Affymetrix GeneChip microarrays.

We used seven estimators (GLASSO, MB, CLIME, LLW, R-GLASSO, R-NADS and R-ACLIME) to reconstruct the regulatory network. The first three estimators are performed after taking the log-transformation of the original data, and the other four estimators are directly applied to the original data. To be more conservative, we only considered the integration by union for the neighborhood selection procedures. We generated $100$ independent Bootstrap samples and computed the frequency of each edge being selected by each estimator. The final model by each method only includes edges selected by at least $80$ times over $100$ Bootstrap samples. We report the number of selected edges by each estimator in Table 7. The rank-based graphical lasso performs similarly to the LLW method. The rank-based adaptive CLIME produces the sparsest graphs. We also compared pairwise intersections of the selected edges among different estimators. More than $70\%$ of the selected edges by GLASSO, MB or CLIME turn out to be validated by both LLW and R-GLASSO, and more than $40\%$ of the selected edges by GLASSO, MB or CLIME are justified by R-NADS and R-ACLIME. The selected models support the biological arguments that the interactions between the pathways do exist although they operate independently under normal conditions [Laule et al. (2003), Rodríguez-Concepción et al. (2004)].

Discussion.

Appendix: Technical proofs

We now prove Lemma 1. First, Spearman’s rank correlation $\hat{r}_{ij}$ can be written in terms of the Hoeffding decomposition [Hoeffding (1948)]

where $d_{ij}=\frac{1}{n(n-1)}\sum_{k\neq l}\operatorname{sign}(x_{ki}-x_{li})\cdot\operatorname{sign}(x_{kj}-x_{lj}),$ and

Applying (20) and (21) yields $\hat{r}_{ij}-E(u_{ij})=u_{ij}-E(u_{ij})+\frac{3}{n+1}d_{ij}-\frac{3}{n+1}u_{ij}.$ Note $|u_{ij}|\leq 3$ and $|d_{ij}|\leq 1$. Hence, $|u_{ij}|\leq\frac{\varepsilon}{4\pi}(n+1)$ and $|d_{ij}|\leq\frac{\varepsilon}{4\pi}(n+1)$ always hold provided that $n>12\pi/\varepsilon$, which are satisfied by the assumption in Lemma 1. For such chosen $n$, we have

Finally, we observe that $u_{ij}$ is a function of independent samples $(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})$. Now we make a claim that if we replace the $t$th sample by some $\hat{\mathbf{x}}_{t}$, the change in $u_{ij}$ will be bounded as

Then we can apply the McDiarmid’s inequality [McDiarmid (1989)] to conclude the desired concentration bound for some absolute constant $c_{0}>0$,

where the third inequality holds if and only if $n_{i}=n_{j}=n-\frac{3}{2}.$

Proof of Theorem 2 For space of consideration, we only show the sketch of the proof, and the detailed proof is relegated to the supplementary file [Xue and Zou (2012)] and also the technical report version of this paper [Xue and Zou (2011b)]. We begin with an important observation that we only need to prove the risk bound for $\hat{\bolds{\Theta}}^{s}_{nd}$ because

where $C_{0}$ is some quantity depending on $b$, $B$ and $M$ only.

Now we can use (23) to further bound $|\hat{\theta}^{s.nd}_{kk}-\theta^{*}_{kk}|$ under the same event. To this end, we first derive an upper bound for $|(\hat{\theta}^{s.nd}_{kk})^{-1}-(\theta_{kk}^{*})^{-1}|$ as

Notice that $|\hat{\theta}^{s.nd}_{kk}-\theta^{*}_{kk}|=|(\hat{\theta}^{s.nd}_{kk})^{-1}-(\theta_{kk}^{*})^{-1}|\cdot|\hat{\theta}^{s.nd}_{kk}|\cdot|\theta^{*}_{kk}|$ and also $|\hat{\theta}^{s.nd}_{kk}|\leq|\hat{\theta}^{s.nd}_{kk}-\theta^{*}_{kk}|+|\theta^{*}_{kk}|$. Then $|\hat{\theta}^{s.nd}_{kk}-\theta^{*}_{kk}|$ can be upper bounded by

Since $d\lambda=o(1)$, we denote the right-hand side as $C_{1}d\lambda$ for some $C_{1}>0$.

Thus we can combine (24) and (Appendix: Technical proofs) to derive the desired upper bound under the same event. This completes the proof of Theorem 2.

Proof of Theorem 3 Throughout the proof, we consider the event

For ease of notation, define $\lambda_{d}=\lambda_{0}$ and $\lambda_{ad}=\lambda_{1}$. We focus on the proof of the sign consistency of $\hat{\bolds{\beta}}^{s.nad}_{k}$ in the sequel.

Under event (26), $\hat{\mathbf{R}}^{s}_{\mathcal{A}_{k}\mathcal{A}_{k}}$ is always positive-definite. To see this, the Weyl’s inequality yields $\lambda_{\min}(\hat{\mathbf{R}}^{s}_{\mathcal{A}_{k}\mathcal{A}_{k}})+\lambda_{\max}(\hat{\mathbf{R}}^{s}_{\mathcal{A}_{k}\mathcal{A}_{k}}-\bolds{\Sigma}^{*}_{\mathcal{A}_{k}\mathcal{A}_{k}})\geq\lambda_{\min}(\bolds{\Sigma}^{*}_{\mathcal{A}_{k}\mathcal{A}_{k}})$, and then we can bound the minimal eigenvalue of $\hat{\mathbf{R}}^{s}_{\mathcal{A}_{k}\mathcal{A}_{k}}$,

where $\circ$ denotes the Hadamard product. Due to the strong duality of linear programming [Boyd and Vandenberghe (2004)], the complementary slackness condition holds for the primal problem with respect to any primal and dual solution pair ($\bolds{\beta},\bolds{\alpha}^{+}_{k},\bolds{\alpha}^{-}_{k}$), which implies that $\alpha^{+}_{j}[(\hat{\mathbf{R}}^{s}_{(k)}\bolds{\beta}-\hat{\mathbf{r}}^{s}_{(k)})_{j}-\lambda_{1}w^{d}_{j}]=0$ and $\alpha^{-}_{j}[-(\hat{\mathbf{R}}^{s}_{(k)}\bolds{\beta}-\hat{\mathbf{r}}^{s}_{(k)})_{j}-\lambda_{1}w^{d}_{j}]=0$ for any $j\neq k$. Observe that only one of $\alpha^{+}_{j}$ and $\alpha^{-}_{j}$ can be zero since only one of $(\hat{\mathbf{R}}^{s}_{(k)}\bolds{\beta}-\hat{\mathbf{r}}^{s}_{(k)})_{j}=\lambda_{1}w^{d}_{j}$ and $(\hat{\mathbf{R}}^{s}_{(k)}\bolds{\beta}-\hat{\mathbf{r}}^{s}_{(k)})_{j}=-\lambda_{1}w^{d}_{j}$ can hold indeed, and thus we can uniquely define $\bolds{\alpha}_{k}=\bolds{\alpha}^{+}_{k}-\bolds{\alpha}^{-}_{k}$. Then we can rewrite the Lagrange dual function as

By the Lagrange duality, the dual problem of (15) is

where (27) and (29) are primal constraints, and (28) and (30) are dual constraints.

Note $K_{1}=(\hat{\mathbf{R}}^{s}_{{\mathcal{A}_{k}\mathcal{A}_{k}}})^{-1}\cdot(\hat{\mathbf{R}}^{s}_{{\mathcal{A}_{k}\mathcal{A}_{k}}}-\bolds{\Sigma}^{*}_{{\mathcal{A}_{k}\mathcal{A}_{k}}})\cdot(\bolds{\Sigma}^{*}_{{\mathcal{A}_{k}\mathcal{A}_{k}}})^{-1}$, and then we have

Some simple calculation shows $K_{1}\leq\frac{d\lambda_{1}G_{k}^{2}}{1-d\lambda_{1}G_{k}}.$ On the other hand,

Under probability event (26), we claim about $\mathbf{w}^{d}_{k}$ that

This claim is very useful to prove the other three optimality conditions (28), (29) and (30), and their proofs will be provided later.

Then we apply the triangle inequality to obtain an upper bound,

Next, we can easily obtain (30) via the triangular inequality

where the last inequality can be easily shown by combining (31) and (32).

Now it remains to prove (29). Using the facts that $\bolds{\theta}^{*}_{{\mathcal{A}_{k}^{c}}}=\mathbf{0}$ and $\bolds{\Sigma}^{*}\bolds{\Theta}^{*}=\mathbf{I}$, simple calculation yields that $\bolds{\Sigma}^{*}_{{\mathcal{A}_{k}^{c}\mathcal{A}_{k}}}(\bolds{\Sigma}^{*}_{{\mathcal{A}_{k}\mathcal{A}_{k}}})^{-1}\bolds{\sigma}^{*}_{\mathcal{A}_{k}}=\bolds{\sigma}^{*}_{\mathcal{A}_{k}^{c}}.$ Then we can rewrite the left-hand side of (29) as

Again we apply the triangle inequality to obtain an upper bound as follows:

where the last inequality is due to (31) and (32).

which immediately yields the desired lower bound by noting that

where both inequalities follow from the proper choices of tuning parameters $\lambda_{0}$ and $\lambda_{1}$ as stated in Theorem 3. On the other hand,

where the last inequality follows from the proper choice of $\lambda_{1}$ as stated in Theorem 3. Likewise we can prove the second claim (32) by noticing that

where we use facts that $\psi_{k}\geq 2C_{0}d\lambda_{0}$ and $\psi^{2}_{k}\geq 8G_{k}\lambda_{1}$. The two claims are proved, which completes the proof of Theorem 3.

Proof of Theorem 5 The techniques we use are similar to these for the proof of Theorem 3. The detailed proof of Theorem 5 is relegated to the supplementary material [Xue and Zou (2012)] and also the technical report version of this paper [Xue and Zou (2011b)] for the sake of space.

Acknowledgments.

We thank the Editor, the Associate Editor and three referees for their helpful comments.

Supplement material for “Regularized rank-based estimation of high-dimensional nonparanormal graphical models” \slink[doi,text=10.1214/12-AOS1041SUPP]10.1214/12-AOS1041SUPP \sdatatype.pdf \sfilenameaos1041_supp.pdf \sdescriptionIn this supplementary note, we give the complete proofs of Theorems 2 and 5.

References