A comparative review of estimates of the proportion unchanged genes and the false discovery rate

The proportion unchanged

In the two-component model for the distribution of the test statistic the mixing parameter p₀, which represents the proportion unchanged genes, is not estimable without strong distributional assumptions, see [1]. Assuming this model the probability density function (pdf) f ^t of a test statistic t may be written as the weighted sum of the null distribution pdf and the alternative distribution pdf

If, on the other hand, we know the value of p₀ we can estimate e.g. through a bootstrap procedure as described in [1], and thus obtain also .

and

One may choose a p-value threshold α _min , which minimises the amount of errors FP(α) + FN(α). Alternatively, one may want to fine-tune the test statistic such that it will minimise the errors at a given threshold. Or, one may try to do both, as suggested in [5], see also [6].

Earlier research providing estimates of p₀ include [1, 3, 7–13]. Articles that compare methods for estimating p₀ and FDR include [12, 14]. I will focus on the FDR as the main use of p₀.

In this article, the formulation of the theory is in terms of p-values rather than in terms of test statistics. Two basic assumptions are made concerning their distribution. First, it is assumed that test statistics corresponding to true null hypotheses will generate p-values that follow a uniform distribution on the unit interval, e.g. [15]. Thus, under the null distribution, the probability that a p-value falls below some cutpoint α equals α. Second, p-values are, unless stated otherwise, assumed to be independent. Empirical investigations will assess the effects of deviations from the second assumption.

The use of p-values means lumping up- and downregulated genes together. However, one may look separately at the two tails of the distribution of the test statistic to assess differential expression corresponding to up- and downregulation.

This article will not concern how p-values are calculated, but rather how they are used to calculate estimates of p₀ and FDR, and draw conclusions based on this evidence. It is assumed that p-values capture the essence of the research problem. Neither does the article treat the choice of an optimal test statistic. For illustration the t-test will be used repeatedly without regard to whether there are better methods or not. By the t-test we mean the unequal variance t-test: for sample means mean₁ and mean₂, sample variances , , and sample sizes n₁ and n₂. We apply the t-test to simulated normally distributed data and a permutation t-test to real data, where normality may be uncertain. All calculations were performed in R. The methods presented are available within packages for the free statistical software R [16, 17] and take a vector of p-values as input and output an estimate of p₀ and of FDR. Emphasis lies on methods available within R packages downloadable from CRAN [18] or Bioconductor [19]. Inevitably any review will exclude interesting work, but time and space limitations will not permit an all comprehensive review. The new and highly interesting concept of a local false discovery rate (LFDR) [1] only receives a cursory treatment.

This article builds on and finds motivation from the experience of the analysis of microarrays, which typically assay the expression of 10,000 or more genes. However, the methods presented apply equally well to other high dimensional technologies, such as fMRI or Mass Spectrometry.

False discovery rate

In the analysis of microarray experiments, the traditional multiple test procedures are often considered too stringent, e.g. [20] and [3]. In the last decade alternatives based on the concept of an FDR have emerged. For more details consult e.g. [2, 3, 9, 11, 21, 22]. There are different definitions proposed, but loosely speaking one would want to measure the proportion of false positive genes among those selected or significant. Loosely put the FDR may be interpreted as the proportion of false positives among those genes judged significantly regulated. Equation (2) is the FDR estimate presented in [3].

Denote by E[X] the expectation (or true mean) of any random variable X. With V the number of false positives given a certain cut-off and R the number of rejected null hypotheses, one may define the FDR as the expectation of the ratio of these quantities, or

,

where care is taken to avoid division by zero.

In [10] and [11] the FDR is estimated as the ratio of the expected proportion of false positives given the cut-off to the expected proportion selected. Viewed as a function of a cut-off α, such that genes g_i with p _i less than α are judged significant in terms of p-values, following the continuous cumulative distribution function (cdf) F, the FDR estimate is

which is nearly equal to (2) with the exception that the P_(L)has been replaced by its upper bound, and the step-wise empirical distribution by a smooth version, either a parametric model or a smoothed version of the empirical distribution. Thus, the FDR is now a continuous function instead of piece-wise continuous with jumps at each observed p-value. This ratio of expected values tries to estimate the expectation of a ratio: In general such an approach will give an overestimation, but in practice this will have little effect, see the Additional file.

