# Transitivity of Wilcoxon signed-rank test

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, and that for $$Y$$ and $$Z$$ is $$p_{Y. Yet I have $$p_{X.

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

• You're doing signed rank tests on three sets of pair-differences? Is this with complete data (no missing values anywhere)? Sep 18 at 0:34
• @Glen_b - yes and yes! Sep 18 at 0:36
• 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. Sep 18 at 0:48
• You do have transitivity. Sep 18 at 0:56
• ... at least in large samples, by the look of it. +1 Sep 18 at 4:18

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

r<-replicate(1000,{x<-rnorm(500)
y<-rnorm(500,0.15,s=0.1)
z<-rnorm(500,0.4)

c(wilcox.test(x,y,paired=TRUE)$$p.value, wilcox.test(y,z,paired=TRUE)$$p.value,
wilcox.test(x,z,paired=TRUE)$p.value)}) table(r[3,]<pmin(r[1,],r[2,])) table(r[3,]<pmax(r[1,],r[2,]))  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


Now:

> 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