13 - Nonparametric Statistics: Wilcoxon and Mann-Whitney

You measured the salaries of seven people in one department and six in another, and one of the values is 9.5 thousand euros while all the others lie between 1.8 and 3.8. With a t-test you would compare two means, and that value moves the mean of the first by about one point. With a non-parametric test you ask something different: whether, picking one person at random from each department, the first earns more than the second roughly one time in two. These tests do not assume that the data follow a particular distribution, such as the normal: they work on the order of the values, not on their size. This article, the thirteenth in the Probability and Statistics series, computes by hand the Mann-Whitney test, the Wilcoxon test for paired data, the Kruskal-Wallis test and rank correlations, compares them with SciPy, and uses simulations to measure when it pays to use them and what they really test. If you do not know the p-value and t-tests, start there.

From values to ranks: the Mann-Whitney U test

The rank of a value is its position in the ordered list of all the data together: the smallest has rank 1, the next 2, and so on; tied values get the average of the positions they would occupy. The Mann-Whitney U statistic (1947) can be read in two equivalent ways: it is the number of pairs (a, b), with a from the first group and b from the second, in which a exceeds b (with half a point for ties), or it is the sum of the ranks of the first group minus the smallest possible value, n₁(n₁ + 1)/2. The ratio U / (n₁ n₂) is the probability that a value drawn from the first group exceeds one drawn from the second.

The listing computes U both ways and then obtains the exact p-value without formulas: under the null hypothesis (the two groups come from the same distribution) all the assignments of the 13 ranks to the two groups are equally likely, and the possible assignments are C(13, 6) = 1716, few enough to list. The p-value is the share of those that produce a U at least as far from the expected value as the observed one. There are two ways of defining «as extreme» and the code prints both, because the distribution of U with ties is not symmetric. Welch's test (a t-test that does not assume equal variances) serves as a comparison.

import itertools
import numpy as np
from scipy import stats

# Net monthly salaries (thousands of euros): department A (7 people, one high value), department B (6 people)
A = np.array([1.8, 2.1, 2.1, 2.4, 2.6, 3.0, 9.5])
B = np.array([2.2, 2.5, 2.9, 3.1, 3.4, 3.8])
n1, n2 = len(A), len(B)

# U = number of pairs (a, b) with a > b, counting 1/2 for each tie
U = sum((a > b) + 0.5 * (a == b) for a in A for b in B)
print(f"U(A) = {U}   out of {n1 * n2} pairs   P(A > B) = {U / (n1 * n2):.4f}")

# Same thing from the ranks: U = R_A - n1 (n1 + 1) / 2, with average ranks for ties
ranks = stats.rankdata(np.concatenate([A, B]))
R_A = ranks[:n1].sum()
print(f"sum of the ranks of A = {R_A}   U from ranks = {R_A - n1 * (n1 + 1) / 2}")

# Exact permutation p: try all C(13,6) = 1716 assignments of the ranks to the two groups
allU = np.array([ranks[list(idx)].sum() - n1 * (n1 + 1) / 2
                 for idx in itertools.combinations(range(n1 + n2), n1)])
mean_U = n1 * n2 / 2
p_distance = np.mean(np.abs(allU - mean_U) >= abs(U - mean_U) - 1e-12)   # "as far from the mean"
p_tails = 2 * min(np.mean(allU <= U + 1e-12), np.mean(allU >= U - 1e-12))  # twice the smaller tail
print(f"exact p: {int(p_distance * len(allU))}/{len(allU)} = {p_distance:.4f} (distance from the mean)   {p_tails:.4f} (twice the tail)")

asym = stats.mannwhitneyu(A, B, alternative="two-sided", method="asymptotic")
perm = stats.mannwhitneyu(A, B, alternative="two-sided", method=stats.PermutationMethod(n_resamples=np.inf))
print(f"scipy asymptotic (tie correction) p = {asym.pvalue:.4f}   scipy exact permutation p = {perm.pvalue:.4f}")
t = stats.ttest_ind(A, B, equal_var=False)
print(f"Welch t p = {t.pvalue:.4f}   means {A.mean():.2f} vs {B.mean():.2f}   medians {np.median(A):.2f} vs {np.median(B):.2f}")

