Effects of dependence in highdimensional multiple testing problems
 Kyung In Kim^{1}Email author and
 Mark A van de Wiel^{2, 3}
DOI: 10.1186/147121059114
© Kim and van de Wiel; licensee BioMed Central Ltd. 2008
Received: 13 August 2007
Accepted: 25 February 2008
Published: 25 February 2008
Abstract
Background
We consider effects of dependence among variables of highdimensional data in multiple hypothesis testing problems, in particular the False Discovery Rate (FDR) control procedures. Recent simulation studies consider only simple correlation structures among variables, which is hardly inspired by real data features. Our aim is to systematically study effects of several network features like sparsity and correlation strength by imposing dependence structures among variables using random correlation matrices.
Results
We study the robustness against dependence of several FDR procedures that are popular in microarray studies, such as BenjaminHochberg FDR, Storey's qvalue, SAM and resampling based FDR procedures. False Nondiscovery Rates and estimates of the number of null hypotheses are computed from those methods and compared. Our simulation study shows that methods such as SAM and the qvalue do not adequately control the FDR to the level claimed under dependence conditions. On the other hand, the adaptive BenjaminiHochberg procedure seems to be most robust while remaining conservative. Finally, the estimates of the number of true null hypotheses under various dependence conditions are variable.
Conclusion
We discuss a new method for efficient guided simulation of dependent data, which satisfy imposed network constraints as conditional independence structures. Our simulation setup allows for a structural study of the effect of dependencies on multiple testing criterions and is useful for testing a potentially new method on π_{0} or FDR estimation in a dependency context.
Background
Scientists regularly face multiple testing of a large number of hypotheses nowadays. Typically in microarray data, one performs hypothesis testing for each gene and the number of genes is usually more than thousands. In this situation, direct application of single hypothesis testing thousands times produces a large number of false discoveries. Hence, alternative testing criterions for controlling errors of false discoveries have been introduced.
It is widely recognized that dependencies are omnipresent in many highthroughput studies. Such dependencies may be regulatory or functional as in gene pathways, but also spatial such as in SNP or DNA copy number arrays because of the genomic order. Although attempts to infer such interactions from data have been made, it is a notoriously difficult problem. Usually solutions focus on some modules with relatively few elements and many samples, in particular for model organisms (see e.g. [1]). With this in mind, one prefers multiple testing methods that are robust to several degrees of dependency in these networktype data. Therefore, we set out to develop a simulation program that allows us testing any multiple testing method for its robustness with respect to dependency parameters using realistic nested network structures.
One of the most widely used multiple testing criterions for controlling errors of false discoveries is False Discovery Rate (FDR). FDR is introduced in Benjamini et al. [2] and is defined as the expected proportion of the number of falsely rejected hypotheses among total number of rejected hypotheses. Since in most cases, underlying distributions of data are unknown, there are several implementations of FDR under different assumptions.
under independence assumption of test statistics. Later, Bejamini and Yekutieli [3] prove the BH procedure controls under positive regression dependency condition and they introduce a modification of the above procedure to control arbitrary dependence circumstances (BY). Storey [4] introduces the positive false discovery rate (pFDR) and the qvalue. pFDR is known to control error rate under the clumpy dependency condition [5]. Significance Analysis of Microarray (SAM) is developed on the purpose of statistical analysis of microarray data [6]. SAM FDR is known to estimate the expected number of false discoveries over the observed number of total rejections under the complete null hypothesis [7].
In (1), there still remains some space of improvement for tighter control if we know π_{0}. Adaptive procedures are developed to gain more power by estimating π_{0} in this sense. To estimate π_{0}, Storey et al. [5] use the fact that independent null pvalues are distributed uniformly on [0, 1] and then plug the estimator ${\widehat{\pi}}_{0}$ into the FDRestimator. Benjamini et al. [8] estimate m_{0} in a twostage adaptive control of FDR (ABH). Under the assumption of independence of test statistics, they show the ABH procedure controls nominal significance level. Careful simulation studies on the comparison of performance of π_{0} estimation methods are done by Black [9] and Langaas et al. [10]. Black [9] notes that the violation of uniformity of pvalues due to the presence of nonnull cases could bias estimates of π_{0} upward.
Recently, several works incorporates correlations among test statistics to estimate FDR. Resampling based approaches are immediate in utilizing sample correlation structure [11]. However, since it is difficult to resample from the true null and the false null distributions separately, it is common to assume the complete null hypothesis and set the number of true discoveries fixed in order to estimate FDR conservatively, as is shown in the resampling based method of Yekutieli and his coworkers [12, 13]. Since the procedures mentioned above are often used, we would like to study validity of those procedures under fairly general dependence circumstances and how correlations among test statistics affect results of FDR multiple testings. Also, we would like to examine effects of violation of independence of pvalues on π_{0} estimations. Hence, designing general dependence conditions is our main concern. In previous works, for convenience of simulations, data are often assumed multivariate normally distributed. Typically in microarray data analysis, equicorrelated normal structures such as single pairwise correlation matrices or block diagonal matrices with a single pairwise correlation in each block are considered [14, 15].
Equicorrelated structures are easy to understand and implement. Moreover, control of dependence conditions is easy by increasing or decreasing single correlations. But they are hardly regarded to represent reality. Random correlation matrices are more realistic candidates, because they reflect heterogeneity between the correlations. However, since the class of random correlation matrices is too large, multiple testing results generated from two arbitrary random correlation matrices are difficult to compare.
Therefore, we propose constrained random correlation matrices to reflect the generality of random correlations and the comparability like equicorrelation models to simulation studies. For simulation studies, we generate a sequence of random correlation matrices and as constraints we impose conditional independence structures on the random correlation matrices in a 'nested' way. Then the sequence of random correlation matrices is ordered in terms of a dependence parameter and we control the strength of dependence by the dependence parameter. An alternative, nonnested, approach is discussed by Jung et al. [16] who simulate multivariate normal test statistics while conserving the correlation structure as present in the data in an asymptotic sense.
In our simulation results, we show how the dependence parameter can serve as a measure of FDR behavior under correlationbased dependence conditions. We prove that this dependence parameter is in fact strongly related to the variance of pairwise correlations. Using this structural simulation setting, we compare the performance of several FDR estimating methods.
Results
We illustrate simulation results. Here, we consider two cases for the proportion of true null hypotheses: π_{0} = 0.8 and π_{0} = 0.95. Both cases show similar results, so we focus on the first case. For π_{0} = 0.95, we refer to Figure S12S14 [see Additional file 1]. We do not take into account for small π_{0}'s because in highdimensional data with thousands hypotheses one is usually interested in the case when only small portions of total hypotheses are truly significant. We generate 16 correlation matrices Σ based on 16 edge densities, which are the proportions of nonzero partial correlations over all possible pairs of partial correlations, 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1. Note for a nested sequence of random correlation matrices, we use one initial Z matrix (see Algorithm 2) for each π_{0}. The total number of hypotheses is set to m = 1000 and sample sizes for X and Y are 10 each. The number of resamples to compute average FDR in (5) are B = 2000. The fixed true difference is chosen to have 80% power for individual two group tstatistic when FDR significance level is 0.1 under independence assumption.
It is hard to decide that which one is recommended in practice when apparent dependence is observed. However, in this simulation, if most weight is given on adhering strict control level and gaining more power is a secondary goal, the ABH seems to be most robust and desirable under dependence cases.
An illustration with real data
In this section, we address an example on how to apply biological information such as pathways using random correlation matrices. Basically, we use estimated marginal mean and variance from data and apply pathway information such as gene sets to correlation structures. Algorithm 3 shows the detailed procedure. It uses Algorithms 1 and 2, which are discussed in the Methods section.
 1.
