16
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ There has been a number of similar questions on this site already. Please look for median test. My own answer/comments is here. $\endgroup$
    – ttnphns
    Jan 10, 2014 at 22:55
  • $\begingroup$ As to comparing the medians of all three simultaneously, see my edit for the slightly modified R code. $\endgroup$ Jan 11, 2014 at 4:09

4 Answers 4

6
$\begingroup$

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

$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$ Apr 26, 2018 at 19:05
12
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – Alexander
    Jan 10, 2014 at 20:54
  • 1
    $\begingroup$ @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. $\endgroup$
    – Greg Snow
    Jan 10, 2014 at 23:00
  • $\begingroup$ +1. If you are speaking about Wilcoxon sum-rank test wouldn't you mind to convert "Wilcox" to that name? $\endgroup$
    – ttnphns
    Jan 10, 2014 at 23:08
  • $\begingroup$ @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? $\endgroup$
    – Glen_b
    Jan 11, 2014 at 2:13
10
$\begingroup$

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
}
$\endgroup$
4
  • 1
    $\begingroup$ +1 for naming the test. I didn't know that median test is also called Mood's test. $\endgroup$
    – ttnphns
    Jan 10, 2014 at 23:00
  • $\begingroup$ +1 Thanks for helping me with this, I really appreciate it! $\endgroup$
    – Alexander
    Jan 12, 2014 at 0:03
  • $\begingroup$ 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. $\endgroup$ Apr 26, 2018 at 18:56
  • $\begingroup$ Does this assume continuous (not discrete) distributions? I ran a two-group version with x=c(2,rep(1:2,50)) and y=rep(2,101) and got a p-value of 5e-19 (basically zero) even though the sample medians are identical (both equal 2), so presumably the p-value should be 1? But I think this failure relies critically on having discrete (not continuous) data. $\endgroup$ Sep 16, 2022 at 16:48
1
$\begingroup$

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)

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Jan 7, 2017 at 18:17
  • 1
    $\begingroup$ 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. $\endgroup$
    – JRF1111
    Jan 8, 2017 at 3:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.