2
$\begingroup$

I am doing a within-subject design experiment and trying to analyze my data. To simplify my problem, let's say that there is condition A and condition B, and I have 24 pariticipants. I want to do permutation statistics in order to test if there is a difference between condition A and B. From what I understood, there is 2^24 possible permutations, as we compute either A(i) - B(i) or B(i) - A(i) for each participant. As this is not possible to do all the permutation, I only compute N=2000 permutations but then, if I run multiple times my code, since there is always a different permutation, sometimes I can get a significant results p < 0.05 and sometimes a non significant results. This is annoying as I don't want to do "p-hacking".

My question is then : How to deal with permutation statistics yielding different results when we only do a part of possible permutations ?

I tried to read different explanation of permutation statistics and some says that 5000 permutations "is enough" but that does not solve the problem of non-determinism...

I hope that I am clear, feel free to ask for more details if needed.

Thank you !

$\endgroup$
4
  • $\begingroup$ Put set.seed(2021)(R) or np.random.seed(2021) (Python) at the top of your code; I'm sure other software packages like SAS, Stata, SPSS, and Matlab have similar commands. Then it will give the same permutations every time. As long as you don't tune the seed to give you the desired result, you are fine. (Perhaps even consider using your birthday or lucky number as the seed, so it doesn't change every year.) $\endgroup$
    – Dave
    Apr 26, 2021 at 13:46
  • $\begingroup$ Thank you for your answer. I was not sure if fixing the seed would be valid for a scientific analysis, as it involves randomization, but now that you say it, it makes sense! $\endgroup$
    – Tehelka
    Apr 26, 2021 at 13:59
  • $\begingroup$ 2000 values is pretty small (for all that it seems to be a common default). It made sense in the late 60s, but not at any time since then. It gives a margin of error of 0.01 (about two figures of accuracy for a p-value near 0.05) -- you can compute an estimate of the standard error of your estimated p-value easily (it's a binomial proportion) and for any application you can decide how much margin of error you can tolerate in your sample-estimate of the p-value. I tend to use 100,000 to 1,000,000 resamples (If I have to wait a minute, so what?) because I seek nearer to three figure accuracy ... $\endgroup$
    – Glen_b
    Dec 17, 2022 at 0:27
  • $\begingroup$ ctd ... If I have one that (for some reason) runs very slow I may do 10,000 ... but I might also then set it going overnight, since I rarely need a more precise result instantly. You can always consider saying that if your p-value is anywhere close to being within the margin of error (say within twice the margin of error) of your significance level, you should simulate much more, as long as the protocol is clearly defined before the first simulation is run so there's no danger of people thinking you had done so because you didn't like the outcome. $\endgroup$
    – Glen_b
    Dec 17, 2022 at 0:27

1 Answer 1

2
$\begingroup$

If you're down to setting seeds to get reproducible results from permutation tests, something may be wrong. Maybe the sample size (e.g., your 24) is too small, you're not using enough permutations (e.g., your 2000), the metric does not take enough uniquely different values, or there is a programming error. Obviously, a simulated permutation test is subject to simulation error, but in my experience, running such a test several times with different seeds (or no seed) gives substantially similar results on each run.

Let's take a look at the kind of simulated permutation test you are using---for fictitious non-normal data---to see how my version of it works in R.

Suppose we have $n = 24$ differences $B - A$ as in my vector dif below: The observed average difference is $1.67.$ Also, descriptive graphics show the differences are right skewed. Thus they are not excellent candidates for testing either with a one-sample t test or a one-sample nonparametric Wilcoxon signed=rank test. (But, for the record, both tests do reject that the differences are centered at $0$, the Wicoxon test with a warning about ties.)

ad.obs = mean(dif); ad.obs
[1] 1.666667

enter image description here

R code for figure:

par(mfrow = c(1,3))
 hist(dif, prob=T, col="skyblue2"); rug(dif)
 boxplot(dif, col="skyblue2")
 qqnorm(dif, pch=19); qqline(dif, col="green2", lwd=2)
par(mfrow = c(1,1))

A one-sided permutation test that averages the differences after randomly switching the 24 signs shows a significant difference with P-value about $1.5\%,$

set.seed(2022)
ad.prm = replicate(5000, mean(dif*sample(c(-1,1),24,rep=T)))
mean(ad.prm > mean(dif))
[1] 0.0142  #  P-value of simulated permutation test

enter image description here

hist(ad.prm, prob=T, br=30, col="skyblue2")
 abline(v=ad.obs, col="red", lwd=2, lty="dotted")
length(unique(ad.prm))
[1] 3838

Two additional runs with different seeds and different numbers of permuted means give essentially the same result.

set.seed(1776)
ad.prm = replicate(10^5, mean(dif*sample(c(-1,1),24,rep=T)))
mean(ad.prm > mean(dif))
[1] 0.01146

length(unique(ad.prm))
[1] 15437

set.seed(14)
ad.prm = replicate(5000, mean(dif*sample(c(-1,1),24,rep=T)))
[1] 0.0122
length(unique(ad.prm))
[1] 3827

Note: Here is R code used to sample the fictitious differences used above.

set.seed(426)
dif = round(rexp(24,1/5)-rexp(24, 1/2), 2)
$\endgroup$
7
  • 1
    $\begingroup$ Thank you for the answer. Maybe I don't understand it fully but isn't there an error in the permutation created ? sample(...) creates a vector of -1 or 1 to say if we invert or not each differences, but you take the mean of this vector of -1/1. For me it should be mean(dif * sample(...)) In this situation, take for instance set.seed(14) and 5000 permutations and you have a case when the tests fails : here is a screenshot Did i miss something ? $\endgroup$
    – Tehelka
    Apr 27, 2021 at 10:27
  • 1
    $\begingroup$ Thanks for comment. Obviously, you have followed my answer carefully! You are right about the typo in the code (now corrected). I got P-values about 0.011 to 0.014 in three runs with different seeds and nrs of replications. (I used your seed and number of replications in the last run now posted.) // But I can't replicate your non-significant result or histogram. Wondering if you may have used different data. $\endgroup$
    – BruceET
    Apr 27, 2021 at 16:03
  • $\begingroup$ Oh indeed, my bad, I was using the seed before generating the data and not changing it after permutation. With your seed for the data (426) I cannot replicate my non-significant result, but if you try with set.seed(14) dif = ... ad.rpm = ... directly, you get my previous histogram. I wonder why is that, but that "proves" that my behavior can still happen, or did I do something wrong? Thanks again ! $\endgroup$
    – Tehelka
    Apr 28, 2021 at 6:19
  • $\begingroup$ Here is my code that do that just in case : here $\endgroup$
    – Tehelka
    Apr 28, 2021 at 6:30
  • 1
    $\begingroup$ If you had normal data with $\sigma=2$ trying to use a t-test to detect a real difference of $\Delta=2,$ with only $n=24$ observations, then you would have power only about 65%. That is, 65% of samples of size 24 will give P-value below 5% and 35% will not. Similarly, if you generate data using my algorithm, a true permutation test (with combinatorics, not simu) will reject sometimes and sometimes not. Sim perm test will act similarly. Not mainly because perm test is (a bit) random, but because expts with 24 obs are underpowered. // BTW, my last msg should have 'perm test', not 'bootstrap'. $\endgroup$
    – BruceET
    Apr 28, 2021 at 13:43

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.