Bootstrap inference for network construction with an application to a breast cancer microarray study\thanksrefT1

# Bootstrap inference for network construction with an application to a breast cancer microarray study\thanksrefT1

[ [    [ [    [ [    [ [ Fred Hutchinson Cancer Research Center, Fred Hutchinson Cancer Research Center, University of California, Davis and Fred Hutchinson Cancer Research Center S. Li
L. Hsu
P. Wang
Fred Hutchinson Cancer Research Center
M2-B500, 1100 Fairview Ave N.
Seattle, Washington 98109
USA
J. Peng
Department of Statistics
University of California, Davis
Mathematical Sciences Building
One Shields Avenue
Davis, California 95616
USA
\smonth11 \syear2011\smonth8 \syear2012
\smonth11 \syear2011\smonth8 \syear2012
\smonth11 \syear2011\smonth8 \syear2012
###### Abstract

Gaussian Graphical Models (GGMs) have been used to construct genetic regulatory networks where regularization techniques are widely used since the network inference usually falls into a high–dimension–low–sample–size scenario. Yet, finding the right amount of regularization can be challenging, especially in an unsupervised setting where traditional methods such as BIC or cross-validation often do not work well. In this paper, we propose a new method—Bootstrap Inference for Network COnstruction (BINCO)—to infer networks by directly controlling the false discovery rates (FDRs) of the selected edges. This method fits a mixture model for the distribution of edge selection frequencies to estimate the FDRs, where the selection frequencies are calculated via model aggregation. This method is applicable to a wide range of applications beyond network construction. When we applied our proposed method to building a gene regulatory network with microarray expression breast cancer data, we were able to identify high-confidence edges and well-connected hub genes that could potentially play important roles in understanding the underlying biological processes of breast cancer.

\kwd
\doi

10.1214/12-AOAS589 \volume7 \issue1 2013 \firstpage391 \lastpage417 \newproclaimremarkRemark

\runtitle

BINCO \thankstextT1Supported by NIH Grants R01GM082802 (SL, JP, PW), P01CA53996 (LH, PW), R01AG014358 (SL, LH), P50CA138293 (LH, PW), U24CA086368 (PW), NSF Grant DBI-08-20854 (JP) and DMS-10-07583 (JP).

{aug}

A]\fnmsShuang \snmLilabel=e1]Shuangli@fhcrc.org, A]\fnmsLi \snmHsulabel=e2]lih@fhcrc.org, B]\fnmsJie \snmPenglabel=e3]jiepeng@ucdavis.edu and A]\fnmsPei \snmWang\correflabel=e4]pwang@fhcrc.org

High dimensional data \kwdGGM \kwdmodel aggregation \kwdmixture model \kwdFDR

## 1 Introduction

The emergence of high-throughput technologies hasmade it feasible to measure molecular signatures of thousands of genes/proteins simultaneously. This provides scientists an opportunity to study the global genetic regulatory networks, shedding light on the functional interconnections among the regulatory genes, and leading to a better understanding of underlying biological processes. In this paper, we propose a network building procedure for learning genetic regulatory networks. Our work is motivated by an expression study of breast cancer (BC) that aims to infer the network structure based on 414 BC tumor samples [Loi et al. (2007)]. The proposed method enables us to detect high-confidence edges and well-connected hub genes that include both those previously implicated in BC and novel ones that may warrant further follow-up.

In practice, dependency structures of molecular activities such as correlation matrix and partial correlation matrix have been used to infer regulatory networks [Pollack et al. (2002), Kim et al. (2006), Varambally et al. (2005), Nie, Wu and Zhang (2006)]. Such dependency structures are often represented by graphical models in which nodes of a graph represent biological components such as genes or proteins, and the edges represent their interactions. These interactions may be indirect (e.g., two genes are co-regulated by a third gene) or direct (e.g., one gene is regulated by another gene). For the latter case, Gaussian Graphical Models (GGMs), which represent dependencies between pairs of nodes conditioning on the remaining of nodes, are often used.

For the data obtained from high-throughput technologies, the number of nodes is typically much larger than the number of samples, which is where the classical GGM theory [Whittaker (1990)] generally fails [Friedman (1989), Schaäfer and Strimmer (2005)]. This large--small- scenario is usually addressed by assuming that the conditional dependency structure is sparse [Dobra et al. (2004), Li and Gui (2006), Meinshausen and Buühlmann (2006), Yuan and Lin (2007), Friedman, Hastie and Tibshirani (2008), Rothman et al. (2008), Peng et al. (2009)]. However, like many high-dimensional regularization problems, finding the appropriate level of sparsity remains a challenge. This is particularly true for network structure learning, since the problem is unsupervised in nature. Traditional methods, such as Bayes information criteria [Schwarz (1978)] and cross-validation, aim to find a model that minimizes prediction error or maximizes a targeted likelihood function. They tend to include many irrelevant features [e.g., Efron (2004b), Efron et al. (2004), Meinshausen and Bühlmann (2006) and Peng et al. (2010)], and thus are not appropriate for learning the interaction structures.

Choosing the amount of regularization by directly controlling the false positive level would be ideal for structure learning. Recently, a few model aggregation methods have been proposed, and some of them provide certain control of false positives. For example, Bach (2008) proposed Bolasso, which chooses variables that are selected by all the lasso models [Tibshirani (1996)] built on bootstrapped data sets. In the context of network reconstruction, Peng et al. (2010) proposed choosing edges that are consistently selected across at least half of the cross-validation folds. More recently, Meinshausen and Bühlmann (2010) proposed the stability selection procedure to choose variables with selection frequencies exceeding a threshold. Under suitable conditions, they derived an upper bound for the expected number of false positives. In the same paper they also proposed the randomized lasso penalty, which aggregates models from perturbing the regularization parameters. Combined with stability selection, randomized lasso achieves model selection consistency without requiring the irrepresentable condition [Zhao and Yu (2006)] that is necessary for lasso to achieve model selection consistency. In another work, Wang et al. (2011) proposed a modified lasso regression—random lasso—by aggregating models based on bootstrap samples and random subsets of variables. All these works have greatly advanced research in model selection in the high-dimensional regime. However, none of these methods provide direct estimation and control of the false discovery rate (FDR).

In this paper, we address the problem of finding the right amount of regularization in the context of high-dimension GGMs learning. In a spirit similar to the aforementioned methods, we first obtain selection frequencies from a collection of models built by perturbing both the data and the regularization parameters. We then model these selection frequencies by a mixture distribution to yield an estimate of FDR on the selected edges, which is then used to determine the cut-off threshold for the selection frequencies. This framework is rather general, as it only depends on the empirical distribution of the selection frequencies. Thus, it can be applied to a wide range of problems beyond GGMs.

The rest of this paper is organized as follows. In Section 2 we describe in detail the proposed method. In Section 3 an extensive simulation study is conducted to compare the method with the stability selection procedure and then evaluate its performance under different settings. In Section 4 the method is illustrated by building a genetic interaction network based on microarray expression data from a BC study. The paper is concluded with some discussion in Section 5.

## 2 Method

### 2.1 Gaussian graphical models

In a Gaussian Graphical Model (GGM) network construction is defined by the conditional dependence relationships among the random variables. Let denote a -dimension random vector following a multivariate normal distribution , where is a positive definite matrix. The conditional dependence structure among is represented by an undirected graph with vertices representing and the edge set defined as

 E={(i,j)\dvtxYi and Yj are dependent given Y−{i,j},1≤i,j≤p},

where . The goal of network construction is to identify the edge set . Under the normality assumption, the conditional independence between and is equivalent to the partial correlation between and given being zero. It is also equivalent to the entry of the concentration matrix () being zero, that is, [Dempster (1972), Cox and Wermuth (1996)], since .

There are two main types of approaches to fitting a GGM. One is the maximum-likelihood-based approach, which estimates the concentration matrix directly. The other is the regression-based approach, which fits the GGM through identifying nonzero regression coefficients of the following regression:

 Yi=∑j≠iβijYj+εi,1≤i≤p,

where is uncorrelated with . The nonzero ’s in the above regression setting correspond to nonzero entries in the concentration matrix since it can be shown that . In both approaches, there are parameters to estimate, which requires proper regularization on the model if is larger than the sample size . This can be achieved by making a sparsity assumption on the network structure, that is, assuming that most pairs of variables are conditionally independent given all other variables. Such an assumption is reasonable for many real life networks, including genetic regulatory networks [Gardner et al. (2003), Jeong et al. (2011), Tegner et al. (2003)]. Methods have been developed along these lines by using regularization. For example, Yuan and Lin (2007) proposed a sparse estimator of the concentration matrix via maximizing the penalized log-likelihood. Efficient algorithms were subsequently developed to fit this model with high-dimensional data [Friedman, Hastie and Tibshirani (2008), Rothman et al. (2008)]. For regression-based approaches, Meinshausen and Bühlmann (2006) considered the neighborhood selection estimator by minimizing individual loss functions

 L(i)(β,Y)=12∥∥∥Yi−∑j:j≠iβijYj∥∥∥2+λ∑j:j≠i|βij|,i=1,…,p, (1)