# Output:
# U(A) = 12.0   out of 42 pairs   P(A > B) = 0.2857
# sum of the ranks of A = 40.0   U from ranks = 12.0
# exact p: 377/1716 = 0.2197 (distance from the mean)   0.2214 (twice the tail)
# scipy asymptotic (tie correction) p = 0.2240   scipy exact permutation p = 0.2214
# Welch t p = 0.7357   means 3.36 vs 2.98   medians 2.40 vs 3.00

U is 12 out of 42 pairs: in 28.6% of comparisons department A earns more than B. The two ways of computing give the same U. The exact p-value is 0.2197 with the «distance from the mean» rule and 0.2214 with the «twice the smaller tail» rule, the latter used by SciPy; the normal approximation with tie correction gives 0.2240. Three different figures and no different conclusion: with 13 people there is no evidence of a difference. Welch's t gives 0.7357, because the mean of department A (3.36) exceeds that of B (2.98) only because of the value 9.5, while the median of A (2.40) is lower than that of B (3.00). The case shows that the outcome of the two tests is not the same when an extreme value drives the mean, but not that the rank test is «more powerful»: here neither rejects. The exact method of mannwhitneyu does not correct for ties, as the documentation states, and that is why the listing uses PermutationMethod.

What Mann-Whitney really tests

In textbooks you often read that Mann-Whitney «compares medians». This is true only under one more assumption: that the two distributions have the same shape and differ by a shift. In general the test detects whether P(X > Y) departs from 1/2, that is, whether one variable tends to exceed the other, and two populations with identical medians can violate this condition. The following simulation shows it with X exponential and Y normal, both with median ln 2 = 0.6931. The exact probability P(X > Y), computed by numerical integration, is 0.5218, not 0.5.

import numpy as np
from scipy import stats, integrate

rng = np.random.default_rng(2027)
alpha = 0.05

# False belief: "Mann-Whitney compares medians". Two populations with the SAME median:
# X ~ exponential (median ln 2), Y ~ normal with mean ln 2 and standard deviation 0.3
med = np.log(2)
X, Y = stats.expon(), stats.norm(med, 0.3)
p_sup = integrate.quad(lambda y: Y.pdf(y) * X.sf(y), -np.inf, np.inf)[0]
print(f"medians: {X.median():.4f} and {Y.median():.4f}   P(X > Y) exact = {p_sup:.4f}")
for m in (50, 200, 1000):
    x = X.rvs((2000, m), random_state=rng); y = Y.rvs((2000, m), random_state=rng)
    rej = np.mean(stats.mannwhitneyu(x, y, axis=1, method="asymptotic").pvalue < alpha)
    print(f"n = {m:4d} per group: rejects H0 in {rej:.3f} of cases")

# Output:
# medians: 0.6931 and 0.6931   P(X > Y) exact = 0.5218
# n =   50 per group: rejects H0 in 0.081 of cases
# n =  200 per group: rejects H0 in 0.122 of cases
# n = 1000 per group: rejects H0 in 0.386 of cases

The medians coincide and the probability of superiority is 0.5218, but the test rejects the null hypothesis with a frequency that grows as the sample grows: 8.1% with 50 cases per group, 12.2% with 200 and 38.6% with 1000 (simulation error about 0.01). This is not a defect of the test, it answers a different question from the one you thought you were asking. If the question is «are the medians equal?», Mann-Whitney is not the right tool without the same-shape assumption; if it is «does one population tend to give higher values than the other?», it is. That is why it pays to report, next to the p-value, the estimate of P(X > Y), that is, U / (n₁ n₂), which also says by how much.

What you lose by giving up the t-test, and when you gain

If the data are truly normal the t-test is the most powerful of those possible; using ranks costs little. The asymptotic relative efficiency (ARE), that is, the ratio between the sample sizes of two tests that reach the same power, is 3/π ≈ 0.955 for the Wilcoxon test against the t-test under normality: you lose about 4.5% of efficiency. Hodges and Lehmann (1956) proved that, for any distribution, the ARE does not fall below 0.864. The gain, on the other hand, can be large when the tails are heavy, because a few extreme values inflate the variance used by the t-test, while ranks bound them.

The listing checks this with 4000 simulations per row and with 30 observations per group: same variance (equal to 1) and same shift of 0.6 standard deviations between the two groups, only the shape changes. It counts the fraction of simulations in which each test rejects the null hypothesis at the 5% level, that is, the power.

import numpy as np
from scipy import stats, integrate

