Non-parametric tests: Wilcoxon and Kruskal-Wallis

What this method is for

“The data do not look normal, so I switched to Wilcoxon” is how a rank test usually enters a paper, and it almost always arrives with an unstated misunderstanding attached: the belief that Wilcoxon compares medians.

It does not. The Wilcoxon rank-sum test (also known as the Mann-Whitney U test) is about stochastic ordering: draw one person at random from group A and one at random from group B, and ask whether the probability that the first exceeds the second is one half. As a formula, the null hypothesis is $P(X > Y) = 1/2$. A median is one summary of a distribution; this probability depends on the whole distribution.

Most of the time the two agree, but they are not the same thing — and when they come apart, the results table carries a pair of numbers that look self-contradictory: two identical medians in Table 1 with a significant p-value underneath. Section two below builds exactly that situation.

This page uses two datasets: MASS::birthwt (189 mothers, outcome birth weight in grams) and datasets::sleep (Student’s own 1908 soporific data). The variable definitions are the same ones used on t-tests and analysis of variance, so numbers on the two pages can be read against each other.

The null hypothesis is stochastic ordering, not equal medians

The data below are constructed on purpose, to demonstrate something that really does happen. The setting is length of stay in days after the same operation at two hospitals:

• Hospital A has no fixed discharge pathway, so each patient goes home when they are ready. That gives a long lower tail (plenty of patients leave the day after surgery) and a ceiling at 9 days.
• Hospital B runs a pathway targeting discharge on day seven, so almost nobody leaves before day 6 and all the remaining variation runs upward — complications keep a few patients for three or four weeks.

Hospital A (30 patients): 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9

Hospital B (30 patients): 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 9, 10, 11, 12, 13, 14, 17, 18, 20, 22, 25

The difference between the medians is 0 days — exactly nothing. The Wilcoxon rank-sum test returns W = 629, p = 0.007.

The middle panel carries the argument. Hospital B’s cumulative distribution function never rises above hospital A’s — wherever you set a threshold, the proportion of hospital B patients above it is never lower. That is the strict meaning of “hospital B tends to keep people longer”, and it is entirely compatible with equal medians: both step functions cross 0.5 on the step at day 7.

For completeness, the difference in means on the same data is 3.87 days, with a Welch t-test p < 0.001. What the mean can see, the median cannot — and the handful of patients in that right tail who stayed three weeks are precisely the ones the clinic needs to see.

Turning stochastic ordering into a number you can report

$P(X > Y)$ does not have to stay abstract; it can be estimated directly. Take every cross-group pair, score one when the first is larger and half each when they are equal, and divide by the number of pairs. Panel C above draws exactly that: 900 cells, one per pair of patients.

For the constructed example the estimate is 0.699: pick one patient from each hospital at random and the one from hospital B stays longer about 70% of the time.

Real data works the same way. For birth weight by maternal smoking in birthwt, the probability is 0.617 — pick one mother-infant pair from each group and the infant of the non-smoker is heavier about 62% of the time.

Its virtue is that you can say it out loud: “pick a patient on the new treatment at random and a patient on standard care at random, and the first recovers better seven times out of ten” carries far more information than a p-value, and it needs no distributional assumption at all.

A p-value is not enough: the Hodges-Lehmann shift estimate

By default the Wilcoxon output is a p-value and nothing else, so the results section becomes “median A vs B (p = …)”. That sentence stitches together two different things: the medians are summaries of each group, the p-value is a test of the comparison, and the effect size — how far apart, and how precisely — is missing from between them.

What fills that gap is the Hodges-Lehmann shift estimate: the median of all cross-group pairwise differences. wilcox.test(conf.int = TRUE) in R returns it with a confidence interval, in the units of the original data.

Taking smoking in birthwt as the example, one comparison yields three different numbers:

All three land near three hundred grams, but they are not one quantity. The middle row is worth reading twice: Hodges-Lehmann estimates “the median of the differences”, not “the difference between the medians” — R’s own documentation labels the latter a common misconception. The two are close under symmetric distributions and are not close under skewed ones.

In practice, write it as “Hodges-Lehmann shift estimate 306.2 g (95% CI 85 to 512; Wilcoxon rank-sum p = 0.007)”, adding the estimate of $P(X > Y)$, 0.617, where it helps. As a check: computing the median of the 8510 pairwise differences directly gives 307 g, very close to the 306.2 g that wilcox.test() reports, though not the same number. The difference is one of algorithm: wilcox.test() inverts the test and solves with uniroot rather than taking the median of those 8510 differences.

Rank test or log transform? Two different estimands

A right-skewed variable — length of stay, cost, biomarker concentration — has a second route available: take logs and run a t-test. Both routes are defensible, but they do not produce the same quantity, and a paper has to say which one it took.

On the same birthwt data the two routes give:

The second row says: the geometric mean birth weight among non-smokers (2950.8 g) is 1.101 times that among smokers (2679.8 g), which as a percentage is 10.1% higher (95% CI 1.4% to 19.6%), p = 0.022.