The related concept of the positive FDR, pFDR = E[V/R|R > 0], the expectation conditional on at least one rejection, appears in [2]. Other forms of FDR have been proposed such as the conditional FDR [2], cFDR, defined as the expected proportion of false positives conditional on the event that R = r rejected have been observed : cFDR(r) = E[V|R = r]/r. This would answer to the question "What proportion of false positives may I expect in my top list of r genes?". Under independence and identical distribution in a Bayesian setting it is proved in [23], that pFDR, cFDR and the marginal FDR, mFDR = E[V]/E[R] [2], all coincide with p₀α/F(α) at the cutpoint α, cf. (6).

Instead of p-values it has been suggested in to calculate q-values that verbally have the following meaning for an individual gene [2, 9]:

The q-value for a particular gene is the minimum false discovery rate that can be attained when calling all genes up through that one significant [9].

These q-values can be used to determine a cut-off similar to the classic 5% cut-off for univariate tests developed in statistics long ago. But in many applications one should not be too rigid about any particular value, since the emphasis often is on discovery rather than hypothesis testing: we generate hypotheses worthy of further investigation. Thus the balance between false positives and false negatives will be crucial: Rather than keeping the risk of erroneously selecting one individual gene at a fixed level, it is the decision involving thousands of genes given the amount of genes we can follow up on that is the focus, and where both types of error must be considered. The q-value does not fully address this problem, but nevertheless represents an improvement over the classical multiple test procedures in these applications.

More mathematically the q-value can be expressed as

Taking minimum in (7) enforces monotonicity in p _i , so that the q-value will be increasing (non-decreasing) in the observed p-value. If the FDR is non-increasing, as it should, then q (p _i ) = FDR (p _i ).

Additionally, the FDR offers a framework for power and sample size calculations, see [24] and the new developments in [25].

Simulated data

Two simulation models were used: one generating values independent between genes and the other generating observations displaying clumpy dependence [14, 30].

To generate dependent data the protocol from [14] was followed. This generates data following clumpy dependence in the sense of [30] such that blocks of genes have dependent expression. First a logarithmic normal distribution with mean 1.5 and standard deviation 0.3 generated a profile for each gene. Denoting by N(μ, σ) a normal distribution with mean μ and standard deviation σ, random errors following a standard normal distribution N(0,1) were added. To create dependencies genes were partitioned into sets of 50 and for each sample the same term from a N(0, 1) distribution was added to the expression of each gene in the set. Finally, genes were randomly assigned to become DEGs with probability 1-p₀ and for each gene a regulation term following either N(0.5, 0.2) or N(0.7, 0.2), with equal probabilities, was added to the expression of one of the groups of samples of size 30. The power to detect either of these two alternatives with a t-test at the 5% significance level equals 31% and 50%, respectively. The procedure generated for each of 400 iterations a set of observations of 10,000 genes. The protocol gives rise to high correlation within the blocks of 50 (on average on the order of 0.5). Results for smaller or weakly dependent datasets appear in the Additional file. Here weakly dependent means that the clumpy dependence term follows a N(0,1/20) distribution (correlation within blocks slumps to 0.003 on an average).

With the simulated independent data all methods for estimating p₀ perform rather well, see Table 2 and 3, and Figure 1.

	BUM	SPLOSH	smoother	bootstrap	SEP	LSL	mgf	PRE
Mean	0.039	0.061	0.038	0.036	0.045	0.18	0.072	0.022
Sd	0.048	0.078	0.032	0.034	0.032	0.12	0.048	0.016

True p 0		0.5	0.6	0.7	0.8	0.9	0.95	0.99
BUM	mean bias	0.013	0.043	0.038	0.028	-0.090	-0.050	-0.010
	Sd	0.0044	0.0074	0.0072	0.0076	0.033	0.0021	0.0010
SPLOSH	mean bias	-0.083	-0.068	-0.043	-0.012	0.018	0.026	0.14
	Sd	0.014	0.020	0.028	0.034	0.040	0.041	0.063
QVALUE	mean bias	-0.066	-0.057	-0.046	-0.034	-0.019	-0.015	0.0012
	Sd	0.022	0.022	0.023	0.027	0.027	0.024	0.0160