rng = np.random.default_rng(2026)
n, R, alpha = 30, 4000, 0.05

def power(draw, shift):
    x = draw((R, n))
    y = draw((R, n)) + shift
    p_t = stats.ttest_ind(x, y, axis=1, equal_var=False).pvalue
    p_u = stats.mannwhitneyu(x, y, axis=1, method="asymptotic").pvalue
    return np.mean(p_t < alpha), np.mean(p_u < alpha)

cases = [("normal",           lambda s: rng.standard_normal(s),          0.6),
         ("Laplace (heavy tails)", lambda s: rng.laplace(size=s) / np.sqrt(2), 0.6),   # variance 1
         ("t with 3 df",      lambda s: rng.standard_t(3, s) / np.sqrt(3),  0.6)]      # variance 1
print("distribution             power t     power Mann-Whitney")
for name, f, d in cases:
    a, b = power(f, d)
    print(f"{name:24s} {a:8.3f}   {b:8.3f}    (error ~0.008)")

# Output:
# distribution             power t     power Mann-Whitney
# normal                      0.628      0.594    (error ~0.008)
# Laplace (heavy tails)       0.636      0.749    (error ~0.008)
# t with 3 df                 0.696      0.856    (error ~0.008)

With normal data the power is 0.628 for the t-test and 0.594 for Mann-Whitney: a loss of 3-4 percentage points, consistent with 0.955. With the Laplace distribution (heavier tails) the roles are reversed, 0.636 against 0.749, and with Student's t with 3 degrees of freedom, which has even heavier tails, the gap reaches 0.696 against 0.856. The simulation error is about 0.008. The practical rule that comes out is asymmetric: the cost of using ranks when they were not needed is small, while the cost of using the t-test on heavy tails reaches 16 percentage points of power (0.696 against 0.856). With few data, moreover, normality cannot be checked well, and the rank test is a prudent choice.

Paired data: Wilcoxon's signed-rank test

When the measurements are paired, such as the same group of pages before and after an intervention, you work on the differences d = before − after. Wilcoxon's test (1945) orders the absolute values |d| and sums their ranks separately for the positive differences (W⁺) and the negative ones (W⁻); their sum is always n(n + 1)/2. Under the null hypothesis, each difference has a positive or negative sign with probability 1/2, so with 8 pairs the possible assignments are 2⁸ = 256 and can be listed. The test assumes that the distribution of the differences is symmetric about zero under the null hypothesis.

import itertools
import numpy as np
from scipy import stats

# Loading time of 8 pages (seconds) before and after an intervention
before = np.array([4.1, 3.8, 5.0, 2.9, 6.2, 3.5, 4.7, 5.5])
after  = np.array([3.5, 3.9, 4.1, 2.6, 5.0, 3.0, 4.5, 4.2])
d = before - after                     # paired differences: the test works on these
print("differences:", np.round(d, 2))

ranks = stats.rankdata(np.abs(d))      # ranks of the absolute values
W_plus = ranks[d > 0].sum()            # sum of the ranks of the positive differences
W_minus = ranks[d < 0].sum()
print("ranks |d| :", ranks, "  W+ =", W_plus, " W- =", W_minus, " (W+ + W- =", ranks.sum(), ")")

# exact p: under H0 each difference has sign + or - with prob. 1/2, hence 2^8 = 256 assignments
sums = np.array([sum(r for r, s in zip(ranks, signs) if s)
                 for signs in itertools.product([False, True], repeat=len(d))])
W = min(W_plus, W_minus)
p_exact = np.mean(sums <= W) * 2       # the distribution is symmetric: double the tail
print(f"exact p by enumeration = {int(np.sum(sums <= W))}/{len(sums)} x 2 = {p_exact:.4f}")

res = stats.wilcoxon(before, after, method="exact")
print(f"scipy.stats.wilcoxon: statistic = {res.statistic}, p = {res.pvalue:.4f}")
print(f"paired t: p = {stats.ttest_rel(before, after).pvalue:.4f}")

# With a zero difference the pair is dropped (zero_method='wilcox'): 7 differences remain
after2 = after.copy(); after2[1] = before[1]
r2 = stats.wilcoxon(before, after2)
print(f"with one zero difference: statistic = {r2.statistic}, p = {r2.pvalue:.4f}")

