Wilcoxon.test in R will not calculate exact distribution due to ties (scipy.stats.wilcoxon will)

Question

I'm running a pre-post hypothesis test on a small dataset, due to it's size n = 12 I am running the test for the exact distribution.

df = 
ID  Period  Varieble    Value
1   0 Month Something   18
1   3 Month Something   26
2   0 Month Something   23
2   3 Month Something   4
3   0 Month Something   24
3   3 Month Something   3
4   0 Month Something   27
4   3 Month Something   26
5   0 Month Something   9
5   3 Month Something   0
6   0 Month Something   40
6   3 Month Something   3
7   0 Month Something   17
7   3 Month Something   10
8   0 Month Something   33
8   3 Month Something   9
9   0 Month Something   6
9   3 Month Something   7
10  0 Month Something   8
10  3 Month Something   1
11  0 Month Something   9
11  3 Month Something   4
13  0 Month Something   26
13  3 Month Something   9

This works fine in python:

dfb = df.query('Period == "0 Month"')

dfa = df.query('Period == "3 Month"')

wilcoxon(dfb.Value, dfa.Value, alternative='two-sided', correction=True, 
          mode ='exact',)

Result = WilcoxonResult(statistic=7.5, pvalue=0.01220703125)

When I try to do the same in R I get a different result as it defaults to the Normal Approximation and provides the warning:

Warning message:
"In wilcox.test.default(dfb$Value, dfa$Value, paired = TRUE, 
         exact = TRUE,  :
  cannot compute exact p-value with ties"

R Code:

dfb <- filter(df, Period == "0 Month")
dfa <- filter(df, Period == "3 Month")

wilcox.test(dfb$Value,dfa$Value,paired=TRUE, exact = FALSE, 
       alternative = 'two.sided')

Result:V = 70.5, p-value = 0.01494

Is there a statistically sound reason why R is defaulting to the normal approximation due to the existing ties? I can't find this in any of the literature.

Conversely, is the method that the scipy library is using valid?

Answer 1

You should note the following paragraph in the scipy manual:

To derive the p-value, the exact distribution (mode == 'exact') can be used for sample sizes of up to 25. The default mode == 'auto' uses the exact distribution if there are at most 25 observations and no ties, otherwise a normal approximation is used (mode == 'approx').

(bolding is my own). The question remains why the results differ.

Answer 2

I don't have enough reputation points to add this as a comment, so although not an answer, hope it does help.

I've adjusted your code a bit for convenience:

import numpy as np
dfa = np.array([26, 4, 3, 26, 0, 3, 10, 9, 7, 1, 4, 9])
dfb = np.array([18, 23, 24, 27, 9, 40, 17, 33, 6, 8, 9, 26])
d = dfb - dfa

One thing to notice is that in the differences a few will have the same absolute value, and will therefore get the same rank.

If I run the following:

wilcoxon(d, alternative='two-sided', correction=False, mode ='exact')

The result of a statistic of 7.5 and a pvalue of 0.012207 is the same as when I run my own exact test but with a statistic of 8. I've made my own using the information from https://www.real-statistics.com/non-parametric-tests/wilcoxon-signed-ranks-test/wilcoxon-signed-ranks-exact-test/

It could be that scipy simply rounds the test statistic, but that is just a guess and it might be using something far more complex (perhaps something from https://www.jstor.org/stable/2284536 but I haven't read the article myself yet).

The result from R can be obtained using:

wilcoxon(d, alternative='two-sided', correction=True, mode ='approx')

So R uses the normal approximation, with a continuity correction (which is different from the correction for ties).

Hope this helps a bit.