Bootstrap LSE	mean bias	-0.065	-0.054	-0.040	-0.025	-0.0067	-0.0037	0.0057
	Sd	0.023	0.023	0.021	0.024	0.026	0.024	0.023
SEP	mean bias	-0.084	-0.072	-0.059	-0.043	-0.028	-0.020	0.0068
	Sd	0.014	0.013	0.014	0.013	0.013	0.012	0.014
LSL	Mean bias	-0.36	-0.31	-0.25	-0.18	-0.096	-0.049	-0.010
	Sd	0.025	0.019	0.014	0.0086	0.0036	0.0022	0.0010
mgf	mean bias	-0.15	-0.12	-0.095	-0.067	-0.040	-0.027	-0.0095
	Sd	0.0052	0.0052	0.0051	0.0056	0.0056	0.01098	0.0018
PRE	mean bias	-0.036	-0.033	-0.028	-0.018	-0.0078	-0.018	-0.0018
	Sd	0.0098	0.010	0.0097	0.012	0.016	0.0080	0.0067

The new methods mgf and PRE were very competitive on these data, and had both low bias and variation, excluding mgf at the 0.5 and 0.6 level. Since mgf tends to overestimate p₀ rather much in the lower range, one may prefer PRE. For practical purposes though overestimation is desirable and enables control of the error rate (exact control in the terminology from [31]).

The smoother and the bootstrap had good and quite similar performance. They give more or less the same variation and bias over the whole range. This variation can be a bit high though, especially when comparing to PRE.

In the higher range BUM gives a crude estimate of the true p₀. In a certain lower range however it underestimates. As we can see in Figure 1, however, the method considerably overestimates p₀ in the higher range, which brings down the power to detect DEGs.

SPLOSH has the advantage of fitting the observed distribution quite well, judging from some tests (data not shown), compare also [11]. This enables a Bayesian analysis as in (3). However, it does the fitting of the p-values close to 0 sometimes at the expense of the accuracy concerning the values at the other end, thus misses the plateau and the minimisation in (10) will give a misleading result. In particular, this tends to happen when there are few DEGs. As we can see in Figure 1 the method underestimates p₀ at the higher range (p₀ ≥ 0.9), which is worrisome and may lead to underestimation of the error rate, which is undesirable for a method of statistical inference.

From Tables 2 and 3 we can see that PRE has the best over-all performance, followed by the smoother and the bootstrap. This does not however imply that the other methods could not be considered. The results vary quite a lot depending on the value of p₀: LSL is quite competitive for p₀ = 0.99, but too conservative for p₀ = 0.5.

With simulated data it is possible to calculate the actual false discovery rate. Figure 3 shows boxplots based on the simulated independent data. Here qvalue and pava FDR stand out as most reliable. Especially BUM FDR has severe problems, often over the whole range of cut-offs, while SPLOSH FDR does rather well.

The results concerning FDR appear in Figure 4. This time the BUM method captures FDR in a competitive way, with the lowest median error but at the same time the second worst mean error. Again the results for pava FDR appear quite stable and accurate, giving the second lowest median error and the lowest mean error.

Results for 300, 5,000 and 10,000 simulated weakly dependent genes appear in the Additional file. Briefly, they resemble those of the independent case. For 300 genes however the variation is such that the value of using these methods seems doubtful.

Real data

Data from Golub et al

These data represent a case where there are many DEGs. The data set concerning two types of leukaemia, ALL and AML, appeared in [32, 33]. Samples of both types were hybridised to 38 Affymetrix HG6800 arrays, representing 27 ALL and 11 AML. In reference [32] 50 genes were identified as DEGs using statistical analysis. The data consisted of Average Difference values and Absolute Calls, giving for each gene (probe set), respectively, the abundance and a categorical assessment of whether the gene was deemed Absent, Marginal or Present. The Average Difference values were pre-processed as in [5], and the proportion of samples where each probe set was scored Present was calculated, giving a present rate. A permutation t-test was used to compare ALL and AML. Figure 5 shows a histogram of t-test statistics revealing a bell-shaped region around the origin and large tails, indicating that a substantial part of the genes are DEGs.

