Appropriate practice for setting seeds in permutation analyses? I am currently using permutation ANCOVA for a study of mine, in R, using the permuco package. Running multiple models for different brain regions of interest (ROIs). What I am finding is that the p-values generated change subtly, depending on the seed I use. While it is not by much, some p-values skirt around the 0.05 cutoff, which reduces my confidence in some of the findings.
So far, I have fixed the seed for every permutation ANCOVA, using a value pulled from random.org, using the "random" package (an integer between 1 and 1,000,000). What I am curious is if it is generally "frowned upon" to set fixed seed values for tests that utilize some form of PRNGs, or is this generally well-accepted, for reproducibility purposes?
 A: Several comments:

*

*There is nothing special about the 5% level (except possibly in the minds
of journal editors). Inaccuracies in simulation on either side of 5% are no
more or less common than for other significance levels.


*If you are doing enough iterations for the permutation test, any seed should
give sufficient accuracy. Maybe more iterations than you usually use: $100\,000$
should give two or maybe three places of accuracy in the P-value; a million should give three.


*If you do two runs with two different seeds and get remarkably different answers,
then you should investigate whether your understanding of that particular permutation test is correct and whether your programming is correct.


*It is best to set a particular seed at the start of any procedure that depends
on the pseudo-random generator, then the same software should always give exactly
the same simulation result. (Keep a record of the seed so you can double-check
the result later if necessary.)


*For permutation tests with small datasets, it is a good idea to check how many
uniquely different values of the metric you have in the simulated permutation
distribution. If the metric involves an order statistic (e.g., median), then the
number of uniquely different results may be small, so that the rejection probability
is more discrete than you may realize.
Example:  Consider a test for the difference in population means of two gamma populations. A Welch t test finds a difference at the 5% level. However, data are far enough from normal that a t test might not be reliable:
set.seed(123)
x1 = rgamma(20, 3, .12)
x2 = rgamma(20, 4, .1)
t.test(x1,x2)

        Welch Two Sample t-test

data:  x1 and x2
t = -2.1976, df = 37.865, p-value = 0.03417
alternative hypothesis: 
  true difference in means is not equal to 0
95 percent confidence interval:
 -17.6586291  -0.7232855
sample estimates:
mean of x mean of y 
 22.76063  31.95159 

A permutation test with difference in sample means as the metric seems safer,
and still finds a difference:
set.seed(2020)
x = c(x1,x2)
da.obs = mean(x2)-mean(x1)
m = 10^6;  da.prm = numeric(m)
for(i in 1:m) {
  x.prm = sample(x)
  da.prm[i] = mean(x.prm[1:20]) - mean(x.prm[21:40]) }
mean(abs(da.prm) >+ abs(da.obs))
[1] 0.034531     # P-value about same as Welch t
2*sd(abs(da.prm) >+ abs(da.obs))/100
[1] 0.003651775  # aprx 95% margin of sim error for P-value

However, using difference in sample medians as metric is not a good choice:
set.seed(2020)
x = c(x1,x2)
dh.obs = median(x2)-median(x1)
m = 10^6;  dh.prm = numeric(m)
for(i in 1:m) {
 x.prm = sample(x)
 dh.prm[i] = median(x.prm[1:20]) - median(x.prm[21:40]) }
mean(abs(dh.prm) >= abs(dh.obs))
[1] 0.090859
2*sd(abs(dh.prm) >= abs(dh.obs))/1000
[1] 0.0005748173
length(unique(dh.prm))
[1] 3612

Not nearly significant at 5% level. Among a million iterations there are
fewer than 4000 uniquely different values of this metric, maybe sparse in the tails. So the P-value
is essentially discrete.
Using a different seed (selected by R), we get about the same P-value, so at
least the non-significant result is stable.
...
mean(abs(dh.prm) >= abs(dh.obs))
[1] 0.090846

Do not mistake a permutation test using difference in medians as metric to be an alternative to a Wilcoxon rank sum test, which does show a significant difference
and seems an appropriate way to see if the two populations have the same
location:
wilcox.test(x1,x2)$p.val
[1] 0.03264336

