Common Implementation Errors In FDR Analysis

Why Implementation Details Matter

False discovery rate procedures turn a collection of p-values into a discovery list. In genomic analyses, the collection may contain thousands of tests, so a small coding error can change which genes or variants appear in the final list. Errors can also make results difficult to reproduce: a reported threshold has little meaning if the number and identity of the tests it was applied to are unclear.

This tutorial focuses on implementation checks for the Benjamini–Hochberg (BH) procedure and Storey-style q-values. The formulas are straightforward; the main risks are mismatching p-values with hypotheses, using the wrong rank or family size, omitting the monotonicity adjustment, treating an estimate of the null proportion as a fixed truth, and leaving the tested family unspecified.

State the Procedure Before Calculating

Suppose a fixed family contains \(m\) hypotheses and has p-values \(p_1,\ldots,p_m\). Under each true null hypothesis, a valid p-value is super-uniform: for every \(t\) between 0 and 1, \(P(p_i\le t)\le t\). The standard BH false discovery rate guarantee applies when the null p-values are valid and the tests are independent or satisfy positive regression dependence on a subset (PRDS) of the true nulls. The dependence material earlier in this unit distinguishes these conditions from arbitrary dependence; implementation accuracy does not itself establish them.

Definition: Write the p-values in nondecreasing order as \(p_{(1)}\le\cdots\le p_{(m)}\), retaining the identity of the hypothesis associated with each value. At target level \(\alpha\), the BH step-up procedure sets \(k=\max\{i:p_{(i)}\le i\alpha/m\}\), with \(k=0\) if no rank satisfies the inequality, and rejects the hypotheses corresponding to ranks \(1,\ldots,k\).

The rank is the position in the sorted list, starting at 1; \(m\) is the total number of tests in the stated family. BH is a step-up rule: after finding the largest qualifying rank \(k\), it rejects every hypothesis with rank at most \(k\). It does not reject only the p-values that individually pass their own rank-specific inequality.

A BH adjusted p-value for rank \(i\) is

p^{\mathrm{BH}}_{(i)}=\min\left(1,\min_{j\ge i}\frac{m p_{(j)}}{j}\right).

The inner minimum enforces monotonicity; the outer minimum caps the value at 1. Rejecting a hypothesis when its BH adjusted p-value is at most \(\alpha\) gives the same rejection set as the step-up rule.

A common Storey-style procedure replaces \(m\) in the adjusted-value calculation by an estimate of the number of true null hypotheses, \(m\hat{\pi}_0\). One basic estimate, for a chosen tuning parameter \(0\le\lambda<1\), is

