# Adjusting P values for multiple comparisons using permutation tests

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:

library(mlbench)

data(Ionosphere)

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

dim(Ionosphere)
# [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(
Ionosphere[-ncol(Ionosphere)],
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))
)

set.seed(123)

system.time(

for (i in 1:n.perm) {

temp <- Ionosphere

temp$$Class <- sample(temp$$Class)

perm.matr[i + 1, ] <- sapply(
temp[-ncol(temp)],
function(x) {
wilcox.test(x ~ temp$$Class)$$p.value
}
)

rm(temp)

}

)
#  user  system elapsed
# 76.14    0.06   76.35

rm(i)

# Calculate permuted P values:

p.perm <- apply(
perm.matr,
2,
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:

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

for (i in 1:length(p)) {

apply(
perm.matr,
1,
min
) <= p[i]
) / nrow(perm.matr)

}

rm(i)

plot(
)


## 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.

data.frame(
P = round(p, 4),
P.permuted = p.perm,