# Output:
# differences: [ 0.6 -0.1  0.9  0.3  1.2  0.5  0.2  1.3]
# ranks |d| : [5. 1. 6. 3. 7. 4. 2. 8.]   W+ = 35.0  W- = 1.0  (W+ + W- = 36.0 )
# exact p by enumeration = 2/256 x 2 = 0.0156
# scipy.stats.wilcoxon: statistic = 1.0, p = 0.0156
# paired t: p = 0.0096
# with one zero difference: statistic = 0.0, p = 0.0156

Of the eight differences, only one is negative (−0.1) and has rank 1, so W⁻ = 1 and W⁺ = 35. Only 2 of the 256 assignments of signs produce a sum of negative ranks equal to or below 1, and the two-sided p-value is 2 · 2/256 = 0.0156, identical to that of scipy.stats.wilcoxon. The paired t gives 0.0096: same conclusion, different p-value. A real trap: two differences that are equal by hand, such as 4.1 − 3.6 and 3.5 − 3.0 (both 0.5), are worth 0.49999999999999956 and 0.5 in floating point, and the two values receive different ranks instead of a common average rank. The SciPy documentation warns that it is better to compute the difference, round it, and pass it to the test as the only argument. Pairs with a zero difference are dropped (zero_method='wilcox', the default): the last line shows it with 7 useful differences.

More than two groups: Kruskal-Wallis and pairwise comparisons

The generalisation of Mann-Whitney to k groups is the Kruskal-Wallis test (1952), the non-parametric analogue of the one-way ANOVA seen in the article on ANOVA and multiple comparisons. You order all N data together and compute H = 12 / (N(N + 1)) · Σ R_i² / n_i − 3(N + 1), where R_i is the sum of the ranks of group i and n_i its size; for samples that are not too small, H is compared with a chi-square distribution with k − 1 degrees of freedom. If the test rejects, it says that at least one group tends to give different values, not which one: to find out you make pairwise comparisons, and here you need to correct the p-values for multiple comparisons, otherwise the probability of at least one false positive exceeds the stated level. Holm's correction orders them from the smallest, multiplies the first by the number of comparisons m, the second by m − 1 and so on, and makes the result non-decreasing.

import numpy as np
from scipy import stats

# Response times (ms) of three versions of a service, 5 measurements each
g = [np.array([212, 230, 245, 251, 300.0]),
     np.array([260, 275, 290, 310, 335.0]),
     np.array([280, 320, 340, 360, 410.0])]
everyone = np.concatenate(g)
N = len(everyone)
ranks = stats.rankdata(everyone)
R = [ranks[:5], ranks[5:10], ranks[10:]]
H = 12 / (N * (N + 1)) * sum(r.sum() ** 2 / len(r) for r in R) - 3 * (N + 1)   # without ties
p_chi2 = stats.chi2.sf(H, df=len(g) - 1)
kw = stats.kruskal(*g)
print(f"rank sums: {[float(r.sum()) for r in R]}   H by hand = {H:.4f}  p (chi2, 2 df) = {p_chi2:.4f}")
print(f"scipy.stats.kruskal: H = {kw.statistic:.4f}  p = {kw.pvalue:.4f}")

# Pairwise comparisons with Holm's correction: sort the p's, multiply by m, m-1, ..., 1 and make monotone
pairs = [(0, 1), (0, 2), (1, 2)]
p = np.array([stats.mannwhitneyu(g[i], g[j], method="exact").pvalue for i, j in pairs])
order = np.argsort(p)
m = len(p)
holm = np.empty(m); current = 0.0
for rank, k in enumerate(order):
    current = max(current, min(1.0, (m - rank) * p[k])); holm[k] = current
for (i, j), pg, ph in zip(pairs, p, holm):
    print(f"version {i + 1} vs {j + 1}: p = {pg:.4f}  corrected p (Holm) = {ph:.4f}")

# Output:
# rank sums: [19.0, 41.0, 60.0]   H by hand = 8.4200  p (chi2, 2 df) = 0.0148
# scipy.stats.kruskal: H = 8.4200  p = 0.0148
# version 1 vs 2: p = 0.0556  corrected p (Holm) = 0.1111
# version 1 vs 3: p = 0.0159  corrected p (Holm) = 0.0476
# version 2 vs 3: p = 0.0952  corrected p (Holm) = 0.1111

