What a rank test actually tests

hypothesis testing

nonparametric statistics

ecology tutorial

The Wilcoxon rank sum test does not compare medians. It compares whole distributions, and with unequal spread it rejects when the medians agree. Simulated in R.

Author

Tidy Ecology

Published

2026-08-05

Modified

2026-10-14

Corrected 14 October 2026: the post said the Wilcoxon failure needs the groups to differ in spread and is worst when the sample sizes are unequal; it now says a difference in shape does the same, and that unequal sizes made it worst only when the smaller group was the more variable one. The paragraph on ranking before an analysis of variance no longer says the four cells occupy nearly separate stretches of the rank scale: the two cells with one factor each share one stretch, and the sum of their mean ranks is what holds the interaction contrast near zero.

Corrected 1 October 2026: the opening said the null hypothesis of the test is that a random value from one group is equally likely to be above or below one from the other; it now says the null is that both groups come from one distribution, which implies that equality but is not implied by it. The calibration condition is now stated as the same distribution rather than the same spread, and a rank comparison is now called a natural summary for ordinal responses, with ordinal regression as a model-based alternative, rather than the only sensible one.

Corrected 30 September 2026: the post said the Wilcoxon failure needs unequal sample sizes as well as unequal spread; it now says unequal sizes make it worse but the equal-size design was still inflated.

Updated 27 September 2026: a new section, When R says it cannot compute an exact p-value, measures the normal approximation that R before 4.6.0 falls back on with tied cover classes against the exact conditional test, and the limits paragraph no longer says the rank statistic ignores ties: it counts a tie as half.

The most common statistical decision in ecology is probably this one: the histogram looks skewed, or the sample is small, or a reviewer once said something about normality, so the t test is replaced by a Wilcoxon test. The replacement is usually described as testing whether the medians differ.

It does not test that. The Wilcoxon rank sum test, which is the same procedure as the Mann-Whitney U test, has a null hypothesis about the two distributions as wholes: that both groups come from one and the same distribution. That makes a value drawn at random from one group equally likely to be above or below a value drawn at random from the other, but the equality alone is not the null. Medians appear nowhere in it. When the two groups have the same centre but different spread, which happens constantly in field data, the test can reject at several times its nominal rate while the medians are identical by construction.

Two samples with the same centre

A disturbed patch and a reference set. The disturbed patch has the same mean and the same median as the reference, and it is more variable, which is what disturbance usually does. Fifteen quadrats in the patch, forty-five in the reference set.

library(ggplot2)

te_paper  <- "#f5f4ee"
te_ink    <- "#16241d"
te_body   <- "#2c3a31"
te_forest <- "#275139"
te_rust   <- "#b5534e"
te_gold   <- "#c9b458"
te_line   <- "#dad9ca"

theme_datasheet <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.background  = element_rect(fill = te_paper, colour = NA),
          panel.background = element_rect(fill = te_paper, colour = NA),
          panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
          panel.grid.minor = element_blank(),
          text             = element_text(colour = te_body),
          plot.title       = element_text(colour = te_ink, face = "bold"),
          axis.text        = element_text(colour = te_body))
}

draw_pair <- function(n1, sd1, n2, sd2, centre = 20) {
  list(a = rnorm(n1, centre, sd1), b = rnorm(n2, centre, sd2))
}

Both groups are normal with mean 20. Their medians are 20 as well, and the probability that a random patch quadrat exceeds a random reference quadrat is exactly one half, because the difference of two independent normals with equal means is symmetric about zero. Every null hypothesis anyone would want to state about location is true here.

rate <- function(n1, sd1, n2, sd2, reps = 5000, seed = 1) {
  set.seed(seed)
  out <- replicate(reps, {
    d <- draw_pair(n1, sd1, n2, sd2)
    c(wilcox = suppressWarnings(wilcox.test(d$a, d$b)$p.value),
      welch  = t.test(d$a, d$b)$p.value,
      pooled = t.test(d$a, d$b, var.equal = TRUE)$p.value)
  })
  round(100 * rowMeans(out < 0.05), 2)
}

small_group_varies <- rate(15, 4, 45, 1, seed = 101)
large_group_varies <- rate(45, 4, 15, 1, seed = 102)
equal_sizes        <- rate(30, 4, 30, 1, seed = 103)

rbind(`n = 15 varies, n = 45 tight` = small_group_varies,
      `n = 45 varies, n = 15 tight` = large_group_varies,
      `equal n, one varies`         = equal_sizes)

                            wilcox welch pooled
n = 15 varies, n = 45 tight  14.04  5.02  23.50
n = 45 varies, n = 15 tight   1.32  4.30   0.06
equal n, one varies           7.58  4.92   5.30

With the smaller group more variable, the Wilcoxon test rejects 14.0 per cent of the time against a nominal five. The pooled t test is worse at 23.5 per cent, which is the familiar Behrens-Fisher problem and the reason nobody should use it. The Welch t test, which is what t.test gives you by default, sits at 5.0 per cent, right where it should be.

Reverse the design so the larger group is the variable one and the Wilcoxon rate drops to 1.3 per cent: now it is conservative rather than liberal. Even with equal sample sizes it is 7.6 per cent rather than five. The distortion is a function of the design, not a constant you could correct for.

rates <- data.frame(
  design = rep(c("smaller group varies", "equal sizes", "larger group varies"),
               each = 3),
  test   = rep(c("Wilcoxon", "Welch t", "pooled t"), 3),
  pct    = c(small_group_varies[c("wilcox", "welch", "pooled")],
             equal_sizes[c("wilcox", "welch", "pooled")],
             large_group_varies[c("wilcox", "welch", "pooled")]))
rates$design <- factor(rates$design, levels = c("smaller group varies",
                                                "equal sizes",
                                                "larger group varies"))
rates$test <- factor(rates$test, levels = c("Wilcoxon", "Welch t", "pooled t"))

ggplot(rates, aes(x = design, y = pct, fill = test)) +
  geom_hline(yintercept = 5, linetype = "dotted", colour = te_ink) +
  geom_col(position = position_dodge(width = 0.75), width = 0.66) +
  geom_text(aes(label = sprintf("%.1f", pct)),
            position = position_dodge(width = 0.75), vjust = -0.45,
            size = 3.3, colour = te_body) +
  annotate("text", x = 3.45, y = 7.6, label = "nominal 5 per cent",
           hjust = 1, size = 3.3, colour = te_ink) +
  scale_fill_manual(values = c("Wilcoxon" = te_rust, "Welch t" = te_forest,
                               "pooled t" = te_gold)) +
  labs(x = NULL, y = "per cent of tests rejecting at 0.05", fill = NULL,
       title = "All three nulls of location are true in every bar") +
  theme_datasheet() +
  theme(legend.position = "top")

