R code to apply bonferroni-holm adjustment for multiple comparisons

bonferroni-holm.R

bonferroni.holm <- function(pValues, alpha=0.05, sidak=FALSE) {
  pOrders <- order(pValues)
  unadjusted.p <- pValues[pOrders]
  undoOrder = order(pOrders)
  m <- rev(rank(unadjusted.p, ties.method = "min"))
  if (sidak) {
    adjusted.alpha <- 1 - (1 - alpha) ^ (1 / m)
    adjusted.p <- 1 - (1 - unadjusted.p) ^ m
  } else {
    adjusted.alpha <- alpha / m
    adjusted.p <- pmin(unadjusted.p * m, 1.0)
  }
  adjusted.p <- Reduce(max, adjusted.p, accumulate = TRUE)

  returnVar <- data.frame(unadjusted.p, adjusted.alpha, adjusted.p)
  returnVar$significant <- returnVar$unadjusted.p < returnVar$adjusted.alpha
  whichNonSig <- which(!returnVar$significant)
  if (length(whichNonSig) > 0) {
    nonSigIndicies <- whichNonSig[1]:nrow(returnVar)
    returnVar$significant[nonSigIndicies] <- FALSE
    returnVar$adjusted.alpha[nonSigIndicies[-1]] <- NA
  }
  return(returnVar[undoOrder,])
}

example-use.R

> source("~/bonferroni-holm.R")
> b <- bonferroni.holm(c(0.01, 0.04, 0.03, 0.005))
> b
  unadjusted.p adjusted.alpha adjusted.p significant
2        0.010     0.01666667       0.03        TRUE
4        0.040             NA       0.06       FALSE
3        0.030     0.02500000       0.06       FALSE
1        0.005     0.01250000       0.02        TRUE
> b <- bonferroni.holm(c(0.01, 0.04, 0.03, 0.005), sidak = TRUE)
> b
  unadjusted.p adjusted.alpha adjusted.p significant
2        0.010     0.01695243  0.0297010        TRUE
4        0.040             NA  0.0591000       FALSE
3        0.030     0.02532057  0.0591000       FALSE
1        0.005     0.01274146  0.0198505        TRUE