# Power analysis for correlations, controlling for multiple comparisons

I want to estimate number of participants required for an experiment. I will be conducting 28 correlation tests. I am seeking 80% power for an effect size of 0.15 and alpha of 0.05 (all two-tailed).

Given the number of tests, I want to apply a multiple test correction and still achieve the desired power.

Doing that for Bonferroni is straightforward enough, I believe: my adjusted p-value would be 0.05/28 = 0.0018. Using this p-value I can now estimate that for the desired paramters above; the sample size I need is 691 (calculated using GPower 3.1).

Now, because the Bonferroni adjustment is quite restrictive, I would prefer to use the Benjamini-Hochberg one, but I am unsure how to apply it. After looking through a few papers with complicated computations, I stumbled upon this R package, pwrFDR (see here), that simplifies the job.

I used the following specification for the main pwrFDR function, attempting to get the sample size.

pwrFDR(effect.size=0.15, average.power=0.80, N.tests=28, r.1 = 10/28,alpha=0.05, FDP.control.method='BHFDR')


The suggest sample size from this is 959. Even if I use the highest (and unrealistic) possible r.1 (proportion of significant tests, if I understand correctly), say 27/28, I still get an estimated sample size of 755, which is higher the one calculated using the Bonferroni procedure.

Given that the Bonferroni method has been shown to have significantly lower power than the Benjamini-Hochberg correction, why do the above estimates show that I need fewer participants with the Bonferroni correction? I think it's likely that I am misunderstanding the pwrFDR function specifically, and power analysis using the BH procedure in general.

• If these correlations concern all possible pairs of eight variables, then neither of your approaches is appropriate due to the strong, complex dependencies among the tests. One reliable way to obtain the power curve is through simulation. Exactly how to do that depends on the specific correlation tests you are conducting.
– whuber
Commented Oct 24, 2021 at 13:23
• @whuber you are correct in that I am correlating all pairs of 8 variables but I cannot figure out why that would be an issue. I appreciate that 1) the amount of tests inflates Type I errors (hence my desire to use an adjustment procedure) and 2) that I might very well be capturing substantial noise, which is a limitation I am willing to endure as this is supposed to be exploratory... And I am not sure I understand why the specific type of correlation test would make a difference: I am still testing a hypothesis (that the correlation is different from 0). Commented Oct 25, 2021 at 5:53
• One simple-minded but effective way to think of the situation is to note that only 8 variables are involved but you estimate 28 correlation coefficients. Those coefficients are thereby pretty strongly interdependent. How, exactly, do you propose to adjust the 28 resulting p-values? Treating them as 28 independent values would give quite the wrong answer--but so would treating them as only 8 independent values. How to adjust them depends on the dependence structure; and that in turn relates to how you estimate the coefficients.
– whuber
Commented Oct 25, 2021 at 13:24
• Correlated tests are a concern. Consider it this way, if correlation was 1 then there would not be 28 multiple tests but 1. For n correlated tests with correlations between -1 and 1, the "effective" number of independent tests is less than n. The point is that adjusting for n simultaneous tests will be too conservative. The latest version of pwrFDR has a provision for correlated variables. In addition to calculating average- or TPX- power under BH-FDR, it now offers these power calculations for control under the methods of Bonferroni, Holm, Hochberg (FWER) and Lehmann-Romano (FDX control). Commented Jul 25 at 9:29

@npetrov937 Sorry it took me three years to respond. I would have seen this much sooner had you submitted an issue.

So it appears that you're comparing average power under Bonferroni with average power under BH-FDR. First let's look at your Bonferonni calculation. I ran the Bonferroni power by hand. At alpha=0.05, and effect size 0.15, the Bonferroni corrected test has average power 79.3% at a sample size of 691

1-pnorm(qnorm(1-0.05/28/2)-(691)^0.5*0.15)


(The average power for any test with a constant threshold is the single test power provided effect sizes are identical)

This is a two sided one sample test assuming normally distributed data. Close enough. But in order to have 80% average power you need 700

1-pnorm(qnorm(1-0.05/28/2)-(700)^0.5*0.15)


Notice that this average power is calculated assuming an identical effect size, e.g. r.1=1.0

Now, some clarifications are in order.

The default settings for pwrFDR are t-distributed data and a 2-group comparison.

The latter of course makes the biggest difference. You can set

control=list(distopt=0)


to get the normally distributed assumption which will be slightly less conservative at 28 tests. You can set r.1=1-1e-8 which may reduce the sample size by 1, BUT the clincher, as you've guessed by now, is the ONE group comparison. You've opted for the default which is a two group comparison (far more common in practice). So, comparing apples to apples, here is the sample size required for average power in a one group comparison at alpha=0.05 under BHFDR, at effect size 0.15 I even put my thumb on the scale in favor of Bonferroni a teeny bit by assuming t-distributed tests with r.1 set as you had in your second try.

pwrFDR(effect.size = 0.15, r.1 = 27/28, alpha = 0.05, groups = 1, N.tests = 28, average.power = 0.8)

n.sample      379
average.power 0.7999168
c.g           2.075293
gamma         0.7727283
err.III       3.221475e-07
se.Rom        0.09279357
se.VoR        0.009078492
se.ToM        0.09017978


At your more realistic guesstimate for proportion of true alternatives 10/28 you get a bit larger required sample size but still roughly 2/3 of that which is required for Bonferroni:

pwrFDR(effect.size = 0.15, r.1 = 10/28, alpha = 0.05, groups = 1,
N.tests = 28, average.power = 0.8)

n.sample      482
average.power 0.7998492
c.g           2.447131
gamma         0.2951473
err.III       5.399081e-09
se.Rom        0.1020887
se.VoR        0.06206888
se.ToM        0.1481223


So would you be interested in helping to make the package more intuitive for people to use? You already have made several suggestions here but I think you could easily provide more. You'll find my contact info in the package maintainer block in pwrFDR.

Best Regards,
Grant Izmirlian
Division of Cancer Prevention
NCI

• It is seldom interesting to test whether a true correlation is zero. But a variety of simulations can be used to determine the adequacy of a sample size for studying a variety of complex questions about correlation matrices as exemplified here. Commented Aug 25 at 12:38
• Simulation is also an option, and in fact available as an option in the package. However, it becomes impractical for deriving sample size in a sequence of runs required for a solver to converge. Additionally the package is a collection of cutting edge sequential procedures with a variety of definitions for power. Not to mention new methodology provide which allows control of the FDX with only slightly more conservativism than BH-FDR. Maybe you should give it a look, Frank. Commented Aug 26 at 17:13
• Not a believer in FDR unless FNR is highly emphasized while using it. Commented Aug 26 at 22:02
• pwrFDR computes average power or TPX power under controlled FDR or FDX. FDX can be controlled via Lehmann-Romano, slightly less conservative than Holm's procedure or my new procedure which when a solution exists, is only slightly more conservative than the BH-FDR procedure. So, you see, it's all about the FNR and not exclusively about about the FDR Commented Aug 27 at 23:20
• For most applications the power of an FDR-controlled analysis will be to low to trust non-discoveries. Commented Aug 28 at 13:03