Search code examples
rstatisticsbonferroni

Why am I getting the same result when using p.adjust for the Holm and Bonferroni methods?


I've run several regressions and I'm trying to see if any of the significant findings hold after a multiple comparison adjustment. I know that Bonferroni adjustment is more conservative than Holm adjustment, so I would expect the outputs to be different, however they weren't. One of the values varied slightly, but not as much as I would have expected given that one of these tests is supposed to be less conservative than the other. I tried another correction method and got similar results. Is there something wrong with my syntax or is this a normal result?

p.vect <- c(.003125, .008947)
p.adjust(p.vect, method = "bonferroni", n=80)
[1] 0.25000 0.71576

p.adjust(p.vect, method = "holm", n=80)
[1] 0.250000 0.706813

p.adjust(p.vect, method = "hochberg", n = 80)
[1] 0.250000 0.706813

Solution

  • Holm and Hochberg just don't differ from each other for length(p)==2.

    Given that lp is length(na.omit(p)) (equal to 2 in this case), and p is the vector of probabilities, here's the code for method=="holm":

    i <- seq_len(lp)   ## (1,2)
    o <- order(p)      ## (1,2)
    ro <- order(o)     ## (1,2) 
    pmin(1, cummax((n + 1L - i) * p[o]))[ro]  ## c(80,79)*p
    

    and the code for method=="hochberg":

    i <- lp:1L                        ## (2,1) 
    o <- order(p, decreasing = TRUE)  ## (2,1)
    ro <- order(o)                    ## (2,1) 
    pmin(1, cummin((n + 1L - i) * p[o]))[ro]  ## c(80,79)*p[c(2,1)][c(2,1)]
    

    If you step through the details you can see how they give the same answer for your case. (Holm uses cummax() on the sorted vector of adjusted probabilities, Hochberg uses cummin() on the reverse-sorted vector, but neither of these changes anything in this case.)

    Bonferroni is pmin(1, n*p). In this case the only difference is a factor of 80/79 on the second element (Hochberg and Holm multiply by (n+1-i) = c(80,79), Bonferroni multiplies by n=79.)

    You can print out the code by typing p.adjust by itself, or in a prettier form here