Visual inspection of Figure 6 showing the p-value distribution also suggests that many genes are altered. The methods tested vindicate this : BUM, SPLOSH, bootstrap, qvalue, mgf, PRE, SEP, LSL give the estimates 0.57, 0.62, 0.65, 0.65, 0.69, 0.60, 0.64 and 0.86, respectively. Thus, roughly 35% of the genes are regarded as changed. The function qvalue finds 873 probe sets with a q-value less than 5%. The estimated of FDR appear in Figures 7 and 8, and show a great deal of concordance between methods. The curves rise to the respective estimate of p₀ and and in doing so slowly diverge. SEP approaches roughly 0.07 in a vicinity of zero, while the other methods continue the decline.

The function locfdr applied to Z = Φ^-1 (F₃₆(t)) with F₃₆ the cdf of the t-distribution with 36 degrees of freedom gives the estimate of p₀ 0.66 and outputs Figure 9 which bears witness of the skewness and of the different local false discovery rates in the two tails of the distribution.

Removing probe sets with less than 20% present rate will leave us with 2999 probe sets, and qvalue indicates that there are 977 of them that are significant with a q-value less than 5% (p₀ = 0.47). In general it is wise to remove probe sets with low presence rate prior to analysis, since this will make the inference more reliable, compare [20]. Doing so will most likely produce more true positives.

Data from Spira et al

Next let us turn to a case where there are rather few DEGs. In [34] the results from a microarray experiment where bronchial epithelial brush biopsies have been hybridised to Affymetrix U133A arrays are presented. The biopsies come from three different subject categories: Current smokers, Never smokers and Former smokers. The cel intensity files were downloaded from the NCBI Gene Expression Omnibus (accession no. GSE994) [35]. The Bioconductor package affy [19] was used to normalise intensities with the quantile method, and to calculate the RMA measure of abundance. The function mas5calls in affy output absolute calls.

Here we will take a brief look at the comparison between Former smokers and Never smokers. The comparison may help identify genes that remain changed after smoke cessation. A fuller analysis would include more analyses, such as the Current smokers vs. Former smokers comparison, and possibly also adjust for the fact that Former smokers tend to be older than Never smokers (Mean Age 45 and 53 years, respectively).

Graphical representations of the results appear in Figures 10 through 14.

Figure 10 suggests that rather few genes are changed. The methods give the following estimates of p₀ : bootstrap 0.92, mgf 0.97, PRE 0.97, LSL 1.00, SEP 0.98, SPLOSH 0.88, BUM 1.00 and smoother 0.92, compare Figure 11. The estimates of FDR appear in Figures 12 and 13, and show some separation between methods. The trio BH, pava and qvalue largely agree. SEP LFDR and FDR stabilise at a value above 0.8 when the cut-off approaches zero. This is consistent: If the local FDR levels off, then the Averaging Theorem, see Methods, implies that so will the FDR. However, this level may seem incompatible with the assessment of SEP that 2% of the genes are truly changed. A rough calculation, replacing f by a histogram, would yield the estimate of LFDR(p) = p₀/f(p) = 0.22 in a vicinity of p = 0.00025. Neither does it agree with locfdr, see Figure 14. Most likely the spline function employed by SEP is playing tricks here. SPLOSH approaches 0.4 in a vicinity of zero. BUM returns the estimate FDR ≡ 1, which is consistent with its estimate of p₀, but does not seem to agree with Figure 11. The method may display numerical problems with these data. The function qvalue finds one gene significant with a q-value less than 5%. However, if we restrict attention to the probe sets with at least 20% present rate (9841 genes) there are fifteen genes fulfilling the criterion.

Figure 14 displays the graphical output from locfdr, where the LFDR seems to approach zero in the tails, contradicting SEP, SPLOSH and BUM, but essentially agreeing with the concordant trio BH, pava and qvalue.

LSL

Let R(α) = # {i : p _i ≤ α}, the number of rejected given the cut-off α. In [36] the approximation

N - R(α) ≈ E[N - R(α)] ≈ N₀(1 - α)