Compute mdimensional sample mean and sample variance vectors, ${\widehat{\mu}}_{X},{\widehat{\mu}}_{Y},{\widehat{\sigma}}_{X}^{2},{\widehat{\sigma}}_{Y}^{2}$ from data ${X}_{1},\mathrm{...},{X}_{{n}_{1}}$ and ${Y}_{1},\mathrm{...},{Y}_{{n}_{2}}$.
 2.
Prepare interested gene sets GS_{1},..., GS_{ k }and make a sequence of nested gene sets N_{1},..., N_{ k }by iterative merging. That is, for each j = 1,..., k, N_{ j }= GS_{1} ∪ ... ∪ GS_{ j }.
 3.
Generate a sequence of binary adjacency matrices A_{1},..., A_{ k }from N_{1},..., N_{ k }. Components of adjacency matrices are encoded as 1 for presence of edge and 0 for absence of edge. For example, [A_{ l }]_{i, j}= 1 means both ith and jth gene are in N_{ l }.
 4.
According to A_{1},..., A_{ k }, generate a sequence of random correlation matrices, R_{1},..., R_{ k }, using Algorithms 1 and 2.
 5.
Generate sample from ${X}_{1}^{\ast b},\mathrm{...},{X}_{{n}_{1}}^{\ast b}~N({\widehat{\mu}}_{X},\text{diag}({\widehat{\sigma}}_{X}){R}_{j}\text{diag}({\widehat{\sigma}}_{X}))$ and ${Y}_{1}^{\ast b},\mathrm{...},{Y}_{{n}_{1}}^{\ast b}~N({\widehat{\mu}}_{Y},\text{diag}({\widehat{\sigma}}_{Y}){R}_{j}\text{diag}({\widehat{\sigma}}_{Y}))$ for b = 1,..., B.
 6.
