I have a number of continuous predictors (biomarker measurements) which I would like to test for association with a binary outcome variable (disease status), adjusting for multiple comparisons. As some of my predictors are correlated, I understand that procedures such as Benjamini-Hochberg might not be valid, so I would like to use permutation testing to adjust for multiple comparisons.

I did not manage to find an R code to do it, so I wrote it myself and used example data (see below). I wonder:

  1. whether the code is correct for calculating P values adjusted for multiple comparisons using permutation testing?
  2. if I did something wrong, as contrary to what I expected, P values adjusted for multiple comparisons using permutation testing are higher than e.g. Benjamini-Hochberg-adjusted P values?
  3. why P values adjusted for multiple comparisons using permutation testing do not always increase monotonically with (unadjusted) permuted P values?

Here is the code:

Prepare example data:



Ionosphere <- Ionosphere[-c(1:2)]  # Remove factor variables.

# [1] 351  33

# Columns 1-32: continuous independent variables.
# Column 33: binary dependent variable ('Class': bad/good).

Calculate original and permuted P values:

# Calculate original P values:

p <- sapply(
  function(x) {
    wilcox.test(x ~ Ionosphere$Class)$p.value

# Permute:

n.perm <- 999

perm.matr <- rbind(
  p,  # Original P values included as one of the permutations.
  matrix(, nrow = n.perm, ncol = length(p))



  for (i in 1:n.perm) {
    temp <- Ionosphere
    temp$Class <- sample(temp$Class)
    perm.matr[i + 1, ] <- sapply(
      function(x) {
        wilcox.test(x ~ temp$Class)$p.value
  #  user  system elapsed 
  # 76.14    0.06   76.35 


# Calculate permuted P values:

p.perm <- apply(
  function(x) {
    sum(x <= x[1]) / length(x)  # x[1] is the original P value

Calculate P values adjusted for multiple comparisons using permutation testing:

Based on this link.

p.adj <- as.numeric(rep(NA, length(p)))

names(p.adj) <- names(p)

for (i in 1:length(p)) {
  p.adj[i] <- sum(
    ) <= p[i]
  ) / nrow(perm.matr)


  p.adj ~ p.adjust(p, method = "BH"),
  main = "Permutation-adjusted vs. Benjamini-Hochberg-adjusted P"

Permutation-adjusted vs. Benjamini-Hochberg-adjusted P

An example where permutation-adjusted P values increase, but (unadjusted) permuted P values decrease:

This happens also with a higher number of permutations, e.g. n.perm = 9999.

This example also shows the discrepancy between the permutation-adjusted P values and Benjamini-Hochberg-adjusted P values.

  P = round(p, 4),
  P.permuted = p.perm,
  P.adjusted.perm = p.adj,
  P.adjusted.BH = round(p.adjust(p, method = "BH"), 4)
)[order(p), ][11:12, ]
#          P P.permuted P.adjusted.perm P.adjusted.BH
# V21 0.0031      0.006           0.075        0.0085
# V4  0.0032      0.003           0.076        0.0085

1 Answer 1


Not a full answer - p values 0.0031 and 0.0032 are very close. Even with 10,000 permutations you are looking at the difference between 31 vs 32 samples passing the observed cutoff. So it may be just by chance that you see an inconsistency between your method and Benjamini-Hochberg's?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.