Assessing Lambda Sensitivity In Storey’s Method

Why Assess Sensitivity Across Lambda?

Storey’s method estimates the proportion of hypotheses that are truly null, denoted by $\pi_0$, from the observed p-values. Its estimate depends on a tuning value $\lambda$: only p-values above $\lambda$ contribute to the estimate. A single estimate can conceal important information. If modest changes to $\lambda$ produce large changes in $\hat{\pi}_0$, then the resulting q-values may also depend substantially on this choice.

The earlier material on choosing lambda introduced the estimator and ways to select a value. Here the emphasis is diagnostic: interpret estimates across a grid, identify why they fluctuate, and determine whether bootstrap selection or smoothing changes the conclusions. Stability is evidence about sensitivity to these choices; it is not proof that the estimator’s assumptions hold or that every reported q-value is accurate.

Definition: For $m$ tested hypotheses with p-values $p_1,\ldots,p_m$, Storey’s estimate at $0\leq\lambda<1$ is $$ \hat{\pi}_0(\lambda)=\frac{\#\{i:p_i>\lambda\}}{m(1-\lambda)}. $$ The numerator counts p-values above $\lambda$, and the denominator is the expected number above $\lambda$ if all $m$ tests were null and their p-values were independent Uniform(0,1). In practice, an estimate used to calculate q-values is often constrained to the interval $[0,1]$.

The estimator relies on null p-values being Uniform(0,1), or at least having the appropriate uniform upper tail. To interpret the estimate in the usual way, the tests must also be defined as a coherent family and the p-values must be valid under their null hypotheses. Independence is useful for the variance calculations and ordinary bootstrap discussion below; it is not a universal prerequisite for the estimator’s definition. Dependence can affect variability and the performance of resampling, so a smooth-looking grid does not establish independence.

What the Lambda Grid Reveals

Under the null, a p-value exceeds $\lambda$ with probability $1-\lambda$. If $m_0$ of the $m$ tests are null, their expected contribution to the numerator is $m_0(1-\lambda)$. If the alternative p-values also have some probability of exceeding $\lambda$, then, writing $P_A$ for a p-value under an alternative,

E\{\hat{\pi}_0(\lambda)\} = \frac{m_0(1-\lambda)+m_1P(P_A>\lambda)}{m(1-\lambda)} = \pi_0+\frac{m_1}{m}\frac{P(P_A>\lambda)}{1-\lambda},

where $m_1=m-m_0$ and $\pi_0=m_0/m$. Thus, alternatives with large p-values can add an upward contribution. At larger $\lambda$, that contribution may shrink if alternative p-values tend to be small. This motivates using a relatively high threshold. But as $\lambda$ approaches 1, fewer p-values remain in the numerator, making the estimate more variable.

For example, if the p-values were independent and all null, the count above $\lambda$ would have a Binomial$(m,1-\lambda)$ distribution. The estimator would then have expectation 1 and variance

\operatorname{Var}\{\hat{\pi}_0(\lambda)\} = \frac{m(1-\lambda)\lambda}{m^2(1-\lambda)^2} = \frac{\lambda}{m(1-\lambda)}.

This variance increases with $\lambda$. When alternatives are present, the count is generally not a single binomial count with that same success probability, but the basic warning remains: high thresholds leave fewer observations to inform the estimate. The grid is therefore a view of a bias–variability tradeoff, not a sequence that must rise or fall monotonically.

Worked Example: Reading a Lambda Grid

Suppose $m=1{,}000$ hypothetical tests produce the following counts. Each entry counts p-values strictly greater than the stated $\lambda$.

$\lambda$	Count above $\lambda$	$\hat{\pi}_0(\lambda)$
0.50	360	0.720
0.60	300	0.750
0.70	225	0.750
0.80	160	0.800
0.90	88	0.880

At $\lambda=0.80$, the estimate is $$ \hat{\pi}_0(0.80)=\frac{160}{1{,}000(1-0.80)} =\frac{160}{200}=0.800. $$ At $\lambda=0.90$, it is $88/[1{,}000(0.10)]=88/100=0.880$. The estimate rises at the high end, but this alone does not show that the true null proportion rises: the same data are being summarized at different thresholds, and the highest threshold uses only 88 p-values. The two equal estimates at 0.60 and 0.70 do not establish that this interval is optimal either.

How Estimate Changes Carry Through to Q-Values

For ordered p-values $p_{(1)}\leq\cdots\leq p_{(m)}$, the Storey q-value at rank $i$, using a selected estimate $\hat{\pi}_0$, is computed by taking the minimum adjusted value over ranks $j\geq i$, and capping the result at 1:

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

Holding the p-values and family fixed, this formula makes every uncapped q-value proportional to $\hat{\pi}_0$. A larger estimate cannot lower a q-value. Yet the practical impact depends on where the q-values lie relative to a chosen reporting threshold: a modest numerical change can leave a discovery set unchanged, while a similar change can move a borderline result across the threshold. These q-values are estimates, so a threshold comparison should not be described as certainty about the truth of an individual hypothesis.

Worked Example: A Lambda Choice Changes a Borderline Result

For a hypothetical family of $m=8$ tests, suppose the ordered p-values are $0.001, 0.004, 0.010, 0.020, 0.040, 0.080, 0.200, 0.700$. Compare estimates $\hat{\pi}_0=0.70$ and $\hat{\pi}_0=0.90$. At rank 5, the smallest $p_{(j)}/j$ for $j\geq5$ is $0.040/5=0.008$: the later ratios are $0.080/6=0.0133$, $0.200/7=0.0286$, and $0.700/8=0.0875$.

Therefore, $$ \hat{q}_{(5)}(0.70)=0.70(8)(0.008)=0.0448, \qquad \hat{q}_{(5)}(0.90)=0.90(8)(0.008)=0.0576. $$ At a q-value cutoff of 0.05, the rank-5 test is below the cutoff under the first estimate and above it under the second. The rank-4 q-value uses the minimum ratio $0.020/4=0.005$ among ranks $j\geq4$, giving $0.70(8)(0.005)=0.0280$ and $0.90(8)(0.005)=0.0360$. Thus the difference here is not a general change to every decision: it affects a result near the chosen threshold.

Bootstrap Selection and Smoothing

A bootstrap selection rule can summarize the grid with a data-dependent choice of $\lambda$. One common version first chooses a grid $\lambda_1,\ldots,\lambda_K$, calculates $\hat{\pi}_0(\lambda_k)$ at each point, and treats the estimate at the largest grid value as a reference. It then resamples the observed p-values with replacement, recomputes the estimates, and selects the grid value with the smallest bootstrap mean squared difference from that reference. The reference is a proxy, not the unknown $\pi_0$; if the largest-grid estimate is noisy or biased, the selection rule inherits that weakness.

Smoothing takes a different approach. A curve is fitted to the grid of estimates, and an estimate is read from the fitted curve, often at or near the upper end of the grid. A fitted curve can reduce jaggedness, but the result depends on choices such as the grid, spline flexibility, boundary behavior, and whether fitted values are constrained to $[0,1]$. Smoothing does not create information absent from a sparse upper tail.

Worked Example: Comparing a Bootstrap Choice with a Smoothed Estimate

Suppose a hypothetical analysis uses the grid $0.50,0.60,0.70,0.80,0.90$. Across 500 bootstrap resamples, the average squared differences from the estimate at $\lambda=0.90$ are $0.0072, 0.0038, 0.0016, 0.0021,$ and $0.0049$, respectively. Under the stated rule, the minimum is $0.0016$, so the selected grid value is $\lambda=0.70$. This selects a value based on estimated resampling error relative to the reference; it does not demonstrate that $0.70$ is closest to the unknown true $\pi_0$.

Now suppose a fitted smoothing curve evaluated at the upper grid endpoint gives $\hat{\pi}_0=0.82$, while the raw estimate at $\lambda=0.90$ is $0.88$. For the rank-5 p-value from the previous example, the smoothed estimate gives $$ \hat{q}_{(5)}=0.82(8)(0.008)=0.0525. $$ This is between the two earlier q-values, $0.0448$ and $0.0576$. It remains just above 0.05. The arithmetic illustrates how a smoothing choice can alter a borderline result; the smoother’s output is not automatically more trustworthy than either raw estimate.

Diagnostics, Comparisons, and Limitations

Assess stability by looking at both the $\hat{\pi}_0(\lambda)$ curve and the q-values or discovery set derived from it. Useful checks include:

Report the grid and the numerator count at each value, so readers can see how much tail information supports each estimate.
Look for abrupt changes, isolated spikes, estimates above 1 before truncation, or a sharp rise near the boundary. These can flag sparse-tail variability or a problematic choice of grid.
Compare raw grid estimates with the selected bootstrap value and the smoothed estimate. Record the bootstrap reference, resampling scheme, and smoothing settings.
Recalculate q-values and the thresholded discovery set under reasonable alternative choices. Focus on whether conclusions change, not only on whether the estimates look smooth.

Under independent sampling of p-values, an ordinary nonparametric bootstrap resamples the observed p-values as if they were exchangeable observations. In genomic applications, p-values can be dependent. Resampling individual p-values then need not reproduce the dependence structure, so bootstrap stability may be misleading. Even without dependence, a bootstrap sample cannot repair a misspecified null distribution, invalid p-values, or a poor choice of testing family.

The raw grid, bootstrap selection, and smoothing are alternatives for handling the tuning choice, not interchangeable guarantees. The raw grid is transparent and shows sensitivity directly, but leaves the user to summarize it. Bootstrap selection supplies a rule but depends on its reference estimate and resampling assumptions. Smoothing can stabilize a noisy curve but depends on model flexibility and boundary behavior. A practical analysis can report the grid and use a prespecified selection or smoothing rule, while showing whether the scientific conclusions are robust to plausible choices.

Common pitfall: “The smoothed curve is flat, so the estimate is stable and the q-values are reliable.” A smooth curve can result from the smoothing settings, and it does not establish valid p-values, a correct null model, or suitable resampling. Check the raw counts and estimates, compare reasonable analysis choices, and examine whether q-values or the discovery set change.

Sensitivity assessment has limits. A stable estimate across a chosen grid does not guarantee that the grid covers informative thresholds, that the null p-values are calibrated, or that the estimated q-values equal the unobservable true error rates. Conversely, visible variation does not by itself prove that any one estimate is wrong. It signals that the analysis should disclose the tuning dependence and be cautious about borderline findings.

Check Your Understanding

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

For $m=500$, 75 p-values exceed $\lambda=0.80$. Calculate $\hat{\pi}_0(0.80)$.
Why can an estimate at $\lambda=0.95$ be more variable than one at $\lambda=0.70$, even though it uses a higher threshold?
If $\hat{\pi}_0$ increases while the p-values and test family stay fixed, can a q-value decrease? Explain.
A bootstrap rule selects the grid point with the smallest mean squared difference from the estimate at the largest grid value. What important limitation does this reference introduce?
A smoothed estimate is stable, but raw estimates vary sharply at the upper end of the grid. Name one diagnostic or sensitivity check that should be reported.

Solutions. (1) $75/[500(1-0.80)]=75/100=0.750$. (2) The higher threshold leaves fewer p-values in the count; a smaller tail count generally gives greater sampling variability. (3) No. Each candidate adjusted value is proportional to $\hat{\pi}_0$, so increasing it cannot lower the minimum or its cap at 1. (4) The largest-grid estimate is only a proxy for the unknown $\pi_0$, and may itself be noisy or biased. (5) For example, report the raw counts and estimates across the grid, compare the smoothed result with the bootstrap-selected value, or recalculate the discovery set under reasonable alternative choices.