This is the code to perform the Bhapkar V test. I’ve rapidly wrote it, in 2 hours. The code is then quite brutal and it could be done better. As soon as possible, I will correct it.

WARNING: it works *ONLY* with 3 groups, for now!

bhapkar.test.3g <- function(data1=list){

sample <- c()
for(i in 1:length(data1)){
sample <- c(sample, rep(i, length(data1[[i]])))
}

obs <- c()
for(i in 1:length(data1)){
obs <- c(obs, data1[[i]])
}
rank <- rank(obs)

cplets <- list()
vec <- c()
for(i in 1:length(data1[[1]])){
vec <- c(vec, (length(data1[[2]][data1[[2]]>data1[[1]][i]]) * length(data1[[3]][data1[[3]]>data1[[1]][i]])))
}
cplets[[1]] <- vec

vec <- c()
for(i in 1:length(data1[[2]])){
vec <- c(vec, (length(data1[[1]][data1[[1]]>data1[[2]][i]]) * length(data1[[3]][data1[[3]]>data1[[2]][i]])))
}
cplets[[2]] <- vec

vec <- c()
for(i in 1:length(data1[[3]])){
vec <- c(vec, (length(data1[[2]][data1[[2]]>data1[[3]][i]]) * length(data1[[1]][data1[[1]]>data1[[3]][i]])))
}
cplets[[3]] <- vec

cplets1 <- c(cplets[[1]], cplets[[2]], cplets[[3]])
mydata <- data.frame(obs=obs, sample=sample, rank=rank, cplets=cplets1)

v1 <- sum(cplets[[1]])
v2 <- sum(cplets[[2]])
v3 <- sum(cplets[[3]])

vtot <- v1+v2+v3
u1 <- v1/vtot
u2 <- v2/vtot
u3 <- v3/vtot
u <- c(u1,u2,u3)

lengths <- c(length(data1[[1]]), length(data1[[2]]), length(data1[[3]]))
N <- sum(lengths)
P <- c(lengths / N)
ngroup <- length(data1)

V <- N * (2*length(data1)-1)* (sum(P*((u-1/ngroup)^2)) - (sum(P*((u-1/ngroup))))^2)

prop <- pchisq(V, df=length(data1)-1)
names(V) = "V = "
method = "Bhapkar V-test"
rval <- list(method = method, statistic = V, p.value = prop)
class(rval) = "htest"
return(rval)



}

An example:

a <- c(42, 46, 48.5, 49, 68, 51)
b <- c(70.5, 54, 60,72)
c <- c(66, 54, 43, 105, 94)

mydata <- list(a,b,c)

bhapkar.test.3g(mydata)


        Bhapkar V-test

data:  
V = 6.7713, p-value = 0.9661

REFERENCES:
Statistical analysis of nonnormal data
By J. V. Deshpande, A. P. Gore, A. Shanubhogue
pag. 61