[SciPy-user] help with scipy.stats.mannwhitneyu

josef.pktd@gmai... josef.pktd@gmai...
Thu Feb 5 11:23:16 CST 2009


On Thu, Feb 5, 2009 at 11:56 AM, Sturla Molden <sturla@molden.no> wrote:
> On 2/5/2009 5:39 PM, josef.pktd@gmail.com wrote:
>
>> According to R:
>> wilcox.test(x,y)
>> Performs one and two sample Wilcoxon tests on vectors of data; the
>> latter is also known as 'Mann-Whitney' test.
>>
>> I tried a normal random variable example ( no ties): the test
>> statistic returned is exactly the same as the one returned by
>> stats.mannwhitneyu(x,y) however the p-values differ. the pvalue in
>> stats is half of the one in R (up to 1e-17) as stated in the
>> docstring:  one-tailed p-value.
>
>
> I believe there is a bug in SciPy:
>
>
> def mannwhitneyu(x, y):
>     """Calculates a Mann-Whitney U statistic on the provided scores and
>     returns the result.  Use only when the n in each condition is < 20 and
>     you have 2 independent samples of ranks.  REMEMBER: Mann-Whitney U is
>     significant if the u-obtained is LESS THAN or equal to the critical
>     value of U.
>
>     Returns: u-statistic, one-tailed p-value (i.e., p(z(U)))
>     """
>     x = asarray(x)
>     y = asarray(y)
>     n1 = len(x)
>     n2 = len(y)
>     ranked = rankdata(np.concatenate((x,y)))
>     rankx = ranked[0:n1]       # get the x-ranks
>     #ranky = ranked[n1:]        # the rest are y-ranks
>     u1 = n1*n2 + (n1*(n1+1))/2.0 - np.sum(rankx,axis=0)  # calc U for x
>     u2 = n1*n2 - u1                            # remainder is U for y
>     bigu = max(u1,u2)
>     smallu = min(u1,u2)
>     T = np.sqrt(tiecorrect(ranked))  # correction factor for tied scores
>     if T == 0:
>         raise ValueError, 'All numbers are identical in amannwhitneyu'
>     sd = np.sqrt(T*n1*n2*(n1+n2+1)/12.0)
>     z = abs((bigu-n1*n2/2.0) / sd)  # normal approximation for prob calc
>     return smallu, 1.0 - zprob(z)
>
>
> Take a look at the last two lines? Do you see something peculiar?
>
> Sturla Molden
>

you mean that it uses bigu for the p-value calculation but reports
smallu as the test-statistic?

I didn't try to figure out what the formula for the p-value actually
is, but I'm pretty happy that we get the same result as R, except for
the times 2.

I looked some more at the R implementation :

the main difference is that R uses by default a continuity correction
"correct a logical indicating whether to apply continuity correction
in the normal approximation for the p-value"

>>> rwilcox=rpy.r('wilcox.test')
>>> stats.mannwhitneyu(rvs1,rvs2)[1]*2 - rwilcox(rvs1,rvs2,correct = False)['p.value']
-1.5265566588595902e-016

The test statistic in R is not symmetric in its argument, although the
p-values are, stats.mannwhitneyu is symmetric in statistic and
p-value.

>>> rresult = rwilcox(rvs2,rvs1)
>>> rresult['statistic']
{'W': 5637.0}
>>> rresult['p.value']
0.11989439052971607

>>> rresult = rwilcox(rvs1,rvs2)
>>> rresult['statistic']
{'W': 4363.0}
>>> rresult['p.value']
0.11989439052971618

So overall stats.mannwhitney, I think, looks pretty good but it could
be expanded to include some of the options that R offers, and I also
think we should multiply the pvalue by 2, so that the reported p-value
actually corresponds to the test.

Josef


More information about the SciPy-user mailing list