 Research Article
 Open access
 Published:
Control procedures and estimators of the false discovery rate and their application in lowdimensional settings: an empirical investigation
BMC Bioinformatics volume 19, Article number: 78 (2018)
Abstract
Background
When many (up to millions) of statistical tests are conducted in discovery set analyses such as genomewide association studies (GWAS), approaches controlling familywise error rate (FWER) or false discovery rate (FDR) are required to reduce the number of false positive decisions. Some methods were specifically developed in the context of highdimensional settings and partially rely on the estimation of the proportion of true null hypotheses. However, these approaches are also applied in lowdimensional settings such as replication set analyses that might be restricted to a small number of specific hypotheses. The aim of this study was to compare different approaches in lowdimensional settings using (a) real data from the CKDGen Consortium and (b) a simulation study.
Results
In both application and simulation FWER approaches were less powerful compared to FDR control methods, whether a larger number of hypotheses were tested or not. Most powerful was the qvalue method. However, the specificity of this method to maintain true null hypotheses was especially decreased when the number of tested hypotheses was small. In this lowdimensional situation, estimation of the proportion of true null hypotheses was biased.
Conclusions
The results highlight the importance of a sizeable data set for a reliable estimation of the proportion of true null hypotheses. Consequently, methods relying on this estimation should only be applied in highdimensional settings. Furthermore, if the focus lies on testing of a small number of hypotheses such as in replication settings, FWER methods rather than FDR methods should be preferred to maintain high specificity.
Background
Advances in molecular biology and laboratory techniques allow for evaluating a multitude of different features in humans on a large scale to elucidate (patho)physiology and risk factors for a specific disease or its progression. In recent studies, up to millions of features are often assessed simultaneously in discovery set analyses such as in genomewide association studies (GWAS) where single nucleotide polymorphisms (SNPs) are evaluated with respect to a single trait or clinical outcome [1]. For reasons of practicability, the usual analysis procedure of such highdimensional data comprises statistical testing of each single feature separately with the outcome of interest [2].
Statistical testing aims to verify a hypothesis, which is either rejected or accepted based on the observed test statistic [3]. Depending on the decision, there are two possible mistakes that can occur: The null hypothesis might be erroneously rejected although it is true (false positive decision, type I error) or failed to reject although it is false (false negative decision, type II error). The type I error can be controlled by defining a significance threshold. For a single hypothesis, a commonly used threshold is α=0.05. However, when testing multiple hypotheses such as in GWAS, the application of a threshold like 0.05 across all tests will result in an unacceptable large number of false positive results. Consequently, other ways to control the type I error are required.
In general, there are different approaches: the control of the familywise error rate (FWER) and the control or the estimation of the false discovery rate (FDR) [4]. FWER methods such as the wellknown Bonferroni correction [5] were already proposed when the number of tested hypotheses was not as large as, for example, in GWAS nowadays. Although often applied, these methods are thought to be too conservative in a highdimensional setting. Alternatively, FDR methods that are less conservative and partially developed in the context of highdimensional data can be used. In addition, there are approaches to estimate a significance measure for each individual hypothesis, such as the local false discovery rate (LFDR) [6] and the qvalue [7].
FDR methods are also used quite frequently nowadays and not only in highdimensional settings but also in situations where the number of assessed features is small such as in a replication set analysis restricted to the significant hypotheses of the discovery set analysis. For a small number of features, however, there are limited data on the performance of FDR methods. The aim of this study was thus to assess FDR methods in lowdimensional data and to compare them to classic FWER methods. For this purpose, we used real data obtained from the CKDGen Consortium [8] to illustrate the different control methods. Moreover, we conducted a simulation study to evaluate different control methods in different settings.
Methods
Control methods
In order to describe different error control and estimation methods, we adopted the notation of Benjamini and Hochberg [9] on test decisions (Table 1). Assuming m hypotheses H_{1}, …, H_{m} were tested leading to the observation of the respective m pvalues p_{1}, …, p_{m}. If the truth would be known, type I errors are described by V and type II errors by T. However, only m and the total number of rejections, R, are observable in practice. The overall significance threshold is called α.
Methods controlling the familywise error rate (FWER)
FWER is defined as the probability of making at least one false positive decision: FWER = Pr(V > 0). The error rate can be controlled by a fixed threshold α. In the following, four well known methods are considered (Table 2a):
The simplest and likely most often applied control method of the FWER is the Bonferroni correction [10]. It compares each individual pvalue p_{1}, …, p_{m} with the fixed threshold\( \frac{\upalpha}{\mathrm{m}} \). Pvalues that are smaller than the threshold lead to the rejection of the respective null hypothesis. The Bonferroni correction guarantees the control of the FWER at level α in a strong sense, which means that the control is ensured for every proportion of true null hypotheses. Bonferroni correction does not demand independent pvalues and hence can be applied to any dependency structures. Nevertheless, Bonferroni can be conservative; true alternatives might therefore be missed.
To reduce the number of missed true alternatives, approaches to adjust Bonferroni correction were proposed that use the number of independent tests (also: effective number) instead of the actual number of conducted tests (e.g. Li et al. [11]). Therefore, these approaches gain in power over the traditional Bonferroni correction. In the specific context of GWAS, for example, an adjusted Bonferroni correction frequently applied was proposed by Pe’er et al. [12] that accounts for correlation between SNPs due to linkage disequilibrium (LD) by estimating the number of independent genomewide loci (n = 1,000,000 in individuals of European ancestry). Instead of using the much larger number of all SNPs tested for association (often several millions), the overall significance threshold such as α=0.05 is divided by the number of independent SNPs to define an adjusted significance threshold. For GWAS on Europeans, for example, the significance threshold becomes \( \frac{0.05}{\mathrm{1,000,000}}=5\times {10}^{8} \). Similarly, the number of independent tests in the field of metabolomics can be estimated with help of principle component analysis to reduce the number of all tests used in Bonferroni correction (e.g. Grams et al. [13]).
The other three FWER control methods considered below are sequential methods for which pvalues need to be ranked in increasing order: p_{(1)} ≤ … ≤ p_{(m)}.
Holm’s stepdown procedure [10] rejects at least as many hypotheses as Bonferroni correction does. The gain in power of Holm’s procedure by defining more features significant is larger with larger number of alternative hypotheses. Like the Bonferroni correction, Holm’s procedure has no restrictions with respect to the dependency structure of pvalues.
Hochberg’s stepup procedure [14] and also Hommel’s procedure [15] make use of the assumption that the pvalues under the true null hypotheses hold a positive regression dependency. Positive dependency structure assumes the probability of a pvalue belonging to the null hypothesis to be increasing with increasing pvalue. In situations of a positive dependency structure, Hochberg’s procedure is more powerful than Holm’s [4]. Hommel’s procedure, however, is the most powerful FWER control procedure of the previously mentioned methods when the assumption holds since it rejects at least as many hypotheses as Hochberg does. One criticism of the method lies in the higher computational load.
Methods controlling the false discovery rate (FDR)
In contrast to FWER, the false discovery rate (FDR) represents the proportion of false positives. This error rate is defined as following: \( FDR=E\left[\frac{V}{R}R>0\right]\Pr \left(R>0\right). \) FDR can be controlled at a fixed significance threshold as well. Furthermore, Benjamini and Hochberg [9] proved that every FWER control method controls the FDR likewise. The three most common FDR control methods that also require ordered pvalues are considered below (Table 2b):
BenjaminiHochberg’s linear stepup procedure [9] controls the FDR at level α assuming positive dependent pvalues (see description above) under the true null hypotheses such as Hommel’s and Hochberg’s FWER procedures. It shows greater power than any of the above mentioned FWER methods.
The twostage linear stepup procedure [16] is an adapted procedure of BenjaminiHochberg’s that takes the estimation of the proportion of the true null hypotheses, π_{0}, into account. The gain in power of the twostage procedure compared to the classical BenjaminiHochberg’s linear stepup procedure is dependent on the proportion of true null hypotheses (π_{0}) [4]. For π_{0} close to 1, the adapted version has low power. The adaptive approach has been proven for independent pvalues only.
Finally, BenjaminiYekutieli’s linear stepup procedure [17] has no restrictions on the dependency structure of pvalues at all. It is more conservative compared to the BenjaminiHochberg’s linear stepup procedure [4] and the twostage linear stepup procedure [16].
Methods estimating the false discovery rate (FDR)
Recent approaches do not control the FDR in the traditional sense, but rather estimate the proportion of false discoveries. In order to estimate the FDR, the estimation of the proportion of the true null hypotheses, π_{0}, is conducted first which can lead to a gain in power compared to the classic FWER and FDR control methods. Two common FDR estimation methods are described in the following:
Storey’s qvalue method [7] uses a Bayesian approach to estimate the socalled positive false discovery rate (pFDR), a modified definition of the false discovery rate which assumes at least one rejection: \( pFDR=E\left[\frac{V}{R}R>0\right] \). The approach is based on the idea of estimating the pFDR for a particular rejection region, γ, to achieve a control of the pFDR. In order to determine a rejection region, the qvalue was introduced as the pFDR analogue of the pvalue. The qvalue provides an error measure for each observed pvalue. It denotes the smallest pFDR that can occur when calling that particular pvalue significant: \( q(p)=\underset{\left\{\gamma \ge p\right\}}{\min } pFDR\left(\gamma \right) \). The approach assumes independent, respectively “weak dependent” pvalues, whose dependency effect becomes negligible for a large number of pvalues [18]. The method provides an improvement in power compared to the classic BenjaminiHochberg’s linear stepup procedure due to its estimation of π_{0} [7].
Likewise, Strimmer [19] proposed an alternative method to estimate qvalues based on pFDR (Strimmer’s qvalue method). In addition, the method provides estimates of the socalled local false discovery rate (LFDR, Strimmer’s LFDR approach) that again present individual significance measures such as the qvalues for each pvalue. It describes the probability that a pvalue leads to a false positive decision given the observed data information. Estimations are based on a Bayesian approach using a modified Grenander density estimator [19].
Software implementation
R packages are available for all described control methods via CRAN [20] or Bioconductor [21]. Specifically, we used the packages multtest [22], qvalue [23] (Bioconductor), mutoss [24] and fdrtool [25] (CRAN) in our study. We applied the methods using default options of the packages. However, Storey’s qvalue application displayed an error whenever the estimated proportion of true null hypotheses (π_{0}) was close to zero, which occurred when all pvalues happened to be (very) small. Therefore, we adjusted the range of input pvalues (“lambda”) in a stepwise manner until the application allowed the estimation of π_{0}. Further details on our Rcode and the stepwise algorithm can be obtained directly from the authors. Statistical significance using either FWER, FDR controlling or FDR estimation methods such as the qvalue methods or LFDR, was defined as a cutoff of 0.05.
Data example
For illustration of the different control methods, we obtained data from the CKDGen Consortium [8]. The aim of this project was to identify genetic variants associated with estimated glomerular filtration rate (eGFR), a measure for kidney function, and chronic kidney disease (CKD). Altogether, 48 study groups provided genomewide summary statistics (GWAS results) from 50 study populations for SNP associations with eGFR based on serum creatinine (eGFRcrea) (2 study groups provided GWAS results for 2 subpopulations separately). The discovery metaanalysis of all GWAS was carried out using an inverse varianceweighted fixed effect model and incorporated data from 133,413 individuals of European ancestry. Genomic control had been applied before and also after metaanalysis to reduce inflation and thus limit the possibility of false positive results. In the metaanalysis, 29 previously identified loci and 34 independent novel loci (pvalue < 10^{−6}) were detected. Novel loci were then verified in an independent replication set (14 studies; N = 42,166). For 16 of the 34 novel loci, replication analysis showed directionconsistent results with pvalue combining discovery and replication < 5×10^{−8} (see Table 1 in Pattaro et al. [8]). For all but 1 SNP (rs6795744), the reported qvalues in the replication study were < 0.05.
The results of the discovery metaanalyses for different traits including eGFRcrea (approximately 2.2 million SNPs) are publicly available [26]. Moreover, we obtained the summary statistics from GWAS results for eGFRcrea of all studies contributing to the discovery (48 studies, 50 result files) for our project. For the illustration of the different control methods in both discovery (highdimensional) setting and replication (lowdimensional) setting, we split the 50 study contributions into two sets taking into account general study characteristics (populationbased study versus diseased cohort) and imputation reference (HapMap versus 1000 Genomes [27]). By conditioning on the presence of at least one study from each of the 4 categories in either setting and on a sample size ratio of 2:1, study contributions were randomly assigned to discovery set or replication set. The final discovery set contained 35 studies with 90,565 individuals (67.9%) and the replication set 15 studies with 42,848 individuals (32.1%).
Based on the same set of SNPs as in the publicly available data set, our discovery set was processed similarly to the original analysis [8] by using an inverse varianceweighted fixed effect model and genomic control before and after that step. For simplicity reasons we considered twosided pvalues in the discovery and replication set analysis. To select independently associated SNPs, SNPs were clustered based on LD pruning using the clump command of Plink v1.90b2 (r^{2}: 0.2, window: 1000 kb, significance threshold for index SNP: 10^{−6}) [28] and data of 1000 Genomes project (phase 3) as the LD reference. SNPs with the lowest pvalue within a specific region were considered as index SNPs. Few SNPs that were either not present in the reference or triallelic were excluded at this point. Using the prepared discovery data, the various FDR and FWER methods were then applied exploratively.
Similar to the published analysis by the CKDGen Consortium (Pattaro et al. [8]), independent index SNPs with pvalue < 10^{−6} were selected from the discovery set to be followed up in the replication set. The various control methods were subsequently applied to the results of the metaanalysis (same model as before but without genomic control) in the replication set to identify significant findings.
Simulation study
In order to assess power and specificity of the described FWER and FDR methods in detail, we conducted a simulation study with varying settings, with special emphasis on situations with a smaller number of tested features. The Rcode of the simulation study can be requested from the author.
For this purpose, test statistics for varying numbers of features (N = 4, 8, 16, 32, 64, 1000) were simulated to generate data sets. Test statistics for single features were simulated by drawing from \( \mathcal{N}\left(\upbeta, 1\right) \) with either β = 0 (null hypothesis) or β ∈ {1.0, 2.5} (alternative or nonnull hypothesis). Depending on the number of features in a given data set, the proportion of the true null hypotheses π_{0} ∈ {25%, 50%, 75%, 100%} was apriori defined. Each scenario defined by the different combinations of parameters was repeated 100 times. In preparation of the subsequent application of control methods, simulated test statistics were transformed into twosided pvalues.
The power of each approach was defined as proportion of correctly rejected hypotheses among all true alternative hypotheses whereas the specificity was defined as the proportion of correctly maintained hypotheses among all true null hypotheses. Furthermore, we evaluated the estimation results of the proportion of true null hypotheses of Storey’s and Strimmer’s qvalue methods within the simulation study.
Results
Data example
For the purpose of illustration, the 50 GWAS summary statistics provided by contributing study groups included in the original CKDGen discovery metaanalysis of eGFRcrea were split into 2 sets resembling a highdimensional discovery set (35 studies, 90,565 individuals) and a lowdimensional replication set (15 studies, 42,848 individuals). Details on the two sets are provided in Additional file 1 and Additional file 2.
Similar to the published analysis by the CKDGen Consortium (Pattaro et al. [8]), the discovery set was processed to select independent variants to be moved forward to a lowdimensional replication analysis. Based on pvalue threshold < 10^{−6} followed by LD pruning, 57 index SNPs from different genomic regions were selected from the discovery set. The replication analysis of the 57 selected index SNPs showed directionconsistent effect estimates for 56 SNPs.
Subsequently, the various control methods were applied to the metaanalysis results of the replication set to identify significant findings. Figure 1 presents the number of significant results of the different control procedures. Since the FWER methods Holm, Hochberg, and Hommel declared the same pvalues as significant, we decided to display the performance of Hommel’s approach only.
In contrast to FDR methods, FWER methods rejected the smallest number of hypotheses with Bonferroni being least powerful. Among the FDR methods, FDR estimating methods by Strimmer and Storey provided more power. Storey’s qvalue method rejected all hypotheses and it was the only approach which declared the directioninconsistent SNP as significant.
As expected, the applied FWER and FDR methods showed a monotone subset behavior related to rejected hypotheses, i.e. that the pvalues declared significant from a more conservative approach were always included in the set of pvalues declared significant from a less conservative method. This is a consequence of the methods’ property that – if a specific pvalue is declared significant – all other smaller pvalues are also declared significant.
Simulation study
Power and specificity of control methods
In a setting where the proportion of true null hypotheses, π_{0}, is 100%, Storey’s and Strimmer’s qvalue methods most often falsely rejected true null hypotheses when the number of tested hypotheses N is small (≤32), while for larger numbers of tested hypotheses and/or other methods the number of erroneous decisions mostly did not exceed 5 (Fig. 2a). BenjaminiYekutieli’s procedure and Strimmer’s LFDR approach performed best with 0 to 3 repetitions of falsely rejected hypotheses for all N. As a remark, Strimmer’s LFDR approach could not provide any results for N = 4. Specificity of methods to correctly maintain hypotheses is similarly good on average; only Storey’s qvalue method showed decreased specificity when the number of tested hypotheses was small.
When the proportion of true null hypotheses was < 100%, the power to correctly reject hypotheses was dependent on π_{0}, the effect size (β) and N. On average, it increased with decreasing π_{0}, increasing β and decreasing N overall. Figure 2b, c and d exemplarily show the average power for varying π_{0} and β_{1} = 2.5 under the alternative hypothesis, in dependence on N. Further figures for an effect size of β_{1}= 1 can be found in the Additional file 3.
As expected, FDR methods, especially the two qvalues methods, were more powerful than FWER methods. In terms of specificity, Storey’s qvalue method followed by Strimmer’s qvalue method showed lower specificity results for small N (≤16) than other methods. We observed similarity in specificities among the other methods. Again, Strimmer’s LFDR approach did not provide results when number of hypotheses were < 8 (Fig. 2b) or < 16 (Fig. 2c and d).
Estimation of proportion of true null hypotheses
LFDR and qvalue methods rely on the estimation of π_{0}. Figure 3 displays its estimations using Storey’s and Strimmer’s qvalue approaches for varying π_{0} and β_{1} = 2.5 under the alternative hypotheses (if present), while remaining figures are in the Additional file 4.
For small N, both estimations showed large variability within repetitions. Throughout all scenarios, Storey’s method showed greater estimation ranges of π_{0} compared to Strimmer’s qvalue approach. Moreover, estimation of π_{0} was often biased. Only when β_{1} = 2.5 and N was larger than 32, bias essentially disappeared. When β_{1} = 1, however, π_{0} was overestimated on average, even for larger N.
Discussion
FDR estimation methods such as Strimmer’s LFDR or Storey’s qvalue method have been mainly developed for highdimensional settings, of which discovery GWAS is one. They provide a less conservative approach compared to standard FWER and FDR control methods. The LFDR as well as the qvalue methods are Bayesian approaches which take the whole information on the data itself into account when estimating the proportion of true null hypotheses, π_{0}. Consequently, for the purposes of FDR estimation, a highdimensional setting is a great advantage allowing reasonable estimation of π_{0}. Though controversial, the qvalue methods as well as other FDR methods have been used in lowdimensional settings as well, such as in the analysis of replication data sets consisting of only limited number of SNPs. We thus aimed to compare various FWER and FDR methods including the qvalue method in order to assess their power and specificity in lowdimensional settings using simulated data and application to real data.
The analysis of our example data from the CKDGen Consortium [8] showed that the FDR estimation methods by Strimmer and Storey declared the largest number of SNPs significant in the lowdimensional replication analysis of 57 SNPs, followed by the FDR control methods of BenjaminiHochberg and BenjaminiYekutieli. As expected, the FWER control methods showed the lowest power by declaring the least number of pvalues significant. Of note, Storey’s qvalue method was the only approach which declared the single SNP (rs10201691) that showed directioninconsistent results between the discovery and replication analyses as significant in the replication analysis.
To deepen the understanding, we conducted a simulation study to systematically assess different scenarios. As one result, the differences between the methods that were seen in the application could be confirmed. For example, Storey’s qvalue method showed the highest power especially for a small number of hypotheses. At the same time, however, the specificity results for Storey’s method were lowest when number of tested hypotheses was small. In the presence of alternative hypotheses (π_{0} < 100%), we also observed that the FDR methods, BenjaminiHochberg and the twostage approach, − although less powerful than both qvalue methods – were more powerful than the FWER control methods of Bonferroni and Hommel, but of similar specificity.
Since both qvalue methods as well as LFDR rely on the estimation of π_{0}, we also investigated its estimation accuracy using the different approaches. For both methods, the estimate of π_{0} was often biased, especially when numbers of tested hypotheses were small. In addition, Storey’s qvalue method showed much higher variance compared to Strimmer’s approach. In summary, the qvalue methods rejected in general the largest number of hypotheses which is especially of advantage if researchers wish to obtain a greater pool of significant features to be followedup in subsequent studies, at the expense of specificity. However, their application should be restricted to highdimensional settings.
The gain in power for both qvalue methods, however, was not observed for LFDR in the simulation study. Strimmer reported the gain in power of the qvalue method compared to the LFDR as well and explained it as the tendency of qvalues being smaller or equal compared to LFDR for a given set of pvalues [19]. In the context of gene expression, Lai [29] mentioned a tendency of the qvalue to underestimate the true FDR leading to a larger number of low qvalues especially when the proportion of differentially expressed genes is small or the overall differential expression signal is weak. We also observed an underestimation in our simulation study, especially for a smaller number of pvalues. To overcome this issue, Lai [29] suggested a conservative adjustment of the estimation of the proportion of true null hypotheses, the pvalues or the number of identified genes.
Moreover, when applying qvalue methods or LFDR, correct interpretation of these estimates is requested that is different for the qvalues and for LFDR. Strimmer [19] highlighted the easier interpretation of the LFDR compared to the qvalue since the LFDR provides point estimates for the proportion of false discoveries for individual hypotheses whereas the qvalue of a pvalue is the expected proportion of false positives when calling that feature significant [18]. In any case, when applying FDR estimation methods, there is a critical need for a sizeable data set [18, 19]. Storey and Tibshirani [18] described their qvalue method as a more explorative tool compared to FWER methods and therefore as wellperforming procedure in highdimensional data. A more recent FDR estimation approach by Stephens [30] provides an alternative to the LFDR, the so called local false sign rate. This empirical Bayes approach describes the probability of making an error in the sign of a certain variant if forced to declare it either as true or false discovery. Simulation studies showed smaller and more accurate estimation of π_{0} by Stephens’ approach compared to Storey’s qvalue method leading to more significant discoveries [30]. However, small sample sizes represent a challenge for this FDR estimation approach as well.
Another observation of our simulation study worth mentioning was that the FDR method by BenjaminiYekutieli for arbitrary dependencies, and thus assumed to be more conservative than the BenjaminiHochberg method, was not only outperformed by this method in terms of power in our application data and simulation, but also less powerful than FWER control methods in some scenarios of our simulation. The latter had already been observed, especially if the expected number of alternative hypotheses is very small [4]. Since BenjaminiHochberg’s approach controls the FDR at level π_{0}α, adaptive FDR control methods such as the twostage approach were developed to control the FDR directly at level α by taking estimated π_{0} into account and thereby gaining power. Especially if π_{0} is substantially smaller than 1, the adaptive approaches might outperform BenjaminiHochberg’s procedure [4].
Before concluding the discussion on results, some limitations of this study warrant mentioning: Although it was important for us to illustrate the effect of the different control methods on the results in real data, observed differences may not be transferrable to every other study setting in general. To overcome this limitation, we conducted a simulation study. Still, the simulation study has limitations of its own: We used a simplified approach to generate data by simulating test statistics rather than analytical data sets to which control methods would have been applied after analysis. Furthermore, we explored a limited set of scenarios and did not consider dependency structures but evaluated pvalues that were derived from independently simulated test statistics. Hence, additional work could add to the current understanding.
In the face of all the different control methods, it is clear that the decision on what method is actually applied in a given setting should be made not only before the analysis is conducted but also on reasonable ground. Among others, aspects to consider include: (a) the amount of tests to be conducted, (b) the general aim of testing, (c) what is known or can be assumed about dependency structure of pvalues under the true null hypothesis and (d) what is the assumed proportion of null hypotheses.
If the general aim of the analysis lies on the specific testing of individual hypotheses, FWER control methods should be preferred to FDR control or estimation methods because they provide higher specificity by correctly maintaining true null hypotheses. Within FWER control methods, the power might differ slightly and is, especially, in dependence of given pvalue structure. If a positive structure can be assumed, Hochberg’s or Hommel’s procedures are preferable to gain power. The computational burden that comes along with Hommel’s procedure should not be a true issue nowadays. Goeman and Solari [4] especially expected a gain in power of Hochberg’s and Hommel’s compared to Bonferroni’s and Holm’s methods if the proportion of alternative hypotheses is rather large. We, however, observed only a rather small gain in power in our simulation study that might be induced by the simulation of independent test statistics.
If researchers, however, wish to identify a promising set of hypotheses for followup rather than specific testing of single hypotheses with high specificity, we agree with Goeman and Solari [4] who recommended the use of FDR control methods. To reach highest power, one may even apply the FDR estimating method of qvalues, when the number of tests is reasonably large.
Conclusions
In summary, our findings highlight the importance of a larger data set for the application of FDR estimation methods in order to guarantee reliable estimation of the proportion of true null hypotheses. The choice of control method mainly depends on the specific setting and the aims of an analysis. For example, when high specificity in testing of a limited number of hypotheses as in a replication study is desired, we recommend to utilize FWER methods rather than FDR methods.
Abbreviations
 FDR:

False discovery rate
 FWER:

Familywise error rate
 GWAS:

Genomewide association study
 LD:

Linkage disequilibrium
 LFDR:

Local false discovery rate
 pFDR:

Positive false discovery rate
 SNP:

Single nucleotide polymorphism
References
Balding DJ. A tutorial on statistical methods for population association studies. Nat Rev Genet. 2006;7(10):781–91.
Zeng P, Zhao Y, Qian C, Zhang L, Zhang R, Gou J, Liu J, Liu L, Chen F. Statistical analysis for genomewide association study. J Biomed Res. 2015;29(4):285–97.
Rothman KJ, Greenland S, Lash TL. Modern Epidemiology. 3rd ed. Philadelphia: Lippincott Williams & Wilkins; 2008. p. 151–6.
Goeman JJ, Solari A. Multiple hypothesis testing in genomics. Stat Med. 2014;33(11):1946–78.
Bland JM, Altman DG. Multiple significance tests: the Bonferroni method. BMJ. 1995;310(6973):170.
Efron B. Empirical Bayes analysis of a microarray experiment. J Am Stat Assoc. 2001;96(456):1151–60.
Storey JD. A direct approach to false discovery rates. J R Stat Soc. 2002;64(3):479–98.
Pattaro C, Teumer A, Gorski M, Chu AY, Li M, Mijatovic V, Garnaas M, Tin A, Sorice R, Li Y, et al. Genetic associations at 53 loci highlight cell types and biological pathways relevant for kidney function. Nat Commun. 2016;7:10023.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995;57(1):289–300.
Holm SA. Simple sequentially Rejective multiple test procedure. Scand J Stat. 1979;6(2):65–70.
Li J, Ji L. Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. Heredity. 2005;95(3):221–7.
Pe'er I, Yelensky R, Altshuler D, Daly MJ. Estimation of the multiple testing burden for genomewide association studies of nearly all common variants. Genet Epidemiol. 2008;32(4):381–5.
Grams ME, Tin A, Rebholz CM, Shafi T, Kottgen A, Perrone RD, Sarnak MJ, Inker LA, Levey AS, Coresh J. Metabolomic alterations associated with cause of CKD. Clin J Am Soc Nephrol. 2017;12(11):1787–94.
Hochberg Y. A sharper Bonferroni procedure for multiple tests of significance. Biometrika. 1988;75(4):800–2.
Hommel G. A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika. 1988;75(2):383–6.
Benjamini Y, Krieger AM, Yekutieli D. Adaptive linear stepup procedures that control the false discovery rate. Biometrika. 2006;93(3):491–507.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001;29(4):1165–88.
Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci U S A. 2003;100(16):9440–5.
Strimmer K. A unified approach to false discovery rate estimation. BMC Bioinform. 2008;9:303.
The Comprehensive R Archive Network. https://cran.rproject.org/. Accessed 17 Oct 2017.
Bioconductor. https://www.bioconductor.org/. Accessed 17 Oct 2017.
Pollard KS, Dudoit S, van der Laan MJ. Multiple testing procedures: the multtest package and applications to genomics. In: Gentleman R., Carey VJ, Huber W, Irizarry RA, Dudoit S, editors. Bioinformatics and Computational Biology Solutions Using R and Bioconductor: Springer, New York; 2005.
Storey JD, Bass AJ, Dabney A, Robinson D. qvalue: Qvalue estimation for false discovery rate control. 2015. R package version 2.10.0. http://github.com/jdstorey/qvalue. Accessed 17 Oct 2017.
MuToss Coding Team (Berlin 2010), Blanchard G, Dickhaus T, Hack N, Konietschke F, Rohmeyer K, Rosenblatt J, Scheer M, Werft W. mutoss: Unified Multiple Testing Procedures. 0.1–10 edn; 2015; The Mutoss package and accompanying mutossGUI package are designed to ease the application and comparison of multiple hypothesis testing procedures. https://CRAN.Rproject.org/package=mutoss. Accessed 17 Oct 2017.
Klaus B, Strimmer S. fdrtool: Estimation of (Local) False Discovery Rates and Higher Criticism. 2015: Estimates both tail areabased false discovery rates (Fdr) as well as local false discovery rates (fdr) for a variety of null models (pvalues, zscores, correlation coefficients, tscores). The proportion of null values and the parameters of the null distribution are adaptively estimated from the data. In addition, the package contains functions for nonparametric density estimation (Grenander estimator), for monotone regression (isotonic regression and antitonic regression with weights), for computing the greatest convex minorant (GCM) and the least concave majorant (LCM), for the halfnormal and correlation distributions, and for computing empirical higher criticism (HC) scores and the corresponding decision threshold. https://CRAN.Rproject.org/package=fdrtool. Accessed 17 Oct 2017.
CKDGen Consortium. MetaAnalysis Data. http://ckdgen.imbi.unifreiburg.de/. Accessed 17 Oct 2017.
Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, Korbel JO, Marchini JL, McCarthy S, McVean GA, Abecasis GRA. Global reference for human genetic variation. Nature. 2015;526(7571):68–74.
Purcell S, Neale B, ToddBrown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al. PLINK: a tool set for wholegenome association and populationbased linkage analyses. Am J Hum Genet. 2007;81(3):559–75.
Lai YA. Statistical method for the conservative adjustment of false discovery rate (qvalue). BMC Bioinformatics. 2017;18(Suppl 3):69.
Stephens M. False discovery rates: a new deal. Biostatistics (Oxford, England). 2017;18(2):275–94.
Acknowledgments
RB worked on this topic in the context of her master thesis at the Faculty of Mathematics and Physics at the University of Freiburg.
Funding
We acknowledge financial support by Deutsche Forschungsgemeinschaft and RuprechtKarlsUniversität Heidelberg within the funding programme Open Access Publishing.
Availability of data and materials
Summary statistics published by the CKDGen Consortium (Pattaro et al. (2016)) are publicly available at [26]. GWAS results from single studies contributing to discovery analysis of that project were directly requested from the CKDGen group and were released for this specific purpose.
Any further details and Rcode can be requested from the corresponding author (brinster@imbi.uniheidelberg.de).
Author information
Authors and Affiliations
Consortia
Contributions
Design: MS, PS. Data acquisition/generation, Analysis: RB, PS. Drafting manuscript: RB, PS. Revision and final approval of manuscript: RB, PS, MS, AK, BOT.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1
CKDGen study contributions assigned to discovery set for illustration of procedures controlling type I error. (DOCX 17 kb)
Additional file 2
CKDGen study contributions assigned to replication set for illustration of procedures controlling type I error. (DOCX 18 kb)
Additional file 3
Simulation – Average power and specificity for β_{1} = 1 and π_{0} = 75% (a), 50% (b), 25% (c). Applied procedures controlling the type I error: Bonferroni correction, Hommel’s procedure, BenjaminiHochberg’s procedure, Twostage procedure, BenjaminiYekutieli’s procedure, Storey’s qvalue method, Strimmer’s qvalue method, Strimmer’s LFDR method. Power is defined as the proportion of correctly rejected hypotheses and specificity as the proportion of correctly maintained hypotheses. Both proportions potentially range from 0 to 1. Simulations for each scenario were repeated 100 times. (PNG 457 kb)
Additional file 4
Simulation – Observed estimations of π_{0} for Storey’s (qv) and Strimmer’s qvalue methods (fdr)for β_{1} = 1 and π_{0} = 75% (a), 50% (b), 25% (c). (PNG 209 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Brinster, R., Köttgen, A., Tayo, B.O. et al. Control procedures and estimators of the false discovery rate and their application in lowdimensional settings: an empirical investigation. BMC Bioinformatics 19, 78 (2018). https://doi.org/10.1186/s128590182081x
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s128590182081x