The difference is more than a change of units:

• An additive shift assumes the groups differ by a fixed amount. Birth weight suits that reading, because the clinical threshold (low birth weight, 2500 g) is additive and “300 g lighter” means roughly the same thing for every infant.
• A ratio assumes the groups differ by a fixed multiple. Quantities spanning orders of magnitude — viral load, C-reactive protein, hospital cost — suit that reading, because “10% higher” holds at both ends while “300 units higher” may be a tripling at the low end and negligible at the high end.
• The log-scale t-test estimates a ratio of geometric means, not of arithmetic means. Back-transforming a mean of logs gives the geometric mean, which is always less than or equal to the arithmetic mean; a claim about arithmetic means needs separate handling.
• Logs require every value to be positive. The usual log(x + 1) fix for zeroes lets that constant influence the answer, and the ratio loses its clean interpretation.

Ties and the exact p-value

The exact p-value of a rank test comes from every way of reallocating the N observations to the two groups. When the data contain ties, the tied observations can only be given the same average rank, which changes the permutation reference distribution.

Birth weight in birthwt is recorded to the gram and still has ties: 131 distinct values across 189 records, 97 of which share a value with somebody else, the largest group holding 5. Clinical data is essentially never tie-free.

Most textbooks say that ties make the exact p-value impossible and the software falls back on the normal approximation. That sentence is no longer true of current R (this page was verified on R 4.6.0): wilcox.test() carries the Streitberg-Röhmel shift algorithm, which under ties switches to the permutation distribution conditional on the observed ranks and still delivers exact inference.

What actually flips the default is sample size, not ties: exact when both samples hold fewer than 50 values, the normal approximation otherwise. So the constructed example on this page (30 per group) receives an exact p-value, while the smoking comparison in birthwt (115 against 74) receives the normal approximation.

Pulling out the rows whose birth weight is unique across the whole dataset (92 records, no ties at all) and then rounding that same subset to the nearest 250 g — the precision a chart review would realistically give you — shows what ties cost:

Same patients, same comparison; measured more coarsely, and both the test statistic and the p-value move. Measurement precision is part of the analysis, not something that stops mattering once data collection ends.

A few practical points:

• Report the software and its version. The same data under older R, SPSS or SAS may return a normal approximation rather than an exact p-value, and the two do not agree exactly.
• The continuity correction only affects the approximation route. The smoking comparison gives p = 0.00677 with the correction and p = 0.00674 without; the gap is usually tiny, but not when a p-value sits against a threshold.
• Do not jitter the data to break ties. That hands the answer to the random seed.
• Past a certain density of ties the fix is a model, not a test: when the outcome is genuinely an ordered category (pain grade, NYHA class), ordinal logistic regression gives an effect size and covariate adjustment at the same time.

The paired version: Wilcoxon signed-rank

The same person measured before and after, the same person’s two eyes, a matched design — a one-to-one link means the paired version, for exactly the reason given for the paired t-test: pairing removes between-person variation entirely.

The Wilcoxon signed-rank test works on each person’s difference: rank the differences by absolute size, then ask whether the sum of the ranks attached to positive differences is off balance. On Student’s soporific data (10 subjects, each given both drugs), the differences are 1.2, 2.4, 1.3, 1.3, 0, 1, 1.8, 0.8, 4.6, 1.4.

Note that 1 of the subjects has a difference of exactly zero. That is common in paired data, and software handles it in more than one way:

The two treatments of the zero difference give the same p-value here, but the point estimate and the interval both differ, and nothing in the output tells you a choice has just been made. When there are zero differences, put the software and its version in the methods.

One more detail that is easy to miss: the achieved coverage of the two rank-based intervals above is 95.6% and 96.5%, not exactly 95%. Rank statistics are discrete, only a finite set of levels is attainable, and R picks the closest and says so.

Three groups or more: Kruskal-Wallis and its post hoc comparisons

The Kruskal-Wallis test is the rank version of ANOVA: rank every observation together, then ask whether the groups have the same mean rank. It is likewise an omnibus test — significance means “not all alike”, not which pair differs.

Splitting birthwt by maternal race gives three groups:

The omnibus result is $\chi^2$ = 8.52 on 2 df, p = 0.014.

There is a common way of running the pairwise comparisons that is wrong: three separate Wilcoxon rank-sum tests. The problem is that each of them re-ranks only the two groups involved, discarding what the third group tells you, and the three reference distributions are mutually inconsistent. The correct approach (Dunn’s test) ranks once, compares the group mean ranks using a variance computed from the whole sample, and shares one tie correction across all comparisons.

This table is the point of the section. The omnibus test (p = 0.014) is significant, yet after the Holm adjustment none of the three pairwise comparisons reaches statistical significance. That is not a contradiction, and it does not mean the groups are alike — it means the evidence in this dataset is not strong enough to name which pair differs while holding the family-wise error rate. The sentence to write is “did not reach statistical significance”, not “there was no difference”.