Grouped bar chart with three designs on the horizontal axis and three tests within each. A dotted line marks the nominal five per cent rate. The Wilcoxon bar sits well above the line when the smaller group is the variable one, below it when the larger group is, and slightly above it at equal sample sizes. The Welch bars sit on the line in all three designs and the pooled t bars swing furthest in both directions. — Figure 1: Rejection rates over 5000 simulations for three designs in which both groups have the same mean and median and differ only in spread.

The null that is actually being tested

The rank sum statistic counts, over all pairs of one observation from each group, how many times the first exceeds the second. Its expectation and variance under the null are worked out on the assumption that the two samples come from the same distribution, so that every ordering of the pooled ranks is equally likely. That is a strong null: identical shape, identical spread, identical everything.

If the distributions differ in any way, the permutation variance no longer matches the sampling variance, and the calibration breaks. That is what the table above is showing. The test does not know it is looking at a spread difference; it only knows that the ranks are arranged more unevenly than exchangeability would allow, and it reports that as a significant result which the reader then attributes to location.

The quantity the test is sensitive to has a name and is worth reporting in its own right. Written as the probability that a random observation from group A exceeds a random observation from group B, it is a whole-distribution summary that needs no normality and no equal-variance assumption. It is the same number as the area under a receiver operating characteristic curve, and it comes straight out of the rank sum statistic.

p_superior <- function(a, b) {
  u <- suppressWarnings(wilcox.test(a, b)$statistic)
  as.numeric(u) / (length(a) * length(b))
}

set.seed(9)
same_spread <- replicate(2000, {
  d <- draw_pair(15, 1, 45, 1); p_superior(d$a, d$b) })
diff_spread <- replicate(2000, {
  d <- draw_pair(15, 4, 45, 1); p_superior(d$a, d$b) })

shifted <- replicate(2000, {
  d <- draw_pair(15, 1, 45, 1); p_superior(d$a + 1.2, d$b) })

round(c(mean_equal_spread = mean(same_spread), sd_equal_spread = sd(same_spread),
        mean_wider_group  = mean(diff_spread), sd_wider_group  = sd(diff_spread),
        mean_shifted      = mean(shifted),     sd_shifted      = sd(shifted),
        shifted_theory    = pnorm(1.2 / sqrt(2))), 3)

mean_equal_spread   sd_equal_spread  mean_wider_group    sd_wider_group 
            0.500             0.087             0.498             0.115 
     mean_shifted        sd_shifted    shifted_theory 
            0.801             0.067             0.802

Averaged over 2000 replicate surveys the estimate sits at 0.500 when the two groups share a spread and 0.498 when the smaller group is four times as variable. Both are the true one half, as they must be. What changes is the scatter: the standard deviation of the estimate goes from 0.087 to 0.115, a factor of 1.3. The rank sum test’s p value is computed from the first of those numbers whatever the data look like, and that single substitution is the whole of the inflation in the table above.

Reported as a summary rather than a test, the quantity behaves well. Move one group up by 1.2 units with equal spread and the estimate averages 0.801 against a true value of 0.802, and it reads directly: a quadrat from the shifted group beats a reference quadrat about 80 times in a hundred. Its scatter across surveys is 0.067, which is the number a reader needs in order to know how seriously to take it. That statement survives any monotone transformation of the response and does not pretend to be about medians.

When a rank test is exactly the right tool

None of this makes rank tests bad. It makes them tests of a different thing, and sometimes that thing is the question. Cover classes, Braun-Blanquet scores, behavioural ranks and any ordinal response have no meaningful arithmetic, so a difference in means is not defined; a rank comparison is a natural summary, and an ordinal regression model, as in Ordinal regression for ordered cover classes, is a model-based alternative. Where the question really is “does a random individual from one group tend to be larger”, the rank sum test is the direct answer, and the probability above is the effect size to report with it.

The failure mode is narrower than the reputation suggests. It is the case where the groups differ in spread or in shape, and where the reader is told the result is about location. Both conditions have to be present. In the spread designs above the inflation was worst when the smaller group was also the more variable one, smaller but still present with equal sizes, and reversed into a conservative test when the larger group was the more variable one.

shape_check <- rate(15, 4, 45, 4, reps = 5000, seed = 44)
shape_check

wilcox  welch pooled 
  5.18   5.12   5.18

Give the two groups the same distribution, however non-normal it might be, and the Wilcoxon test is calibrated again: 5.2 per cent here with both groups at the wider standard deviation and the same unequal sample sizes. Equal distributions under the null is the assumption that matters, and normality is not.

Ranking first, then modelling, is a separate trap

The other place ranks turn up is as a pre-treatment: replace the response by its rank and then run an ordinary analysis of variance, on the grounds that this makes the procedure distribution-free. For a one-way comparison that is close to the Kruskal-Wallis test and behaves reasonably. For a factorial design it does not work, because ranking is a nonlinear transformation applied to the pooled response, so what the interaction term sees depends on how far apart the main effects have pushed the cells.

inter_rate <- function(reps = 600, main = 3, inter = 0, seed = 55) {
  set.seed(seed)
  raw <- rankt <- numeric(reps)
  for (i in seq_len(reps)) {
    g <- expand.grid(A = c(0, 1), B = c(0, 1), rep = 1:12)
    g$y <- main * g$A + main * g$B + inter * g$A * g$B + rnorm(nrow(g))
    raw[i]   <- anova(lm(y ~ A * B, data = g))["A:B", "Pr(>F)"]
    g$r      <- rank(g$y)
    rankt[i] <- anova(lm(r ~ A * B, data = g))["A:B", "Pr(>F)"]
  }
  round(100 * c(raw = mean(raw < 0.05), ranked = mean(rankt < 0.05)), 1)
}

rt <- rbind(`no interaction, strong main effects`   = inter_rate(inter = 0),
            `real interaction, strong main effects` = inter_rate(inter = 1.5),
            `real interaction, weak main effects`   = inter_rate(main = 1, inter = 1.5))
rt

                                       raw ranked
no interaction, strong main effects    3.5    0.0
real interaction, strong main effects 70.2    0.0
real interaction, weak main effects   70.2   39.3

The first row is the calibration check, and the rank version passes it in the least reassuring way possible: it rejects 0.0 per cent of the time where the raw analysis rejects 3.5. The second row explains why. With an interaction genuinely present, the raw analysis detects it 70 per cent of the time and the rank analysis detects it 0.0 per cent. The interaction has not been mis-estimated; it has been deleted.