while Peng et al. (2009) proposed the space algorithm by minimizing the joint loss

 L(Y,θ)=12{p∑i=1∥∥∥Yi−∑j:j≠i√σjjσiiρijYj∥∥∥2}+λ∑1≤i

From objective functions (1) and (2), it is clear that the selected edge set depends on the regularization parameter . Since the goal here is to recover the true edge set, ideally should be determined based on considerations such as FDR and power with respect to edge selection. Moreover, when the sample size is limited, a model-aggregation-based strategy can improve the selection result compared to simply tuning the regularization parameter. Thus, in the following section, we introduce a new model-aggregation-based procedure that selects edges based on directly controlling the FDRs.

Throughout the rest of this paper, we refer to the set of all pairs of variables as the candidate edge set (denoted by ), the subset of those edges in the true model as the true edge set (denoted by ) and the rest as the null edge set (denoted by ). We denote the size of a set of edges by . Note that and the total number of edges in is .

### 2.2 Model aggregation

Consider a good network construction procedure, where good is in the sense that the true edges are stochastically more likely to be selected than the null edges. Then it would be reasonable to choose edges with high selection probabilities. In practice, these selection probabilities can be estimated by the selection frequencies over networks constructed based on perturbed data sets. In the following, we formalize this idea.

Let be an edge selection procedure with a regularization parameter and be the set of selected edges by applying to data . The selection probability of edge is defined as

 pij=E(I{(i,j)∈Sλ(Y)}),

where is the indicator function. Let be the space of resamples from (e.g., through bootstrapping or subsampling). For a random resample from , we define

In many cases (see Section C in the supplemental article [Li et al. (2013)]), ’s and ’s are close. For these cases, we can estimate by the selection frequency , which is the proportion of resamples in which the edge is selected:

 Xλij ≡ Xij(A(λ);Y1,…,YB∈R(Y)) = 1BB∑k=1I{(i,j)∈Sλ(Yk)},1≤i

The aggregation-based procedures for choosing edges of large selection frequencies can be represented as

 Sλc={(i,j)\dvtxXλij≥c}%forc∈(0,1].

is reasonable as long as most true edges have selection frequencies greater than or equal to and most null edges have selection frequencies less than . Ideally, we want to find a threshold satisfying

 (4)

so that the corresponding procedure is consistent, that is, . In fact, if is selection consistent and , then

 Pr({⋂(i,j)∈E{Xλij=1}}∩{⋂(i,j)∈Ec{Xλij=0}})→1as n→∞, (5)

and thus any satisfies (4). Note that (4) is in general a much weaker condition than (5), which suggests that we might find a consistent even when is not consistent.

For the finite data case, an aggregation-based procedure could also perform better than the original procedure, as illustrated by the following simulation example (the simulation setup is provided in Section 3). Figure 1(a) shows the empirical distribution of selection frequencies based on a simulated data set and Figure 1(b) shows the empirical distributions of true edges (green triangles) and null edges (red crosses). Note that most null edges have low selection frequencies , while most true edges have large selection frequencies . This suggests that with a properly chosen (say, ), will select mostly true edges and only a small number of null edges. In fact, by simply choosing the cutoff , outperforms in both FDR and power (Figure 2).

### 2.3 Modeling selection frequency

Now we introduce a mixture model, similar in spirit to Efron (2004a), for estimating the FDR of an aggregation-based procedure . We will use this estimate to choose the optimal and by controlling FDR while maximizing power. Assume that the selection frequencies , generated from resamples, fall into two categories, “true” or “null,” depending on whether is a true edge or a null edge. Let be the proportion of the true edges. We also assume that has density or if it belongs to the “true” or the “null” categories, respectively. Note that both and depend on the sample size , but such dependence is not explicitly expressed in order to keep the notation simple. The mixture density for can be written as

 fλ(x)=(1−π)fλ0(x)+πfλ1(x),x∈{0,1/B,2/B,…,1}. (6)

Based on this mixture model, the (positive) FDR [Storey (2003)] of the aggregation-based procedure is

 FDR(Sλc)=Pr((i,j)∈Ec|(i,j)∈Sλc)=∑x≥c(1−π)fλ0(x)∑x≥cfλ(x). (7)

Given an estimate (which will be discussed below) from (7), the number of true edges in can be estimated by

 ˆNE(Sλc)=∣∣Sλc∣∣(1−ˆFDR(Sλc)), (8)

which can be used to compare the power of across various choices of and , as the total number of true edges is a constant. Consequently, for a given targeted FDR level , we first seek for the optimal threshold for each ,

 c∗(λ)=min{c\dvtxˆFDR(Sλc)≤α}, (9)

and then we find the optimal regularization parameter

 λ∗=argmaxλ∈ΛˆNE(Sλc∗(λ)), (10)

such that the corresponding procedure achieves the largest power among all competitors with estimated FDR not exceeding .

The above procedure depends on a good FDR estimate, which in turn requires good estimates of the mixture density and its null-edge contribution . A natural estimator of is simply the empirical selection frequencies, that is,

 ^fλ(kB)=nλkNΩ,k=0,1,…,B,

where is the total number of candidate edges and is the number of edges with selection frequencies equal to .

Before describing an approach to estimating and , we note two observations from Figure 1(b). First, the contribution from the true edges to the mixture density is small in the range where the selection frequencies are small. Second, the empirical distribution of is monotonically decreasing. These can be formally summarized as the following condition.

{pc*}

There exist and , , such that as : {longlist}[(C1)]

on ;

is monotonically decreasing on .

This proper condition is satisfied by a class of procedures as described in the lemma below (the proof is provided in the Appendix).

###### Lemma 1

A selection procedure satisfies the proper condition if, as the sample size increases, tends to one uniformly for all true edges and has a limit superior strictly less than one for all null edges.

{remark}

It is easy to verify that all consistent procedures applied to subsampling resamples satisfy the condition in Lemma 1. Other examples are procedures that use randomized lasso penalties [Meinshausen and Bühlmann (2010)]. See Section 2.5 for more details.

The proper condition motivates us to estimate and by fitting a parametric model for in the region and then extrapolating the fit to the region . This is because if C1 is satisfied, then can be well approximated based on the empirical mixture density from the region . If C2 is also satisfied, the extrapolation of will be a good approximation to on for a reasonably chosen family of .

We choose the parametric family as follows. Given , it is natural to model the selection frequency by a (rescaled) binomial distribution, denoted by , due to the independent and identical nature of resampling conditional on the original data. Moreover, we use a powered beta distribution [i.e., the distribution of where beta, ] as the prior for ’s, denoted by with . This is motivated by the fact that the beta family is a commonly used conjugate prior for the binomial family, and the additional power parameter  simply provides more flexibility in fitting. Thus, the distribution of selection frequencies of null edges is modeled as

 hθ(x)=∫10b1(x|τ)b2(τ|θ)dτ.

The null-edge contribution can be estimated by fitting to the empirical mixture density in the fitting range , which, in practice, is determined based on the shape of (details are given in Section 2.4). Specifically, we estimate and by and , via

 (^π,^θ)=argminπ,θL(^fλ(⋅),(1−π)hθ(⋅)), (11)

where , which amounts to the Kullback–Leibler distance.

### 2.4 Proper regularization range

Following what we propose in Section 2.3, we can evaluate the aggregation-based procedure for different choices of with regard to model–selection–based criteria: the FDR and the number of selected true edges. For the range of , we consider those that yield “U-shaped” empirical distributions of selection frequencies, that is, decreases in the small-selection-frequency range and then increases in the large-selection-frequency range [see Figure 1(a) and Figure 3 for examples of “U-shaped” distribution]. The decreasing trend is needed for the proper condition to hold, while the increasing trend helps to control the FDR, since an with implies, by (7), that

 ∑x≥cfλ(x)≥(1−π)∑x≥cfλ0(x)α. (12)
U-shape detection procedure
1. INPUT , the empirical density of selection frequencies. Set (the U-shape indicator).
2. Check U-shape.
2.1. Check valley point.
2.1.1. Calculate , the valley point position,
where is a smooth curve fitted based on . (We use
the R-function , where the degree of
freedom parameter is determined such that the derivative
of has only one sign change.)
2.1.2. IF
Set , GOTO Step 3.
END IF
2.2. Calculate , the peak before .
2.3. Check if is “roughly” decreasing on .
2.3.1. Calculate , and
.
2.3.2. IF
Set , GOTO Step 3.
END IF
2.4. Check if is “roughly” increasing on .
2.4.1. Calculate , and
.
2.4.2. IF
Set , GOTO Step 3.
END IF
3. RETURN .