Do multiple testing B times and estimate average FDR from (5).
We applied the above algorithm to the "Two Class" example data of Excel addin of SAM which consists of 1000 genes with 10 control and 10 treatment experiments. Along with the data provided, we used gene sets file, "c2.v2.symbols.gmt" for pathway information from MSigDB [21]. There are 1687 gene sets in the file and we chose those 10 gene sets (Gene Set 291, 698, 861, 885, 1069, 1177, 1179, 1237, 1345, 1453) which overlap more than 50 genes with the gene list of the "Two Class" data.
Median number of detected genes under increasing edge densities and the corresponding correlation matrices
R _{1}  R _{2}  R _{3}  R _{4}  R _{5}  R _{6}  R _{7}  R _{8}  R _{9}  R _{10}  

edge density  0.003  0.012  0.022  0.037  0.067  0.089  0.107  0.140  0.169  0.182 
#total discoveries  110  110  110  110  110  109  108  107  107  106 
16 genes showing different significance feature under nested 10 correlation matrices
R _{1}  R _{2}  R _{3}  R _{4}  R _{5}  R _{6}  R _{7}  R _{8}  R _{9}  R _{10}  

PRKCZ  1  1  1  1  1  1  1  0  1  1 
HSPA4  1  1  0  0  0  1  0  0  0  1 
SIAT7B  0  0  1  1  0  0  0  0  1  0 
40222_s_at  1  1  1  1  1  0  1  0  1  1 
36374_at  0  0  0  0  0  0  0  1  1  0 
1627_at  0  0  0  0  0  1  0  0  0  0 
SSR1  1  1  1  1  1  0  1  1  1  0 
SEDLP  1  1  1  1  1  0  1  0  0  0 
VG5Q  0  0  0  0  1  1  0  0  0  0 
MAN2B1  1  1  1  1  1  1  1  1  0  1 
NDUFS1  0  0  0  1  1  0  0  0  0  0 
AMT  1  1  1  0  1  1  0  1  1  1 
STX3A  1  1  1  1  1  1  1  1  0  1 
AP3S2  1  1  1  1  1  1  1  1  0  0 
SLC35A2  0  0  0  0  1  0  0  0  0  0 
METTL3  1  1  1  1  0  1  1  1  1  0 
Interpretation on the results of Table 2 depends on the specific correlation structures given in R_{1},..., R_{10} and there does not seem clear trends in rejections for 16 genes. Since marginal distributions of single genes do not change when we apply various correlation structures to correlation matrices of the multivariate normal distribution, the result that almost all detected genes were the same confirms our expectation.
Discussion and Conclusion
We considered effects of dependence on FDR multiple testing results using multivariate normal samples. We found that in all our simulations, the simple adaptive BenjaminiHochberg procedure [8] is most optimal under dependence, since it achieves relatively high power while remaining conservative. By definition, FDR is the expected value of a nonlinear function of indicator random variables of rejection. Hence, for computations of FDR, we need to take into account of the joint distribution of the indicator random variables. To focus on joint distributional properties of FDR, we have designed to observe variations of FDR in terms of variations of correlation structures and we have fixed other parameters such as marginal distributions and probabilities of rejections for true null and false null hypotheses. Therefore, our results could be additional useful guideline to FDR estimation methods which have been developed based on marginal distributional assumptions.
Nowadays, explaining highdimensional data with conditional independence structures is quite active especially in microarray data analysis [1, 22–24]. Such methods focus on testing on partial correlation coefficients. The necessary and sufficient condition of zero partial correlation is the same as (2). The results of testings on partial correlations is a network which can be used directly in our simulation framework when, for example, testing on difference of means between two groups of samples. Then, our simulation setup can be regarded as a dataguided simulation to study whether a particular multiple testing method is useful for the data at hand. As a dataguided simulation using known gene sets [21, 25], we introduce an algorithm for using real data in the Results section. Although a very slight downward trend for the number of discoveries with respect to increasing edge density (dependence) is found, we observe that the BH FDR method is very robust in this setting as well.
In our simulation study, we did not categorize test statistics. Most of the FDR methods in the Results section are based on simple gene specific tstatistic, while SAM uses its own statistic using the fudge factor which stabilize estimates of genewise variances. The effects of using such modified tstatistic are not clear but we can reflect those effects from the viewpoint of sample sizes. As sample sizes increase, the fudge factor of SAM shows a convergence feature, although it does not improve SAM's anticonservative bias under dependence conditions. As an alternative to the fudge factor, the random variance model (RVM) by [26] can be used and simple replacement of the pooled variance of tstatistic by the RVM variance results in close control of the FDR to the nominal level under dependence in moderate to large sample size conditions. For the effects of various sample sizes on the fudge factor and π_{0} estimates of SAM and the RVM FDR, see Figure S1S4 of Additional file 1.
Effect size may be another important factor in evaluating FDR methods. We consider the cases for multiple small effect sizes or very small proportion of fixed effect size, for example π_{0} = 0.99. In both cases, we observe overall similar patterns of the FDR estimates shown in the Results section [see Figure S8S11 of the Additional file 1].
Generally in highdimensional situation, we doubt that the permutational based approach to estimate joint distributional properties of test statistics always give a correct answer. In a further simulation study, the estimated spread of ordered SAM statistics under permutational null hypothesis shows to be narrower than that of the true distribution. Note that the difference becomes wider as edge density increases. This seems to cause the anticonservative feature of SAM under dependence. For more detail on the effect of sample size and the performance of SAM and RVM, see Appendix 2 of Additional file 1.
Efron [20] notices that variance of pairwise correlations plays an important role in characterizing FDR, defined somewhat differently as the expected number of false rejections over the observed number of rejections, E(V)/R. We confirm this finding, but in our networkbased simulation setup, we found it natural to characterize FDR using two parameters: first, edge density to decide the proportion of interactions present and second the variance of pairwise correlations. This allows to study multiple testing methods for a given prior network.
Other interesting works on the effects of dependence on the number of false discoveries rather than FDR are Owen [27] and Korn et al. [15] who discuss that large positive correlations may result in high variation on the number of false discoveries. Under simple correlation structures, Qiu et al. [14] investigate the vulnerability of application of empirical Bayes theory under strong correlations.
One can extend our simulation framework by considering the distribution of the Z matrix. Until now, we have considered the constrained random correlation matrices depending on the fixed single Z matrix and given nested structures. Taking into account the distributional properties of Z as a prior, one may attain explicit posterior distribution of covariance matrices Σ_{1},...,Σ_{ d }. A family of covariance matrices as a Gaussian ensemble can also be considered as described in [28]. However, both approaches require very complicated mathematical computations so we remain these as future works.
Our simulation setup is also useful for testing a potentially new method on π_{0} or FDR estimation in a dependency context. One may not have designed the procedure for the multivariate normal setting in particular; however, it seems reasonable that the new method should perform well under these conditions to be of general use. Or one may at least sketch the boundaries of the usefulness of the method in terms of type of network, edge density, and correlation strength.
Methods
In this section, firstly, we introduce the property of conditional independence in multivariate normal distributions and its implications as graphical representations. Secondly, we introduce how to incorporate conditional independence structures to random correlation matrices and how to generate constrained random correlation matrices in a 'nested' way. Thirdly, we introduce FDR methods and π_{0} estimation methods used in this simulation study.
Conditional independence correlation structures
Here, "h" represents independence between random variables.
Given m dimensional multivariate normal distribution, however, there are ${2}^{\left(\begin{array}{c}m\\ 2\end{array}\right)}$ different conditional independence structures or graphs. Comparing every pair of structures for large m is infeasible. But note that the class of structures is a partially ordered set by inclusions or exclusions of edges. In the partially ordered set, the minimal element is the totally independence structure corresponding to identity variancecovariance matrix in a matrix form. Maximal elements are completely dependent structures without any conditional independence constraints, that is every entry of inverse of the variancecovariance matrix is nonzero. Hence, it is meaningful to regard a partially ordered path as a sequence of dependence conditions as in the single correlation structures. Comparisons through a partially ordered path give insights on dependence effects. Then, it is natural to regard the proportions of edges in a path as a dependence parameter. Figure 7 shows such an instance of the partially ordered path. In following sections, we use the term 'edge density' as proportion of edges and by a 'nested' sequence we mean a partially ordered path of conditional independence structures.
Generating constrained random correlation matrices
Unconstrained random correlation matrices are generated simply by products of matrix transposes and its standardizations [31]. Let Z be an M × m matrix whose entries are generated from independent standard normal distributions. If M is greater than m, then the matrix W = (Z^{ T }Z)^{1} is a symmetric positive definite matrix with probability one. M will be used as a parameter to control the strength of the correlations. After standardizing W, we obtainΣ = diag(W)^{1/2} W diag(W)^{1/2}.
Then Σ is an unconstrained random positive definite correlation matrix.
To incorporate conditional independence structures into the process (3), we need to transform the Z matrix into a matrix $\tilde{Z}$ such that $\tilde{Z}$ bears the information on the structures. These transformations are basically based on successive orthogonal projections. For a simple example, consider imposing the simple constraint X_{1} ╨X_{2}{rest} on Σ in (3). We carry out the following steps. First, we generate the Z = [z_{1},..., z_{ m }] matrix with m column vectors. Second, we let ${\tilde{z}}_{2}={z}_{2}{z}_{1}{({z}_{1}^{T}{z}_{1})}^{1}{z}_{1}^{T}{z}_{2}$, then ${\tilde{z}}_{2}$ is the residual vector of z_{2} projected onto the linear space spanned by z_{1}. Finally, if we replace matrix W in (3) by ${({\tilde{Z}}^{T}\tilde{Z})}^{1}$ where $\tilde{Z}=[{z}_{1},{\tilde{z}}_{2},{z}_{3},\mathrm{...},{z}_{m}]$, then Σ is a random correlation matrix satisfying the constraint [Σ^{1}]_{12} = 0 by construction.
Constraint matrices corresponding to the graphs in Figure 7
(a) e_{1} = 0/6 $\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& 0\\ 1& 0\\ 1\end{array}\right]$  (b) e_{2} = 1/6 $\left[\begin{array}{cccc}1& 1& 0& 0\\ 1& 0& 0\\ 1& 0\\ 1\end{array}\right]$  (c) e_{3} = 2/6 $\left[\begin{array}{cccc}1& 1& 0& 0\\ 1& 1& 0\\ 1& 0\\ 1\end{array}\right]$  (d) e_{4} = 3/6 $\left[\begin{array}{cccc}1& 1& 0& 0\\ 1& 1& 0\\ 1& 1\\ 1\end{array}\right]$ 

