Mann-Whitney U Test vs. simple Bayesian inference

Specific Domains Statistics
,
3

Your proposed approach might work in cases with few ties, but it will be downwardly biased in the presence of many ties. For example, you wouldn’t want to conclude that the distributions are different if PR(X>Y) = 0.01, PR(X<Y) = 0.01, and PR(X=Y) = 0.98.

When the data are discrete and bounded, a straight-forward way to conduct such a test is to model the two vectors as categorical. Then, you can use the posterior draws to compare PR(X>Y) to PR(X<Y) without PR(X=Y) contaminating the results. This also motivates the Bayesian version, since you can apply whatever prior you want to the probabilities, including one that shrinks the two sets towards each other. You could add a smoother to the probabilities. etc.

using Plots
using StatsPlots
using Turing
using HypothesisTests

# Generate data; y is slightly higher than x
x = rand([1,1,2,3,4], 20)
y = rand([1,2,3,4,4], 20)

# Typical Mann-Whitney U test
MannWhitneyUTest(x, y)

# Bayesian model for category probabilities
# Assumes only categorical/multinomial
# distribution
@model function mann_whitney(x,y)
    px ~ Dirichlet([1,1,1,1])
    py ~ Dirichlet([1,1,1,1])

    x ~ Categorical(px)
    y ~ Categorical(py)
end

mw_mod = mann_whitney(x, y)
mw_results = sample(mw_mod, NUTS(1_000, 0.85), 1_000)

# Extract the category probabilities
px1 = Array(mw_results[["px[1]"]])
px2 = Array(mw_results[["px[2]"]])
px3 = Array(mw_results[["px[3]"]])
px4 = Array(mw_results[["px[4]"]])
py1 = Array(mw_results[["py[1]"]])
py2 = Array(mw_results[["py[2]"]])
py3 = Array(mw_results[["py[3]"]])
py4 = Array(mw_results[["py[4]"]])

# Probability that X and Y are at the same level (ties)
px_equal = px1 .* py1 +
           px2 .* py2 +
           px3 .* py3 +
           px4 .* py4
# Probability that X is higher than Y
px_higher = px2 .* py1 +
            px3 .* (py1 + py2) +
            px4 .* (py1 + py2 + py3)
# Probability that X is lower than Y
px_lower = px1 .* (py2 + py3 + py4) +
           px2 .* (py3 + py4) +
           py3 .* py4

# Summarize: PR(X>Y) - PR(X<Y)
higher_lower = px_higher - px_lower
mean(higher_lower)
std(higher_lower)
quantile(higher_lower, [0.025, 0.5, 0.975])
mean(higher_lower .> 0)

trace_plot = plot(higher_lower)
density_plot = plot(higher_lower, seriestype = "density")
plot(trace_plot, density_plot, layout = (1, 2), legend = false)
1 Like