A few practical points:

• Fix the post hoc method and the adjustment in the protocol. “All pairs” and “each group against the control” need different adjustments, and the arithmetic of accumulating type I error is on t-tests and analysis of variance.
• Holm is almost always preferable to Bonferroni. Both control the family-wise error rate; Holm is the stepwise version, strictly more powerful under the same assumptions.
• Kruskal-Wallis does not compare medians either. Reading it as a shift requires the extra assumption that the group distributions have a similar shape, exactly as in the two-group case.
• Once covariates need adjusting, a rank test is not enough. Rank tests have no adjustment mechanism; move to ordinal logistic regression or to linear regression on a transformed variable.

Run it yourself

library(MASS)
data(birthwt, package = "MASS")

birthwt$smoke_f <- factor(birthwt$smoke, levels = c(0, 1),
                          labels = c("Non-smoker", "Smoker"))
birthwt$race_f  <- factor(birthwt$race, levels = 1:3,
                          labels = c("White", "Black", "Other"))

# Two independent samples: conf.int = TRUE is what yields Hodges-Lehmann
wilcox.test(bwt ~ smoke_f, data = birthwt, conf.int = TRUE)

# P(X > Y): count the pairs directly, ties scoring half
ns <- birthwt$bwt[birthwt$smoke_f == "Non-smoker"]
sm <- birthwt$bwt[birthwt$smoke_f == "Smoker"]
mean(outer(ns, sm, ">")) + 0.5 * mean(outer(ns, sm, "=="))

# The other route: a t-test on logs, giving a ratio of geometric means
fit <- t.test(log(ns), log(sm))
unname(exp(c(diff(rev(fit$estimate)), fit$conf.int)))

# Paired. ⚠️ Since R 4.6 the formula interface rejects paired = TRUE outright,
# so pass two vectors or use the Pair() form.
data(sleep)
d1 <- sleep$extra[sleep$group == 1]
d2 <- sleep$extra[sleep$group == 2]
wilcox.test(d2, d1, paired = TRUE, conf.int = TRUE)
wilcox.test(Pair(d2, d1) ~ 1, conf.int = TRUE)   # equivalent

# Three groups or more
kruskal.test(bwt ~ race_f, data = birthwt)

# Dunn post hoc: rank once, use the whole-sample variance, then adjust with Holm
r <- rank(birthwt$bwt); N <- length(r); tie <- table(birthwt$bwt)
s2 <- N * (N + 1) / 12 - sum(tie^3 - tie) / (12 * (N - 1))
nn <- tapply(r, birthwt$race_f, length); rb <- tapply(r, birthwt$race_f, mean)
cb <- combn(levels(birthwt$race_f), 2)
z  <- apply(cb, 2, function(k) (rb[[k[1]]] - rb[[k[2]]]) /
                               sqrt(s2 * (1 / nn[[k[1]]] + 1 / nn[[k[2]]])))
data.frame(pair = apply(cb, 2, paste, collapse = " vs "),
           z = round(z, 2),
           p_holm = p.adjust(2 * pnorm(-abs(z)), method = "holm"))

import numpy as np
import statsmodels.api as sm
from scipy import stats

bw = sm.datasets.get_rdataset("birthwt", "MASS").data
ns = bw.loc[bw["smoke"] == 0, "bwt"].to_numpy()
smk = bw.loc[bw["smoke"] == 1, "bwt"].to_numpy()

# The test itself: method="exact" / "asymptotic" mirrors R's exact argument
stats.mannwhitneyu(ns, smk, alternative="two-sided")

# ⚠️ scipy has no Hodges-Lehmann estimate and no interval for it.
# The point estimate is easy (it is the median of the pairwise differences);
# an interval has to be bootstrapped, which is NOT the same object as R's
# test-inversion interval, so do not report the two as if they were.
diffs = np.subtract.outer(ns, smk).ravel()
hodges_lehmann = np.median(diffs)

# P(X > Y) likewise has to be counted by hand
p_superior = (np.mean(np.subtract.outer(ns, smk) > 0)
              + 0.5 * np.mean(np.subtract.outer(ns, smk) == 0))

sleep = sm.datasets.get_rdataset("sleep").data
d1 = sleep.loc[sleep["group"] == 1, "extra"].to_numpy()
d2 = sleep.loc[sleep["group"] == 2, "extra"].to_numpy()
stats.wilcoxon(d2, d1)          # zero_method differs from R; set it explicitly

groups = [g["bwt"].to_numpy() for _, g in bw.groupby("race")]
stats.kruskal(*groups)
# ⚠️ Dunn's test is not in scipy (scikit-posthocs has it);
# to stay inside the standard library, code the formula from the R block above.

Common misuses

Regenerating every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B1-04-nonparametric.R

Watch next

▶醫學統計 EP10 t 檢定與非參數法

▶醫學統計 EP11 ANOVA

Report a problem on this page