(e) e_{5} = 4/6 $\left[\begin{array}{cccc}1& 1& 0& 1\\ 1& 1& 0\\ 1& 1\\ 1\end{array}\right]$  (f) e_{6} = 5/6 $\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 0\\ 1& 1\\ 1\end{array}\right]$  (g) e_{7} = 6/6 $\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 1\\ 1& 1\\ 1\end{array}\right]$ 
Now, we provide two algorithms used for our simulation studies. Basically, we apply the second algorithm and the first one is included in the second one.
 1.
Generate an Z = [z_{1},..., z_{ m }] matrix from standard normal distributions.
 2.
Let I_{ l }= {k : [J]_{ kl }= 0 for k = 0,..., l  1} for l = 1,..., m and ${z}_{{I}_{l}}$ be the matrix consisting of column vectors of Z with indices in I_{ l }.
 3.
Let ${\tilde{z}}_{1}$ = z_{1}.
 4.
For each l = 2,..., m, ${\tilde{z}}_{l}$ = z_{ l } P_{ l }z_{ l }where ${P}_{l}={\tilde{z}}_{{I}_{l}}{({\tilde{z}}_{{I}_{l}}^{T}{\tilde{z}}_{{I}_{l}})}^{1}{\tilde{z}}_{{I}_{l}}^{T}$, i.e. the projection of z_{ l }onto the space spanned by {${\tilde{z}}_{i}$ : i ∈ I_{ l }}.
 5.
