Calculate the Mann-Whitney test statistic in bootstrap samples using Hodges-Lehmann shift estimator

bootstraps_mann_whitney.R

bootstraps_mann_whitney <- function(bootstraps_number, sample_1, sample_2) {
  differences <- numeric()
  for (i in sample_1) {
    for (j in sample_2) {
      differences[1 + length(differences)] <- i - j
    }
  }
  bootstraps_statistic <- numeric(bootstraps_number)
  hodges_lehmann_estimator <- median(differences)
  sample_1_delta <- sample_1 - hodges_lehmann_estimator
  for (i in 1:bootstraps_number) {
    sample_1_bootstraps <- sample(sample_1_delta, length(sample_1_delta), 
                                  replace = TRUE)
    sample_2_bootstraps <- sample(sample_2, length(sample_2), 
                                replace = TRUE)
    mann_whitney_test <- wilcox.test(sample_1_bootstraps, 
                                     sample_2_bootstraps, 
                                     paired = FALSE, 
                                     alternative = "two.sided")
    bootstraps_statistic[i] <- mann_whitney_test$statistic
  }
  bootstraps_statistic
}

foo <- c(1,3,5,2,4,2,1,4,1,4,6)
bar <- c(5,4,5,3,4,3,5,4,5,6,4)

test_statistic <- wilcox.test(foo, bar, paired = FALSE, alternative = "two.sided")$statistic
bootstrap_test_statistic <- bootstraps_mann_whitney(5000, foo, bar)

tail <- mean(test_statistic >= bootstrap_test_statistic)
pval <- 2 * min(tail, 1 - tail)

# Reference:
# Boos, Dennis D., and Cavell Brownie. "Bootstrap p-values for tests of nonparametric hypotheses." (1988).