The mechanism is in the third row. When the main effects are large, the cell with neither factor takes nearly all of the lowest ranks and the cell with both nearly all of the highest, while the two cells with one factor each, which have the same expected mean in this design, share the ranks in between. The interaction contrast depends on those two cells only through the sum of their mean ranks, and that sum is nearly the same in every survey, so the contrast is held near zero no matter what the raw means do. Shrink the main effects so all four cells overlap and the rank analysis recovers 39 per cent power against the same interaction. The failure is caused by the terms the analyst was not interested in.

set.seed(56)
g <- expand.grid(A = c(0, 1), B = c(0, 1), rep = 1:200)
g$y <- 3 * g$A + 3 * g$B + 1.5 * g$A * g$B + rnorm(nrow(g))
g$r <- rank(g$y)

panel_means <- function(v, name) {
  m <- tapply(v, list(g$A, g$B), mean)
  contrast <- (m["1", "1"] - m["0", "1"]) - (m["1", "0"] - m["0", "0"])
  list(tab = data.frame(scale = name,
                        A = rep(c(0, 1), 2), B = rep(c(0, 1), each = 2),
                        value = as.numeric(m)),
       share = contrast / diff(range(m)))
}
raw_p <- panel_means(g$y, "raw response")
rnk_p <- panel_means(g$r, "ranked response")

cells <- rbind(raw_p$tab, rnk_p$tab)
cells$scale <- factor(cells$scale, levels = c("raw response", "ranked response"))
cells$B <- factor(cells$B, labels = c("B absent", "B present"))

tags <- data.frame(
  scale = factor(c("raw response", "ranked response"),
                 levels = levels(cells$scale)),
  A = 0.55, B = factor("B absent", levels = levels(cells$B)),
  value = c(max(raw_p$tab$value), max(rnk_p$tab$value)),
  lab = sprintf("lines diverge by %.0f per cent\nof the panel height",
                100 * c(raw_p$share, rnk_p$share)))

ggplot(cells, aes(x = factor(A, labels = c("A absent", "A present")),
                  y = value, colour = B, group = B, shape = B)) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 3) +
  geom_text(data = tags, aes(x = A, y = value, label = lab), inherit.aes = FALSE,
            size = 3.3, colour = te_body, hjust = 0, vjust = 1.05) +
  facet_wrap(~ scale, scales = "free_y") +
  scale_colour_manual(values = c("B absent" = te_gold, "B present" = te_forest)) +
  labs(x = NULL, y = "cell mean", colour = NULL, shape = NULL,
       title = "Ranking flattens a real interaction") +
  theme_datasheet() +
  theme(legend.position = "top",
        strip.text = element_text(colour = te_ink, face = "bold"))

Two panels of interaction plots, each annotated with how far the two lines diverge as a share of the panel height. On the raw scale the lines clearly separate further from left to right. After ranking, the same two lines are close to parallel and the annotated divergence is near zero. — Figure 2: Cell means on the raw scale and after ranking, for a design with strong additive main effects and a real interaction of 1.5.

When R says it cannot compute an exact p-value

Cover classes are one of the responses for which, as above, a rank comparison is a natural summary, and they are full of ties: six classes shared by twelve plots cannot all differ. Versions of R before 4.6.0 print the warning in the title for a default wilcox.test() on such data; from 4.6.0 they do not. Here cover is a latent normal score (centre 2, standard deviation 1.4, from the same draw_pair()) cut into classes 0 to 5, and group A may sit one latent unit higher, as it does in the example. The first chunk runs the default call and records what the running version of R does with it.

class_breaks <- c(0.5, 1.5, 2.5, 3.5, 4.5)
cover_pair <- function(n1, n2, shift = 0) {
  d <- draw_pair(n1, 1.4, n2, 1.4, centre = 2)
  list(a = findInterval(d$a + shift, class_breaks), b = findInterval(d$b, class_breaks))
}
set.seed(12)
plots <- cover_pair(6, 6, shift = 1)
default_warnings <- character(0)
default_fit <- withCallingHandlers(wilcox.test(plots$a, plots$b), warning = function(w) {
  default_warnings <<- c(default_warnings, conditionMessage(w))
  invokeRestart("muffleWarning") })
p_normal <- wilcox.test(plots$a, plots$b, exact = FALSE)$p.value
if (getRversion() < "4.6.0") {
  stopifnot(length(default_warnings) > 0, !grepl("exact", default_fit$method),
            isTRUE(all.equal(default_fit$p.value, p_normal)))
} else stopifnot(length(default_warnings) == 0, grepl("exact", default_fit$method))
stopifnot(isTRUE(all.equal(p_superior(plots$a, plots$b), mean(outer(plots$a, plots$b, ">")) +
                             0.5 * mean(outer(plots$a, plots$b, "==")))))
list(R = R.version.string, cover_A = plots$a, cover_B = plots$b, warnings = default_warnings,
     method = default_fit$method, p_value = round(default_fit$p.value, 4))

$R
[1] "R version 4.5.3 (2026-03-11)"

$cover_A
[1] 1 5 2 2 0 3

$cover_B
[1] 2 1 2 3 1 0

$warnings
[1] "cannot compute exact p-value with ties"

$method
[1] "Wilcoxon rank sum test with continuity correction"

$p_value
[1] 0.5641

The output starts with the version of R that knitted this page, because the answer depends on it. Before R 4.6.0 the default call warns cannot compute exact p-value with ties and falls back on the normal approximation with a tie-corrected variance and a continuity correction; exact = FALSE gives that approximation in every version. From R 4.6.0 the release notes say that wilcox.test() “can now perform exact (conditional) inference in case of ties”: when both groups have fewer than 50 values the same call computes the exact conditional p-value itself, and the warning is gone. The stopifnot() lines test whichever branch the running R is on. The last one checks the rule inside p_superior(): a tied pair counts as half a win.

exact_tied_p <- function(a, b) {  # counts the ways each rank sum can arise
  n1 <- length(a); score <- as.integer(round(2 * rank(c(a, b))))
  top <- sum(sort(score, decreasing = TRUE)[seq_len(n1)])
  ways <- matrix(0, top + 1, n1 + 1); ways[1, 1] <- 1
  for (s in score) {
    rows <- (s + 1):(top + 1)
    ways[rows, -1] <- ways[rows, -1] + ways[rows - s, -(n1 + 1)]
  }
  w <- (0:top) / 2 - n1 * (n1 + 1) / 2; prob <- ways[, n1 + 1] / choose(length(score), n1)
  w_obs <- sum(score[seq_len(n1)]) / 2 - n1 * (n1 + 1) / 2
  min(1, 2 * sum(prob[w <= w_obs + 1e-8]), 2 * sum(prob[w >= w_obs - 1e-8]))
}
r12 <- rank(c(plots$a, plots$b))
w_here <- sum(r12[1:6]) - 21
splits <- combn(12, 6)  # every way to pick six of the twelve plots
stopifnot(choose(12, 6) == ncol(splits))
w_all <- colSums(matrix(r12[splits], 6)) - 21
p_listed <- min(1, 2 * mean(w_all <= w_here + 1e-8), 2 * mean(w_all >= w_here - 1e-8))
stopifnot(isTRUE(all.equal(p_listed, exact_tied_p(plots$a, plots$b))))
if (getRversion() >= "4.6.0") stopifnot(isTRUE(all.equal(default_fit$p.value, p_listed)))
p_three <- c(normal = wilcox.test(c(2, 2, 2), c(0, 0, 0), exact = FALSE)$p.value,
             exact = exact_tied_p(c(2, 2, 2), c(0, 0, 0)))  # three plots per group