\hat{\pi}_0(\lambda)=\min\left(1,\frac{\#\{i:p_i>\lambda\}}{m(1-\lambda)}\right).

A corresponding Storey-style q-value is

q_{(i)}=\min\left(1,\min_{j\ge i}\frac{m\hat{\pi}_0 p_{(j)}}{j}\right).

This formula describes one common calculation, not every software implementation. Smoothing, bootstrap selection of \(\lambda\), or other stabilization choices can produce different estimates. Report the method used. The earlier material on sensitivity to \(\lambda\) explains why the estimate can vary with that choice.

Worked BH Calculation and Label Tracking

Worked Example: Sorting Without Losing the Hypothesis Labels

Suppose four hypothetical genomic tests have these p-values: gene A, 0.018; gene B, 0.004; gene C, 0.200; and gene D, 0.013. Apply BH at \(\alpha=0.05\).

First sort the values while carrying the gene identifiers with them: rank 1 is B, 0.004; rank 2 is D, 0.013; rank 3 is A, 0.018; and rank 4 is C, 0.200. Here \(m=4\). The rank thresholds are \(i\alpha/m\): 0.0125, 0.0250, 0.0375, and 0.0500.

Compare each sorted p-value with its threshold. At rank 1, \(0.004\le0.0125\); at rank 2, \(0.013\le0.0250\); at rank 3, \(0.018\le0.0375\); and at rank 4, \(0.200>0.0500\). Thus the largest qualifying rank is \(k=3\), and the BH discovery list is genes B, D, and A—not the first three labels in the original input order.

As a check, the unadjusted rank calculations \(mp_{(i)}/i\) are 0.016, 0.026, 0.024, and 0.200. Taking suffix minima gives BH adjusted values 0.016, 0.024, 0.024, and 0.200 for ranks 1 through 4. The first three are at most 0.05, agreeing with the step-up decision.

Why Adjusted Values Need a Monotone Pass

The expression \(mp_{(i)}/i\) is not, by itself, the BH adjusted p-value. It may decrease as rank increases, even though a hypothesis with a larger p-value should not end up with a smaller adjusted value than a preceding hypothesis in the sorted list. The suffix minimum in the definition corrects this.

Justification: Let \(r_i=mp_{(i)}/i\) and \(a_i=\min_{j\ge i}r_j\). Then \(a_i=\min(r_i,a_{i+1})\) for \(i<m\), so \(a_i\le a_{i+1}\). Also, \(a_i\le r_i\), and no value \(r_j\) with \(j\ge i\) is smaller than \(a_i\), by the definition of a minimum. Thus the suffix minimum is the smallest rank-based value available from rank \(i\) onward and produces a nondecreasing adjusted-value sequence.

Worked Example: Correcting Nonmonotone Raw Values

For a hypothetical family of \(m=4\) tests, let the sorted p-values be 0.001, 0.006, 0.007, and 0.040. The raw values \(mp_{(i)}/i\) are, in rank order, \(4(0.001)/1=0.004\), \(4(0.006)/2=0.012\), \(4(0.007)/3=0.00933\), and \(4(0.040)/4=0.040\).

These raw values decrease from 0.012 at rank 2 to about 0.00933 at rank 3. Reporting them directly would therefore give a larger adjusted value to the rank-2 p-value than to the larger rank-3 p-value. Apply the suffix minimum: rank 4 is 0.040; rank 3 is \(\min(0.00933,0.040)=0.00933\); rank 2 is \(\min(0.012,0.00933,0.040)=0.00933\); and rank 1 is \(\min(0.004,0.012,0.00933,0.040)=0.004\). The correct adjusted values are 0.004, 0.00933, 0.00933, and 0.040.

At \(\alpha=0.01\), the first three adjusted values are at most 0.01, so the first three hypotheses are rejected. The raw rank-2 value of 0.012 would incorrectly suggest otherwise if it were mistaken for the adjusted value.

Estimating \(\pi_0\) Without Treating It as Known

The Storey estimate uses the share of p-values above \(\lambda\) to estimate the null proportion. The intuition is that p-values from true nulls tend to spread across the interval from 0 to 1, whereas some nonnull p-values concentrate near 0. The estimate is not a count of confirmed null hypotheses: it depends on the observed p-values and the chosen \(\lambda\). A small family can give a coarse or unstable estimate, and a change in \(\lambda\) can alter the resulting q-values.

Worked Example: Calculating a Storey-Style Estimate

Consider five hypothetical p-values, sorted as 0.01, 0.04, 0.20, 0.65, and 0.90. Choose \(\lambda=0.50\). There are two values greater than 0.50, so

\hat{\pi}_0=\min\left(1,\frac{2}{5(1-0.50)}\right) =\min(1,0.80)=0.80.

Using this estimate, the raw Storey-style rank values are \(5(0.80)(0.01)/1=0.040\), \(5(0.80)(0.04)/2=0.080\), \(5(0.80)(0.20)/3=0.2667\), \(5(0.80)(0.65)/4=0.650\), and \(5(0.80)(0.90)/5=0.720\). They are already nondecreasing, so the suffix-minimum pass leaves them unchanged. At a q-value threshold of 0.05, only the first test is selected.

For comparison, setting \(\pi_0=1\) would give a first adjusted value of \(5(0.01)/1=0.050\), rather than 0.040. This illustrates the effect of the estimate; it does not establish that 0.80 is the true null proportion. A different \(\lambda\) or estimation method could yield a different result.

Define the Family Before Interpreting the List

The family is the set of hypotheses to which the multiple-testing procedure is applied. Its size \(m\) belongs in the calculation, so family definition is not just a reporting detail. For example, an analysis might control FDR separately within each phenotype, or across all gene-by-phenotype tests together. These choices use different families and can yield different adjusted values. The earlier material on FDR across phenotypes and annotations considers the implications of such choices.

Worked Example: Separate or Combined Hypothesis Families

Suppose two phenotypes each have five tests. For a hypothetical p-value of 0.012 at rank 1 within its phenotype, BH at \(\alpha=0.05\) gives a rank-1 threshold of \(0.05/5=0.010\) when the family consists of that phenotype’s five tests. The p-value does not pass this rank-1 threshold.

If, instead, all ten gene-by-phenotype tests form one family and this p-value is rank 1 in the combined sorted list, its rank-1 threshold is \(0.05/10=0.005\), so it still does not pass. These rank-1 comparisons alone do not determine the full step-up decision: the largest qualifying rank must be checked in the complete sorted family. The calculation does show why stating whether the family has five or ten hypotheses is essential.

A defensible report names the tested units, the number of hypotheses included, whether analyses were combined or conducted separately, how any excluded tests were handled, the procedure and tuning choices, and the target level. Filtering or removing tests after seeing their p-values can affect validity; it should not be treated as an innocuous change to \(m\).

Diagnostics, Alternatives, and Common Errors

Use implementation checks that preserve the full path from input to decision. Confirm that every p-value is finite and lies in \([0,1]\), that each identifier remains attached to its p-value after sorting, and that the number sorted equals the declared \(m\). Check the sorted sequence is nondecreasing, the BH thresholds use ranks \(1,\ldots,m\), and the adjusted values are nondecreasing and bounded by 1. For Storey-style results, record \(\lambda\), the count above it, the resulting \(\hat{\pi}_0\), and any stabilization or truncation. Recompute a small test case by hand or compare the implementation with a trusted calculation.

Do not silently drop missing p-values and then leave the original \(m\) in the formula, or replace \(m\) without documenting which tests were removed. Either choice changes the calculation. Likewise, using a different family for the BH threshold and the reported discovery list makes the result uninterpretable. A software result is not a substitute for specifying the family and the procedure.

Wrong: Sort the p-values but keep the hypothesis names in their original order. Correction: Sort pairs of identifiers and p-values together, then map discoveries back to the original records.
Wrong: Use \(p_{(i)}\le\alpha/(m-i+1)\), or start ranks at zero. Correction: For BH, use \(p_{(i)}\le i\alpha/m\), with ranks starting at 1.
Wrong: Call \(mp_{(i)}/i\) the adjusted value at every rank. Correction: Take the minimum over that rank and all larger ranks, then cap at 1.
Wrong: Interpret \(\hat{\pi}_0\) as the known fraction of true null hypotheses. Correction: It is an estimate whose value depends on the method and, for this estimator, on \(\lambda\).
Wrong: Say only that “FDR was controlled at 0.05.” Correction: State the family, procedure, relevant dependence assumptions or guarantee, and how the discovery list was formed.

BH is a standard choice under independence or PRDS conditions for valid null p-values; Benjamini–Yekutieli is a more conservative alternative designed for arbitrary dependence. Storey-style q-values can use an estimated null proportion and may be more powerful when that estimate is appropriate, but add estimation choices that must be reported. These procedures are not interchangeable just because they share the phrase “FDR control.” Their guarantees, assumptions, and calculations differ.

The central limitation is that flawless arithmetic cannot repair invalid p-values, an unsuitable family definition, or dependence conditions that do not support the claimed guarantee. Nor does an FDR target mean that each selected finding has that probability of being false. It describes an expected proportion of false discoveries among rejections, with the procedure’s stated assumptions and scope.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For \(m=5\), \(\alpha=0.05\), and sorted p-values 0.003, 0.012, 0.025, 0.040, 0.30, find the largest BH qualifying rank and state the discoveries.
For \(m=3\), sorted p-values 0.002, 0.010, 0.020, calculate the BH adjusted values.
With \(m=8\), \(\lambda=0.75\), and two p-values greater than \(\lambda\), calculate the basic Storey estimate \(\hat{\pi}_0\).
What error occurs if p-values are sorted but the gene labels are not moved with them?
Name two details about the hypothesis family that should appear in a report.

Solutions. (1) The thresholds are 0.010, 0.020, 0.030, 0.040, and 0.050. Ranks 1 through 4 qualify, while rank 5 does not, so \(k=4\) and the first four hypotheses in sorted order are rejected. (2) The raw values \(3p_{(i)}/i\) are 0.006, 0.015, and 0.020; they are nondecreasing, so the adjusted values are 0.006, 0.015, and 0.020. (3) \(\hat{\pi}_0=\min(1,2/[8(1-0.75)])=\min(1,1)=1\). (4) The rejection decision can be assigned to the wrong genes, because the sorted ranks no longer identify the hypotheses that produced those p-values. (5) For example, report which tests were included and the total number \(m\); also state whether tests across phenotypes or annotations were combined or analyzed as separate families.