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?
 A: 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))
Instead try:
set.seed(910401)
distr_data <- rgamma(100,0.1,1000)
wilcox.test(distr_data, conf.int = T, exact = T)
# 6.788116e-06
all_pairs <- rbind(t(combn(distr_data, 2)),cbind(distr_data,distr_data)) 
all_pair_means <- (all_pairs[,1] + all_pairs[,2]) / 2
median(all_pair_means)
# 6.788116e-06

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.