stopifnot(isTRUE(all.equal(suppressWarnings(wilcox.test(c(2, 2, 2), c(0, 0, 0)))$p.value,
                           p_three[[if (getRversion() < "4.6.0") "normal" else "exact"]])))
tie_rates <- function(n1, n2, shift, reps) {
  hit <- replicate(reps, {
    d <- cover_pair(n1, n2, shift)
    c(normal = wilcox.test(d$a, d$b, exact = FALSE)$p.value, exact = exact_tied_p(d$a, d$b))
  }) < 0.05
  c(normal = mean(hit["normal", ]), exact = mean(hit["exact", ]),
    disagree = mean(hit["normal", ] != hit["exact", ]))
}
set.seed(35)
tie_table <- rbind(`6 + 6, no shift`   = tie_rates(6, 6, 0, 4000),
                   `6 + 6, shift 1`    = tie_rates(6, 6, 1, 4000),
                   `15 + 45, no shift` = tie_rates(15, 45, 0, 1000),
                   `15 + 45, shift 1`  = tie_rates(15, 45, 1, 1000))
round(c(p_exact = p_listed, p_normal = p_normal, distinct_W = length(unique(w_all))), 4)

   p_exact   p_normal distinct_W 
    0.5952     0.5641    47.0000

round(tie_table[1:2, ], 4)

                normal  exact disagree
6 + 6, no shift 0.0415 0.0357   0.0122
6 + 6, shift 1  0.1585 0.1485   0.0255

The exact conditional test keeps the tied midranks, finds the share of all relabellings of the plots that give a rank sum at least as far out on the observed side, and doubles it, which is R’s two-sided rule: a permutation test of the midranks, as in Permutation tests from scratch. At six plots per group the 924 possible splits are few enough to list, and the listing checks exact_tied_p(), which counts the ways each rank sum can arise instead of listing them and so also handles 15 + 45. For the example the exact p-value is 0.595 against 0.564 from the normal approximation; on R 4.6.0 and later the chunk also checks that R’s own p-value equals the listed one. Over 4000 pairs of six plots with no difference between the groups, the normal approximation rejected at 0.05 in 4.2 per cent and the exact test in 3.6 per cent, and the two p-values fell on opposite sides of 0.05 in 1.2 per cent of the surveys (Monte Carlo standard error of a rate near five per cent: 0.34 points). With group A one latent unit higher they rejected in 15.8 and 14.8 per cent. At the 15 + 45 design of the first section they disagreed in 1 of 1000 null surveys and 4 of 1000 shifted ones.

So the approximation that the old warning announces was not the problem at these sizes; the discreteness of the conditional null is. The example’s 924 splits give 47 distinct rank sums, so the attainable p-values come in steps, and the exact test rejected in 3.6 per cent of null surveys at 6 + 6, below the nominal five. At three plots per group there are only 20 splits of six plots, so no exact two-sided p-value can fall below 0.10 whatever the data. The normal approximation can: classes 2, 2, 2 against 0, 0, 0 give 0.047 from it and 0.10 from the exact test, so on those data the default call rejects at 0.05 before R 4.6.0 and does not from 4.6.0. The rates here are for six or more plots per group. Likert items: when group means lie runs the same test on five-box answers at 100 per group, where R uses the normal approximation in every version.

What to do instead

For two independent samples with possibly unequal spread, the Welch t test is the default in R for a reason, and the simulations above show it holding its rate in every design. If the response is heavily skewed or bounded, the better move is usually a model that respects the response: a generalised linear model for counts or proportions, a bootstrap interval for an awkward statistic, or a permutation test that permutes within the structure of the design.

If a rank comparison is what the question calls for, run it and report the probability of superiority alongside the p value, so the reader knows which quantity was compared.

Honest limits

The simulations here use normal samples throughout, deliberately: the point is that the calibration problem comes from unequal distributions and not from non-normality, and normal samples make that impossible to confuse. With skewed distributions the same mechanism operates, and the direction of the error still depends on which group is larger.

The probability of superiority computed from the rank sum statistic counts a tied pair as half a win, so with coarse ordinal data or counts with many repeated values it estimates the probability that A exceeds B plus half the probability of a tie, not the first alone. What ties do to the p value is measured in the section When R says it cannot compute an exact p-value above.

The interaction result belongs to one balanced two-by-two design at one noise level, and the rank transform literature reports failures in both directions depending on the configuration: inflated rejection rates are also documented, in designs where the cells overlap more than they do here. The transferable claim is that the rank transform does not preserve the meaning of an interaction term, not the particular percentages above. A method built for the purpose, such as an aligned rank transform, adjusts the response separately for each term before ranking, precisely to break the dependence on the main effects.

Finally, none of the tests here addresses the assumption that usually matters more than any of this: independence. A Wilcoxon test on pseudoreplicated quadrats is wrong for a reason no amount of distribution-free machinery can fix.

References

Wilcoxon F 1945 Biometrics Bulletin 1(6):80-83 (10.2307/3001968)

Mann HB, Whitney DR 1947 Annals of Mathematical Statistics 18(1):50-60 (10.1214/aoms/1177730491)

Fagerland MW, Sandvik L 2009 Statistics in Medicine 28(10):1487-1497 (10.1002/sim.3561)

Divine GW, Norton HJ, Baron AE, Juarez-Colunga E 2018 The American Statistician 72(3):278-286 (10.1080/00031305.2017.1305291)