for small α and N₀ = Np₀ the number of true null hypotheses appears. Consequently, (N - R(p_(i))/(1 - p_(i)) = (N - i)/(1 - p_(i)) will approximate N₀, which lead the pioneering authors to consider plotting 1 - p_(i)against N - i, thus giving them an estimate of N₀. In [26] the Lowest SLope estimator (LSL) of N₀ based on the slopes S _i = (1 - p_(i))/(N - i + 1) is presented. Starting from i = 1, the procedure stops at the first occurrence of Si₀ <Si₀-1, and outputs the estimate

In [12] the two above estimates are presented, derived and compared together with a method called Mean of Differences Method (MD). MD and LSL are motivated by assuming independence and approximating the gaps d_(i)= p_(i)- p_(i-1)(define p₍₀₎ = 0 and p_(N+1))≡ 1) with a Beta(1, N₀) distribution, which has expectation 1/(N₀ + 1). This expectation may be estimated by the inverse of a mean of the form

MD proceeds downward and chooses i₀ equal to the first j satisfying

Of these three methods LSL and MD give very similar results, and outperform their predecessor [12].

LSL is available as function fdr.estimate.eta0 in package GeneTS [18] with the option method= "adaptive".

The smoother

A method here referred to as the smoother appeared in [9]. This method, like all presented, is based on a comparison of the empirical p-value distribution to that of the uniform distribution. There will likely be fewer p-values close to 1 in the empirical than in the null distribution, which is a uniform. The ratio of the proportion of p-values greater than some η to the expected proportion under the uniform distribution, 1-η, will give a measure of the thinning of observed p-values compared to the null distribution. Thus, with F _e denoting the empirical distribution, the ratio {1-F _e (η)}/{1-η} will often be a good estimate of p₀ for an astutely chosen threshold η. A spline is fitted to the function p₀(η) = {1-F _e (η)}/{1-η}, and the resulting function is evaluated at η = 1, yielding the estimate

(tilde '~' above p₀ denoting the spline-smoothed version of p₀(η)) which goes into the calculation of the FDR in (2), see Figure 15. Note the relationship with LSL, p₀(p_(i)) = (N-i)/{(N-i+1)S _i N}. It can be shown that for fixed η this method offers a conservative point estimate of the FDR:

The q-value is estimated by combining (2), (7) and (11), and an implementation is provided as the function qvalue in package qvalue available on CRAN [18].

In [30] the authors go to great lengths to prove that for fixed η, as above, the conservativeness remains under various forms of dependency, such as clumpy dependence.

The Bootstrap LSE

Another approach pioneered by Storey in [3] is to use a bootstrap least squares estimate (LSE), which chooses a value of η in p₀(η), that minimises the variation of the estimate for random samples of the original p-values. The bootstrap standard reference [37] provides more theoretical background. Generate B new samples p-values p* ^b (b = 1,..., B) by sampling with replacement from the observed ones, calculate a measure of the Mean Squared Error (MSE) of the corresponding estimates p* ^b ₀(η) for a lattice of values of η and choose the value minimising the MSE. More formally, the optimal η is obtained through

The version of the bootstrap used in this article uses more samples B than the version available in qvalue (B = 500 instead of B = 100), and seems to perform better (data not shown).

Available in functions qvalue [18] and p0.mom (in package SAGx) [18, 38].

SPLOSH

In [11] a spline function estimates the log-transformed pdf log[f(x)] using a complex algorithm involving splines called spacings LOESS histogram (SPLOSH). To obtain a stable estimate of FDR near zero a technique from mathematical analysis called l'Hospital's rule is used to approximate the ratio in (5) and to yield

,

where the numerator has been estimated as in (4). An R package with the same name is available [39].

The FDR estimate (6) is used with F obtained by the non-parametric estimate of the pdf.

Note that we can now calculate the posterior probability given its p-value that a gene is a DEG as p₁(x) = 1 - p₀/f(x), compare (3).

The method is available in R function splosh [40].

BUM

In [10] the authors assume a beta-uniform (BUM) distribution, i.e. in (1) they replace f₁ by a beta distribution,

f(x) = λ + (1 - λ)ax^a-1    (12)

where in addition to λ which corresponds to p₀ the shape parameter a has to be estimated. Thanks to the simple form of the distribution it is possible to estimate parameters through the maximum likelihood principle, i.e. by choosing values that maximise the likelihood of observing the p-values that were actually observed. However, due to problem in identifying p₀ with λ, the authors instead use an upper bound

