Hypothesis test for difference in medians among more than two samples

Question

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.

Answer 1

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).

Answer 2

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.

Answer 3

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
}

Answer 4

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)