Difference of 'centers' of 2 non-normal samples with Mann-Whitney test

I have 2 non-normally distributed samples of different sizes (N1~=N2).

To evaluate whether there is a significant difference between these samples, I used the Mann Whitney U test (ranksum in MATLAB).

Now I want to evaluate by how much the populations differ. With normally distributed data, I would just use the difference of the means with a SEM confidence interval.

However, for nonparamteric data, my Wikipedia research suggested using the Hodges–Lehmann estimator. My questions are:

1. Is the Hodges–Lehmann estimator indeed correct?
2. How do I calculate confidence intervals for the HL estimator (like the SEM)
3. How can I implement this in MATLAB?

Thanks

• It would be surprising to me if it wasn't already implemented. It's in R, for example - the Mann-Whitney is done with wilcox.test and you just add the argument conf.int=TRUE and you get the HL estimate and an interval (which you can change the coverage for via another argument). Commented Nov 20, 2013 at 1:07
• I am just as shocked as you are, but its not. Heres the MATLAB documentation: mathworks.com/help/stats/ranksum.html. Do you know how to implement it? Commented Nov 20, 2013 at 1:11
• I've added information on computing the intervals Commented Nov 21, 2013 at 5:05

Yes, if your alternative is a shift alternative, the corresponding estimate is the Hodges-Lehmann.

The estimate is easy enough.

The Hodges Lehmann estimate is the sample median of the cross-sample pairwise differences.

For really large samples, there are efficient methods (for example, ones that sort the two samples and then use the information in that to speed up the search for the median).

I'm going to assume your sample is not so large that its worth the effort of doing that, in which case, you simply want to find all pairwise differences and find the median of the result.

For example, if your x-sample is:

16.388  6.775 18.270 17.034 18.825 16.197 16.709 18.188  9.999


 9.305  9.799 11.103


then the matrix of differences, x-y is

   7.083  6.589  5.285
-2.530 -3.024 -4.328
8.965  8.471  7.167
7.729  7.235  5.931
9.520  9.026  7.722
6.892  6.398  5.094
7.404  6.910  5.606
8.883  8.389  7.085
0.694  0.200 -1.104


(here columns correspond to particular y's and rows to particular y's). You then take the median of the whole set of differences.

(In R this is just a matter of:

> median(outer(x,y,"-"))
[1] 6.91


It's doubtless just about as easy in Matlab; would something like [X,Y] = meshgrid(x,y); Z = X - Y; work to generate the matrix of differences?)

If you have samples so large that this isn't feasible, you should indicate this.

I don't know why it didn't click before, but the basic concept of finding an interval for something like the Hodges-Lehmann shift parameter in the Wilcoxon-Mann-Whitney (WMW) test is pretty obvious; the same idea is widely used for other permutation tests.

The underlying idea is this: you shift the second sample left and right until you hit the $$\alpha/2$$ and $$1-\alpha/2$$ quantiles of the WMW statistic (being discrete, these will cut off less than $$\alpha/2$$ in the tails, but you can legitimately have one a bit larger than $$\alpha/2$$ if the sum of tail areas is less than $$\alpha$$; however, most people will prefer symmetry even at the cost of a more-conservative-than-necessary interval).

This approach can be easily implemented in software via use of a canned root-finding algorithm that doesn't require derivatives. Bounds are obvious enough (if you haven't hit it by the time you get complete separation of samples, you never will, so check those 'maximum shift' values first and then use them as bounds; similarly, the HL estimate itself is the bound for the other end of each search).

There are other approaches to doing this that will probably be more efficient, but this should work well enough with a good root-finder (perhaps involving only a dozen or so evaluations of statistics at different shift values).

In response to the comments about the issue being a change of scale, not location:

