Deriving Storey Q-Values From A P-Value List

From a P-Value List to q-Values

A genomic analysis may produce thousands of p-values, but the calculation of Storey q-values can be understood from a short list. The essential steps are to estimate the proportion of true null hypotheses, calculate an estimated false discovery rate at each relevant cutoff, and assign each test the smallest such estimate among cutoffs that would include it. The final minimum step matters: estimated false discovery rates across cutoffs need not increase smoothly.

The q-value calculation operates on one defined family of \(m\) tests. The p-values must be valid for their respective null hypotheses; in particular, under a true null, a p-value should be uniformly distributed on \([0,1]\) when the null model is continuous, or at least superuniform in the sense that \(P(p_i\leq t)\leq t\). The Storey estimate additionally uses the approximate uniformity of null p-values in the upper tail. Calculating q-values does not itself establish that these assumptions hold or guarantee FDR control under every dependence structure.

Definition: For a chosen tuning value \(\lambda\) between 0 and 1, Storey’s estimate of the null proportion is commonly calculated as \(\hat{\pi}_0(\lambda)=\min\{1,\#(p_i>\lambda)/[m(1-\lambda)]\}\). At a cutoff \(t\) with \(R(t)=\#(p_i\leq t)>0\), the estimated false discovery rate is \(\widehat{\mathrm{FDR}}(t)=\hat{\pi}_0 m t/R(t)\), conventionally capped at 1. The q-value for a test with p-value \(p_i\) is the minimum estimated FDR over cutoffs \(t\geq p_i\), also capped at 1.

Here \(R(t)\) is the number of rejections that would be made by rejecting every test with p-value at most \(t\). In the sorted list \(p_{(1)}\leq\cdots\leq p_{(m)}\), a convenient equivalent calculation is

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

The index \(j\) represents the rejection count at an observed cutoff. With ties, all p-values equal to a cutoff are included in \(R(t)\); the rank at the end of a tied group gives that count. The reverse minimum formula handles this automatically.

Why These Calculations Are Used

If there are \(m_0\) true null hypotheses and their p-values are uniform, the expected number of true-null p-values exceeding \(\lambda\) is \(m_0(1-\lambda)\). Dividing the observed number above \(\lambda\) by \(m(1-\lambda)\) therefore estimates \(m_0/m\), the null proportion. This is the motivation for \(\hat{\pi}_0\). In finite data, the ratio varies randomly; moreover, non-null p-values can also exceed \(\lambda\), so the estimate is not automatically exact.

At a cutoff \(t\), the expected number of true-null p-values at or below \(t\) is \(m_0t\) under uniform null p-values. Replacing the unknown \(m_0\) by \(\hat{\pi}_0m\), then dividing by the observed rejection count \(R(t)\), gives the estimated FDR formula. It is an estimate, not an observed fraction of false discoveries: the truth status of each rejected hypothesis is generally unknown.

A test with p-value \(p_{(i)}\) is included at every cutoff \(t\geq p_{(i)}\). Its q-value uses the smallest estimated FDR across those eligible cutoffs. It is therefore possible for the raw estimate at its own cutoff to exceed its q-value: a more inclusive cutoff may have a smaller estimated FDR. For a fixed list, q-values are nondecreasing in p-value rank.

Worked Example: Calculating q-Values for Eight Tests

Suppose a hypothetical genomic analysis has \(m=8\) tests with sorted p-values \(0.002,\ 0.008,\ 0.021,\ 0.049,\ 0.120,\ 0.310,\ 0.620,\ 0.840\). For this calculation, take \(\lambda=0.50\).

Step 1: Estimate \(\pi_0\). Two p-values, 0.620 and 0.840, exceed 0.50. Thus

\hat{\pi}_0=\frac{2}{8(1-0.50)}=\frac{2}{4}=0.5.

The estimate is \(0.5\), or 50.0%, for this calculation. The multiplier \(\hat{\pi}_0m\) is \(0.5(8)=4\).

Step 2: Calculate the cutoff estimates. At the cutoff \(p_{(j)}\), there are \(j\) rejections because all eight p-values are distinct. The estimated FDR is \(4p_{(j)}/j\).

Rank \(j\)	Cutoff \(p_{(j)}\)	\(4p_{(j)}/j\)
1	0.002	0.008
2	0.008	0.016
3	0.021	0.028
4	0.049	0.049
5	0.120	0.096
6	0.310	0.207
7	0.620	0.354
8	0.840	0.420

For example, at cutoff 0.120, \(R(t)=5\), so the estimated FDR is \(0.5(8)(0.120)/5=0.096\). The displayed estimates are already nondecreasing, so the reverse cumulative minimum leaves them unchanged. The q-values in rank order are therefore 0.008, 0.016, 0.028, 0.049, 0.096, 0.207, 0.354, and 0.420.

The first test has q-value 0.008: among the listed cutoffs that include its p-value, the smallest estimated FDR is 0.008. This does not mean there is a 0.8% probability that this particular null hypothesis is true.

The Monotone Minimum: A Necessary Calculation

The raw cutoff estimates \(a_j=\hat{\pi}_0mp_{(j)}/j\) need not increase with \(j\). But a q-value is a minimum over the available cutoffs. For rank \(i\), every later rank \(j\geq i\) gives a cutoff that includes test \(i\), so its q-value must be the smallest \(a_j\) to its right, subject to the cap at 1. This directly proves the reverse cumulative minimum calculation. It also proves monotonicity: when moving from rank \(i\) to \(i+1\), the set of eligible cutoffs becomes smaller, so its minimum cannot decrease.

\hat q_{(i)}=\min\{1,a_i,a_{i+1},\ldots,a_m\}, \qquad \hat q_{(i)}\leq \hat q_{(i+1)}.

The following example shows why the minimum cannot be skipped.

Worked Example: Applying the Reverse Cumulative Minimum

Consider six hypothetical sorted p-values \(0.010,\ 0.030,\ 0.040,\ 0.200,\ 0.700,\ 0.900\), with \(\lambda=0.50\). Two of six p-values exceed \(\lambda\), giving

\hat{\pi}_0=\frac{2}{6(1-0.50)}=\frac{2}{3}=0.6667, \qquad \hat{\pi}_0m=4.

The raw estimates \(4p_{(j)}/j\), rounded to four significant figures, are \(0.04000,\ 0.06000,\ 0.05333,\ 0.2000,\ 0.5600,\ 0.6000\). They decrease from rank 2 to rank 3, so they cannot themselves be the q-values.

Start at the last rank and take running minima from right to left. Rank 6 has q-value 0.6000. Rank 5 has \(\min(0.5600,0.6000)=0.5600\). Rank 4 has \(\min(0.2000,0.5600,0.6000)=0.2000\). Rank 3 has \(\min(0.05333,0.2000,0.5600,0.6000)=0.05333\). Rank 2 has \(\min(0.06000,0.05333,0.2000,0.5600,0.6000)=0.05333\), and rank 1 has q-value 0.04000.

In rank order, the q-values are \(0.04000,\ 0.05333,\ 0.05333,\ 0.2000,\ 0.5600,\ 0.6000\). For the second test, the estimate at its own cutoff is 0.06000, but the cutoff 0.040 includes it and has estimated FDR \(4(0.040)/3=0.05333\), which is smaller. That later cutoff determines its q-value.

Tied Values, Checks, and Interpretation

Worked Example: Counting a Tied Cutoff Correctly

Suppose a hypothetical list of six p-values is \(0.010,\ 0.020,\ 0.020,\ 0.080,\ 0.750,\ 0.900\), again with \(\lambda=0.50\). The estimate is \(\hat{\pi}_0=2/[6(0.50)]=2/3\), so \(\hat{\pi}_0m=4\). At cutoff \(t=0.020\), three p-values are at or below the cutoff, including both ties. Thus

\widehat{\mathrm{FDR}}(0.020)=\frac{4(0.020)}{3}=0.02667.

The tied observations occupy ranks 2 and 3. At the cutoff 0.020 their common q-value is 0.02667, using the rejection count of 3, not 2. The raw rank-2 expression \(4(0.020)/2=0.040\) is not the cutoff estimate for that tied group; the rank-3 expression correctly accounts for all tests included at the cutoff. The full q-value sequence is \(0.02667,\ 0.02667,\ 0.02667,\ 0.08000,\ 0.6000,\ 0.6000\).

Several checks help catch errors. First, \(\hat{\pi}_0\) and the q-values should be between 0 and 1 under the stated caps. Second, q-values in sorted p-value order should be nondecreasing. Third, for each cutoff, verify \(R(t)\) by counting every p-value less than or equal to \(t\), especially at ties. Finally, keep the p-value family fixed: changing which tests are included changes \(m\), the estimate, and the resulting q-values.

Storey’s calculation is closely related to the Benjamini-Hochberg adjusted p-value calculation: both use a rank-based quantity proportional to \(mp_{(j)}/j\) and a reverse minimum. The Storey version scales that quantity by \(\hat{\pi}_0\), rather than using 1 for the null proportion. It can therefore produce smaller q-values when \(\hat{\pi}_0<1\). The estimate of \(\pi_0\) is itself uncertain, and the resulting q-values are not posterior probabilities that individual hypotheses are null. They are minimum estimated FDRs over rejection cutoffs that include the test.

A q-value threshold is useful for organizing a discovery list, but it should not be described as certainty about each finding. For example, reporting tests with q-values at most 0.05 means those tests lie at cutoffs with estimated FDR no greater than 0.05 according to this calculation; it does not prove that exactly 5% of those particular findings are false. The distinction between an estimated error rate and a realized, unknown error proportion remains important in genomic interpretation.

Common Mistakes and Limitations

Using the wrong denominator for \(\hat{\pi}_0\). The denominator is \(m(1-\lambda)\), not \(m\lambda\). For example, with \(m=8\), \(\lambda=0.50\), and two p-values above \(\lambda\), the estimate is \(2/[8(0.50)]=0.5\).
Forgetting the reverse minimum. The raw values \( \hat{\pi}_0mp_{(j)}/j\) can fall at a later rank. Taking them directly as q-values can violate the definition and produce a nonmonotone sequence.
Counting only one observation in a tie. A cutoff includes all p-values equal to that cutoff, so \(R(t)\) must include the full tied group.
Calling a q-value a probability that a null is true. It is a minimum estimated FDR for eligible cutoffs, not a probability assigned to an individual hypothesis.
Treating the estimate as exact. The tail count is random, and alternatives can contribute p-values above \(\lambda\). The calculation does not incorporate uncertainty in \(\hat{\pi}_0\).

The calculation is only as useful as the p-values, family definition, and estimate behind it. If null p-values are poorly calibrated, or the fixed \(\lambda\) yields an unstable tail count, the resulting q-values may be misleading. A numerical q-value list should therefore be accompanied by a clear description of the tested family and the chosen \(\hat{\pi}_0\) calculation.

Check Your Understanding

Use the definitions above. The solutions are included so that each calculation can be checked.

For \(m=10\), \(\lambda=0.60\), and three p-values above \(\lambda\), calculate \(\hat{\pi}_0\).
Suppose \(\hat{\pi}_0m=3\), and at a cutoff \(t=0.04\) there are \(R(t)=4\) rejections. Calculate the estimated FDR at that cutoff.
For a sorted list, the raw cutoff estimates at ranks 1 through 4 are \(0.02,\ 0.07,\ 0.05,\ 0.30\). Find the q-values.
At a cutoff equal to a tied p-value, there are two p-values below it and three equal to it. What rejection count belongs in the denominator?

Solutions. (1) \(\hat{\pi}_0=3/[10(1-0.60)]=3/4=0.75\). (2) \(\widehat{\mathrm{FDR}}(0.04)=3(0.04)/4=0.030\). (3) Take running minima from right to left: rank 4 is 0.30; rank 3 is \(\min(0.05,0.30)=0.05\); rank 2 is \(\min(0.07,0.05,0.30)=0.05\); rank 1 is \(\min(0.02,0.07,0.05,0.30)=0.02\). The q-values are \(0.02,\ 0.05,\ 0.05,\ 0.30\). (4) Use \(R(t)=5\), because the cutoff includes both smaller p-values and all three tied values.