Why does R say 'cannot compute exact p-values with ties' when I can do it with pen and paper?

Question

Suppose I have two sets of three numbers: $x_1, x_2, x_3$ and $y_1, y_2, y_3$ and I want to test the Null hypothesis that they are drawn from the same distribution using the Wilcoxon-Mann-Whitney test.

So what I would do with pen and paper is the following:

1. Compute the test statistic $W_{observed}$ for my sample
2. Make a list of all 20 ways of assigning the label '$x$' to three of my 6 numbers and '$y$' to the three others
3. Remember that under the Null hypothesis all 20 assignments in my list are equally likely
4. Compute the statistic $W$ for all of them.
5. Count how many of these $W$'s are at least as extreme (defined as 'far from the median of my list of 20 $W$-values') as the value $W_{observed}$ of step 1
6. Divide this number by 20
7. Interpret (with step 3 in mind) the outcome of step 6 as the probability under the Null-hypothesis (conditioned on observing these 6 numbers) that I would have observed a $W$-value at least as extreme as the one I did observe.
8. Remember that the probability described in step 7 is the definition of $p$-value
9. Conclude that, in step 6, I computed the $p$-value of my test
10. Feel pleased with myself because I did this computation 'exactly'

Alternatively, I can feed my six numbers in to R as wilcox.test(c(x_1, x_2, x_3), c(y_1, y_2, y_3))

Non-surprisingly, if all 6 numbers are distinct, both methods give the same answer, with R being slightly faster.

However, when some $x_i$ equals some $y_j$ and even when some $x_i$ equals some $x_j$ (this seems even more absurd to me, although I cannot quite pinpoint why) R will complain that it 'cannot compute exact $p$-values with ties.'

My question is why?

Obviously the problem cannot be in step 2 through 7 so it must be in step 1 or 8. Indeed I found an internet-source claiming that the test statistic does not exist in case of ties (supporting the thesis that the problem lies in step 1): because the statistic is defined as a sum of ranks, it becomes undefined when some of the number are equal, as well-defined ranks do no longer exist in that case.

However this argument is invalid because, and here I quote the R helpfile of wilcox test:

R's value can also be computed as the number of all pairs (x[i], y[j]) for which y[j] is not greater than x[i], the most common definition of the Mann-Whitney test.

Indeed, with this definition we can compute $W$ whether there are ties or not. So why does R still have trouble with computing $p$-values?

(Side note: there is a strange irony in the above quote form the help file, in that, in case of ties, R does not use above definition, but instead computes $W$ as the number of all pairs (x[i], y[j]) for which y[j] is strictly smaller than x[i] plus one half of all pairs (x[i], y[j]) for which y[j] is equal to x[i]. However that seems irrelevant to my question since as long as we have a way of computing $W$ we can run through the 10 steps I wrote down to obtain an exact $p$-value.)

Answer 1

The wilcox.test from the standard stats library is limited to cases without ties because it uses an algorithm from the function pwilcox that assumes that there are no ties.

This algorithm is not your pen and paper solution which would become computation intensive for larger sample sizes.

The algorithm in pwilcox is not computing all possibilities, and instead it has a function that counts the number of possibilities for the tail. E.g. if the statistic $U$ which has the range $0$ to $n_1n_2$, had a value close to $0$, say $U = 4$ then we only compute the number of permutation that result in $U = 1$, $U = 2$, $U = 3$ and $U = 4$. The computations are done with iterations. This is faster than making all potential permutations and compute the statistic for each of those permutations.

The documentation of the wilcox.test function mentions that there are functions in other packages that can compute the exact test for the cases with ties.

'wilcox_test' in package coin for exact, asymptotic and Monte Carlo conditional p-values, including in the presence of ties.

The c-code of the algorithm behind the wilcox.test can be viewed on GitHub here: https://github.com/SurajGupta/r-source/blob/master/src/nmath/wilcox.c

An interesting related article might be (although I haven't read it myself because it is behind a paywall):

Bergmann, Reinhard, John Ludbrook, and Will PJM Spooren. "Different outcomes of the Wilcoxon—Mann—Whitney test from different statistics packages." The American Statistician 54.1 (2000): 72-77.