Comparison Brunner-Munzel and Wilcoxon rank sum test (my first try)

try_brunnermunzel.R

library(brunnermunzel)
#https://cran.r-project.org/web/packages/brunnermunzel/vignettes/usage.html

## median of the difference: point a s.t. P(X-Y<a)+.5*P(X-Y=a) = 0.5
conv <- function(x){
  #numerical convolution
  integrate(function(y){dnorm(x+y,mu)*dexp(y)},-Inf,0)$value + 
    integrate(function(y){dnorm(x+y,mu)*dexp(y)},0,Inf)$value
}
conv <- Vectorize(conv)
m_root <- uniroot(function(x)integrate(conv,-Inf,x)$value-0.5, c(-1,0))
print(m_root$root)
#-0.1825963

mu <- qexp(0.5)
print(qnorm(0.5,mu,1))

set.seed(123)
Y1 = rnorm(10, mean = mu)
Y2 = rexp(10)

pvfun_b = function(loc,x,y){
  sapply(loc, function(m)brunnermunzel.test(x-m, y)$p.value)
}

pvfun_w = function(loc,x,y){
  sapply(loc, function(m)wilcox.test(x, y, mu = m)$p.value)
}

# png("pvalue_fun.png")
curve(pvfun_b(x,Y1,Y2), -2, 2, ylim=c(0,1), n=501, col="firebrick", ylab="p-value", lwd=1.5)
curve(pvfun_w(x,Y1,Y2), -2, 2, add = TRUE, n=501, col="steelblue", lty=4, lwd=1.5)
legend("topleft", c("brunnermunzel.test", "wilcox.test"), lty=c(1,4), col=c("firebrick","steelblue"), lwd=2)
points(m_root$root, 1,  pch=4, lwd=2)
# dev.off()

####

pv_b = numeric(10000)
pv_w = numeric(10000)
for(i in 1:10000){
  Y1 = rnorm(100, mean = mu, sd = 1)
  Y2 = rexp(100)
  pv_b[i] = brunnermunzel.test(Y1-m_root$root, Y2)$p.value
  pv_w[i] = wilcox.test(Y1, Y2, mu = m_root$root)$p.value
}

# png("dist_pvalue.png")
plot(ecdf(pv_b),col="firebrick", ylab="actual", xlab = "nominal", main="")
curve(ecdf(pv_w)(x), add=TRUE,col="steelblue", lty=4, lwd=1.5)
legend("topleft", c("brunnermunzel.test", "wilcox.test"), lty=c(1,4), col=c("firebrick","steelblue"), lwd=2)
# dev.off()