Simulation Study Of BH And Storey Procedures

Why Simulate Multiple-Testing Procedures?

A procedure’s theoretical guarantee is conditional on assumptions about the p-values, their dependence, and the null hypotheses. Genomic tests can be dependent, and a guarantee under one dependence condition does not automatically extend to every correlation pattern. A simulation makes the design assumptions explicit and shows how procedures behave in repeated data sets with known truth. It is a comparison under chosen conditions, not a proof of performance in every genomic study.

The central quantities are calculated separately in each simulated data set. Their averages over many replicates describe performance under the specified data-generating mechanism. In particular, the false discovery rate is an expectation of the false discovery proportion; it is not generally the proportion of all tests that are null or the probability that any particular finding is false.

Specify the Simulation Before Comparing Procedures

A useful design records the true null status of every hypothesis, the distribution of each test statistic under the null and alternative, and the dependence structure among statistics. The null p-values should be marginally valid: for each true null and every \(u\) between 0 and 1, \(P(p_i\leq u)\leq u\). Under a continuous exact null, this probability equals \(u\). Otherwise, an apparent excess of false discoveries may reflect invalid p-values rather than a failure of the multiple-testing rule.

For a basic mixture design, choose \(m\) tests and a known null proportion \(\pi_0=m_0/m\), where \(m_0\) is the number of true nulls. For example, null p-values can be uniform and alternative p-values can follow a distribution concentrated near zero, such as a Beta distribution with parameters \(a<1\) and 1. This specifies the marginal distributions, but not their dependence: that must be specified separately.

Use the same simulated p-values for all procedures in each replicate. This paired comparison prevents differences between procedures from being confounded with differences in the generated data. The procedures being compared are:

• BH: sort the p-values \(p_{(1)}\leq\cdots\leq p_{(m)}\), and find the largest \(k\) such that \(p_{(k)}\leq \alpha k/m\). Reject the hypotheses with ranks 1 through \(k\); if no such \(k\) exists, make no rejections.
• Storey q-values: estimate \(\pi_0\) using a preselected tuning value \(\lambda\), then calculate monotone q-values from the ordered p-values. In this tutorial, \(\hat{\pi}_0(\lambda)=\min\{1,\#\{p_i>\lambda\}/[m(1-\lambda)]\}\), and the ordered q-value is \(q_{(i)}=\min\{1,\min_{j\geq i}[\hat{\pi}_0m p_{(j)}/j]\}\). Reject when \(q_i\leq\alpha\). This is one explicit estimator and implementation choice; results can depend on that choice.
• BY: let \(H_m=\sum_{j=1}^m 1/j\). Find the largest \(k\) such that \(p_{(k)}\leq \alpha k/(mH_m)\), and reject ranks 1 through \(k\). This more conservative threshold is designed to control FDR under arbitrary dependence when the p-values are valid.

BH has FDR control under independence and under the PRDS conditions covered in the earlier material on dependence and FDR control. BY offers a guarantee under arbitrary dependence, at the cost of a more stringent threshold. Storey’s method estimates the null proportion and can gain power when \(\pi_0\) is below 1, but its performance depends on the estimator and assumptions. A simulation can compare observed behavior; it should not be used to claim a guarantee outside a procedure’s conditions.

Dependence Patterns and a Fully Specified Design

The simulation should vary dependence while retaining valid marginal null p-values and the same underlying truth pattern. One practical design generates one-sided Gaussian test statistics. Under a true null, \(Z_i\sim N(0,1)\) and \(p_i=1-\Phi(Z_i)\), where \(\Phi\) is the standard normal cumulative distribution function. Under an alternative, use \(Z_i\sim N(\delta,1)\) with \(\delta>0\); larger values tend to produce smaller p-values. Fix which tests are null and which are alternative in advance, or randomly assign statuses according to a stated mixture rule and retain each replicate’s truth labels for scoring.

For independent statistics, generate each \(Z_i\) independently. For positively dependent blocks, generate each block from a multivariate normal distribution with unit marginal variances and positive within-block correlations. For a negative-dependence comparison, pair statistics with correlation \(-0.60\). A two-variable correlation matrix with this value has eigenvalues \(1-0.60=0.40\) and \(1+0.60=1.60\), both positive, so it is a valid covariance matrix. These are illustrative dependence patterns, not models that represent every linkage structure.

Each null statistic remains standard normal in these designs, so its one-sided p-value is uniform marginally, even when statistics are correlated. That checks marginal null validity; it does not establish PRDS. In particular, do not infer that a positive correlation matrix automatically establishes the exact dependence conditions required for a theoretical guarantee in every test construction.

Summarize, Diagnose, and Compare

A clear report includes the number of tests, null proportion, alternative effect distribution, dependence structure, significance level, number of replicates, Storey estimator and tuning value, and the precise definition of power. Report mean FDP as the empirical FDR estimate and mean power, together with Monte Carlo standard errors or intervals. If the number of rejections matters for the scientific use case, also summarize its distribution; equal mean FDPs can conceal substantially different discovery counts.

The following hypothetical results illustrate how dependence comparisons might be presented. Each cell reports empirical mean FDP / mean power for a specified design. They are not universal expected values.

In this hypothetical comparison, the first two designs show similar empirical FDP for BH and Storey, while the negative-pair design shows higher empirical FDP for both. BY is more conservative in all three. Interpretation depends on Monte Carlo error: differences similar in size to their uncertainty should not be described as clear evidence of different performance.

Investigate unexpected results by checking the generated marginal null p-value distributions, the realized null proportion, the correlation matrices, and the code that assigns truth labels. Plot or summarize replicate-level FDP and rejection counts, not only their means; a small number of replicates with large FDP can affect the average. Run enough replicates for the desired Monte Carlo precision, and repeat key comparisons across plausible effect sizes, null proportions, and dependence strengths.

Common alternatives include comparing procedures using analytic guarantees, or using permutations and resampling schemes to model dependence from observed data. Guarantees answer what can be established under assumptions; simulations answer what happens under the specific model implemented. A model-based simulation is easier to control and interpret, but may not capture features such as heterogeneous effects, imperfect calibration, or complex genomic dependence. A favorable simulated result cannot repair invalid marginal p-values, and a single simulated dependence pattern cannot stand in for all arbitrary dependence.