Let $\tilde{Z}=[{\tilde{z}}_{1},\mathrm{...},{\tilde{z}}_{m}]$. Then Σ with $W={({\tilde{Z}}^{T}\tilde{Z})}^{1}$ is a random correlation matrix with constraint matrix J.
 1.
Generate a Z matrix from standard normal distributions.
 2.
Generate a sequence of edge densities, e_{1},..., e_{ d }with 0 = e_{1} < ⋯ <e_{ d }= 1.
 3.
Generate a nested sequence of random constraint matrices J_{1},..., J_{ d }according to edge densities, e_{1},..., e_{ d }. Note the proportion of nonzero offdiagonal elements in J_{ i }is e_{ i }.
 4.
Given the Z matrix, generate Σ_{1},...,Σ_{ d }according constraint matrices J_{1},..., J_{ d }by Algorithm 1. Then Σ_{1} = I and the sequence of random correlation matrices has nested conditional independence structures.
Variance of correlations and the choice of M
when Z is generated from standard normal distributions [see Appendix 1 of Additional file 1]. Hence for large M, we expect average pairwise correlations are close to zero and the effect of dependence when M is large is almost ignorable.
Average variances of offdiagonal entries in (3) decrease very quickly to zero as M increases. Hence, in this paper, when m = 1000, we focus on FDR results for M = 1001 since this value illustrates the effects of dependence in the most unrestricted way. For large M  m, variances of pairwise correlations are close to zero and the effects are almost negligible. In Figure 6, we show the FDR results for such a case.
Simulation details
 1.