One possibility - you could consider looking at estimating the location shift in the logs, which should correspond to the log of the scale change (ratio of scales) in the untransformed values. That is if $$X∼F(x)$$ and $$Y∼F(y/a)$$ then if $$X^∗=\ln X∼G(x), Y^∗=\ln Y∼G(y−\ln a)$$, where $$G(x)= F(\exp(x))$$.

Intervals, at least should carry over nicely.

Here's an example of that in use:

The data are as follows. The x sample is 60 observations from a distribution with scale parameter $$\mu=18$$, while the y sample is 50 observations with scale parameter $$\mu=25$$. They have the same distribution, aside from a shift in scale.

The increase in scale is a factor of $$25/18 \approx 1.3889$$.

Here's a boxplot of the data on the original scale and on the log-scale:

Here we do a Wilcoxon-Mann-Whitney on the log-data:

> lx=log(x)
> ly=log(y)
> wilcox.test(ly,lx,conf.int=TRUE)

Wilcoxon rank sum test with continuity correction

data:  ly and lx
W = 1743, p-value = 0.1455
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-0.09044406  0.68813695
sample estimates:
difference in location
0.2917924


If we exponentiate the interval, it should contain the true ratio of scales 95% of the time. It does contain the population ratio on this occasion:

> exp(wilcox.test(ly,lx,conf.int=TRUE)$conf.int) [1] 0.9135254 1.9900046 attr(,"conf.level") [1] 0.95  If we exponentiate the estimated shift in the logs, it will be a (somewhat) downward-biased estimate of the scale-ratio, but it's fairly good: > exp(wilcox.test(ly,lx,conf.int=TRUE)$estimate)
difference in location
1.338825


[The data used in the example were generated as follows:

> set.seed(94394381)
> x=rexp(60,1/18)
> y=rexp(50,1/25)


... just in case anyone wants to try it]

There are other ways to estimate the scale change; one could, for example, divide the y-sample by a constant, $$c$$, and try different values of $$c$$ until the WMW statistic was equal to its expected value under the null hypothesis. One could also try different values until the statistic was equal to the $$\alpha/2$$ and $$1-\alpha/2$$ quantiles of the null distribution of the WMW. As above, this will be the point estimate of and interval for the ratio of scales.

(That would have the advantage of not introducing a small bias into the central estimate.)

• Hi Glen. Thank you for the answer! My data is not paired and the sample sizes are different. How would I get the estimate then? Also, you mentioned "shift alternative". I guess I am actually looking for the scale estimate, ie r=A/B not s=A-B. Would that be a significantly different proceedure? (ie is r=1+s/median(B) ) Thanks! Commented Nov 20, 2013 at 21:18
• I'll expand on what's involved then. I may have linked to the wrong calculation (sorry, it's trivial for me to write in R, but I don't remember Much matlab and I don't have Matlab handy to try it with). $\quad$ If you're after a scale-shift estimate, it's a bit trickier. Are your data all positive? Commented Nov 20, 2013 at 21:27
• One possibility - you could consider looking at estimating the location shift in the logs, which should correspond to the log of the scale shift in the untransformed values. That is if $X\sim F(x)$ and $Y\sim F(y/a)$ then if $X^*=\ln X\sim G(x)$, $Y^*=\ln Y\sim G(y-\ln a)$. Intervals at least should carry over nicely. Commented Nov 20, 2013 at 22:00
• Glen, I tried the method you advocate and I feel like its not accurate just by visual inspection. To confirm im doing it right, I (1) took the log of the data (2) estimated the 'shift' (3) scale=exp(shift). My background isnt statistics, but I wonder how innacurate it would be to just get the 'shift' & estimate the 'scale' as 'scale'='shift'/(average of the 2 medians), or some other ad-hoc method like that Commented Nov 21, 2013 at 8:01
• DankMasterDan - my apologies; I didn't make my meaning sufficiently clear. I have added an example that should help, and then discuss an alternative suggestion for estimating the ratio of scales. If you have changes in both location and scale this approach won't work. Commented Nov 21, 2013 at 10:24