[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