# What is pseudomedian in R function wilcox.test?

On my journey to understand what pseudomedian is I tried to calculate it manually in R using definition from wikipedia. Calculated on my data I found that reported value of wilcox.test was quite different than value I calculated as median of pairwise means of the data (0.275 vs 0.33 while the confidence interval is 0.330; 0.345). I wanted to ask what am I missing, so I tried to create a reproducible example for this Q and I found that the reported value were almost the same for the data I have generated :

set.seed(910401)
distr_data <- rgamma(100,0.1,1)

wilcox.test(distr_data, conf.int = T)
# 0.006787143

all_pairs <- t(combn(distr_data, 2))
all_pair_means <- (all_pairs[,1] + all_pairs[,2]) / 2
median(all_pair_means)
# 0.006970087


So I played with gamma distribution a bit after trying different sample sizes and parameters I eventually found one that also generates quite different results :

distr_data <- rgamma(100,0.1,1000)

wilcox.test(distr_data, conf.int = T)
# 6.533661e-05
# CI 6.049175e-05 8.016369e-05

all_pairs <- t(combn(distr_data, 2))
all_pair_means <- (all_pairs[,1] + all_pairs[,2]) / 2
median(all_pair_means)
# 9.181717e-06


So here comes my question. What exactly is a pseudomedian in wilcox.test? What does it mean when wilcox.test and wikipedia pseudomedians disagree?

• Although this is framed in terms of R, I am voting to leave it open since it is actually about statistics. – Peter Flom Mar 22 '18 at 11:58

A good clue is to look at the actual code of wilcox.test:

https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/wilcox.test.R

Specifically the bit concerning the estimate of the pseudo median for a 1 sample test is line 91-122:

x <- x + mu             # we want a conf.int for the median
alpha <- 1 - conf.level
diffs <- outer(x, x, "+")
diffs <- sort(diffs[!lower.tri(diffs)]) / 2
...
ESTIMATE <- c("(pseudo)median" = median(diffs))


Your code is typically missing the pairs of twice the same elements when you do all_pairs <- t(combn(distr_data, 2))

set.seed(910401)

Notice that I also added exact==T which you would not want to do with large dataset but actually matters in this case, otherwise your estimate is slightly off.