The ranks of the three groups sum to 19, 41 and 60, and H is 8.42 both with the formula by hand and with SciPy (which applies the tie correction, absent here), with p = 0.0148 under the chi-square approximation. In the pairwise comparisons only version 1 against version 3 stays below 0.05 after Holm's correction (0.0159 becomes 0.0476), while 1 against 2 goes from 0.0556 to 0.1111. The global test and the pairwise comparisons do not have to agree: Kruskal-Wallis uses all the data together, the pairwise comparisons use part of them and pay for the correction. There is also Dunn's post-hoc test, which uses the global ranks, but it is not in SciPy's standard library and is not run here.

Rank correlations: Spearman and Kendall

Pearson's correlation measures linear association and is sensitive to outliers. Rank correlations measure monotone association, that is, how far one variable grows together with the other without requiring that it do so along a straight line. Spearman's ρ is Pearson's correlation computed on the ranks. Kendall's τ (1938) is (C − D) / (number of pairs), where C counts the concordant pairs (the two variables move in the same direction) and D the discordant ones.

import numpy as np
from scipy import stats

# Rank correlations on a monotone but non-linear relationship, with an outlier
x = np.arange(1.0, 11.0)
y = np.exp(0.5 * x)                          # increasing, curved
y_out = y.copy(); y_out[-1] = 1.0            # the last point collapses
def kendall(a, b):                           # (concordant - discordant) / pairs, no ties
    c = d = 0
    for i in range(len(a)):
        for j in range(i + 1, len(a)):
            s = (a[i] - a[j]) * (b[i] - b[j]); c += s > 0; d += s < 0
    return (c - d) / (len(a) * (len(a) - 1) / 2)
for name, yy in (("exponential curve", y), ("with outlier", y_out)):
    sp = np.corrcoef(stats.rankdata(x), stats.rankdata(yy))[0, 1]
    print(f"{name:17s} Pearson {stats.pearsonr(x, yy)[0]:.4f}  Spearman {sp:.4f} (scipy {stats.spearmanr(x, yy)[0]:.4f})  "
          f"Kendall {kendall(x, yy):.4f} (scipy {stats.kendalltau(x, yy)[0]:.4f})")

# Output:
# exponential curve Pearson 0.8654  Spearman 1.0000 (scipy 1.0000)  Kendall 1.0000 (scipy 1.0000)
# with outlier      Pearson 0.5920  Spearman 0.4545 (scipy 0.4545)  Kendall 0.6000 (scipy 0.6000)

On an increasing exponential curve Pearson is 0.8654 because the relationship is not a straight line, while Spearman and Kendall are 1.0000: the order is perfectly respected. If the last point collapses to 1, Pearson drops to 0.5920, Spearman to 0.4545 and Kendall to 0.6000: all three are affected by the outlier, but to different degrees, and no coefficient is «robust» in absolute terms. The calculation by hand coincides with SciPy up to the fourth digit. Which to choose? For data without ties and with few observations τ has a direct interpretation as the difference between the probability of concordance and of discordance of a pair; ρ is used more out of tradition.

A worked exercise and a test to do on your own

Five pages, time differences (seconds) +1.5; −0.5; +2.0; +0.8; +3.0. The sorted absolute values are 0.5 / 0.8 / 1.5 / 2.0 / 3.0 with ranks 1-5, so W⁻ = 1 (the rank of 0.5). Of the 2⁵ = 32 assignments of signs, two have negative sum ≤ 1 (no negative sign, or only rank 1), so the exact two-sided p-value is 2 · 2/32 = 0.125. With five pairs the smallest possible value is 2/32 = 0.0625: no result can reject at 5%, however sharp the difference. The test was checked with scipy.stats.wilcoxon(..., method='exact'), which returns 0.125.

What you can do now

You can compute ranks and the U statistic, obtain the exact p-value by enumeration, and read U as the probability that one group exceeds the other. You know that Mann-Whitney tests stochastic superiority and not medians, that it costs about 4.5% of efficiency under normality and can gain a lot with heavy tails, and you can apply Wilcoxon, Kruskal-Wallis with Holm's correction, and rank correlations. The next article, the bootstrap, builds intervals and tests with no assumption about the distribution, by resampling the data. The chi-square for frequency tables is in article 9; the Markov chain of article 12 comes before this one.

The same series also includes Exploratory Data Analysis (EDA) with Python.