I have performed Wilcoxon signed-rank tests for variables $X$, $Y$, and $Z$, where it appears that (speaking informally) $X < Y < Z$. Indeed, the p-value for the one-sided Wilcoxon signed-rank test for $X$ vs $Y$ is $p_{X<Y} = 8 \cdot 10^{-3}$, and that for $Y$ and $Z$ is $p_{Y<Z} = 4 \cdot 10^{-6}$. Yet I have $p_{X<Z} = 3 \cdot 10^{-4}$.

Clearly $p_{X < Z} \nless \min(p_{X<Y},p_{Y<Z})$, but what sort of relationship is guaranteed to hold? I would think that $p_{X < Z} \le \max(p_{X<Y},p_{Y<Z})$. Is this in fact the case? And is there something stronger that one can say?

  • $\begingroup$ You're doing signed rank tests on three sets of pair-differences? Is this with complete data (no missing values anywhere)? $\endgroup$
    – Glen_b
    Commented Sep 18, 2023 at 0:34
  • $\begingroup$ @Glen_b - yes and yes! $\endgroup$
    – calmcc
    Commented Sep 18, 2023 at 0:36
  • $\begingroup$ I'm not even sure you necessarily have transitivity in the population version of the one-sample Hodges-Lehmann statistic (median Walsh-average) on pair-differences across three populations. If you don't then p-values might end up anywhere. $\endgroup$
    – Glen_b
    Commented Sep 18, 2023 at 0:48
  • 1
    $\begingroup$ You do have transitivity. $\endgroup$ Commented Sep 18, 2023 at 0:56
  • $\begingroup$ ... at least in large samples, by the look of it. +1 $\endgroup$
    – Glen_b
    Commented Sep 18, 2023 at 4:18

2 Answers 2


Update: I'm no longer convinced that the ordering is transitive, but it is still described by the pseudomedians. It's just that when you don't have a location shift the pseudomedians need not be additive: you don't necessarily have $m_{XZ}=m_X-m_Z$, so $m_{XZ}=m_{ZY}-m_{YZ}$. I think everything else holds -- because most of the results below are negative rather than positive.

We can think in terms of the estimated pseudo-medians and their standard errors.

If the distributions are continuous and differ only by location shift, not in shape or scale, then the large-sample distributions of the pseudo-medians will be $N(m_i,\sigma^2/n)$ for $m_X$, $m_Y$ and $m_Z$ pseudo-medians and the same $\sigma^2$. Since $\hat m_Z-\hat m_X=(\hat m_Z-\hat m_Y)+(\hat m_Y-\hat m_X)$, the $X:Z$ $p$-value will be smaller than the minimum of the other two if the other two are both less than 0.5 (so that the pseudomedians are ordered the right way).

When the distributions differ in scale or shape it's more complicated, because the $\sigma^2$s will differ. For example, if $Y$ is less variable than $Z$ and $X$ it's a lot easier to get the $p$-values not to be ordered the way you want. For example, consider




You get quite a few examples where $p_{XY}>p_{XZ}>p_{YZ}$ even with all the $p$-values less than 0.5.

However, if $\hat m_X\sim N(m_X, \sigma^2_X/n)$ in large samples (and similarly for $Y$ and $Z$), we know that the variance of $m_X-m_Y$ is $(\sigma^2_X+\sigma^2_Y)/n$ (and similarly).

Now, $\sigma^2_X+\sigma^2_Z$ is not necessarily less than $\max\{\sigma^2_X+\sigma^2_Y, \sigma^2_Y+\sigma^2_Z\}$, so we might not expect even to have $p_{XZ}<\max \{ p_{XY}, p_{YZ}\}$, but it's going to be a pain to prove.

Simulation to the rescue:

> set.seed(2023-9-18)
> r<-replicate(10000,{x<-rnorm(500)
+ y<-rnorm(500,0.1,s=0.1)
+ z<-rnorm(500,0.2)
+ c(wilcox.test(x,y,paired=TRUE, alternative="less")$p.value,
+ wilcox.test(y,z,paired=TRUE, alternative="less")$p.value,
+ wilcox.test(x,z,paired=TRUE, alternative="less")$p.value)})
> which(!(r[3,]<pmax(r[1,],r[2,])))
[1] 1420 1745 2226 3959 5440 6158 9463 9844


> r[,1420]
[1] 0.4756392 0.5053687 0.6037722

Nope: one is bigger than 0.5

> r[,1745]
[1] 0.6671444 0.6444861 0.6835040

Nope: they are all bigger than 0.5

> r[,2226]
[1] 0.2439124 0.2258590 0.2484954

Success: it is possible to have $p_{XZ}<\max \{ p_{XY}, p_{YZ}\}$. The constraint does seem to be almost true, but not quite.

I think any sharp results will have to be on the Normal-approximation scale rather than on the $p$-value scale.


Does this meet the conditions of a counterexample? Or have I misunderstood something somewhere?

x <- c(3.51, 3.75, 3.914, 4.118, 4.77, 6.31, 6.53, 9.59, 15.117, 14.301)
y <- c(3.29, 3.32, 3.413, 3.66, 3.78, 10.143, 10.22, 10.31, 11.88, 12.4)
z <- c(0.01, 0.55, 5.22, 6.11, 7.115, 7.94, 8.48, 9.51, 10.503, 11.604)

wilcox.test(x,y,"less",paired=TRUE) # p-value = 0.6523
wilcox.test(y,z,"less",paired=TRUE) # p-value = 0.7217
wilcox.test(x,z,"less",paired=TRUE) # p-value = 0.7842

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.