Inconsistent result in checking the approximate exactness of Mood’s median test I am running a simulation study to check the approximate exactness of Mood’s median test with the following R code (using the RVAideMemoire package):
set.seed(1)
typeIErrors <- 0
power <- 0
for (i in 1:reps){
group1 <- rnorm(n1,sd=2)
group2 <- rnorm(n2,sd=2)
bothGroups <- c(group1,group2)
groups <- c(rep(1,length(group1)),rep(2,length(group2)))
tmp <- mood.medtest(bothGroups,groups)
if (tmp$p.value<0.05){
  typeIErrors <- typeIErrors+1
  }
}
print(typeIErrors/reps)

The proportion of type 1 errors seems to go up and down when I increase the samples sizes stepwise (keeping n1 and n2 equal) instead of converging to 0.05, how is this possible?

I am doing Mood’s tests 1000 times where I keep n1 and n2 equal, and both sampled from a normal distribution with a sd of 2. Then out of those thousand tests, I get the number of times the p-value of the test was below .05 (so the type 1 errors) and then I divide this number by 1000 to get the type 1 error proportion. I start with n1=n2=12, then 25, 50 and 75 and 100, which gives me type 1 error rates of 0.037, 0.024, 0.023, 0.052 and 0.043 respectively. Until sample sizes of 200 the default is a Fisher exact test.
 A: You don't need to simulate. Let $n$ be the no. observations in each group; by definition $n$ observations overall exceed the median of the pooled observations, & $n$ fall short of it. The data therefore constitute a contingency table with fixed row margins $(n, n)$ & fixed column margins $(n, n)$; the count in any one cell determines the counts in the other three cells, & follows a hypergeometric distribution under the null hypothesis that all observations are identically & independently distributed, regardless of the group they belong to. So, the probability pertaining to the rejection region of any test performed on that contingency table can be calculated with the hypergeometric mass function.
If you perform an exact test (one based on the hypergeometric distribution), the size of that test, the probability of a Type I error, is less than or equal to the significance level by construction; when the size fails to attain the level it's owing to the discreteness of the sample space. (If you're using an approximate test the true size may exceed the nominal level.) As $n$ increases, the sample space becomes more fine-grained, & the general trend is for the size to approach the level. But the approach needn't be smooth: at each increment of $n$ both the elements of the sample space within the rejection region & the probability of those elements can change. A brief illustration:—When $n=3$, the rejection region of a one-tailed test at the $5\%$ level is $\{0\}$, & the size precisely $5\%$. When $n=4$ & $n=5$, the rejection region is still $\{0\}$, but the size has dropped to $\sim 1.43\%$ & $\sim 0.40\%$ respectively. When $n=6$, the rejection region is now $\{0,1\}$, & the size has gone up to $\sim 4.00\%$. 
Here is R code to calculate Type I errors for sample sizes up to 1000:
n <- 1:1000 # sample sizes
level <- 0.05 # significance level
eps <- 1e-7 # tolerance for equality comparison
size <- numeric(length(n)) # initialized vector of sizes
for (i in 1:length(n)){
  x <- 0:n[i]
  tsp <- phyper(x, n[i], n[i], n[i]) # hypergeometric CDF
  size[i] <- max(tsp[tsp <= (level + eps)])
}
size[size == -Inf] <- 0
plot(n, size, type="l", xlab="n", ylab="Type I error")
abline(h=level, lty="dashed")


In fine, you should expect Type I error to be a non-monotone function of sample size; but though it continues to bump around as the sample size increases, the bumps get smaller & it approaches the nominal significance level.