which corresponds to f(1).

The FDR estimate (5) is used with F the cdf corresponding to (12).

The authors provide R code for the application of their method [39].

A more intricate Hierarchical Bayes model based on the beta-uniform concept allowing for different parameter values in different intervals appears in [41]. The R function localFDR provides an implementation of the method [42].

Poisson regression

In [28] it is suggested to estimate any empirical distribution by dividing the real axis into intervals and regarding the number of hits in each interval as the result of an inhomogeneous Poisson process, much like counting the number of cars arriving at a crossing during different time intervals. This method was used in [27] to model the distribution of a transformed test statistic, it also appears in function locfdr which estimates a local FDR as a function of a test statistic. In our case, of course, the support of the distribution is the unit interval [0,1]. Then the expected number of hits in each subinterval of [0,1] can be modelled as a polynomial in the midpoints of the subintervals by a technique called Poisson regression (PRE). The approach taken here is to choose a polynomial of low degree so that the plateau representing the uniform distribution is well captured. In doing so the ability the capture the distribution at low p-values is sacrificed.

A more mathematical description now follows. The PRE method assumes that the counts S _k follow a Poisson distribution whose intensity is determined by the midpoint t _k of the interval I _k , see [28]. To be specific: in the current application it is assumed that the expected frequency of observations in an interval is given by

where μ ^o _k are the smoothed observed frequencies in each interval I _k . In statistical jargon this is a Poisson regression model with μ ^o _k as offset. This assumes independence between counts in different intervals. Although this does not hold true the model succeeds to capture the essential features of distributions. Standard functions in e.g. R can fit this model. Normalising the curve by the total number of p-values we get an estimate of the pdf. Finally, smooth the pdf f(x) with a spline to obtain a more stable result, and use the estimate (10). An implementation of PRE is provided through R function p0.mom in package SAGx [18, 38].

SEP

The Successive Elimination Procedure (SEP) excludes and includes p _i successively such that the similarity of the distribution of the included tends to behave increasingly like a uniform [13]. Finally, an index set J _final will map to a set of p-values that represent the true null hypotheses. This yields the point estimate

with N _J = # J, the cardinality of the set J. The identification of J _final proceeds by an intricate optimisation algorithm where the objective function consists of two terms : one Kolmogorov-Smirnov score

for the empirical cdf F _J (based on J), to measure the distance to a uniform, and one penalty term to guard against overfitting

for some tuning parameter λ.

A local FDR is obtained from smoothed histogram estimates based on equidistant bins

,

where the function h₀ refers to the J _final set and h to the total set of p-values.

The function twilight in package twilight provides an implementation of SEP [19].

Moment generating function approach

The next approach is based on the moment generating function (mgf), which is a transform of a random distribution, which yields a function M(s) characteristic of the distribution, cf. Fourier or Laplace transforms, e.g. [43]. Knowing the transform means knowing the distribution. It is defined as the expectation (or the true mean) of the antilog transform of s times a random variable X, i.e. the expectation of e ^sX or in mathematical notation:

M(s) = ∫e ^sx f(x)dx.

To calculate the mgf for p-values, we use the fact that the pdf is a mixture of pdf's (8). This yields the weighted sum of two transformed distributions:

,

where we have used the fact that the mgf of a uniform distribution over [0,1] equals g(s) = (e ^s - 1)/s. Denoting the second transform by M₁(s) we finally have

M(s) = p₀g(s) + (1 - p₀)M₁(s).     (13)

Now, the idea is to estimate these mgf's and to solve for p₀. In the above equation M(s) can be estimated based on an observed vector of p-values and g(s) can be calculated exactly, respectively, while p₀ and M₁(s) cannot be estimated independently. The estimable transform is, given the observed p-values p = p₁,..., p _n , estimated by

Then, one can solve (13) for p₀:

Let us do so for s _n > s_n-1, equate the two ratios defined by the right hand side in (14) and solve for M₁(s _n ). This gives the recursion

