# Hypothesis test for difference in medians among more than two samples

### Question

The test scores of three groups of people are saved as separate vectors in R.

set.seed(1)
group1 <- rnorm(100, mean = 75, sd = 10)
group2 <- rnorm(100, mean = 85, sd = 10)
group3 <- rnorm(100, mean = 95, sd = 10)


I want to know if there is a significant difference in the medians between these groups. I know that I could test group 1 versus group 2 using the Wilcoxon test, like so.

wilcox.test(group1, group2)


However, this compares only two groups at a time, and I would like to compare all three simultaneously. I would like a statistical test that yields a p value at the 0.05 significance level. Could someone please help?

### Edit #1 - Mood's median test

Following user Hibernating's suggested answer, I tried Mood's median test.

median.test <- function(x, y){
z <- c(x, y)
g <- rep(1:2, c(length(x), length(y)))
m <- median(z)
fisher.test(z < m, g)$p.value } median.test(group1, group2)  However, this approach allows me to test for a significant difference between the medians of only two groups at a time. I am not sure of how to use it to compare the medians of all three simultaneously. ### Edit #2 - Kruskal-Wallis test User dmartin's suggested answer appears to be more or less what I need, and allows me to test all three groups simultaneously. kruskal.test(list(group1, group2, group3))  ### Edit #3 User Greg Snow helpfully notes in his answer that the Kruskal-Wallis test is appropriate as long as it makes strict assumptions that make it also a test of means. • There has been a number of similar questions on this site already. Please look for median test. My own answer/comments is here. – ttnphns Jan 10 '14 at 22:55 • As to comparing the medians of all three simultaneously, see my edit for the slightly modified R code. – Hibernating Jan 11 '14 at 4:09 ## 4 Answers The Kruskal-Wallis test could also be used, as it's a non-parametric ANOVA. Additionally, it is often considered to be more powerful than Mood's median test. It can be implemented in R using the kruskal.test function in the stats package in R. To respond to your edit, interpreting K-W is similar to a one-way ANOVA. A significant p-value corresponds to rejected the null that all three means are equal. You must use a follow-up test (again, just like an ANOVA), to answer questions about specific groups. This typically follows specific research questions you may have. Just by looking at the parameters of the simulation, all three groups should be significantly different from one another if you do a follow-up test (as they're all 1 SD apart with N = 100). • To clarify a couple of things. 1) Kruskal-Wallis is not a test of medians, unless the distributions of observations in the groups meet certain assumptions. If you are really looking to compare medians it may not be the appropriate test. It is best to choose a test that actually tests the hypothesis you are interested in testing. 2) Kruskal-Wallis is not an "ANOVA". That is, it is not an analysis of variance. 3) The mention of "means" in this answer is incorrect. – Sal Mangiafico Apr 26 '18 at 19:05 First, the Wilcoxon test (or Mann-Whitney test) is not a test of medians (unless you make very strict assumptions that also make it a test of means). And for comparing more than 2 groups the Wilcoxon test can lead to some paradoxical results (see Efron's Dice). Since the Wilcoxon test is just a special case of a permutation test and you are specifically interested in the medians, I would suggest a permutation test on the medians. First choose a measure of the difference, something like the largest of the 3 medians minus the smallest of the 3 (or the variance of the 3 medians, or the MAD, etc.). Now compute your stat for the original data. pool all the data in one set then randomly partition the values into 3 groups of the same sizes as the original and compute the same statistic. repeat many times (like 9998) Compare how the statistic from the real data compares to the distribution of all the stats for your test. • Let's say that I am willing to make the strict assumptions necessary for the Wilcox test that would also make it a test of means. Would that require changing the R code I have written above? Could this also be done for the Kruskal-Wallis test? – Alexander Jan 10 '14 at 20:54 • @Alexander, If you are willing to make those assumptions then the R code is fine and Kruskal Wallis would also be fine. But then if you are willing to make those assumptions then t.test and aov would probably be fine as well. – Greg Snow Jan 10 '14 at 23:00 • +1. If you are speaking about Wilcoxon sum-rank test wouldn't you mind to convert "Wilcox" to that name? – ttnphns Jan 10 '14 at 23:08 • @GregSnow +1 for the points made ... but I assume by 'Wilcox' you mean the test named after Frank Wilcoxon. (This confusion is unfortunately compounded by R, which - misleadingly - calls the corresponding test wilcox.test). Could you edit? – Glen_b -Reinstate Monica Jan 11 '14 at 2:13 The Mood’s median test is a nonparametric test that is used to test the equality of medians from two or more populations. See here for the R part of your question. See also a related question here. Also from here: Mood's median test is the easiest one to do by hand: work out the overall median (of all the data), and count how many values are above and below the median in each group. If the groups are all about the same, the observations should be about 50-50 above and below the overall median in each group... The counts of below-median and above-median... form a two-way table, which is then analyzed using a chi-squared test. Mood's median test is a lot like the sign test generalized to two or more groups. Edit: For three groups, you may consider this simple generalisation of the R code I linked to: median.test2 <- function(x, y, z) { a <- c(x, y, z) g <- rep(1:3, c(length(x), length(y), length(z))) m <- median(a) fisher.test(a < m, g)$p.value
}