Therefore, if is not sufficiently large at the tail, may not be achieved for a small value of . The increasing trend also helps to obtain decent power since it guarantees a substantial size of . Based on our experience, the values chosen based on (9) and (10) indeed always corresponds to a “U-shaped” empirical selection frequency distribution.

Thus, we propose the following simple procedure for identifying “U-shaped” ’s to determine the proper regularization range in practice. An illustration for this procedure is given in Figure 3.

{remark}

Step 2.1 is based on our extensive simulation where we find that a large value of often corresponds to a too-small , yielding too many null edges with high selection frequencies, which makes (12) difficult to hold for reasonably small FDR levels (see Section D1 in the supplemental article [Li et al. (2013)]). If is not recognized as “U-shaped” for a large range of ’s, we would consider the data as lack of signals where a powerful is not attainable. One example is the empty network (see Section 3.2 and Figure S-1 in the supplemental article [Li et al. (2013)]).

Sections 2.22.4 provide a procedure for network inference based on directly estimating FDR. We name the procedure as BINCO—Bootstrap Inference for Network COnstruction, as we suggest to use bootstrap resamples. The main steps are summarized below.

### 2.5 Randomized lasso

For an regularized procedure , the proper condition (Section 2.3) is satisfied if is selection consistent, which usually requires strong conditions, for instance, the well-known irrepresentable condition under the lasso regression setting [Zhao and Yu (2006), Zou (2006), Yuan and Lin (2007), Wainwright (2009)] or the so-called neighborhood stability condition under the GGM setting [Meinshausen and Bühlmann (2006), Peng et al. (2009)]. Recently, Meinshausen and Bühlmann (2010) proposed the randomized lasso, which is a procedure based on randomly sampled regularization parameters. For example, the randomized lasso version of space would be

 L(Y,θ,W)=12{p∑i=1∥∥∥Yi−∑j:j≠i√σjjσiiρijYj∥∥∥2}+λ∑1≤i
BINCO procedure
1. INPUT the initial range of regularization parameter values; the dataset; and the desired FDR level.
2. FOR TO
2.1.
2.2. Generate the empirical density of selection frequencies.
2.3. Check whether is U-shaped based on the output from
the “U-Shape Detection Procedure.”
2.4. IF is U-shaped (i.e., )
2.4.1. Obtain the null density estimate by (11).
2.4.2. Find the optimal threshold by (9), where the FDR
is estimated based on (7) with and replaced by
and ,respectively.
2.4.3. Obtain and calculate , the estimated
number of true edges being selected, based on (8).
END IF
ELSE , .
2.5. OUTPUT and .
NEXT
3. Determine the optimal regularization through (10). The optimal selection is .

where ’s are randomly sampled from a probability distribution supported on for some (note that corresponds to the ordinary penalty). The advantage of this randomized lasso procedure is that, by perturbing the regularization parameters, the irrelevant features may be decorrelated from the true features in some configurations of randomly sampled weights such that the irrepresentable condition is satisfied. Therefore, it selects all true features with probability tending to 1 and any irrelevant feature with a limiting probability strictly less than 1. As a result, a consistent aggregation-based procedure can be achieved under conditions “typically much weaker than the standard assumption of the irrepresentable condition” [Meinshausen and Bühlmann (2010), Theorem 2]. For this case, based on Lemma 1, the proper condition is also satisfied.

If (13) is used as the original (nonaggregated) procedure, an additional parameter , which controls the amount of perturbation of the regularization parameter, needs to be chosen. A small guards better against false positives but damages power, while a large may result in a liberal procedure. Here we provide a two-step data-driven procedure for choosing an appropriate in BINCO. We first fix , that is, the ordinary penalty, to find a proper range for that corresponds to the “U-shaped” empirical mixtures. Then for each , we consider a set of pairs such that , that is, keeping the average amount of regularization unchanged. For example, in the simulation study, we use . We then pick the pair such that is the smallest among those ’s that yield U-shaped empirical mixture distributions. Our simulation shows that such a choice of ensures good power for BINCO while controlling FDR in a slightly conservative fashion.

## 3 Simulation

In this section we first compare the performance of BINCO with stability selection [Meinshausen and Bühlmann (2010)], and then investigate the performance of BINCO with respect to various factors, including the network structure, dimensionality, signal strength and sample size.

We use space [Peng et al. (2009)] coupled with randomized lasso (13) as the original nonaggregate procedure, where the random weights ’s are generated from the uniform distribution for . The selection frequencies are obtained based on resamples. Since subsampling of size is proposed for stability selection, we use subsampling to generate resamples when comparing BINCO and stability selection. For investigating BINCO’s performance, we use bootstrap resamples because it yields slightly better performance (see Remark 3.2).

The performance of both methods are evaluated by true FDRs and power, since for simulations we know whether an edge is true or null. In addition, we define ideal power, which is the best power one can achieve for given the true (in simulation we consider and ). Based on ideal power, we can evaluate the efficiency of the methods under different settings. For each simulation setting, results are based on 20 independent simulation runs.

### 3.1 Comparison between BINCO and stability selection

Stability selection procedure selects , a set of edges with the maximum selection frequency over a prespecified regularization set exceeding a threshold . Assuming an exchangeability condition upon the irrelevant variables (here the null edges), Meinshausen and Bühlmann [(2010), Theorem 1] derived an upper bound for the expected number of falsely selected variables for each choice of . Specifically, under suitable conditions, the expected number of null edges selected by the set , denoted by , satisfies

 E(V)≤q2Λ(2t−1)NΩ, (14)

where is the total number of candidate edges and is the expected number of edges selected under at least one . In practice, can be estimated by . Dividing both sides of (14) by , we obtain

 E(V)|SΛstable(t)|≤q2Λ(2t−1)NΩ⋅|SΛstable(t)|. (15)

Although stability selection is intended to control , for an easier comparison with BINCO, we use to approximate FDR and obtain the optimal by finding the smallest threshold such that the upper bound on the right-hand side of (15) is less than or equal to .

For data generation, we first consider a power-law network with nodes whose degree (i.e., the number of connected edges for each node) distribution follows . The scaling exponent is set to be 2.3, which is consistent with the findings in the literature for biological networks [Newman (2003)]. There are in total 495 true edges in this network and its topology is illustrated in Figure 5(a). The sample size is . Two settings with different signal strengths are considered: (1) strong signal, the mean and standard deviation (SD) of nonzero ’s are 0.34 and 0.13, respectively; (2) weak signal, the mean and SD of nonzero ’s are 0.25 and 0.09, respectively. Note both positive and negative correlations are allowed in this network.

We compare the performance of BINCO and stability selection at a targeted FDR level of 0.05. For BINCO, we consider as the initial range for and then obtain the optimal final selection following the steps at the end of Section 2.4. For stability selection, since no specific guidance was provided for choosing and (the randomized lasso regularization perturbation parameter), we consider three different values for and a collection of intervals with varying from 40 to 100 and . This choice of is due to the fact that the upper bound in (15) cannot be controlled at 0.05 for any for , and the performance of stability selection is largely invariant for .

When the signals are strong, BINCO gives a conservative but still maintains good power [Figures 4(a) and 4(c)]. The performance of stability selection varies for different choices of and . The FDRs are larger than the targeted level 0.05 for some ’s when and for all ’s when . For other cases (some ’s when and all ’s when ), the FDR control is very conservative and the corresponding power is consistently lower than BINCO. When the signals are weak, stability selection is much more conservative than BINCO and results in much lower power [Figures 4(b) and 4(d)]. In Table 1 we report the ideal power, the power for BINCO and the best power for stability selection (among different choices of ) under , 0.8 and 1. We also calculate the power efficiency as the ratio of the power for the method over the ideal power, for BINCO and stability selection, respectively. It can be seen that the power of BINCO is close to the ideal power for both levels of signal strength, while stability selection is too conservative when the signal strength is weak. For more detailed results, see Section A1 in the supplemental article [Li et al. (2013)].

{remark}

In some cases we find that stability selection fails to control FDR. We suspect this may be due to the violation of the exchangeability assumption in Theorem 1 of Meinshausen and Bühlmann (2010). We examine the impact of the exchangeability assumption by simulation and find that when it is violated, the theoretical upper bound in (14) for may not hold (see Section D2 in the supplemental article [Li et al. (2013)] for further details).

### 3.2 Further investigation of BINCO

Now we investigate the effects of the network structure, dimensionality, signal strength and sample size on the performance of BINCO.

Network structure.

We consider four different network topologies: empty network, power-law network, empirical network and hub network. In each network there are five disconnected components with 100 nodes each. Below is a brief description of the network topologies:

{longlist}

[(1)]

Empty network: there is no edge connecting any pair of nodes.

