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!!
 A: 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.
