4
$\begingroup$

I originally posted this on Stack Overflow (https://stackoverflow.com/questions/72649894/determining-direction-of-effect-in-wilcox-test-and-interpreting-pseudomedian), but I realized it fits better here. Apologies if reposting like this is a no-no.

I pretty frequently run Wilcoxon tests in R, but one question I've not been able to find a satisfactory answer to is how to (easily) determine the direction of effect, particularly in a paired test. In other words, if I run a paired test between time1 and time2, and the result is significant, how can I figure out whether there was a net increase or net decrease in the values between time1 and time2? Which one is "higher"?

I thought I had found the answer to this in the form of conf.int = TRUE in the call to wilcox.test, since this returns an "estimate" value that is either the "difference in location" for two sample/unpaired tests, or the "pseudomedian" for paired tests. I will admit I don't fully understand how this value is calculated, but I deduced that a positive value indicates the first sample was bigger (or that time1 - time2 was, on average, positive), and a negative value indicates the second sample was bigger (or that time1 - time2 was, on average, negative).

That is, until I encountered a weird edge case where, when running the paired test, the pseudomedian is positive no matter which way you input the vectors. This confuses me, because I would have thought these would necessarily have opposite signs. I apologize, my MRE is going to be messy because this doesn't seem to be a common issue, but kind of specific to these values:

time1 <- c(3.668e-05, 1.861e-05, 1.839e-05, 2.498e-05, 2.749e-05, 4.54e-05, 3.136e-05, 2.133e-05, 2.658e-05, 2.001e-05, 
                     1.729e-05, 2.648e-05, 2.957e-05, 4.523e-05, 4.783e-05, 0.0001078, 6.264e-05, 2.755e-05, 3.004e-05, 2.812e-05, 
                     3.212e-05, 3.77e-05, 4.832e-05, 6.537e-05, 3.566e-05, 2.286e-05, 2.366e-05, 3.65e-05, 4.119e-05, 2.058e-05, 
                     6.643e-06, 5.849e-05, 2.303e-05, 1.825e-05, 0.0001005, 2.885e-05, 1.82e-05, 1.59e-05, 3.474e-05, 2.058e-05, 
                     2.222e-05, 3.5e-05, 4.341e-06, 2.962e-05, 8.032e-05, 2.575e-05, 2.193e-05, 4.051e-05, 1.707e-05, 4.305e-05, 
                     1.722e-05, 4.717e-05, 1.832e-05, 3.919e-05, 6.009e-05, 2.579e-05, 2.629e-05, 5.719e-06, 1.083e-05, 7.241e-05, 
                     1.342e-05, 3.37e-05, 2.047e-05, 4.44e-05, 4.557e-05)
time2 <- c(3.352e-05, 2.525e-05, 1.942e-05, 1.488e-05, 3.992e-05, 1.855e-05, 3.64e-05, 1.442e-05, 2.41e-05, 6.004e-06, 
                     1.074e-05, 1.201e-05, 5.339e-06, 1.726e-05, 4.42e-05, 8.217e-05, 3.194e-05, 1.856e-05, 1.611e-05, 1.395e-05, 
                     2.696e-05, 3.224e-05, 1.652e-05, 7.501e-05, 1.251e-05, 1.349e-05, 1.277e-05, 2.187e-05, 2.674e-05, 2.048e-05, 
                     1.541e-05, 3.773e-05, 2.235e-05, 0, 4.773e-05, 1.969e-05, 1.176e-05, 1.166e-05, 2.499e-05, 1.127e-05, 2.188e-05, 
                     1.88e-05, 1.86e-05, 1.387e-05, 5.087e-05, 2.192e-05, 1.792e-05, 2.019e-05, 1.042e-05, 5.499e-06, 5.226e-06, 
                     3.641e-05, 1.705e-05, 1.334e-05, 3.281e-05, 3.213e-05, 1.66e-05, 0, 1.225e-05, 2.818e-05, 2.61e-05, 2.173e-05, 
                     1.486e-05, 7.171e-06, 4.991e-05)

> summary(time1)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
4.341e-06 2.058e-05 2.812e-05 3.358e-05 4.119e-05 1.078e-04 
> summary(time2)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 1.334e-05 1.860e-05 2.246e-05 2.696e-05 8.217e-05 
> summary(time1 - time2)
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-1.426e-05  3.160e-06  9.370e-06  1.113e-05  1.825e-05  5.277e-05 

So time1 has the larger median AND the difference is, on average, positive. As I would expect, then, running the paired test with time1 as the first sample gives a positive pseudomedian:

> wilcox.test(time1, time2, paired = TRUE, alternative = "two.sided", mu = 0, conf.int = TRUE)

    Wilcoxon signed rank test with continuity correction

data:  time1 and time2
V = 1892, p-value = 8.693e-08
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -1.4259e-05  5.2770e-05
sample estimates:
(pseudo)median 
     5.277e-05 

But then here's my problem: switching the order still gives a positive (albeit lower) pseudomedian:

> wilcox.test(time2, time1, paired = TRUE, alternative = "two.sided", mu = 0, conf.int = TRUE)

    Wilcoxon signed rank test with continuity correction

data:  time2 and time1
V = 253, p-value = 8.693e-08
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -5.2770e-05  1.4259e-05
sample estimates:
(pseudo)median 
    1.4259e-05

If you change to an unpaired test, the difference in location becomes negative as I would expect:

> wilcox.test(time2, time1, paired = FALSE, alternative = "two.sided", mu = 0, conf.int = TRUE)

    Wilcoxon rank sum test with continuity correction

data:  time2 and time1
W = 1242, p-value = 5.1e-05
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -3.349214e-05  3.521518e-06
sample estimates:
difference in location 
         -1.498531e-05 

It looks like, in the paired case, the pseudomedian is the maximum of the confidence interval (which seems to also just correspond to the min and max of the differences?). Shouldn't it be the middle/median/average of this interval? So my questions are:

  1. What exactly does the pseudomedian mean in the paired/one sample case? My understanding is, in the two sample case, it is the median of all possible pairwise differences, but I'm not sure how that translates here.
  2. Is there any way to determine the direction of effect from the output of wilcox.test? It seems a bit odd that a stats test can tell you that two samples are different, but not give an indication as to how they are different, since that is often just as important in my line of work.
  3. If the answer to #2 is "no", is it correct to take the sign of the median paired difference as an indication of the direction of the trend?
$\endgroup$
6
  • $\begingroup$ "Apologies if reposting like this is a no-no." The ideal behavior is to identify one good place for your post. If that's not where you originally posted it, you flag the original and ask to have it moved. $\endgroup$
    – Glen_b
    Jun 17, 2022 at 3:57
  • $\begingroup$ Oddly, when I multiply your data by 10^5, I get quite different answers (mutatis mutandis) for the pseudomedian (the p-value and test statistic are unchanged). This looks like a potential bug (perhaps!) in the calculation of the pseudomedian. While the estimate is not exactly flipped by reversing the two sets of times after that scaling, it is pretty close - at least accurate to the 4 significant figures in your post). I wonder if there's some hardcoded number in the code that should be relative rather than absolute, maybe. $\endgroup$
    – Glen_b
    Jun 17, 2022 at 4:23
  • $\begingroup$ The relative difference in the two pseudomedians (ignoring sign) gets smaller still when multiplying by $10^6$. Try this: L=10^6;wilcox.test(time1*L,time2*L,paired=TRUE,conf.int=TRUE); wilcox.test(time2*L,time1*L,paired=TRUE,conf.int=TRUE) $\endgroup$
    – Glen_b
    Jun 17, 2022 at 4:33
  • $\begingroup$ That is indeed odd and fascinating - so it seems like the data being such small numbers is the "cause" of this? I tried your suggestion, and then upped L to 10^12, and the pseudomedians / confidence intervals became perfectly symmetric, and the p-value is unchanged regardless of the scaling (good). I think I'm still confused if this is a rounding issue or other bug, or if this is the intended and "correct" behavior for such small values. $\endgroup$
    – tlbello
    Jun 17, 2022 at 19:28
  • 1
    $\begingroup$ Oh my goodness, I feel silly. I found the cause/solution: the tol.root arugment in wilcox.test defaults to 1e-4. If I set it to 1e-8, the issue goes away entirely and all the outputs are as expected. $\endgroup$
    – tlbello
    Jun 17, 2022 at 19:48

1 Answer 1

2
$\begingroup$

Answering my own question for any future searchers: the issue was a rounding/tolerance issue caused by the tol.root option in wilcox.test defaulting to 1e-4, which was larger than my data. Setting tol.root to a lower value fixed the issue:

> wilcox.test(time1, time2, paired = TRUE, alternative = "two.sided", mu = 0, conf.int = TRUE, tol.root = 1e-10)

    Wilcoxon signed rank test with continuity correction

data:  time1 and time2
V = 1892, p-value = 8.693e-08
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 7.135051e-06 1.396800e-05
sample estimates:
(pseudo)median 
  1.010505e-05 

> wilcox.test(time2, time1, paired = TRUE, alternative = "two.sided", mu = 0, conf.int = TRUE, tol.root = 1e-10)

    Wilcoxon signed rank test with continuity correction

data:  time2 and time1
V = 253, p-value = 8.693e-08
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -1.396800e-05 -7.135051e-06
sample estimates:
(pseudo)median 
 -1.010508e-05 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.