Find c_{ γ }satisfying FDR(c_{ γ }) = γ under independence assumption.
 2.
Generate random correlation matrices Σ_{1},..., Σ_{ d }from given structures in Algorithm 2.
 3.
For each Σ_{ j }, ${X}_{1},\mathrm{...},{X}_{{n}_{1}}$ ~ N_{ m }(μ_{ X }, Σ_{ j }) and ${Y}_{1},\mathrm{...},{Y}_{{n}_{2}}$ ~ N_{ m }(μ_{ Y }, Σ_{ j }).
 4.
Apply various multiple testing procedures to these data and compare the corresponding results of FDR, FNR and π_{0} estimates.
FDR procedures, π_{0}estimation methods and software used in the simulations
We briefly introduce the FDR implementations used in the simulation studies. Most of them are regularly used and all of them are developed in R software packages [32].

BenjaminiHochberg procedure (BH): Implemented FDR control by a linear stepup procedure [2]. From ordered observed pvalues p_{(1)},..., p_{(m)}, it finds maximal k such that ${p}_{(k)}\le \gamma \frac{k}{m}$ given significance level γ. It is known to control FDR at level $\gamma \frac{{m}_{0}}{m}$ under independence assumption of test statistics. π_{0} estimation procedure is not implemented, hence π_{0} is assumed to be 1. We use R package multtest for this procedure.

BenjaminiYekutieli procedure (BY): Benjamini et al. [3] extends the BH procedure to control FDR at level γ under arbitrary dependence conditions. It uses the linear stepup procedure, and it finds maximal k such that ${p}_{(k)}\le \gamma \frac{k}{m}{({\displaystyle {\sum}_{i=1}^{m}{i}^{1}})}^{1}$. We use R package multtest for this procedure.

Adaptive BenjaminiHochberg procedure (ABH): The ABH procedure improves the BH procedure by estimating m_{0}. Given significance level γ, the twostage ABH procedure first performs the linear stepup BH procedure to find r_{1}, the number of rejected hypotheses at level γ* = γ/(1 + γ). It estimates ${\widehat{m}}_{0}$ as m  r_{1} and then applies ${\gamma}^{\ast}\frac{m}{{\widehat{m}}_{0}}$ as a new significance level in the second step. Under the independence assumption of test statistics, ABH is known to control FDR at level γ[8]. We use R package FDRAME for this procedure.

Significance Analysis of Microarray (SAM): Based on [6], the SAM procedure is developed. For twoclass, unpaired data, it uses a tstatistic combined with a fudge factor which makes test statistics more stable when sample variance is very small. Using permutations and a userspecified cutoff, it produces asymmetric testing results. To apply the same significance level γ as other FDR procedures, we set median FDR level to be γ instead of using the userspecified cutoff. We use R package samr with internal permutation number 200.

Qvalue: Storey [18] proposes a new multiple testing criterion qvalue based on pFDR. pFDR is defined as the expected proportion of the number of false rejections over the number of rejections given the number of rejections is at least one. qvalue is the minimum pFDR where the statistic is declared significant. The estimates of qvalues are computed from a function of individual pvalues while preserving the order of pvalues. We use R package qvalue and choose the default option "smoother" as "pi0.method".

