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:
- 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.
- 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. - 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?
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$