# 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. • Welcome to Cross Validated! Interesting question, but lacks detail that would facilitate/motivate answers. Can you explain your code please? - not everyone reads R. What test for homogeneity are you performing after splitting by the overall median? - exact or approximate? And show your results - how many tests are you performing for each sample size? - what values of$n_1=n_2\$ are you considering? Feb 20, 2020 at 14:08
• 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. But I see your answer below and now it makes sense, thank you! Feb 21, 2020 at 16:30
• You're welcome. I added some code & a plot of results. For your simulation work note: (1) 1000 isn't very many repetitions (if you're expecting around 5% rejection rate, the std error of your Type I error estimate will be about 0.7%), & (2) you should be counting p-values less than or equal to 5% (which gives you the problem of testing equality of floating-point representions of numbers). Feb 21, 2020 at 17:37
• I see, I now increased my repetitions and the maximum of sample sizes with setting my p-value less than or equal to 5% and indeed I get the same results with the test as your plot. Thank you again! Feb 24, 2020 at 8:18
• You're still welcome. I've edited the information you provided into the question - the idea is that other people with the same question, or potential answerers can easily get what they need without wading through a comment thread (which can eventually be deleted). If my answer solved your problem you can indicate that to others by "accepting" it - clicking on the tick next to it. (That won't prevent others answering as well, if they've something to add or criticize, & you can change your mind.) Feb 24, 2020 at 9:09

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.