Resampling based FDR adjustments (RBH): Resampling based FDR estimation is based on the resampling distribution of pvalues under the complete null hypothesis. Basically, it is defined as E_{R(γ)*}[R(γ)*/(R(γ)* + $\widehat{s}$(γ))] where R*(γ) is the number of resamplingbased pvalues less than γ. Two estimators conditioned on $\widehat{s}$(γ) are proposed. The point RBH estimator is based on $\widehat{s}$(γ) = r(γ)  mγ and the upper limit RBH estimator is based on $\widehat{s}(\gamma )=r(\gamma ){r}_{\beta}^{\ast}(\gamma )$ where ${r}_{\beta}^{\ast}(\gamma )$ is 1  β quantile of R*(γ) [12]. We use R package FDRAME for this procedure.
ABH, SAM and Qvalue contain internal π_{0} estimation. Recently, another π_{0} estimation method is introduced by Langaas et al. [10]. Here, pvalues are modeled as f(p) = π_{0} + (1  π_{0})h(p) where h(p) is a convex decreasing density of false null hypotheses with h(1) = 0. In this setup, nonparametric maximum likelihood estimation is employed to compute estimate of π_{0}. For the case of nonconvexity of f, the authors advise to restrict the domain to the convex part of f, but this is not implemented yet. We use the convest function in the limma R packages in the default option.
Simulation program
We developed R code [32] for this simulation studies. The code can also be used in a supervised sense, using known gene sets. Please contact the authors for obtaining the R program.
Declarations
Acknowledgements
We thank the referees for their stimulating remarks on earlier versions of this paper.
Authors’ Affiliations
References
 Wille A, Zimmermann P, Vranova E, Furholz A, Laule O, Bleuler S, Hennig L, Prelic A, von Rohr P, Thiele L, Zitzler E, Gruissem W, Buhlmann P: Sparse graphical Gaussian modeling of the isoprenoid gene network in Arabidopsis thaliana. Genome Biol 2004, 5(11):R92. 10.1186/gb2004511r92PubMed CentralView ArticlePubMedGoogle Scholar
 Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Statist Soc Ser B 1995, 57: 289–300.Google Scholar
 Benjamini Y, Yekutieli D: The control of the false discovery rate in multiple testing under dependency. Ann Statist 2001, 29(4):1165–1188. 10.1214/aos/1013699998View ArticleGoogle Scholar
 Storey JD: The positive false discovery rate: a Bayesian interpretation and the q value. Ann Statist 2003, 31(6):2013–2035. 10.1214/aos/1074290335View ArticleGoogle Scholar
 Storey J, Tibshirani R: Estimating false discovery rates under dependence, with applications to DNA microarrays. Tech Rep 2001–12 Stanford University; 2001. [http://wwwstat.stanford.edu/reports/papers2001.html]Google Scholar
 Tusher V, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci 2001, 98: 5116–5121. 10.1073/pnas.091062498PubMed CentralView ArticlePubMedGoogle Scholar
 Dudoit S, Shaffer JP, Boldrick JC: Multiple hypothesis testing in microarray experiments. Statist Sci 2003, 18: 71–103. 10.1214/ss/1056397487View ArticleGoogle Scholar
 Benjamini Y, Krieger AM, Yekutieli D: Adaptive linear stepup procedures that control the false discovery rate. Biometrika 2006., 93(3):
 Black MA: A note on the adaptive control of false discovery rates. J R Stat Soc Ser B Stat Methodol 2004, 66(2):297–304. 10.1111/j.13697412.2003.05527.xView ArticleGoogle Scholar
 Langaas M, Lindqvist BH, Ferkingstad E: Estimating the proportion of true null hypotheses, with application to DNA microarray data. J R Stat Soc Ser B Stat Methodol 2005, 67(4):555–572. 10.1111/j.14679868.2005.00515.xView ArticleGoogle Scholar
 Westfall P, Young S: Resamplingbased multiple testing: examples and methods for pvalue adjustment. Wiley, New York; 1993.Google Scholar
 Yekutieli D, Benjamini Y: Resamplingbased false discovery rate controlling multiple test procedures for correlated test statistics. J Statist Plann Inference 1999, 82(1–2):171–196. 10.1016/S03783758(99)000415View ArticleGoogle Scholar
 Reiner A, Yekutieli D, Benjamini Y: Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics 2003, 19(3):368–375. [http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve\&db=pubmed\&dopt=Abstract\&list_uids=12584122] 10.1093/bioinformatics/btf877View ArticlePubMedGoogle Scholar
 Qiu X, Klebanov L, Yakovlev A: Correlation Between Gene Expression Levels and Limitations of the Empirical Bayes Methodology for Finding Differentially Expressed Genes. Statistical Applications in Genetics and Molecular Biology 2005., 4(34): Epub 2005 Nov 22.Google Scholar
 Korn EL, Troendle JF, McShane LM, Simon R: Controlling the number of false discoveries: application to highdimensional genomic data. J Statist Plann Inference 2004, 124(2):379–398. 10.1016/S03783758(03)002118View ArticleGoogle Scholar
 Jung SH, Jang W: How accurately can we control the FDR in analyzing microarray data? Bioinformatics 2006, 22(14):1730–1736. 10.1093/bioinformatics/btl161View ArticlePubMedGoogle Scholar
 Genovese C, Wasserman L: Operating characteristics and extensions of the false discovery rate procedure. J R Stat Soc Ser B Stat Methodol 2002, 64(3):499–517. 10.1111/14679868.00347View ArticleGoogle Scholar
 Storey JD: A direct approach to false discovery rates. J R Stat Soc Ser B Stat Methodol 2002, 64(3):479–498. 10.1111/14679868.00346View ArticleGoogle Scholar
 Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci USA 2003, 100(16):9440–9445. 10.1073/pnas.1530509100PubMed CentralView ArticlePubMedGoogle Scholar
 Efron B: Correlation and LargeScale Simultaneous Significance Testing.2006. [http://wwwstat.stanford.edu/~brad/papers/]Google Scholar
 Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: From the Cover: Gene set enrichment analysis: A knowledgebased approach for interpreting genomewide expression profiles. PNAS 2005, 102(43):15545–15550. [http://dx.doi.org/10.1073/pnas.0506580102] 10.1073/pnas.0506580102PubMed CentralView ArticlePubMedGoogle Scholar
 Schafer J, Strimmer K: An empirical Bayes approach to inferring largescale gene association networks. Bioinformatics 2005, 21(6):754–764. [http://dx.doi.org/10.1093/bioinformatics/bti062] 10.1093/bioinformatics/bti062View ArticlePubMedGoogle Scholar
 Dobra A, Hans C, Jones B, Nevins JR, Yao G, West M: Sparse graphical models for exploring gene expression data. J Multivariate Anal 2004, 90: 196–212. 10.1016/j.jmva.2004.02.009View ArticleGoogle Scholar
 Jones B, West M: Covariance decomposition in undirected Gaussian graphical models. Biometrika 2005, 92(4):779–786. 10.1093/biomet/92.4.779View ArticleGoogle Scholar
 Efron B, Tibshirani R: On Testing the Significance of Sets of Genes. Ann Appl Statist 2007, 1: 107–129. 10.1214/07AOAS101View ArticleGoogle Scholar
 Wright G, Simon R: A random variance model for detection of differential gene expression in small microarray experiments. Bioinformatics 2003, 19: 2448–55. 10.1093/bioinformatics/btg345View ArticlePubMedGoogle Scholar
 Owen AB: Variance of the number of false discoveries. J R Stat Soc Ser B Stat Methodol 2005, 67(3):411–426. 10.1111/j.14679868.2005.00509.xView ArticleGoogle Scholar
 Wagner GP: On the eigenvalue distribution of genetic and phenotypic dispersion matrices: evidence for a nonrandom organization of quantitative character variation. J Math Biol 1984, 21: 77–95.View ArticleGoogle Scholar
 Lauritzen SL: Graphical models, of Oxford Statistical Science Series. Volume 17. New York: The Clarendon Press Oxford University Press; 1996. [, Oxford Science Publications]Google Scholar
 Whittaker J: Graphical models in applied multivariate statistics. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics, Chichester: John Wiley & Sons Ltd; 1990.Google Scholar
 Marsaglia G, Olkin I: Generating correlation matrices. SIAM J Sci Statist Comput 1984, 5(2):470–475. 10.1137/0905034View ArticleGoogle Scholar
 Ihaka R, Gentleman R: R: A Language for Data Analysis and Graphics. Journal of Computational and Graphical Statistics 1996, 5(3):299–314. [http://www.amstat.org/publications/jcgs/] 10.2307/1390807Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.