1

I asked this question on cross-validate, but was asked to move it here because it involved coding more than a statistics question.

I'm teaching a statistics class about family-wise error rate, and have created a graph to illustrate how, when we increase the # of comparisons, the rate of type I errors increases rapidly (Family-Wise Error Rate, FWER). I then show how the bonferroni correction keeps the FWER near 0.05, or 5%.

I would like to also add the Holm correction to the graph, but am struggling at figuring out the code in R to do this. I see this example here, but don't quite get how the calculation works for the graph. In case this helps, Bonferroni corrects all significance thresholds by the same value (the number of tests, hence my code BON <- a/pwhere a=alpha=threshopld, p=# of comparisons to make), while Holm's method modifies the threshold sequentially from lowest (raw) p-value to highest. I cannot figure this part (adding Holm) of the code out.

Adding other correction methods is also open (and encouraged)! The more, the better, that way I could make a true all-inclusive graph (this would be really cool). But adding Holm correction would already be really great.

Here is what I have done so far (assuming we have a maximum of 100 comparisons to make)

par(mfrow = c(1, 1), cex=1, pch=10)

p <- seq(1,100)
a <- 0.05 # alpha
FWER <- 1-(1-a)^p
BON <- a/p

plot(p,FWER, xlim = c(1,100), ylim = c(0,1), type ='n', main = '', ylab="Family Wise Error Rate [FWER]", xlab='Number of Comparisons')
points(p,BON,pch = 20,col="red")
points(p,FWER,pch=20, col="blue")
legend(-.1, .95, legend=c("No Correction", "Bonferroni Correction"),
       col=c("blue", "red"), pch=20, cex=2)

enter image description here

The next point to add, if any, would be holm, and would something like this

points(p, HOLM, pch=20, col="green")
Andy
  • 413
  • 2
  • 15
  • Would something like `points(p, p.adjust(p, "holm"), pch=20, col="green")` work for you? – Vincent Guillemot Apr 07 '22 at 06:11
  • hello @VincentGuillemot, I tried this, and that would be cool, but the holm line appears at 1 (soo 100% error rate) from 1:100. So, no it does not, unfortunately, but it is a good start, and I could play with it to see if I could get it to work. – Andy Apr 07 '22 at 06:19
  • I guess this is more of a stats question, but does it even make sense to calculate the FWER **after** performing a multiple-testing correction? – Dan Adams Apr 07 '22 at 06:30
  • Hi @DanAdams : yes and no, the point I am trying to graphically illustrate is that the FWER stays low thanks to the correction. It's easy to say this, but when the students see it, it really helps get the point across. – Andy Apr 07 '22 at 06:38
  • I totally get why you want to show it and how powerful it would be. I'm just not sure it makes sense computationally. I.e. you might have to simulate the reduced error rate some other way because `p.adjust()` expects p values, but I think `p` in your example is the number of comparisons. – Dan Adams Apr 07 '22 at 06:40
  • ```p```is the number of comparisons. I just tried this, and it seems to work actually, but I thought the curve should go up sharply as we reach a higher number of comparisons. First, create vector of p values ```h <- seq(0.01:1, by=0.01) # add sequence of p-values, from 0.01 to 1``` Fill in holm equation: Holm correction equation alpha = alpha / m-j+ 1 , where m=comparisons, j = ordered p values, and + 1 ```scaled_pvalues <- p-h+1 holm <- a/scaled_pvalues ``` This gives me a line that matches, almost 100% the bonferroni curve – Andy Apr 07 '22 at 06:49
  • If you solved it, you can add an answer to help future readers and accept it. – Dan Adams Apr 07 '22 at 06:53
  • 1
    @DanAdams I added an answer, but again, would be interesting if others validate (or not) the graph :) – Andy Apr 07 '22 at 07:00
  • OK, my bad, I read too fast: indeed `p.adjust` will not help you much – Vincent Guillemot Apr 07 '22 at 14:31

1 Answers1

1

I tried this, but want to hear from others if this is truly correct. I have seen Holm curves where the holm correction sharply increases the FWER rate when high number of comparisons are made, but my calculations are not coming up with this answer. Either this is correct (please validate if so!) & the others I've seen are not correct, or vice versa

par(mfrow = c(1, 1), cex=1, pch=10) # set plot window to desired conditions

p <- seq(1,100)
a <- 0.05 # alpha
h <- seq(0.01:1, by=0.01) # vector of scaled p values
FWER <- 1-(1-a)^p
BON <- a/p

#  Holm correction equation alpha = alpha / m-j+ 1 , where m=comparisons, j = ordered p values, and + 1

scaled_pvalues <- p-h+1

holm <- a/scaled_pvalues


plot(p,FWER, xlim = c(1,100), ylim = c(0,1), type ='n', main = '', ylab="Family Wise Error Rate [FWER]", xlab='Number of Comparisons')
points(p,BON,pch = 20,col="red", cex=2)
points(p,FWER,pch=20, col="blue")
points(p, holm, pch=20, col="green",cex=1)
legend(-.1, .95, legend=c("No Correction", "Bonferroni Correction", "Holm Correction"),
       col=c("blue", "red", "Green"), pch=20, cex=2)

enter image description here

Andy
  • 413
  • 2
  • 15