Power-law network: the degree follows a power-law distribution with parameter as described in Section 3.1 [Figure 5(a)].

Empirical network: the topology is simulated according to an empirical degree distribution of one genetic regulatory network [Schadt et al. (2005)] [Figure 5(b)].

Hub network: three nodes per component have a large number of connecting edges (15) and all other nodes have a small number of connecting edges (5) [Figure 5(c)].

We set the sample size . The signal strength for all networks except for the empty network is fixed at the strong level as in Section 3.1.

For the empty network, the empirical mixture distributions of selection frequencies monotonically decrease on a wide range of (Figure S-1) and are not recognized by BINCO as “U-shaped.” Thus, we reach the correct conclusion that there is no signal in this case. In contrast, data sets from the other three networks produce the desired “U-shaped” mixture distributions for some (Figure S-2).

We compare BINCO results across networks 2-4 with FDR targeted at level and 0.1. BINCO gives slightly conservative control on FDR and achieves reasonable power for all three networks (Table 2). The comparison to the ideal power shows that the network topologies investigated here have only a small effect on BINCO’s efficiency (Table 3).

Dimensionality. We investigate the impact of dimensionality on the performance of BINCO. We consider the power-law network and let the number of nodes vary from 500, 800 to 1000. To keep the complexity of each component the same across different choices of , we set the component size constant, being 100, and the number of components . Again the sample size is used for all three cases and the signal strength is fixed at the strong level as in Section 3.1.

For all three choices of , BINCO performs similarly (Table 4), with slightly conservative FDR and power around 0.8. The dimensionality does not demonstrate a significant impact on BINCO. BINCO’s result is also largely invariant when we compare networks of differing numbers of components with fixed (such that component size varies, see Section A3 in the supplemental article [Li et al. (2013)]).

Signal strength. We consider three levels of signal strength: strong, weak and very weak. The corresponding means and SDs of nonzero ’s are , and , respectively. The network is the power-law network with and sample size is for all settings.

BINCO provides good control on FDR, however, the power decreases from 0.8 to 0.3 as the signal weakens (Table 5). Comparing the power of BINCO with the ideal power (Table 6), we see that BINCO remains efficient and the loss in power is largely due to reduction of signal strength.

Sample size. Now we consider the impact of sample size by varying it from 200, 500 and 1000, while keeping the signal strength at the “very weak” level as in the previous simulation. The network structure is again the power-law network with .

With an increased sample size, the power of BINCO is significantly improved from 0.3 to nearly 0.9 while the FDRs are well controlled (Table 7).

In summary, BINCO has good control for FDR under a wide range of scenarios. Its performance is shown to be robust for networks with different topologies and dimensionalities, and its efficiency is not influenced much even when the signal strength is weak. As the sample size increases, the power of BINCO is improved significantly.

{remark}

We propose to use bootstrap over subsampling, as the former appears to give slightly better power. Intuitively, bootstrap contains more distinct samples [0.632, Pathak (1962)] than subsampling (0.5). However, the difference we have observed is rather small. For example, we compare the power over 20 independent samples between bootstrap and subsampling under the power-law network setting. For , the power is 0.810 for bootstrap and 0.801 for subsampling (compare Tables 3 and 1); while for , the power is 0.845 for bootstrap and 0.835 for subsampling [compare Tables 3 and S-7 from Li et al. (2013)]. This observation is in agreement with the conclusions of several others [Menshausen and Bühlmann (2010), Freedman (1977), Bühlmann and Yu (2002)].

## 4 A real data application

