Example 2014.6: Comparing medians and the Wilcoxon rank-sum test

A colleague recently contacted us with the following question: "My outcome is skewed-- how can I compare medians across multiple categories?"

What they were asking for was a generalization of the Wilcoxon rank-sum test (also known as the Mann-Whitney-Wilcoxon test, among other monikers) to more than two groups. For the record, the answer is that the Kruskal-Wallis test is the generalization of the Wilcoxon, the one-way ANOVA to the Wilcoxon's t-test.

But this question is based on a false premise: that the the Wilcoxon rank-sum test is used to compare medians. The premise is based on a misunderstanding of the null hypothesis of the test. The actual null hypothesis is that there is a 50% probability that a random value from one population exceeds an random value from the other population. The practical value of this is hard to see, and thus in many places, including textbooks, the null hypothesis is presented as "the two populations have equal medians". The actual null hypothesis can be expressed as the latter median hypothesis, but only under the additional assumption that the shapes of the distributions are identical in each group.

In other words, our interpretation of the test as comparing the medians of the distributions requires the location-shift-only alternative to be the case. Since this is rarely true, and never assessed, we should probably use extreme caution in using, and especially in interpreting, the Wilcoxon rank-sum test.

To demonstrate this issue, we present a simple simulation, showing a case of two samples with equal medians but very different shapes. In one group, the values are exponential with mean = 1 and therefore median log 2, in the other they are normal with mean and median = log 2 and variance = 1. We generate 10,000 observations and show that the Wilcoxon rank-sum test rejects the null hypothesis. If we interpret the null incorrectly as applying to the medians, we will be misled. If our interest actually centered on the medians for some reasons, an appropriate test that would not be sensitive to the shape of the distribution could be found in a quantile regression. Another, of course, would be the median test. We show that this test does not reject the null hypothesis of equal medians, even with the large sample size.

We leave aside the deeper questions of whether comparing medians is a useful substitute for comparing means, or whether means should not be compared when distributions are skewed.

R
In R, we simulate two separate vectors of data, then feed them directly to the wilcox.test() function (section 2.4.2).

y1 = rexp(10000)
y2 = rnorm(10000) + log(2)

wilcox.test(y1,y2)

This shows a very small p-value, denoting the fact not that the medians are unequal but that one or the other of these distributions generally has larger values. Hint: it might be the one with the long tail and no negative values.

 Wilcoxon rank sum test with continuity correction

data:  y1 and y2
W = 55318328, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0

To set up the quantile regression, we put the observations into a single vector and make a group indicator vector. Then we load the quantreg package and use the rq() function (section 4.4.4) to fit the median regression.

y = c(y1,y2)
c = rep(1:2, each = 10000)

library(quantreg)
summary(rq(y~as.factor(c)))

This shows we would fail to reject the null of equal medians.

Call: rq(formula = y ~ as.factor(c))

tau: [1] 0.5

Coefficients:
              Value    Std. Error t value  Pr(>|t|)
(Intercept)    0.68840  0.01067   64.53047  0.00000
as.factor(c)2 -0.00167  0.01692   -0.09844  0.92159
Warning message:
In rq.fit.br(c, y, tau = tau, ...) : Solution may be nonunique

SAS
In SAS, we make a data set with a group indicator and use it to generate data conditionally.

data wtest;
do i = 1 to 20000;
  c = (i gt 10000);
  if c eq 0 then y = ranexp(0);
    else y = normal(0) + log(2);
  output;
  end;
run;

The Wilcoxon rank-sum test is in proc npar1way (section 2.4.2).

proc npar1way data = wtest wilcoxon;
class c;
var y;
run;

As with R, we find a very small p-value.

                       The NPAR1WAY Procedure

             Wilcoxon Scores (Rank Sums) for Variable y
                      Classified by Variable c

                     Sum of     Expected      Std Dev         Mean
  c          N       Scores     Under H0     Under H0        Score

  0      10000    106061416    100005000   408258.497   10606.1416
  1      10000     93948584    100005000   408258.497    9394.8584


                      Wilcoxon Two-Sample Test

                Statistic             106061416.0000

                Normal Approximation
                Z                            14.8348
                One-Sided Pr >  Z             <.0001
                Two-Sided Pr > |Z|            <.0001

                t Approximation
                One-Sided Pr >  Z             <.0001
                Two-Sided Pr > |Z|            <.0001

             Z includes a continuity correction of 0.5.


                        Kruskal-Wallis Test

                Chi-Square                  220.0700
                DF                                 1
                Pr > Chi-Square               <.0001

The median regression is in proc quantreg (section 4.4.4). As in R, we fail to reject the null of equal medians.

proc quantreg data =wtest;
class c;
model y = c;
run;

                       The QUANTREG Procedure

                  Quantile and Objective Function

              Predicted Value at Mean          0.6909


                        Parameter Estimates

                            Standard    95% Confidence
    Parameter   DF Estimate    Error        Limits       t Value

    Intercept    1   0.6909   0.0125    0.6663    0.7155   55.06
    c         0  1   0.0165   0.0160   -0.0148    0.0479    1.03
    c         1  0   0.0000   0.0000    0.0000    0.0000     .

                        Parameter Estimates

                        Parameter   Pr > |t|

                        Intercept     <.0001
                        c         0   0.3016
                        c         1    .

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.