R removing zeros for pseudomedian and its confidence interval in wilcox.test?

It makes sense to me that we need to remove zeros to calculate the p-value in a Wilcoxon signed rank test. What confuses me is that R seems to leave the zeros deleted when computing the pseudomedian and its accompanying confidence interval (and I have checked this by stepping through the code). The estimated pseudomedian below feels to me like it makes no sense for the data that goes into it. Is it logical for R to do this? If so, is this just a reflection of this data being too far from continuous for Wilcoxon to be useful in this case? And if so, what is the right way to compute a confidence interval for the difference between two ordinal variables (or, in general, discrete variables with very few levels)? The diffs in this case were post-pre diffs on an ordinal scale in some real world data.

> diffs <- c(-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1)

> table(diffs)
diffs
-2 -1  0  1
1 13 16  2

> wilcox.test(diffs,conf.int=TRUE)

Wilcoxon signed rank test with continuity correction

data:  diffs
V = 16, p-value = 0.00314
alternative hypothesis: true location is not equal to 0
95 percent confidence interval:
-1.0000148 -0.9999398
sample estimates:
(pseudo)median
-0.9999517


Any help in sorting this out in my head much appreciated!!

If you want to report the median of differences and the confidence interval for this statistic, then that's what you should do. For data that are discrete with few levels, there may not great method for doing this, but I'll present a couple of methods in R below. For these data, the median is 0 and reasonable 95% confidence limits might be -1 and 0.

Note that if you truly have (merely-) ordinal data, that you can't have differences in values † . That is, if you have ("worst", "bad", "neutral", "good", "best"), it's not the case that "bad" – "worst" is one "goodness unit".

The first step in the Wilcoxon signed-rank test is taking the difference in values. So at that step, you've already made the data interval (more than ordinal, if you will). If you want to treat the data as simply ordinal, you could do the analysis with ordinal regression, for example with the ordinal package in R. On the other hand, if you are treating your data as interval, you could probably report the mean of the differences and construct a confidence interval for this statistic.

Edit: The following is my code to mirror the method to calculate a confidence interval for the median from Conover, Practical Nonparametric Statistics, 3rd, based on the binomial distribution. I think because it's based on the binomial distribution, it should be valid for discrete values, and those with few levels. I'll add this function to the next release of the rcompanion package.

source("http://rcompanion.org/r_script/quantileCI.r")

diffs <- c(-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1)

quantileCI(diffs, tau=0.50, level=0.95)

### tau  n Quantile Nominal.level Actual.level Lower.ci Upper.ci
### 0.5 32        0          0.95        0.965       -1        0


It also works for ordinal data, where the data are coded as an ordered factor.

set.seed(12345)
Pool = factor(c("smallest", "small", "medium", "large", "largest"),
ordered=TRUE,
levels=c("smallest", "small", "medium", "large", "largest"))
Sample = sample(Pool, 24, replace=TRUE)

quantileCI(Sample)

### tau  n Quantile Nominal.level Actual.level Lower.ci Upper.ci
### 0.5 24   medium          0.95        0.957    small    large


Another method would be to use bootstrap. This method is probably not particularly valid with discrete values with few levels, but in practice, the percentile method is similar to the method above. The following uses the percentile method.

require(boot)

Function = function(input, index){
Input = input[index]
return(median(Input))}

Boot = boot(diffs, Function, R=5000)

hist(Boot\$t[,1], col="darkgray")

boot.ci(Boot, conf = 0.95, type = "perc")

### BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
### Based on 5000 bootstrap replicates
###
### Intervals :
### Level     Percentile
### 95%   (-1,  0 )


† Well, if you knew the relative differences in the ordinal values, you could have some way of ordering these differences, even if you couldn't give them interval values, but this would be unusual.

• Thanks a lot! The MedianCI is a new function to me, so that is really interesting. However, I note that you say that neither is ideal for discrete data with few levels, and that is where I really get stuck... Is there an approach that is appropriate for discrete data with few levels? That is perhaps a better description than "ordinal data" -- it's the lack of levels that troubles me with all these approaches and I'm not sure how really to proceed (while, at the same time, this is a very similar statistical question I see many colleagues frequently trying to answer). Thank you! – justme Jan 29 at 16:44
• (I've updated the question to reflect that my true worry is about how to deal with there being so few levels. I do take your point about ordinal data not suiting differences though! That is an extra layer or wrongness in the analysis...) – justme Jan 29 at 16:47
• I don't know if there is any recommended method in this case. I just don't know. I guess it depends what you what expect the confidence interval to tell you. With your data, no matter how you approach it (except the wilcox.test approach !), the median is 0, and the only 95% confidence interval that makes sense is -1 to 0. To me, the bootstrap approach with the percentile method is the most understandable. You resample the data a bunch, and 95% of the medians fall within a range. But I think this approach is not recommended if the distribution of the statistic isn't fairly continuous. – Sal Mangiafico Jan 29 at 18:27
• But this is one of the things that baffles me. I can't see how the non-zero median popping out of the Wilcox approach can be anything other than a failure of the method (or a bug in the code)? – justme Jan 29 at 19:41
• My understanding is that the pseudo-median returned by wilcox.test is the Hodges–Lehmann estimator. Looking for information about that may be helpful. – Sal Mangiafico Jan 29 at 20:39