• +1 for naming the test. I didn't know that median test is also called Mood's test. – ttnphns Jan 10 '14 at 23:00
• +1 Thanks for helping me with this, I really appreciate it! – Alexander Jan 12 '14 at 0:03
• I know of a couple of implementations in R. mood.medtest in the RVAideMemoire package appears to be the usual test except that it uses Fisher exact test by default for smaller sample sizes. The median_test function in the coin package can provide an asymptotic test or use Monte Carlo. – Sal Mangiafico Apr 26 '18 at 18:56

I know this is way late, but I couldn't find a good package for Mood's median test either, so I took it upon myself to make a function in R that seems to do the trick.

#Mood's median test for a data frame with one column containing data (d),
#and another containing a factor/grouping variable (f)

moods.median = function(d,f) {

#make a new matrix data frame
m = cbind(f,d)
colnames(m) = c("group", "value")

#get the names of the factors/groups
facs = unique(f)

#count the number of factors/groups
factorN = length(unique(f))

#Make a 2 by K table that will be saved to the global environment by using "<<-":
#2 rows (number of values > overall median & number of values <= overall median)
#K-many columns for each level of the factor
MoodsMedianTable <<- matrix(NA, nrow = 2, ncol = factorN)

rownames(MoodsMedianTable) <<- c("> overall median", "<= overall median")
colnames(MoodsMedianTable) <<- c(facs[1:factorN])
colnames(MoodsMedianTable) <<- paste("Factor: ",colnames(MoodsMedianTable))

#get the overall median
overallmedian = median(d)

#put the following into the 2 by K table:
for(j in 1:factorN){ #for each factor level

g = facs[j] #assign a temporary "group name"

#count the number of observations in the factor that are greater than
#the overall median and save it to the table
MoodsMedianTable[1,j] <<- sum(m[,2][ which(m[,1]==g)] > overallmedian)

#count the number of observations in the factor that are less than
# or equal to the overall median and save it to the table
MoodsMedianTable[2,j] <<- sum(m[,2][ which(m[,1]==g)] <= overallmedian)

}

#percent of cells with expected values less than 5
percLT5 = ((sum(chisq.test(MoodsMedianTable)$expected < 5)) / (length(chisq.test(MoodsMedianTable)$expected)))

#if >20% of cells have expected values less than 5
#then give chi-squared stat, df, and Fisher's exact p.value
if (percLT5 > 0.2) {
return(list(
"Chi-squared" = chisq.test(MoodsMedianTable)$statistic, "df" = chisq.test(MoodsMedianTable)$parameter,
"Fisher's exact p.value" = fisher.test(MoodsMedianTable)$p.value)) } #if <= 20% of cells have expected values less than 5 #then give chi-squared stat, df, and chi-squared p.value if (percLT5 <= 0.2) { return(list( "Chi-squared" = chisq.test(MoodsMedianTable)$statistic,
"df" = chisq.test(MoodsMedianTable)$parameter, "Chi-squared p.value" = chisq.test(MoodsMedianTable)$p.value))

}

}


For the OP's question, you would first run this to make a new data frame to hold the values from your three group vectors with a matched "group" variable.

require(reshape2)
df = cbind(group1, group2, group3)
df = melt(df)
colnames(df) = c("observation", "group", "value")


and run the function for Mood's median test with moods.median(df$value, df$group)

• It seems that the Kruskal-Wallis test was the answer. The OP needed a solution with 3 groups. It appears that ttnphns has already provided R code for the Mood test. – Michael R. Chernick Jan 7 '17 at 18:17
• The code ttnphns gave only provides a p value, the one I wrote also gives the chi squared stat and df, and it works for any number of groups. I mostly just posted here since this post is the first to come up when searching for how to do Mood's median median test in R. – JRF1111 Jan 8 '17 at 3:41