# Permutation tests for paired data with several variables

I am trying to answer my own previous question, with a permutation method.

I have a reference method $$M_{ref}$$ that I would like to compare with other methods $$M_1$$, $$M_2$$, $$M_3$$ when taking measurements (denoted scores from now on) from a data sample of size $$N (N > 30)$$. However, when comparing $$M_{ref}$$ with a method $$M_i$$, i do not have one kind of continuous score $$S_j$$, but several, say $$S_1$$, $$S_2$$, $$S_3$$, and $$S_4$$. So for each method, I have $$N$$ scores of type $$S1$$, $$N$$ scores of type $$S2$$, etc.

My question (still) is: how to obtain a non-parametric test p-value for each comparison $$C_{i,j} = _{score=C_j}$$. Please note I am not interested in any comparison between methods other than the reference method (e.g. $$M_{2}$$ Vs $$M_3$$ is of no interest), and comparisons between different kind of scores would not make sense.

So now i am trying to think of a way to test all the scores at the same time with permutation tests for a given method comparison, e.g. $$M_{ref}$$ Vs $$M1$$, and then simply use Bonferonni correction by multiplying by the number of times I compare methods this way (so by 4, because there are 4 methods comparison to do: $$M_{ref}$$ Vs $$M1$$, $$M_{ref}$$ Vs $$M2$$, etc.).

Now I am not sure of the proper way to test all the scores at the same time with permutation tests. Would the following (inspired from here) work?

    # Number of random reps
nreps <- 10000
# All these operations are realized across score columns,
# but randomly assigning sign the same way for all columns
# [is that the main trick?]
# Now we will randomly assign the sign of the difference
Mref.vs.M1.S1.to.S4 <- Mref.S1.to.S4 - M1.S1.to.S4
sample.means <- mean(Mref.vs.M1.S1.to.S4)
cat("Mean differences for S1, S2, S3 and S4 are", sample.means, "\n")
unsigned.diff <- abs(Mref.vs.M1.S1.to.S4)
n <- length(unsigned.diff)
for (i in 1:nreps) {
# Draw N obs randomly from 1, -1
signs <- sample(c(1, -1), size = n, replace = TRUE)
signed.diff <- unsigned.diff*signs
# Calculate the mean with that arrangement of signs
# for S1, S2, S3 and S4
means.random[i] <- mean(signed.diff)
}

# The question is what percentage of those differences
# are greater than or equal to plus or minus sample.means
# of S1, S2, S3 and S4
percent.above <- length(means.random[means.random >= abs(sample.measn)])/nreps
percent.below <- length(means.random[means.random <= abs(sample.means)*(-1)])/nreps
percent.extreme <- percent.above + percent.below

cat("The randomization test gives a probability under the null")
cat("for S1, S2, S3 and S4 [4 p-values] of", percent.extreme, "\n\n")


Does this work or is it just equivalent to doing resampling multiple times for each score type, and the resulting p-values still need to be Bonferonni-corrected by multiplying by the number of score types (4 in this case) ?