When R's Mann-Whitney U test returns extremely distorted p-values

The Mann–Whitney U test (also known as the Wilcoxon rank-sum test) is one of the most popular nonparametric statistical tests. In R, it can be accessed using the wilcox.test function, which has been available since R 1.0.0 (February 2000). With its extensive adoption and long-standing presence in R, the wilcox.test function has become a trusted tool for many researchers. But is it truly reliable, and to what extent can we rely on its accuracy by default?

In my work, I often encounter the task of comparing a large sample (e.g., of size 50+) with a small sample (e.g., of size 5). In some cases, the ranges of these samples do not overlap with each other, which is the extreme case of the Mann–Whitney U test: it gives the minimum possible p-value. In one of the previous posts, I presented the exact equation for such a p-value. If we compare two samples of sizes \(n\) and \(m\), the minimum p-value we can observe with the one-tailed Mann–Whitney U test is \(1/C_{n+m}^n\). For example, if \(n=50\) and \(m=5\), we get \(1/C_{55}^5 \approx 0.0000002874587\). Let’s check these calculations using R:

> wilcox.test(101:105, 1:50, alternative = "greater")$p.value
[1] 0.0001337028

The obtained p-value is \(\approx 0.0001337028\), which is \(\approx 465\) times larger than we expected! Have we discovered a critical bug in wilcox.test? Can we now trust this function? Let’s find out!

In fact, the wilcox.test has two implementations of the Mann–Whitney U test: the exact one and the approximated one. The choice of a particular strategy is controlled by the exact parameter of this function. Here is an example from the official documentation:

> x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
> y <- c(1.15, 0.88, 0.90, 0.74, 1.21)
> wilcox.test(x, y, alternative = "greater", exact = TRUE)$p.value
[1] 0.1272061
> wilcox.test(x, y, alternative = "greater", exact = FALSE)$p.value
[1] 0.1223118

In this example, the approximation looks accurate, with a difference ofhis makes it an impra only about \(5\%\). When the exact parameter is not specified, it is determined as follows:

if(is.null(exact))
    exact <- (n.x < 50) && (n.y < 50)

Before we start exploring different extreme cases, let us briefly review both methods to estimate the p-value within this test.

The first step of the Mann–Whitney U test is calculating the U statistic:

\[U(x, y) = \sum_{i=1}^n \sum_{j=1}^m S(x_i, y_j),\quad S(a,b) = \begin{cases} 1, & \text{if } a > b, \\ 0.5, & \text{if } a = b, \\ 0, & \text{if } a < b. \end{cases} \]

It is a simple rank-based statistic, which varies from \(0\) to \(nm\). Once we get the value of \(U\), we should match it with the distribution of \(U\) statistics under the valid null hypothesis. The distribution can be obtained in two different ways:

• The exact implementation
  The straightforward implementation uses dynamic programming to obtain the exact probability of obtaining the given \(U\) statistic for the given \(n\) and \(m\). The only drawback of this solution is its computational complexity (which is \(\mathcal{O}(n^2 m^2)\) for both time and memory). This makes this approach an impractical solution for large samples.
• The approximated implementation
  Fortunately, we have a method of large sample approximation (see [Hollander1973], pages 68–75). Long story short, we approximate the distribution of \(U\) statistics with the following normal distribution (in the case of tied values, a small correction is needed):

\[\mathcal{N}\left( \frac{nm}{2},\, \frac{nm(n+m+1)}{12} \right) \]

For the middle values of \(U\) (around \(nm/2\)), this approximation provides reasonably accurate results. However, it performs extremely poorly at the tails. To illustrate the problems, we perform a simple numerical study:

• Enumerate \(n\) from \(1\) to \(50\), set \(m\) to \(50\) (since \(50\) is a critical value, after which wilcox.text switches to the approximated algorithm);
• For each pair of \(n,\, m\), we consider two non-overlapped samples \(x = \{ 101, 102, \ldots, 100 + n \}\) and \(y = \{ 1, 2, \ldots, m \}\);
• For each pair of samples \(x\) and \(y\), we estimate the exact and approximated p-values using wilcox.test (R 4.3.0). We also evaluate the ratio between two obtained p-values.

This simulation can be performed using the following simple script:

options(scipen = 999)
df <- expand.grid(n = 1:50, m = 50)
df$approx <- sapply(1:nrow(df), function(i) wilcox.test(101:(100+df$n[i]), 1:df$m[i], "g", exact = F)$p.value)
df$exact <-  sapply(1:nrow(df), function(i) wilcox.test(101:(100+df$n[i]), 1:df$m[i], "g", exact = T)$p.value)
df$ratio <- df$approx / df$exact
df

Here are the results:

Thus, for \(n=m=50\), the ratio between the exact (the correct one) and approximated (the default one) solutions is about \(3.5 \cdot 10^{11}\). Of course, both p-values are pretty small in this case. If you are using the notorious \(\alpha = 0.05\), you will most likely avoid such a problem (however, for \(n=1,\, m = 50\), the exact true p-value is \(0.0196\), while the approximated one is \(0.0480\), whish is alarmingly close to \(0.05\)).

If you prefer using extremely low statistical significance levels and there are non-zero chances of getting extreme values of the \(U\) statistic, you may consider making exact = TRUE your default way to call wilcox.test. Such a simple trick may save you from misleading research results.

References

• [Hollander1973]
  Myles Hollander and Douglas A. Wolfe (1973) Nonparametric Statistical Methods. New York: John Wiley & Sons. DOI:10.1002/9781119196037