If we can find a suitable start for this recursion we should be in a position to approximate the increasing function M₁(s) for s = s₁ <s₂ < ... <s _m in (0, 1]. Now, note that 1 ≤ M(s), for any mgf, with close to equality for small values of s. It makes sense to start the recursion with some M₁(s₁) in I = [1, M(s₁)]. In general, it will hold true that 1 ≤ M₁(s) <M(s) <g(s), since f₁ puts weight to the lower range of the p-values at the expense of the higher range, the uniform puts equal weight, and f being a mixture lies somewhere in between. We can calculate g, M and M₁ for an increasing series of values in [0,1], e.g. for s = (0.01, 0.0101, 0.0102, ..., 1). The output from one data set appears in Figure 16. Since all ratios (14) should be equal, a good choice of M₁(s₁) will be one that minimises the variation of the ratios. Standard one-dimensional optimisation will find the value in I that minimises the coefficient of variation (CV, standard deviation divided by mean)

where s = (s₁, ..., s _n ). The CV will in contrast to the variance put both small and high values of the ratios on an equal footing and enable comparison.

Finally, these essentially equal ratios provide an estimate of p₀.

A heuristic convexity argument suggests that mgf over-estimates p₀, see the Additional file. Furthermore, the bias seems to decrease as p₀ grows.

An implementation of mgf appears as function p0.mom in package SAGx [18, 38].

Local FDR and FDR

The concept of a local false discovery rate originates from [1]. Let the (true) local FDR at t be defined as the probability that gene i is unchanged conditional upon that its p-value equals t, or in formulas : LFDR(t) = Pr(gene i unchanged | p _i = t) = p₀/f(t). The Averaging Theorem of [4] states that integrating the local FDR over the rejection region R, such as R = [0, 0.01], yields the FDR : FDR(R) = E[LFDR(y) | y ∈ R]. In [29] it is noted that the estimated q-value equals the mean of a local FDR

where the local FDR at the i ^th ordered p-value p_(i)equals

,

where N denotes the total number of genes. This rephrases the theorem in terms of estimators. The local FDR is meant as an approximation of the probability that an individual gene i is a DEG. As remarked in [29] the q-value does not estimate the probability that a gene is a false positive. Indeed, the theorem shows that it is the mean of that probability for all genes at least as extreme as i. Thus the q-value will tend to give a lower value than LFDR(i).

Under a wide range of models, where f(x) is non-increasing, e.g. the BUM model, the expected local LFDR(i) will be non-increasing, and hence the differences above should tend to increase, see the Additional file. Hence there is a need for enforcing monotonicity as in (7). One tool for enforcing monotonicity is the Pooling of Adjacent Violators (PAVA) algorithm [44]. This algorithm has an intuitive appeal, is less ad-hoc than the local spline approach presented in [29], and is the non-parametric maximum likelihood estimate under the assumption of monotonicity. As an example of how it works, consider the series (1,2,4,3,5), which PAVA turns into the non-decreasing series (1, 2, 3.5, 3.5, 5) by pooling the violators of monotonicity (4, 3) and replacing them by their mean. Though not equivalent to the q-value from (2) and (7), the results from applying PAVA to the terms in (15) agreed rather well with the values obtained from function qvalue. In the Results section this approach combined with the PRE estimate of p₀ is referred to as pava FDR. We could have used mgf for calculating FDR, but it was excluded due to the better over-all performance of PRE.

The bootstrap LSE gives a very similar result to the smoother and thus was excluded in comparison of FDR estimates.

The R function pava.fdr in package SAGx provides an implementation of pava FDR, and returns a list of estimates of FDR, LFDR and p₀ [18, 38].

The reference [25] presents the theory behind the estimation of local false discovery rates provided by R package locfdr [18, 25]. The method procedes by transforming the test statistic t into Z = Φ^-1(m(t)), where Φ is the cdf corresponding to N(0,1) and m is a transform that for the uninteresting/null class renders the distribution of Z into a N(0,1). Typically, m could equal the cdf of the t-distribution. Assuming the model (1) the method estimates the mixture density f ^t using Poisson regression, and fits either a theoretical null sub density (from now on suppressing superscript t and denoting the pdf corresponding to Φ by φ) f⁺₀(z) = p₀ φ(z) around z = 0, or fits an empirical null distribution. Then the procedure estimates the local false discovery rate fdr(z) = f⁺₀(z)/f(z).