We apply the BINCO method to a microarray expression data set of breast cancer (BC) [Loi et al. (2007)] to build a gene expression network related to the disease. The data (http://www.ncbi.nlm. nih.gov/geo/query/acc.cgi?acc=GSE6532) contains measurements of expression levels of 44,928 probes in tumor tissue samples from 414 BC patients based on the Affymetrix Human Genome U133A, U133B and U133 plus 2.0 Microarray platforms.

We preprocess the data as follows. First, a global normalization is applied by centering the median of each array to zero and scaling the Median Absolute Deviation (MAD) to one. Probes with standard deviation (SD) greater than the -trimmed mean of all SDs are selected. We further focus on a subset of probes for genes from cell cycle and DNA-repair related pathways (http://peiwang.fhcrc.org/internal/papers/DNArepair_CellCircle_related.csv/view), as these pathways have been shown to play significant roles in BC tumor initiation and development. Clinical information including age, tumor size, ER-status (positive or negative) and treatment status (tamoxifen treated or not) is incorporated in the analysis as “fake genes” since we are also interested in investigating whether gene expressions are associated with these clinical characteristics. Finally, we standardize each expression level to have mean zero and SD one. The resulting data set has genes/probes (including four clinical variables) and tumor samples.

We generate selection frequencies by applying the space algorithm with randomized lasso regularization to bootstrap resamples. The initial range of the tuning parameter is set to be . We then apply the BINCO procedure and find that the optimal values for the regularization parameters are and . The empirical distribution of selection frequencies of all edges and the null density estimation are given in Figure 6. When the estimated FDR is controlled at 0.05, 0.1 and 0.2, BINCO identifies 125, 222 and 338 edges, respectively. The estimated network for is shown in Figure 7. In this figure, two components of a large connectivity structure are observed. They contain most of the genes that are connected by a large number of high-selection-frequency edges. This constructed network can help to generate a useful biological hypothesis and to design follow-up experiments to better understand the underlying mechanism in BC. For example, BINCO suggests with high confidence for the association between MAP3K4 and STAT3. MAP3K4 plays a role in the signal transduction pathways of BC cell proliferation, survival and apoptosis [Bild and Johnson (2001)], and the constitutive activation of STAT3 is also frequently detected in BC tissues and cell lines [Hsieh, Cheng and Lin (2005)]. Interestingly, both MAP3K4 and STAT3 play roles in the regulation of c-Jun, a novel candidate oncogene whose aberrant expression contributes to the progression of breast and other human cancers [Tront, Hoffman and Liebermann (2006); Shackleford et al. (2011)]. The association between MAP3K4 and STAT3 detected by BINCO suggests their potential cooperative roles in BC. It is also worth noting that for the four clinical variables, the only edge with high selection frequency is the one between age and ER-status (selection frequency ). All edges between clinical variables and the genes/probes are insignificant (selection frequencies ).

Networks built on perturbed data sets can also be used to detect hub genes (i.e., highly connected genes), which are often of great interest due to the central role these genes may play in genetic regulatory networks. The idea is to look for genes that show consistent high connection in estimated networks across perturbed data sets. Here, we propose to detect hub genes by the ranks of their degrees based on the estimated networks using and . The ten genes with the largest means and the smallest SDs of the degree rank across 100 bootstrap resamples (see Figure S-3, in black dots) are MBD4, TARDBP, DDB2, MAP3K4, ORC3L, CDKN1B, REL, ATR, LGMN and CDKN3. Nine out of these ten genes have been reported relevant to BC, while the remaining one (TARDBP) is newly discovered to be related to cancer [Postel-Vinay et al. (2012)], although its role in BC is not clear at present. The neighborhood topologies of these hub genes in the network estimated by BINCO are illustrated in Figure 7. More details of these hub genes are given in the supplemental article [Li et al. (2013)], Section B.

## 5 Discussion

In this paper we propose the BINCO procedure to conduct high-dimensional network inference. BINCO employs model aggregation strategies and selects edges by directly controlling the FDR. This is achieved by modeling the selection frequencies of edges with a two-component mixture model, where a flexible parametric distribution is used to model the density for the null edges. By doing this, BINCO is able to provide a good estimate of FDR and hence properly controls the FDR. To ensure BINCO works, we propose a set of screening rules to identify the U-shape characteristic of empirical selection frequency distributions. Based on our experience, a U-shape corresponds to a proper amount of regularization such that the FDR is well controlled and the power is reasonable. Extensive simulation results show that BINCO performs well under a wide range of scenarios, indicating that it can be used as a practical tool for network inference. Although we focus on the GGM construction problem in this paper, BINCO is applicable to a wide range of problems where model selection is needed because it provides a general approach to modeling the selection frequencies.

We use a mixture distribution with two components, one corresponding to true edges and the other corresponding to null edges, to model the selection frequency distribution. This two-component mixture model is adequate as long as the distribution of the null component is identifiable and can be reasonably estimated, as formalized in the proper condition. Note that the proper condition holds for a wide range of commonly used (nonaggregation) selection procedures (Lemma 1, Remark 2.3). To further ensure the FDR can be controlled at a reasonable level, we propose a U-shape detection procedure and only apply BINCO if the empirical distribution of selection frequencies passes the detection. These rules for U-shape detection are empirical but appear to work very well based on our extensive simulations.

BINCO works well despite the presence of correlations between edges (see Section D1 in the supplemental article [Li et al. (2013)]), because we use the independence of edges only as a working assumption . It is well known that if the marginal distribution is correctly specified, the parameter estimates are consistent even in the presence of correlation. This is similar to the generalized estimating equations, where if the mean function is correctly specified, the parameters will be consistently estimated [Liang and Zeger (1986)]. Toward this end, we use the three-parameter power beta distribution to allow for adequate flexibility in modeling the marginal distribution of selection frequencies.

BINCO is computationally feasible for high-dimensional data. The major computational cost lies in generating the selection frequencies via resampling. For each resample, the computational cost is determined by that of the nonaggregated procedure BINCO coupled with. In terms of space, it is . The processing time for a data set with , , under a given and 100 bootstrap samples to generate selection frequencies is about 20 minutes on a PC with Pentium dual-core CPU at 2.8 GHz and 1 G ram. These selection frequencies can be simultaneously generated through parallel computing for different ’s and weights. Fitting the mixture model takes much less time, which is about 2 minutes for the above example on the same computer.

Although we use GGM as our motivating example, BINCO works well even if the multivariate normality assumption does not hold. Note that the multivariate normality assumption only concerns the interpretation of the edges. Under GGM, the presence of an edge means conditional dependency of the corresponding nodes given all other nodes. Without the normality assumption, one can only conclude nonzero partial correlation between the two nodes given the rest of the nodes. The space method used in this paper is to estimate the concentration network (where an edge is drawn between two nodes if the corresponding partial correlation is nonzero) and has been shown to work well under nonnormal cases such as multivariate- distributions [Peng et al. (2009)]. We also generate data from nonnormal distributions and found that BINCO works well in this situation (see Section D4 in the supplemental article [Li et al. (2013)]).

BINCO is an aggregation-based procedure. In principle, it can be coupled with any selection procedure. In this sense, it has a wide range of applications as long as the features are defined (e.g., edges as in this paper, variables or canonical correlations as in the example below) and the selection procedure is reasonably good, for example, producing probabilities that satisfy the condition in Lemma 1. One application beyond GGM could be on the multi-attribute network construction where the links/edges are defined based on canonical correlations [Waaijenborg, Verselewel de Witt Hamer and Zwinderman (2008), Katenka and Kolaczyk (2012), Witten, Tibshirani and Hastie (2009)]. Another interesting extension may be on the time-varying network construction [Kolar et al. (2010)] where appropriate incorporation of the time-domain structure across aggregated models will be important. These are beyond the scope of this paper and will be pursued in future research.

The R package BINCO is available through CRAN.

## Appendix: Proof of Lemma 1

{pf}

Suppose as the sample size increases, an edge selection procedure gives selection probabilities (with respect to resample space) which uniformly satisfy

 ~p(n)ij→1if (i,j)∈E (.1)

and

 limsup~p(n)ij≤M<1if (i,j)∈Ec. (.2)

Suppose is large such that . Let be a random sample from the set of selection frequencies generated by applying on resamples, that is, , . Also suppose has density if . Then the mixture model (6) becomes

 fλ(x) = (1−π)fλ0(x)+πfλ1(x) = ∑(i,j)∈Ec1NΩfλij(x)+∑(i,j)∈E1NΩfλij(x)

with and .

Because of the i.i.d. nature of resamples given the data, is a binomial density with as the probability of success, that is, for . This binomial density is monotone decreasing for greater than its mode or . By (.2), given and such that , such that for all and hence for any null edge , is monotone decreasing on , which implies C2 since . Also, (.1) implies, for , uniformly for , which implies C1 for any . Taking such that satisfies the proper condition and completes the proof.

## Acknowledgments

We thank anonymous reviewers and editors for helpful comments that significantly improved this paper. We also thank Ms. Noelle Noble for technical editing.

{supplement}

[id=suppA] \stitleSupplement to “Bootstrap inference for network construction with an application to a breast cancer microarray study”
\slink[doi]10.1214/12-AOAS589SUPP \sdatatype.pdf \sfilenameaoas589_supp.pdf \sdescriptionThis supplement contains additional simulation results, details of the hub genes detected by BINCO on the breast cancer data, and examples of and being close.

## References

• Bach (2008) {binproceedings}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmBach, \bfnmF.\binitsF. (\byear2008). \btitleBolasso: Model consistent lasso estimation through the bootstrap. In \bbooktitleProceedings of the 25th International Conference on Machine \bpages33–40. \bpublisherACM, \blocationNew York. \bptokimsref \endbibitem
• Bild and Johnson (2001) {bmisc}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmBild, \bfnmA.\binitsA. \AND\bauthor\bsnmJohnson, \bfnmG.\binitsG. (\byear2001). \bhowpublishedSignaling by erbB receptors in breast cancer: Regulation by compartmentalization of heterodimetric receptor complexes. Annual summary report. Available at http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA400019. \bptokimsref \endbibitem
• Bühlmann and Yu (2002) {barticle}[mr] \bauthor\bsnmBühlmann, \bfnmPeter\binitsP. \AND\bauthor\bsnmYu, \bfnmBin\binitsB. (\byear2002). \btitleAnalyzing bagging. \bjournalAnn. Statist. \bvolume30 \bpages927–961. \biddoi=10.1214/aos/1031689014, issn=0090-5364, mr=1926165 \bptokimsref \endbibitem
• Cox and Wermuth (1996) {bbook}[mr] \bauthor\bsnmCox, \bfnmD. R.\binitsD. R. \AND\bauthor\bsnmWermuth, \bfnmNanny\binitsN. (\byear1996). \btitleMultivariate Dependencies: Models, Analysis and Interpretation. \bseriesMonographs on Statistics and Applied Probability \bvolume67. \bpublisherChapman & Hall, \blocationLondon. \bidmr=1456990 \bptokimsref \endbibitem
• Dempster (1972) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmDempster, \bfnmA. P.\binitsA. P. (\byear1972). \btitleCovariance selection. \bjournalBiometrika \bvolume32 \bpages95–108. \bptokimsref \endbibitem
• Dobra et al. (2004) {barticle}[mr] \bauthor\bsnmDobra, \bfnmAdrian\binitsA., \bauthor\bsnmHans, \bfnmChris\binitsC., \bauthor\bsnmJones, \bfnmBeatrix\binitsB., \bauthor\bsnmNevins, \bfnmJoseph R.\binitsJ. R., \bauthor\bsnmYao, \bfnmGuang\binitsG. \AND\bauthor\bsnmWest, \bfnmMike\binitsM. (\byear2004). \btitleSparse graphical models for exploring gene expression data. \bjournalJ. Multivariate Anal. \bvolume90 \bpages196–212. \biddoi=10.1016/j.jmva.2004.02.009, issn=0047-259X, mr=2064941 \bptokimsref \endbibitem
• Efron (2004a) {barticle}[mr] \bauthor\bsnmEfron, \bfnmBradley\binitsB. (\byear2004a). \btitleLarge-scale simultaneous hypothesis testing: The choice of a null hypothesis. \bjournalJ. Amer. Statist. Assoc. \bvolume99 \bpages96–104. \biddoi=10.1198/016214504000000089, issn=0162-1459, mr=2054289 \bptokimsref \endbibitem
• Efron (2004b) {barticle}[mr] \bauthor\bsnmEfron, \bfnmBradley\binitsB. (\byear2004b). \btitleThe estimation of prediction error: Covariance penalties and cross-validation. \bjournalJ. Amer. Statist. Assoc. \bvolume99 \bpages619–642. \biddoi=10.1198/016214504000000692, issn=0162-1459, mr=2090899 \bptnotecheck related\bptokimsref \endbibitem
• Efron et al. (2004) {barticle}[mr] \bauthor\bsnmEfron, \bfnmBradley\binitsB., \bauthor\bsnmHastie, \bfnmTrevor\binitsT., \bauthor\bsnmJohnstone, \bfnmIain\binitsI. \AND\bauthor\bsnmTibshirani, \bfnmRobert\binitsR. (\byear2004). \btitleLeast angle regression. \bjournalAnn. Statist. \bvolume32 \bpages407–499. \biddoi=10.1214/009053604000000067, issn=0090-5364, mr=2060166 \bptnotecheck related\bptokimsref \endbibitem
• Freedman (1977) {barticle}[mr] \bauthor\bsnmFreedman, \bfnmDavid\binitsD. (\byear1977). \btitleA remark on the difference between sampling with and without replacement. \bjournalJ. Amer. Statist. Assoc. \bvolume72 \bpages681. \bidissn=0162-1459, mr=0445667 \bptokimsref \endbibitem
• Friedman (1989) {barticle}[mr] \bauthor\bsnmFriedman, \bfnmJerome H.\binitsJ. H. (\byear1989). \btitleRegularized discriminant analysis. \bjournalJ. Amer. Statist. Assoc. \bvolume84 \bpages165–175. \bidissn=0162-1459, mr=0999675 \bptokimsref \endbibitem
• Friedman, Hastie and Tibshirani (2008) {barticle}[pbm] \bauthor\bsnmFriedman, \bfnmJerome\binitsJ., \bauthor\bsnmHastie, \bfnmTrevor\binitsT. \AND\bauthor\bsnmTibshirani, \bfnmRobert\binitsR. (\byear2008). \btitleSparse inverse covariance estimation with the graphical lasso. \bjournalBiostatistics \bvolume9 \bpages432–441. \biddoi=10.1093/biostatistics/kxm045, issn=1468-4357, mid=NIHMS248717, pii=kxm045, pmcid=3019769, pmid=18079126 \bptnotecheck year\bptokimsref \endbibitem
• Gardner et al. (2003) {barticle}[pbm] \bauthor\bsnmGardner, \bfnmTimothy S.\binitsT. S., \bauthor\bparticledi \bsnmBernardo, \bfnmDiego\binitsD., \bauthor\bsnmLorenz, \bfnmDavid\binitsD. \AND\bauthor\bsnmCollins, \bfnmJames J.\binitsJ. J. (\byear2003). \btitleInferring genetic networks and identifying compound mode of action via expression profiling. \bjournalScience \bvolume301 \bpages102–105. \biddoi=10.1126/science.1081900, issn=1095-9203, pii=301/5629/102, pmid=12843395 \bptokimsref \endbibitem
• Hsieh, Cheng and Lin (2005) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmHsieh, \bfnmF. C.\binitsF. C., \bauthor\bsnmCheng, \bfnmG.\binitsG. \AND\bauthor\bsnmLin, \bfnmJ.\binitsJ. (\byear2005). \btitleEvaluation of potential Stat3-regulated genes in human breast cancer. \bjournalBiochem Biophys Res. Commun. \bvolume335 \bpages292–299. \bptokimsref \endbibitem
• Jeong et al. (2011) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmJeong, \bfnmH.\binitsH., \bauthor\bsnmMasion, \bfnmS.\binitsS., \bauthor\bsnmBarabasi, \bfnmA.\binitsA. \AND\bauthor\bsnmOltvai, \bfnmZ.\binitsZ. (\byear2011). \btitleLethality and centrality in protein networks. \bjournalNature \bvolume411 \bpages41–42. \bptokimsref \endbibitem
• Katenka and Kolaczyk (2012) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmKatenka, \bfnmN.\binitsN. \AND\bauthor\bsnmKolaczyk, \bfnmE.\binitsE. (\byear2012). \btitleInference and characterization of multi-attribute networks with application to computational biology. \bjournalAnn. Appl. Stat. \bvolume6 \bpages1068–1094. \bptokimsref \endbibitem
• Kim et al. (2006) {barticle}[pbm] \bauthor\bsnmKim, \bfnmY. H.\binitsY. H., \bauthor\bsnmGirard, \bfnmL.\binitsL., \bauthor\bsnmGiacomini, \bfnmC. P.\binitsC. P., \bauthor\bsnmWang, \bfnmP.\binitsP., \bauthor\bsnmHernandez-Boussard, \bfnmT.\binitsT., \bauthor\bsnmTibshirani, \bfnmR.\binitsR., \bauthor\bsnmMinna, \bfnmJ. D.\binitsJ. D. \AND\bauthor\bsnmPollack, \bfnmJ. R.\binitsJ. R. (\byear2006). \btitleCombined microarray analysis of small cell lung cancer reveals altered apoptotic balance and distinct expression signatures of MYC family gene amplification. \bjournalOncogene \bvolume25 \bpages130–138. \biddoi=10.1038/sj.onc.1208997, issn=0950-9232, pii=1208997, pmid=16116477 \bptokimsref \endbibitem
• Kolar et al. (2010) {barticle}[mr] \bauthor\bsnmKolar, \bfnmMladen\binitsM., \bauthor\bsnmSong, \bfnmLe\binitsL., \bauthor\bsnmAhmed, \bfnmAmr\binitsA. \AND\bauthor\bsnmXing, \bfnmEric P.\binitsE. P. (\byear2010). \btitleEstimating time-varying networks. \bjournalAnn. Appl. Stat. \bvolume4 \bpages94–123. \biddoi=10.1214/09-AOAS308, issn=1932-6157, mr=2758086 \bptokimsref \endbibitem
• Li and Gui (2006) {barticle}[pbm] \bauthor\bsnmLi, \bfnmHongzhe\binitsH. \AND\bauthor\bsnmGui, \bfnmJiang\binitsJ. (\byear2006). \btitleGradient directed regularization for sparse Gaussian concentration graphs, with applications to inference of genetic networks. \bjournalBiostatistics \bvolume7 \bpages302–317. \biddoi=10.1093/biostatistics/kxj008, issn=1465-4644, pii=kxj008, pmid=16326758 \bptokimsref \endbibitem
• Li et al. (2013) {bmisc}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmLi, \bfnmS.\binitsS., \bauthor\bsnmHsu, \bfnmL.\binitsL., \bauthor\bsnmPeng, \bfnmJ.\binitsJ. \AND\bauthor\bsnmWang, \bfnmP.\binitsP. (\byear2013). \bhowpublishedSupplement to “Bootstrap inference for network construction with an application to a breast cancer microarray study.” DOI:\doiurl10.1214/12-AOAS589SUPP. \bptokimsref \endbibitem
• Liang and Zeger (1986) {barticle}[auto] \bauthor\bsnmLiang, \bfnmKung-yee\binitsK.-y. \AND\bauthor\bsnmZeger, \bfnmScott L.\binitsS. L. (\byear1986). \btitleLongitudinal data analysis using generalized linear models. \bjournalBiometrika \bvolume73 \bpages13–22. \bptokimsref \endbibitem
• Loi et al. (2007) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmLoi, \bfnmS.\binitsS., \bauthor\bsnmHaibe-Kains, \bfnmH.\binitsH., \bauthor\bsnmDesmedt, \bfnmC.\binitsC., \bauthor\bsnmLallemand, \bfnmF.\binitsF., \bauthor\bsnmTutt, \bfnmA.\binitsA., \bauthor\bsnmGillet, \bfnmC.\binitsC., \bauthor\bsnmEllis, \bfnmP.\binitsP., \bauthor\bsnmHarris, \bfnmA.\binitsA., \bauthor\bsnmBergh, \bfnmJ.\binitsJ., \bauthor\bsnmFoekens, \bfnmJ.\binitsJ., \bauthor\bsnmKlijn, \bfnmJ.\binitsJ., \bauthor\bsnmLarsimont, \bfnmD.\binitsD., \bauthor\bsnmBuyse, \bfnmM.\binitsM., \bauthor\bsnmBontempi, \bfnmG.\binitsG., \bauthor\bsnmDelorenzi, \bfnmM.\binitsM., \bauthor\bsnmPiccart, \bfnmM.\binitsM. \AND\bauthor\bsnmSotiriou, \bfnmC.\binitsC. (\byear2007). \btitleDefinition of clinically distinct molecular subtypes in estrogen receptor–positive breast carcinomas through genomic grade. \bjournalJournal of Clinical Oncology \bvolume25 \bpages1239–1246. \bptokimsref \endbibitem
• Meinshausen and Bühlmann (2006) {barticle}[mr] \bauthor\bsnmMeinshausen, \bfnmNicolai\binitsN. \AND\bauthor\bsnmBühlmann, \bfnmPeter\binitsP. (\byear2006). \btitleHigh-dimensional graphs and variable selection with the lasso. \bjournalAnn. Statist. \bvolume34 \bpages1436–1462. \biddoi=10.1214/009053606000000281, issn=0090-5364, mr=2278363 \bptokimsref \endbibitem
• Meinshausen and Bühlmann (2010) {barticle}[mr] \bauthor\bsnmMeinshausen, \bfnmNicolai\binitsN. \AND\bauthor\bsnmBühlmann, \bfnmPeter\binitsP. (\byear2010). \btitleStability selection. \bjournalJ. R. Stat. Soc. Ser. B Stat. Methodol. \bvolume72 \bpages417–473. \biddoi=10.1111/j.1467-9868.2010.00740.x, issn=1369-7412, mr=2758523 \bptokimsref \endbibitem
• Newman (2003) {barticle}[mr] \bauthor\bsnmNewman, \bfnmM. E. J.\binitsM. E. J. (\byear2003). \btitleThe structure and function of complex networks. \bjournalSIAM Rev. \bvolume45 \bpages167–256 (electronic). \biddoi=10.1137/S003614450342480, issn=0036-1445, mr=2010377 \bptokimsref \endbibitem
• Nie, Wu and Zhang (2006) {barticle}[pbm] \bauthor\bsnmNie, \bfnmLei\binitsL., \bauthor\bsnmWu, \bfnmGang\binitsG. \AND\bauthor\bsnmZhang, \bfnmWeiwen\binitsW. (\byear2006). \btitleCorrelation between mRNA and protein abundance in Desulfovibrio vulgaris: A multiple regression to identify sources of variations. \bjournalBiochem. Biophys. Res. Commun. \bvolume339 \bpages603–610. \biddoi=10.1016/j.bbrc.2005.11.055, issn=0006-291X, pii=S0006-291X(05)02584-2, pmid=16310166 \bptokimsref \endbibitem
• Pathak (1962) {barticle}[mr] \bauthor\bsnmPathak, \bfnmP. K.\binitsP. K. (\byear1962). \btitleOn simple random sampling with replacement. \bjournalSankhyā Ser. A \bvolume24 \bpages287–302. \bidissn=0581-572X, mr=0169350 \bptokimsref \endbibitem
• Peng et al. (2009) {barticle}[mr] \bauthor\bsnmPeng, \bfnmJie\binitsJ., \bauthor\bsnmWang, \bfnmPei\binitsP., \bauthor\bsnmZhou, \bfnmNengfeng\binitsN. \AND\bauthor\bsnmZhu, \bfnmJi\binitsJ. (\byear2009). \btitlePartial correlation estimation by joint sparse regression models. \bjournalJ. Amer. Statist. Assoc. \bvolume104 \bpages735–746. \biddoi=10.1198/jasa.2009.0126, issn=0162-1459, mr=2541591 \bptokimsref \endbibitem
• Peng et al. (2010) {barticle}[mr] \bauthor\bsnmPeng, \bfnmJie\binitsJ., \bauthor\bsnmZhu, \bfnmJi\binitsJ., \bauthor\bsnmBergamaschi, \bfnmAnna\binitsA., \bauthor\bsnmHan, \bfnmWonshik\binitsW., \bauthor\bsnmNoh, \bfnmDong-Young\binitsD.-Y., \bauthor\bsnmPollack, \bfnmJonathan R.\binitsJ. R. \AND\bauthor\bsnmWang, \bfnmPei\binitsP. (\byear2010). \btitleRegularized multivariate regression for identifying master predictors with application to integrative genomics study of breast cancer. \bjournalAnn. Appl. Stat. \bvolume4 \bpages53–77. \biddoi=10.1214/09-AOAS271, issn=1932-6157, mr=2758084 \bptokimsref \endbibitem
• Pollack et al. (2002) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmPollack, \bfnmJ.\binitsJ., \bauthor\bsnmSrlie, \bfnmT.\binitsT., \bauthor\bsnmPerou, \bfnmC.\binitsC., \bauthor\bsnmRees, \bfnmC.\binitsC., \bauthor\bsnmJeffrey, \bfnmS.\binitsS., \bauthor\bsnmLonning, \bfnmP.\binitsP., \bauthor\bsnmTibshirani, \bfnmR.\binitsR., \bauthor\bsnmBotstein, \bfnmD.\binitsD., \bauthor\bsnmBrresen-Dale, \bfnmA.\binitsA. \AND\bauthor\bsnmBrown, \bfnmP.\binitsP. (\byear2002). \btitleMicroarray analysis reveals a major direct role of DNA copy number alteration in the transcriptional program of human breast tumors. \bjournalProc. Natl. Acad. Sci. USA \bvolume99 \bpages12963–12968. \bptokimsref \endbibitem
• Postel-Vinay et al. (2012) {barticle}[auto] \bauthor\bsnmPostel-Vinay, \bfnmS.\binitsS., \bauthor\bsnmVéron, \bfnmA. S.\binitsA. S., \bauthor\bsnmTirode, \bfnmF.\binitsF., \bauthor\bsnmPierron, \bfnmG.\binitsG., \bauthor\bsnmReynaud, \bfnmS.\binitsS., \bauthor\bsnmKovar, \bfnmH.\binitsH., \bauthor\bsnmOberlin, \bfnmO.\binitsO., \bauthor\bsnmLapouble, \bfnmE.\binitsE., \bauthor\bsnmBallet, \bfnmS.\binitsS., \bauthor\bsnmLucchesi, \bfnmC.\binitsC., \bauthor\bsnmKontny, \bfnmU.\binitsU., \bauthor\bsnmGonzález-Neira, \bfnmA.\binitsA., \bauthor\bsnmPicci, \bfnmP.\binitsP., \bauthor\bsnmAlonso, \bfnmJ.\binitsJ., \bauthor\bsnmPatino-Garcia, \bfnmA.\binitsA., \bauthor\bsnmde Paillerets, \bfnmB. B.\binitsB. B., \bauthor\bsnmLaud, \bfnmK.\binitsK., \bauthor\bsnmDina, \bfnmC.\binitsC., \bauthor\bsnmFroguel, \bfnmP.\binitsP., \bauthor\bsnmClavel-Chapelon, \bfnmF.\binitsF., \bauthor\bsnmDoz, \bfnmF.\binitsF., \bauthor\bsnmMichon, \bfnmJ.\binitsJ., \bauthor\bsnmChanock, \bfnmS. J.\binitsS. J., \bauthor\bsnmThomas, \bfnmG.\binitsG., \bauthor\bsnmCox, \bfnmD. G.\binitsD. G. \AND\bauthor\bsnmDelattre, \bfnmO.\binitsO. (\byear2012). \btitleCommon variants near TARDBP and EGR2 are associated with susceptibility to Ewing sarcoma. \bjournalNat. Genet. \bvolume44 \bpages323–327. \bptokimsref \endbibitem
• Rothman et al. (2008) {barticle}[mr] \bauthor\bsnmRothman, \bfnmAdam J.\binitsA. J., \bauthor\bsnmBickel, \bfnmPeter J.\binitsP. J., \bauthor\bsnmLevina, \bfnmElizaveta\binitsE. \AND\bauthor\bsnmZhu, \bfnmJi\binitsJ. (\byear2008). \btitleSparse permutation invariant covariance estimation. \bjournalElectron. J. Stat. \bvolume2 \bpages494–515. \biddoi=10.1214/08-EJS176, issn=1935-7524, mr=2417391 \bptokimsref \endbibitem
• Schadt et al. (2005) {barticle}[pbm] \bauthor\bsnmSchadt, \bfnmEric E.\binitsE. E., \bauthor\bsnmLamb, \bfnmJohn\binitsJ., \bauthor\bsnmYang, \bfnmXia\binitsX., \bauthor\bsnmZhu, \bfnmJun\binitsJ., \bauthor\bsnmEdwards, \bfnmSteve\binitsS., \bauthor\bsnmGuhathakurta, \bfnmDebraj\binitsD., \bauthor\bsnmSieberts, \bfnmSolveig K.\binitsS. K., \bauthor\bsnmMonks, \bfnmStephanie\binitsS., \bauthor\bsnmReitman, \bfnmMarc\binitsM., \bauthor\bsnmZhang, \bfnmChunsheng\binitsC., \bauthor\bsnmLum, \bfnmPek Yee\binitsP. Y., \bauthor\bsnmLeonardson, \bfnmAmy\binitsA., \bauthor\bsnmThieringer, \bfnmRolf\binitsR., \bauthor\bsnmMetzger, \bfnmJoseph M.\binitsJ. M., \bauthor\bsnmYang, \bfnmLiming\binitsL., \bauthor\bsnmCastle, \bfnmJohn\binitsJ., \bauthor\bsnmZhu, \bfnmHaoyuan\binitsH., \bauthor\bsnmKash, \bfnmShera F.\binitsS. F., \bauthor\bsnmDrake, \bfnmThomas A.\binitsT. A., \bauthor\bsnmSachs, \bfnmAlan\binitsA. \AND\bauthor\bsnmLusis, \bfnmAldons J.\binitsA. J. (\byear2005). \btitleAn integrative genomics approach to infer causal associations between gene expression and disease. \bjournalNat. Genet. \bvolume37 \bpages710–717. \biddoi=10.1038/ng1589, issn=1061-4036, mid=NIHMS175481, pii=ng1589, pmcid=2841396, pmid=15965475 \bptokimsref \endbibitem
• Schäfer and Strimmer (2005) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmSchäfer, \bfnmJ.\binitsJ. \AND\bauthor\bsnmStrimmer, \bfnmK.\binitsK. (\byear2005). \btitleLearning large-scale graphical Gaussian models from genomic data. \bjournalAIP Conf. Proc. \bvolume776 \bpages263–276. \bptokimsref \endbibitem
• Schwarz (1978) {barticle}[mr] \bauthor\bsnmSchwarz, \bfnmGideon\binitsG. (\byear1978). \btitleEstimating the dimension of a model. \bjournalAnn. Statist. \bvolume6 \bpages461–464. \bidissn=0090-5364, mr=0468014 \bptokimsref \endbibitem
• Shackleford et al. (2011) {barticle}[pbm] \bauthor\bsnmShackleford, \bfnmTerry J.\binitsT. J., \bauthor\bsnmZhang, \bfnmQingxiu\binitsQ., \bauthor\bsnmTian, \bfnmLing\binitsL., \bauthor\bsnmVu, \bfnmThuy T.\binitsT. T., \bauthor\bsnmKorapati, \bfnmAnita L.\binitsA. L., \bauthor\bsnmBaumgartner, \bfnmAngela M.\binitsA. M., \bauthor\bsnmLe, \bfnmXiao-Feng\binitsX.-F., \bauthor\bsnmLiao, \bfnmWarren S.\binitsW. S. \AND\bauthor\bsnmClaret, \bfnmFrancois X.\binitsF. X. (\byear2011). \btitleStat3 and CCAAT/enhancer binding protein beta (C/EBP-beta) regulate Jab1/CSN5 expression in mammary carcinoma cells. \bjournalBreast Cancer Res. \bvolume13 \bpagesR65. \biddoi=10.1186/bcr2902, issn=1465-542X, pii=bcr2902, pmcid=3218954, pmid=21689417 \bptokimsref \endbibitem
• Storey (2003) {barticle}[mr] \bauthor\bsnmStorey, \bfnmJohn D.\binitsJ. D. (\byear2003). \btitleThe positive false discovery rate: A Bayesian interpretation and the -value. \bjournalAnn. Statist. \bvolume31 \bpages2013–2035. \biddoi=10.1214/aos/1074290335, issn=0090-5364, mr=2036398 \bptokimsref \endbibitem
• Tegner et al. (2003) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmTegner, \bfnmJ.\binitsJ., \bauthor\bsnmYeung, \bfnmM.\binitsM., \bauthor\bsnmHasty, \bfnmJ.\binitsJ. \AND\bauthor\bsnmCollins, \bfnmJ.\binitsJ. (\byear2003). \btitleReverse engineering gene networks: Integrating genetic perturbations with dynamical modeling. \bjournalProc. Natl. Acad. Sci. USA \bvolume100 \bpages5944–5949. \bptokimsref \endbibitem
• Tibshirani (1996) {barticle}[mr] \bauthor\bsnmTibshirani, \bfnmRobert\binitsR. (\byear1996). \btitleRegression shrinkage and selection via the lasso. \bjournalJ. Roy. Statist. Soc. Ser. B \bvolume58 \bpages267–288. \bidissn=0035-9246, mr=1379242 \bptokimsref \endbibitem
• Tront, Hoffman and Liebermann (2006) {barticle}[pbm] \bauthor\bsnmTront, \bfnmJennifer S.\binitsJ. S., \bauthor\bsnmHoffman, \bfnmBarbara\binitsB. \AND\bauthor\bsnmLiebermann, \bfnmDan A.\binitsD. A. (\byear2006). \btitleGadd45a suppresses Ras-driven mammary tumorigenesis by activation of c-Jun NH2-terminal kinase and p38 stress signaling resulting in apoptosis and senescence. \bjournalCancer Res. \bvolume66 \bpages8448–8454. \biddoi=10.1158/0008-5472.CAN-06-2013, issn=1538-7445, pii=66/17/8448, pmid=16951155 \bptokimsref \endbibitem
• Varambally et al. (2005) {barticle}[pbm] \bauthor\bsnmVarambally, \bfnmSooryanarayana\binitsS., \bauthor\bsnmYu, \bfnmJianjun\binitsJ., \bauthor\bsnmLaxman, \bfnmBharathi\binitsB., \bauthor\bsnmRhodes, \bfnmDaniel R.\binitsD. R., \bauthor\bsnmMehra, \bfnmRohit\binitsR., \bauthor\bsnmTomlins, \bfnmScott A.\binitsS. A., \bauthor\bsnmShah, \bfnmRajal B.\binitsR. B., \bauthor\bsnmChandran, \bfnmUma\binitsU., \bauthor\bsnmMonzon, \bfnmFederico A.\binitsF. A., \bauthor\bsnmBecich, \bfnmMichael J.\binitsM. J., \bauthor\bsnmWei, \bfnmJohn T.\binitsJ. T., \bauthor\bsnmPienta, \bfnmKenneth J.\binitsK. J., \bauthor\bsnmGhosh, \bfnmDebashis\binitsD., \bauthor\bsnmRubin, \bfnmMark A.\binitsM. A. \AND\bauthor\bsnmChinnaiyan, \bfnmArul M.\binitsA. M. (\byear2005). \btitleIntegrative genomic and proteomic analysis of prostate cancer reveals signatures of metastatic progression. \bjournalCancer Cell \bvolume8 \bpages393–406. \biddoi=10.1016/j.ccr.2005.10.001, issn=1535-6108, pii=S1535-6108(05)00305-3, pmid=16286247 \bptokimsref \endbibitem
• Waaijenborg, Verselewel de Witt Hamer and Zwinderman (2008) {barticle}[auto:STB—2012/09/13—11:39:40] \bauthor\bsnmWaaijenborg, \bfnmSandra\binitsS., \bauthor\bsnmVerselewel de Witt Hamer, \bfnmPhilip\binitsP. \AND\bauthor\bsnmZwinderman, \bfnmAeilko\binitsA. (\byear2008). \btitleQuantifying the association between gene expressions and DNA-markers by penalized canonical correlation analysis. \bjournalStat. Appl. Genet. Mol. Biol. \bvolume7 \bpages1329. \bptokimsref \endbibitem
• Wainwright (2009) {barticle}[mr] \bauthor\bsnmWainwright, \bfnmMartin J.\binitsM. J. (\byear2009). \btitleSharp thresholds for high-dimensional and noisy sparsity recovery using -constrained quadratic programming (Lasso). \bjournalIEEE Trans. Inform. Theory \bvolume55 \bpages2183–2202. \biddoi=10.1109/TIT.2009.2016018, issn=0018-9448, mr=2729873 \bptokimsref \endbibitem
• Wang et al. (2011) {barticle}[mr] \bauthor\bsnmWang, \bfnmSijian\binitsS., \bauthor\bsnmNan, \bfnmBin\binitsB., \bauthor\bsnmRosset, \bfnmSaharon\binitsS. \AND\bauthor\bsnmZhu, \bfnmJ.\binitsJ. (\byear2011). \btitleRandom Lasso. \bjournalAnn. Appl. Stat. \bvolume5 \bpages468–485. \biddoi=10.1214/10-AOAS377, issn=1932-6157, mr=2810406 \bptokimsref \endbibitem
• Whittaker (1990) {bbook}[mr] \bauthor\bsnmWhittaker, \bfnmJoe\binitsJ. (\byear1990). \btitleGraphical Models in Applied Multivariate Statistics. \bpublisherWiley, \blocationChichester. \bidmr=1112133 \bptokimsref \endbibitem
• Witten, Tibshirani and Hastie (2009) {barticle}[pbm] \bauthor\bsnmWitten, \bfnmDaniela M.\binitsD. M., \bauthor\bsnmTibshirani, \bfnmRobert\binitsR. \AND\bauthor\bsnmHastie, \bfnmTrevor\binitsT. (\byear2009). \btitleA penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. \bjournalBiostatistics \bvolume10 \bpages515–534. \biddoi=10.1093/biostatistics/kxp008, issn=1468-4357, pii=kxp008, pmcid=2697346, pmid=19377034 \bptokimsref \endbibitem
• Yuan and Lin (2007) {barticle}[mr] \bauthor\bsnmYuan, \bfnmMing\binitsM. \AND\bauthor\bsnmLin, \bfnmYi\binitsY. (\byear2007). \btitleModel selection and estimation in the Gaussian graphical model. \bjournalBiometrika \bvolume94 \bpages19–35. \biddoi=10.1093/biomet/asm018, issn=0006-3444, mr=2367824 \bptokimsref \endbibitem
• Zhao and Yu (2006) {barticle}[mr] \bauthor\bsnmZhao, \bfnmPeng\binitsP. \AND\bauthor\bsnmYu, \bfnmBin\binitsB. (\byear2006). \btitleOn model selection consistency of Lasso. \bjournalJ. Mach. Learn. Res. \bvolume7 \bpages2541–2563. \bidissn=1532-4435, mr=2274449 \bptokimsref \endbibitem
• Zou (2006) {barticle}[mr] \bauthor\bsnmZou, \bfnmHui\binitsH. (\byear2006). \btitleThe adaptive lasso and its oracle properties. \bjournalJ. Amer. Statist. Assoc. \bvolume101 \bpages1418–1429. \biddoi=10.1198/016214506000000735, issn=0162-1459, mr=2279469 \bptokimsref \endbibitem
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters