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
  • $\begingroup$ Don't proceed without doing a sample size calculation. I fear that your sample is far too small to get reliable results, since the outcome is a minimum-information binary response. And note that large P-values mean nothing more than "get more data". Also I would take a different approach. Compute the first 10 principal components and get a perfectly multiplicity-adjusted "chunk" test with 10 d.f. relating the 10 PC scores to the outcome. $\endgroup$ Jan 6 at 9:26

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 and acknowledge you have read our privacy policy.

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