6
$\begingroup$

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?

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

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.