--- title: "What a rank test actually tests" description: "The Wilcoxon rank sum test does not compare medians. It compares whole distributions, and with unequal spread it rejects when the medians agree. Simulated in R." date: "2026-08-05 15:00" date-modified: "2026-10-14" categories: [R, hypothesis testing, nonparametric statistics, ecology tutorial] image: thumbnail.png image-alt: "Two overlapping density curves with the same mean and the same median but very different spread, with a dashed line marking their shared centre." --- *Corrected 14 October 2026: the post said the Wilcoxon failure needs the groups to differ in spread and is worst when the sample sizes are unequal; it now says a difference in shape does the same, and that unequal sizes made it worst only when the smaller group was the more variable one. The paragraph on ranking before an analysis of variance no longer says the four cells occupy nearly separate stretches of the rank scale: the two cells with one factor each share one stretch, and the sum of their mean ranks is what holds the interaction contrast near zero.* *Corrected 1 October 2026: the opening said the null hypothesis of the test is that a random value from one group is equally likely to be above or below one from the other; it now says the null is that both groups come from one distribution, which implies that equality but is not implied by it. The calibration condition is now stated as the same distribution rather than the same spread, and a rank comparison is now called a natural summary for ordinal responses, with ordinal regression as a model-based alternative, rather than the only sensible one.* *Corrected 30 September 2026: the post said the Wilcoxon failure needs unequal sample sizes as well as unequal spread; it now says unequal sizes make it worse but the equal-size design was still inflated.* *Updated 27 September 2026: a new section, When R says it cannot compute an exact p-value, measures the normal approximation that R before 4.6.0 falls back on with tied cover classes against the exact conditional test, and the limits paragraph no longer says the rank statistic ignores ties: it counts a tie as half.* The most common statistical decision in ecology is probably this one: the histogram looks skewed, or the sample is small, or a reviewer once said something about normality, so the t test is replaced by a Wilcoxon test. The replacement is usually described as testing whether the medians differ. It does not test that. The Wilcoxon rank sum test, which is the same procedure as the Mann-Whitney U test, has a null hypothesis about the two distributions as wholes: that both groups come from one and the same distribution. That makes a value drawn at random from one group equally likely to be above or below a value drawn at random from the other, but the equality alone is not the null. Medians appear nowhere in it. When the two groups have the same centre but different spread, which happens constantly in field data, the test can reject at several times its nominal rate while the medians are identical by construction. ## Two samples with the same centre A disturbed patch and a reference set. The disturbed patch has the same mean and the same median as the reference, and it is more variable, which is what disturbance usually does. Fifteen quadrats in the patch, forty-five in the reference set. ```{r setup} #| message: false #| warning: false library(ggplot2) te_paper <- "#f5f4ee" te_ink <- "#16241d" te_body <- "#2c3a31" te_forest <- "#275139" te_rust <- "#b5534e" te_gold <- "#c9b458" te_line <- "#dad9ca" theme_datasheet <- function() { theme_minimal(base_size = 12) + theme(plot.background = element_rect(fill = te_paper, colour = NA), panel.background = element_rect(fill = te_paper, colour = NA), panel.grid.major = element_line(colour = te_line, linewidth = 0.3), panel.grid.minor = element_blank(), text = element_text(colour = te_body), plot.title = element_text(colour = te_ink, face = "bold"), axis.text = element_text(colour = te_body)) } draw_pair <- function(n1, sd1, n2, sd2, centre = 20) { list(a = rnorm(n1, centre, sd1), b = rnorm(n2, centre, sd2)) } ``` Both groups are normal with mean 20. Their medians are 20 as well, and the probability that a random patch quadrat exceeds a random reference quadrat is exactly one half, because the difference of two independent normals with equal means is symmetric about zero. Every null hypothesis anyone would want to state about location is true here. ```{r rates} rate <- function(n1, sd1, n2, sd2, reps = 5000, seed = 1) { set.seed(seed) out <- replicate(reps, { d <- draw_pair(n1, sd1, n2, sd2) c(wilcox = suppressWarnings(wilcox.test(d$a, d$b)$p.value), welch = t.test(d$a, d$b)$p.value, pooled = t.test(d$a, d$b, var.equal = TRUE)$p.value) }) round(100 * rowMeans(out < 0.05), 2) } small_group_varies <- rate(15, 4, 45, 1, seed = 101) large_group_varies <- rate(45, 4, 15, 1, seed = 102) equal_sizes <- rate(30, 4, 30, 1, seed = 103) rbind(`n = 15 varies, n = 45 tight` = small_group_varies, `n = 45 varies, n = 15 tight` = large_group_varies, `equal n, one varies` = equal_sizes) ``` With the smaller group more variable, the Wilcoxon test rejects `r sprintf("%.1f", small_group_varies[["wilcox"]])` per cent of the time against a nominal five. The pooled t test is worse at `r sprintf("%.1f", small_group_varies[["pooled"]])` per cent, which is the familiar Behrens-Fisher problem and the reason nobody should use it. The Welch t test, which is what `t.test` gives you by default, sits at `r sprintf("%.1f", small_group_varies[["welch"]])` per cent, right where it should be. Reverse the design so the larger group is the variable one and the Wilcoxon rate drops to `r sprintf("%.1f", large_group_varies[["wilcox"]])` per cent: now it is conservative rather than liberal. Even with equal sample sizes it is `r sprintf("%.1f", equal_sizes[["wilcox"]])` per cent rather than five. The distortion is a function of the design, not a constant you could correct for. ```{r fig-rates} #| fig-cap: "Rejection rates over 5000 simulations for three designs in which both groups have the same mean and median and differ only in spread." #| fig-alt: "Grouped bar chart with three designs on the horizontal axis and three tests within each. A dotted line marks the nominal five per cent rate. The Wilcoxon bar sits well above the line when the smaller group is the variable one, below it when the larger group is, and slightly above it at equal sample sizes. The Welch bars sit on the line in all three designs and the pooled t bars swing furthest in both directions." rates <- data.frame( design = rep(c("smaller group varies", "equal sizes", "larger group varies"), each = 3), test = rep(c("Wilcoxon", "Welch t", "pooled t"), 3), pct = c(small_group_varies[c("wilcox", "welch", "pooled")], equal_sizes[c("wilcox", "welch", "pooled")], large_group_varies[c("wilcox", "welch", "pooled")])) rates$design <- factor(rates$design, levels = c("smaller group varies", "equal sizes", "larger group varies")) rates$test <- factor(rates$test, levels = c("Wilcoxon", "Welch t", "pooled t")) ggplot(rates, aes(x = design, y = pct, fill = test)) + geom_hline(yintercept = 5, linetype = "dotted", colour = te_ink) + geom_col(position = position_dodge(width = 0.75), width = 0.66) + geom_text(aes(label = sprintf("%.1f", pct)), position = position_dodge(width = 0.75), vjust = -0.45, size = 3.3, colour = te_body) + annotate("text", x = 3.45, y = 7.6, label = "nominal 5 per cent", hjust = 1, size = 3.3, colour = te_ink) + scale_fill_manual(values = c("Wilcoxon" = te_rust, "Welch t" = te_forest, "pooled t" = te_gold)) + labs(x = NULL, y = "per cent of tests rejecting at 0.05", fill = NULL, title = "All three nulls of location are true in every bar") + theme_datasheet() + theme(legend.position = "top") ``` ## The null that is actually being tested The rank sum statistic counts, over all pairs of one observation from each group, how many times the first exceeds the second. Its expectation and variance under the null are worked out on the assumption that the two samples come from the *same* distribution, so that every ordering of the pooled ranks is equally likely. That is a strong null: identical shape, identical spread, identical everything. If the distributions differ in any way, the permutation variance no longer matches the sampling variance, and the calibration breaks. That is what the table above is showing. The test does not know it is looking at a spread difference; it only knows that the ranks are arranged more unevenly than exchangeability would allow, and it reports that as a significant result which the reader then attributes to location. The quantity the test is sensitive to has a name and is worth reporting in its own right. Written as the probability that a random observation from group A exceeds a random observation from group B, it is a whole-distribution summary that needs no normality and no equal-variance assumption. It is the same number as the area under a receiver operating characteristic curve, and it comes straight out of the rank sum statistic. ```{r estimand} p_superior <- function(a, b) { u <- suppressWarnings(wilcox.test(a, b)$statistic) as.numeric(u) / (length(a) * length(b)) } set.seed(9) same_spread <- replicate(2000, { d <- draw_pair(15, 1, 45, 1); p_superior(d$a, d$b) }) diff_spread <- replicate(2000, { d <- draw_pair(15, 4, 45, 1); p_superior(d$a, d$b) }) shifted <- replicate(2000, { d <- draw_pair(15, 1, 45, 1); p_superior(d$a + 1.2, d$b) }) round(c(mean_equal_spread = mean(same_spread), sd_equal_spread = sd(same_spread), mean_wider_group = mean(diff_spread), sd_wider_group = sd(diff_spread), mean_shifted = mean(shifted), sd_shifted = sd(shifted), shifted_theory = pnorm(1.2 / sqrt(2))), 3) ``` Averaged over 2000 replicate surveys the estimate sits at `r sprintf("%.3f", mean(same_spread))` when the two groups share a spread and `r sprintf("%.3f", mean(diff_spread))` when the smaller group is four times as variable. Both are the true one half, as they must be. What changes is the scatter: the standard deviation of the estimate goes from `r sprintf("%.3f", sd(same_spread))` to `r sprintf("%.3f", sd(diff_spread))`, a factor of `r sprintf("%.1f", sd(diff_spread) / sd(same_spread))`. The rank sum test's p value is computed from the first of those numbers whatever the data look like, and that single substitution is the whole of the inflation in the table above. Reported as a summary rather than a test, the quantity behaves well. Move one group up by 1.2 units with equal spread and the estimate averages `r sprintf("%.3f", mean(shifted))` against a true value of `r sprintf("%.3f", pnorm(1.2 / sqrt(2)))`, and it reads directly: a quadrat from the shifted group beats a reference quadrat about `r sprintf("%.0f", 100 * mean(shifted))` times in a hundred. Its scatter across surveys is `r sprintf("%.3f", sd(shifted))`, which is the number a reader needs in order to know how seriously to take it. That statement survives any monotone transformation of the response and does not pretend to be about medians. ## When a rank test is exactly the right tool None of this makes rank tests bad. It makes them tests of a different thing, and sometimes that thing is the question. Cover classes, Braun-Blanquet scores, behavioural ranks and any ordinal response have no meaningful arithmetic, so a difference in means is not defined; a rank comparison is a natural summary, and an ordinal regression model, as in [Ordinal regression for ordered cover classes](../ordinal-regression-cover-classes/), is a model-based alternative. Where the question really is "does a random individual from one group tend to be larger", the rank sum test is the direct answer, and the probability above is the effect size to report with it. The failure mode is narrower than the reputation suggests. It is the case where the groups differ in spread or in shape, and where the reader is told the result is about location. Both conditions have to be present. In the spread designs above the inflation was worst when the smaller group was also the more variable one, smaller but still present with equal sizes, and reversed into a conservative test when the larger group was the more variable one. ```{r shapes} shape_check <- rate(15, 4, 45, 4, reps = 5000, seed = 44) shape_check ``` Give the two groups the same distribution, however non-normal it might be, and the Wilcoxon test is calibrated again: `r sprintf("%.1f", shape_check[["wilcox"]])` per cent here with both groups at the wider standard deviation and the same unequal sample sizes. Equal distributions under the null is the assumption that matters, and normality is not. ## Ranking first, then modelling, is a separate trap The other place ranks turn up is as a pre-treatment: replace the response by its rank and then run an ordinary analysis of variance, on the grounds that this makes the procedure distribution-free. For a one-way comparison that is close to the Kruskal-Wallis test and behaves reasonably. For a factorial design it does not work, because ranking is a nonlinear transformation applied to the pooled response, so what the interaction term sees depends on how far apart the main effects have pushed the cells. ```{r rank-anova} inter_rate <- function(reps = 600, main = 3, inter = 0, seed = 55) { set.seed(seed) raw <- rankt <- numeric(reps) for (i in seq_len(reps)) { g <- expand.grid(A = c(0, 1), B = c(0, 1), rep = 1:12) g$y <- main * g$A + main * g$B + inter * g$A * g$B + rnorm(nrow(g)) raw[i] <- anova(lm(y ~ A * B, data = g))["A:B", "Pr(>F)"] g$r <- rank(g$y) rankt[i] <- anova(lm(r ~ A * B, data = g))["A:B", "Pr(>F)"] } round(100 * c(raw = mean(raw < 0.05), ranked = mean(rankt < 0.05)), 1) } rt <- rbind(`no interaction, strong main effects` = inter_rate(inter = 0), `real interaction, strong main effects` = inter_rate(inter = 1.5), `real interaction, weak main effects` = inter_rate(main = 1, inter = 1.5)) rt ``` The first row is the calibration check, and the rank version passes it in the least reassuring way possible: it rejects `r sprintf("%.1f", rt[1, "ranked"])` per cent of the time where the raw analysis rejects `r sprintf("%.1f", rt[1, "raw"])`. The second row explains why. With an interaction genuinely present, the raw analysis detects it `r sprintf("%.0f", rt[2, "raw"])` per cent of the time and the rank analysis detects it `r sprintf("%.1f", rt[2, "ranked"])` per cent. The interaction has not been mis-estimated; it has been deleted. The mechanism is in the third row. When the main effects are large, the cell with neither factor takes nearly all of the lowest ranks and the cell with both nearly all of the highest, while the two cells with one factor each, which have the same expected mean in this design, share the ranks in between. The interaction contrast depends on those two cells only through the sum of their mean ranks, and that sum is nearly the same in every survey, so the contrast is held near zero no matter what the raw means do. Shrink the main effects so all four cells overlap and the rank analysis recovers `r sprintf("%.0f", rt[3, "ranked"])` per cent power against the same interaction. The failure is caused by the terms the analyst was not interested in. ```{r fig-ranks} #| fig-cap: "Cell means on the raw scale and after ranking, for a design with strong additive main effects and a real interaction of 1.5." #| fig-alt: "Two panels of interaction plots, each annotated with how far the two lines diverge as a share of the panel height. On the raw scale the lines clearly separate further from left to right. After ranking, the same two lines are close to parallel and the annotated divergence is near zero." set.seed(56) g <- expand.grid(A = c(0, 1), B = c(0, 1), rep = 1:200) g$y <- 3 * g$A + 3 * g$B + 1.5 * g$A * g$B + rnorm(nrow(g)) g$r <- rank(g$y) panel_means <- function(v, name) { m <- tapply(v, list(g$A, g$B), mean) contrast <- (m["1", "1"] - m["0", "1"]) - (m["1", "0"] - m["0", "0"]) list(tab = data.frame(scale = name, A = rep(c(0, 1), 2), B = rep(c(0, 1), each = 2), value = as.numeric(m)), share = contrast / diff(range(m))) } raw_p <- panel_means(g$y, "raw response") rnk_p <- panel_means(g$r, "ranked response") cells <- rbind(raw_p$tab, rnk_p$tab) cells$scale <- factor(cells$scale, levels = c("raw response", "ranked response")) cells$B <- factor(cells$B, labels = c("B absent", "B present")) tags <- data.frame( scale = factor(c("raw response", "ranked response"), levels = levels(cells$scale)), A = 0.55, B = factor("B absent", levels = levels(cells$B)), value = c(max(raw_p$tab$value), max(rnk_p$tab$value)), lab = sprintf("lines diverge by %.0f per cent\nof the panel height", 100 * c(raw_p$share, rnk_p$share))) ggplot(cells, aes(x = factor(A, labels = c("A absent", "A present")), y = value, colour = B, group = B, shape = B)) + geom_line(linewidth = 0.9) + geom_point(size = 3) + geom_text(data = tags, aes(x = A, y = value, label = lab), inherit.aes = FALSE, size = 3.3, colour = te_body, hjust = 0, vjust = 1.05) + facet_wrap(~ scale, scales = "free_y") + scale_colour_manual(values = c("B absent" = te_gold, "B present" = te_forest)) + labs(x = NULL, y = "cell mean", colour = NULL, shape = NULL, title = "Ranking flattens a real interaction") + theme_datasheet() + theme(legend.position = "top", strip.text = element_text(colour = te_ink, face = "bold")) ``` ## When R says it cannot compute an exact p-value Cover classes are one of the responses for which, as above, a rank comparison is a natural summary, and they are full of ties: six classes shared by twelve plots cannot all differ. Versions of R before 4.6.0 print the warning in the title for a default `wilcox.test()` on such data; from 4.6.0 they do not. Here cover is a latent normal score (centre 2, standard deviation 1.4, from the same `draw_pair()`) cut into classes 0 to 5, and group A may sit one latent unit higher, as it does in the example. The first chunk runs the default call and records what the running version of R does with it. ```{r ties-default} class_breaks <- c(0.5, 1.5, 2.5, 3.5, 4.5) cover_pair <- function(n1, n2, shift = 0) { d <- draw_pair(n1, 1.4, n2, 1.4, centre = 2) list(a = findInterval(d$a + shift, class_breaks), b = findInterval(d$b, class_breaks)) } set.seed(12) plots <- cover_pair(6, 6, shift = 1) default_warnings <- character(0) default_fit <- withCallingHandlers(wilcox.test(plots$a, plots$b), warning = function(w) { default_warnings <<- c(default_warnings, conditionMessage(w)) invokeRestart("muffleWarning") }) p_normal <- wilcox.test(plots$a, plots$b, exact = FALSE)$p.value if (getRversion() < "4.6.0") { stopifnot(length(default_warnings) > 0, !grepl("exact", default_fit$method), isTRUE(all.equal(default_fit$p.value, p_normal))) } else stopifnot(length(default_warnings) == 0, grepl("exact", default_fit$method)) stopifnot(isTRUE(all.equal(p_superior(plots$a, plots$b), mean(outer(plots$a, plots$b, ">")) + 0.5 * mean(outer(plots$a, plots$b, "=="))))) list(R = R.version.string, cover_A = plots$a, cover_B = plots$b, warnings = default_warnings, method = default_fit$method, p_value = round(default_fit$p.value, 4)) ``` The output starts with the version of R that knitted this page, because the answer depends on it. Before R 4.6.0 the default call warns `cannot compute exact p-value with ties` and falls back on the normal approximation with a tie-corrected variance and a continuity correction; `exact = FALSE` gives that approximation in every version. From R 4.6.0 the release notes say that `wilcox.test()` "can now perform exact (conditional) inference in case of ties": when both groups have fewer than 50 values the same call computes the exact conditional p-value itself, and the warning is gone. The `stopifnot()` lines test whichever branch the running R is on. The last one checks the rule inside `p_superior()`: a tied pair counts as half a win. ```{r ties-exact} exact_tied_p <- function(a, b) { # counts the ways each rank sum can arise n1 <- length(a); score <- as.integer(round(2 * rank(c(a, b)))) top <- sum(sort(score, decreasing = TRUE)[seq_len(n1)]) ways <- matrix(0, top + 1, n1 + 1); ways[1, 1] <- 1 for (s in score) { rows <- (s + 1):(top + 1) ways[rows, -1] <- ways[rows, -1] + ways[rows - s, -(n1 + 1)] } w <- (0:top) / 2 - n1 * (n1 + 1) / 2; prob <- ways[, n1 + 1] / choose(length(score), n1) w_obs <- sum(score[seq_len(n1)]) / 2 - n1 * (n1 + 1) / 2 min(1, 2 * sum(prob[w <= w_obs + 1e-8]), 2 * sum(prob[w >= w_obs - 1e-8])) } r12 <- rank(c(plots$a, plots$b)) w_here <- sum(r12[1:6]) - 21 splits <- combn(12, 6) # every way to pick six of the twelve plots stopifnot(choose(12, 6) == ncol(splits)) w_all <- colSums(matrix(r12[splits], 6)) - 21 p_listed <- min(1, 2 * mean(w_all <= w_here + 1e-8), 2 * mean(w_all >= w_here - 1e-8)) stopifnot(isTRUE(all.equal(p_listed, exact_tied_p(plots$a, plots$b)))) if (getRversion() >= "4.6.0") stopifnot(isTRUE(all.equal(default_fit$p.value, p_listed))) p_three <- c(normal = wilcox.test(c(2, 2, 2), c(0, 0, 0), exact = FALSE)$p.value, exact = exact_tied_p(c(2, 2, 2), c(0, 0, 0))) # three plots per group stopifnot(isTRUE(all.equal(suppressWarnings(wilcox.test(c(2, 2, 2), c(0, 0, 0)))$p.value, p_three[[if (getRversion() < "4.6.0") "normal" else "exact"]]))) tie_rates <- function(n1, n2, shift, reps) { hit <- replicate(reps, { d <- cover_pair(n1, n2, shift) c(normal = wilcox.test(d$a, d$b, exact = FALSE)$p.value, exact = exact_tied_p(d$a, d$b)) }) < 0.05 c(normal = mean(hit["normal", ]), exact = mean(hit["exact", ]), disagree = mean(hit["normal", ] != hit["exact", ])) } set.seed(35) tie_table <- rbind(`6 + 6, no shift` = tie_rates(6, 6, 0, 4000), `6 + 6, shift 1` = tie_rates(6, 6, 1, 4000), `15 + 45, no shift` = tie_rates(15, 45, 0, 1000), `15 + 45, shift 1` = tie_rates(15, 45, 1, 1000)) round(c(p_exact = p_listed, p_normal = p_normal, distinct_W = length(unique(w_all))), 4) round(tie_table[1:2, ], 4) ``` The exact conditional test keeps the tied midranks, finds the share of all relabellings of the plots that give a rank sum at least as far out on the observed side, and doubles it, which is R's two-sided rule: a permutation test of the midranks, as in [Permutation tests from scratch](../permutation-tests-from-scratch/). At six plots per group the `r sprintf("%d", ncol(splits))` possible splits are few enough to list, and the listing checks `exact_tied_p()`, which counts the ways each rank sum can arise instead of listing them and so also handles 15 + 45. For the example the exact p-value is `r sprintf("%.3f", p_listed)` against `r sprintf("%.3f", p_normal)` from the normal approximation; on R 4.6.0 and later the chunk also checks that R's own p-value equals the listed one. Over 4000 pairs of six plots with no difference between the groups, the normal approximation rejected at 0.05 in `r sprintf("%.1f", 100 * tie_table["6 + 6, no shift", "normal"])` per cent and the exact test in `r sprintf("%.1f", 100 * tie_table["6 + 6, no shift", "exact"])` per cent, and the two p-values fell on opposite sides of 0.05 in `r sprintf("%.1f", 100 * tie_table["6 + 6, no shift", "disagree"])` per cent of the surveys (Monte Carlo standard error of a rate near five per cent: `r sprintf("%.2f", 100 * sqrt(0.05 * 0.95 / 4000))` points). With group A one latent unit higher they rejected in `r sprintf("%.1f", 100 * tie_table["6 + 6, shift 1", "normal"])` and `r sprintf("%.1f", 100 * tie_table["6 + 6, shift 1", "exact"])` per cent. At the 15 + 45 design of the first section they disagreed in `r sprintf("%d", as.integer(round(1000 * tie_table["15 + 45, no shift", "disagree"])))` of 1000 null surveys and `r sprintf("%d", as.integer(round(1000 * tie_table["15 + 45, shift 1", "disagree"])))` of 1000 shifted ones. So the approximation that the old warning announces was not the problem at these sizes; the discreteness of the conditional null is. The example's `r sprintf("%d", ncol(splits))` splits give `r sprintf("%d", length(unique(w_all)))` distinct rank sums, so the attainable p-values come in steps, and the exact test rejected in `r sprintf("%.1f", 100 * tie_table["6 + 6, no shift", "exact"])` per cent of null surveys at 6 + 6, below the nominal five. At three plots per group there are only `r sprintf("%d", as.integer(choose(6, 3)))` splits of six plots, so no exact two-sided p-value can fall below `r sprintf("%.2f", 2 / choose(6, 3))` whatever the data. The normal approximation can: classes 2, 2, 2 against 0, 0, 0 give `r sprintf("%.3f", p_three[["normal"]])` from it and `r sprintf("%.2f", p_three[["exact"]])` from the exact test, so on those data the default call rejects at 0.05 before R 4.6.0 and does not from 4.6.0. The rates here are for six or more plots per group. [Likert items: when group means lie](../likert-means-and-latent-variance/) runs the same test on five-box answers at 100 per group, where R uses the normal approximation in every version. ## What to do instead For two independent samples with possibly unequal spread, the Welch t test is the default in R for a reason, and the simulations above show it holding its rate in every design. If the response is heavily skewed or bounded, the better move is usually a model that respects the response: a generalised linear model for counts or proportions, a bootstrap interval for an awkward statistic, or a permutation test that permutes within the structure of the design. If a rank comparison is what the question calls for, run it and report the probability of superiority alongside the p value, so the reader knows which quantity was compared. ## Honest limits The simulations here use normal samples throughout, deliberately: the point is that the calibration problem comes from unequal distributions and not from non-normality, and normal samples make that impossible to confuse. With skewed distributions the same mechanism operates, and the direction of the error still depends on which group is larger. The probability of superiority computed from the rank sum statistic counts a tied pair as half a win, so with coarse ordinal data or counts with many repeated values it estimates the probability that A exceeds B plus half the probability of a tie, not the first alone. What ties do to the p value is measured in the section [When R says it cannot compute an exact p-value](#when-r-says-it-cannot-compute-an-exact-p-value) above. The interaction result belongs to one balanced two-by-two design at one noise level, and the rank transform literature reports failures in both directions depending on the configuration: inflated rejection rates are also documented, in designs where the cells overlap more than they do here. The transferable claim is that the rank transform does not preserve the meaning of an interaction term, not the particular percentages above. A method built for the purpose, such as an aligned rank transform, adjusts the response separately for each term before ranking, precisely to break the dependence on the main effects. Finally, none of the tests here addresses the assumption that usually matters more than any of this: independence. A Wilcoxon test on pseudoreplicated quadrats is wrong for a reason no amount of distribution-free machinery can fix. ## References Wilcoxon F 1945 Biometrics Bulletin 1(6):80-83 (10.2307/3001968) Mann HB, Whitney DR 1947 Annals of Mathematical Statistics 18(1):50-60 (10.1214/aoms/1177730491) Fagerland MW, Sandvik L 2009 Statistics in Medicine 28(10):1487-1497 (10.1002/sim.3561) Divine GW, Norton HJ, Baron AE, Juarez-Colunga E 2018 The American Statistician 72(3):278-286 (10.1080/00031305.2017.1305291) ## Related tutorials - [Classical tests as linear models](../linear-models-t-test-anova/) - [Pretesting variances before a t test](../pretesting-assumptions-distorts-the-test/) - [Likert items: when group means lie](../likert-means-and-latent-variance/) - [Permutation tests from scratch](../permutation-tests-from-scratch/)