Kruskal-Wallis test in R, with Dunn's follow-up

To run a Kruskal-Wallis test in R, call kruskal.test(outcome ~ group, data = ...) from the base stats package. It tells you whether at least one group tends to have higher or lower values than the others, but not which groups differ or by how much. A complete analysis therefore adds two more calls: an effect size with its interval, and Dunn’s pairwise test with a named p-value adjustment.

The worked example below runs all three on one public dataset in one script, reads every line of the output, and ends with the sentence you would put in a paper. The numbers come from a deterministic Celeus engine run of the same analysis, and the R script it generated is below the article.

The data and the question

The example uses InsectSprays, a data frame that ships with R. It records the number of insects counted on 72 field plots, 12 plots for each of six insecticides labelled A to F. Each plot received exactly one spray, so the six groups are independent, and the outcome is a count. The data come from Beall (1942), a paper about transforming insect counts so that a classic analysis of variance could be applied to them.

The question is: do insect counts per plot differ between the six sprays, and if they do, which pairs of sprays differ?

The design matters more than the function. The Kruskal-Wallis test needs groups that do not share units, an outcome that can be put in order, and a question about whether values tend to be larger in some groups than in others. A count of insects meets all three. Plots measured repeatedly under different sprays would not, and the FAQ below covers that case.

The whole script

The whole analysis is three calls, shown here as a hand-written script. kruskal.test() is in base R; the effect size needs the effectsize package and Dunn’s test needs PMCMRplus, both on CRAN.

library(effectsize)
library(PMCMRplus)

data <- InsectSprays
data$spray <- factor(data$spray)   # already a factor here; needed if read from a CSV

# 1. The omnibus test
kw <- kruskal.test(count ~ spray, data = data)
kw

# 2. Effect size: epsilon squared for ranks, two-sided percentile-bootstrap interval
set.seed(20261101)
rank_epsilon_squared(data$count, data$spray, ci = 0.95,
                     alternative = "two.sided", iterations = 1000)

# 3. Follow-up: Dunn's test on the joint ranking, Holm-adjusted p-values
dunn <- kwAllPairsDunnTest(count ~ spray, data = data, p.adjust.method = "holm")
dunn
summary(dunn)   # one row per pair, with z and the adjusted p-value

The engine’s own script, under the article, does the same three steps on a copy of the data read from a CSV file, which is why it converts the group to a factor and drops incomplete rows explicitly. On the same package versions it gives the same numbers.

Three ways to call kruskal.test()

kruskal.test() accepts the data in three shapes, and all three give the same answer. In a check on this dataset under R 4.6.1, the formula form, the vector form and the list form returned the identical statistic.

kruskal.test(count ~ spray, data = InsectSprays)            # formula
kruskal.test(InsectSprays$count, InsectSprays$spray)        # x and g
kruskal.test(split(InsectSprays$count, InsectSprays$spray)) # a list of groups

The formula form is the one to use. It takes a subset argument for running the test on part of the data, and it names the columns in the output’s data: line without the data$ prefix, so a reader can see at a glance which columns were compared. The list form is handy when your groups already sit in separate vectors, for example one per file.

Reading the kruskal.test() output line by line

The first step printed the following when the engine’s script ran:

	Kruskal-Wallis rank sum test

data:  count by spray
Kruskal-Wallis chi-squared = 54.691, df = 5, p-value = 1.511e-10

The title says which test ran. The name “rank sum” is a reminder that the test does not use the counts themselves. It pools all 72 observations, ranks them from smallest to largest, gives tied values the average of the ranks they span, and then compares the sums of the ranks in each group.

The data line repeats your outcome and group. Check it. A swapped formula, such as spray ~ count, does not stop with an error. In a check on this data it ranked the spray codes, treated each distinct count as a group and printed data: spray by count with far more degrees of freedom than there are sprays, so this line and the df are where the mistake shows.

Kruskal-Wallis chi-squared is the Kruskal-Wallis statistic H, here 54.69. R labels it chi-squared because H is referred to a chi-squared distribution. The value printed already includes the correction for ties: the R source divides H by a tie factor on every call, and there is no argument to turn it off. Insect counts are whole numbers and repeat often, so the correction is at work in this example. Report the printed value as H.

df is the number of groups minus one, here 5 for six sprays. If it is lower than you expect, a group has no rows. kruskal.test() drops empty factor levels before counting groups: on a three-spray subset of this data that still carried all six levels, it counted only the three sprays present.

p-value is the probability, if all six sprays had the same distribution of counts, of an H at least this large, taken from the chi-squared approximation. It is far below any usual threshold here, so the data are hard to reconcile with six identical distributions. The p-value is not the probability that the sprays are alike, and a small one says nothing about how large the differences are. That is the effect size’s job.

What the output does not show is the sample size. kruskal.test() silently drops any row with a missing outcome or group, so compare sum(complete.cases(data[, c("count", "spray")])) with the row count you expect. The engine’s record shows n = 72 with no rows dropped.

Why kruskal.test() here and not aov()

A one-way ANOVA compares means. The Kruskal-Wallis test compares the groups’ distributions through their ranks, and it needs only that larger counts mean more insects. Beall argued that counts like these do not suit a classic ANOVA as they stand and should be transformed first; a rank-based test sidesteps the transformation. Because the outcome is a count, there is also a third route: model the counts directly, for example with Poisson or negative binomial regression, which answers a question about rates rather than ranks. This post takes the rank-based route because its question is which sprays tend to leave more insects, not by what factor the counts change.

The engine’s assumption output shows the same thing from another angle. It ran a Shapiro-Wilk normality test in each group as context only, because the Kruskal-Wallis test does not assume normality. Two of the six groups showed evidence of non-normality, spray C (p = 0.0476) and spray D (p = 0.0027). The only check that is not context-only is the sample size, which passed.

Do not use those normality tests to choose between aov() and kruskal.test() after looking at the data. A test chosen by a preliminary test has distorted error rates conditional on that choice (Rochon and colleagues, 2012, show this for the two-sample case, although the overall error rate of the two-stage procedure stayed acceptable in their simulations), and with small groups the normality test has little power anyway. Choose the test from the design and the question, before you see the results, and say in the analysis plan which one is primary. The parametric test assumptions post goes through the checks and their fallbacks, and the guide to one-way ANOVA, Welch’s ANOVA and Kruskal-Wallis compares the three tests in detail.

What the result does and does not say

A significant Kruskal-Wallis result says that at least one group tends to produce larger or smaller values than the others. The test is often described as comparing medians, but that reading needs an extra assumption: that the groups’ distributions have the same shape and spread and differ only by a shift. Counts often spread more where they are higher, and nothing in this analysis checks that the six sprays share one shape, so the safe wording is about one group tending to have more insects than another, not about medians.

The Kruskal-Wallis test also says nothing about which groups differ. With six groups there are 15 pairs, and the omnibus result is compatible with one spray differing from the other five, or with two sets of sprays, or with every pair differing. That is why the next two steps exist.

Effect size: epsilon squared with an interval

Epsilon squared for ranks (Tomczak and Tomczak, 2014) is H divided by n minus 1. It reads as the share of the variation in the ranks that the grouping accounts for, from 0 to 1. For this example the engine reported 0.7703 with a 95% percentile-bootstrap interval from 0.7302 to 0.8336. Dividing the H above by 71 gives the same point estimate, a quick check you can run by hand.

The script’s second step printed it at two decimals:

Epsilon2 (rank) |       95% CI
------------------------------
0.77            | [0.73, 0.83]

The effect is large: the grouping accounts for most of the variation in the ranks. That is a statement about ranks, not about the raw counts, and it should be worded that way.

Two settings in that call change the interval, so set both on purpose. In effectsize 1.0.2, rank_epsilon_squared() defaults to alternative = "greater", a one-sided interval whose upper limit is fixed at 1, and to 200 bootstrap resamples. The engine asks for a two-sided interval, so it has a real upper limit, and for 1000 resamples, so the limits rest on more draws in each tail.

Because the interval comes from resampling, it also changes with the seed (setting seeds and pinning package versions in R shows by how much). set.seed() immediately before the call makes it reproducible; a different seed gives slightly different limits, never a different point estimate. The post on Dunn’s test ran the same call on the same data with another seed, which is why its interval differs a little from this one around the same 0.7703.

Some tutorials report a different rank-based effect size, eta squared based on H, which is H minus the number of groups plus one, divided by n minus the number of groups, and comes out a little smaller than epsilon squared on this data. The two are not interchangeable, so name the one you report. The post on effect sizes covers how rank-based and mean-based effect sizes relate.

Dunn’s follow-up in R

Dunn’s test (Dunn, 1964) compares each pair of groups using the mean ranks from the same joint ranking the Kruskal-Wallis test built, with a standard error that uses the spread of all the ranks. PMCMRplus::kwAllPairsDunnTest() runs it and, with no p.adjust.method given, applies the Holm adjustment (Holm, 1979). The engine named Holm explicitly. With Holm, the family-wise error rate across the 15 pairs is controlled only approximately, because Dunn’s z statistics rest on a large-sample normal approximation. When the script ran, print() showed a lower-triangle matrix of adjusted p-values, and R raised a warning. The printout below is trimmed: the header lines above the matrix are left out and the warning’s call is shortened.

Warning message:
In kwAllPairsDunnTest.default(...) :
  Ties are present. z-quantiles were corrected for ties.

  A       B       C       D       E
B 1.00000 -       -       -       -
C 2.3e-05 5.1e-06 -       -       -
D 0.01276 0.00482 0.58571 -       -
E 0.00118 0.00035 1.00000 1.00000 -
F 1.00000 1.00000 3.3e-06 0.00384 0.00025

The warning is not an error. It says the counts contain ties and that the z statistics were corrected for them, the same situation kruskal.test() handled silently.

Each cell is one pair’s adjusted p-value. Read across the row and up the column: the cell in row D, column A is A vs D. The engine recorded the same values at more digits. Nine pairs fall below 0.05: each of A, B and F against each of C, D and E. For example, A vs D had an adjusted p-value of 0.01276, B vs D 0.004823, A vs E 0.001179 and D vs F 0.003838. The other six pairs, the ones within A, B and F and within C, D and E, had adjusted p-values from 0.5857 to 1.

The engine condenses this into letters: A, B and F share the letter a, and C, D and E share the letter b. Groups sharing a letter were not separated at the adjusted threshold, which is not the same as being shown equal.

Direction comes from the effect size, not the p-value. The engine added a rank-biserial correlation for each pair, from effectsize::rank_biserial(), with an unadjusted 95% interval. For A vs D it was 0.9514, from 0.8798 to 0.9808: plots under A almost always had more insects than plots under D. For A vs B it was -0.1389, from -0.545 to 0.32, an interval that includes zero.

Two pairs are worth a second look. For C vs D the rank-biserial correlation was -0.7222, from -0.8818 to -0.4143, and for C vs E it was -0.4653, from -0.7511 to -0.03252. Both intervals exclude zero, while Dunn’s adjusted p-values were 0.5857 and 1. The two answer different questions: the correlation re-ranks only the 24 plots of the two sprays in the pair, while Dunn’s test works on all 72 ranks. The post on Dunn’s test explains this disagreement in detail and puts the Holm and Bonferroni adjustments side by side on the same 15 pairs.

Pitfalls and silent behaviours

A numeric group column. If your groups are coded 1 to 6, kruskal.test() still works: it converts the group to a factor internally, and in a check on this data a numeric version of the spray column gave the same statistic and degrees of freedom. rank_epsilon_squared() does not. It stops with “groups cannot be numeric”, so wrap the column in factor() first, as the script does.

Missing values. Rows with a missing outcome or group are dropped without a message by kruskal.test() and kwAllPairsDunnTest(); rank_epsilon_squared() at least warns that it dropped them. Count the complete rows yourself and report that n.

pairwise.wilcox.test() is a different test. Base R’s pairwise.wilcox.test() runs a separate Wilcoxon rank-sum test on each pair, re-ranking the two groups each time. That is a valid procedure with its own adjustment, but it is not Dunn’s test and it no longer uses the ranking the omnibus test used. Say which one you ran.

Not every package defaults to Holm. kwAllPairsDunnTest() does, and so does FSA::dunnTest(). Other implementations, in R and elsewhere, may adjust nothing unless asked. Pass the adjustment explicitly so the code says what the report says. The guide to multiple comparisons covers when Holm, Bonferroni or a false discovery rate method fits.

Gating Dunn’s test on the omnibus result. Nothing in R stops you running kwAllPairsDunnTest() after a non-significant kruskal.test(), or forces you to. Whether the follow-up waits for a significant omnibus result is a decision for the analysis plan, and the post on Dunn’s test weighs the two rules.

Exact p-values. kruskal.test() always uses the chi-squared approximation. With very small groups, an exact or permutation version from another package may be preferable; with 12 plots per group, as here, the approximation is the usual choice.

How to report it

A Kruskal-Wallis report names the test, the statistic with its degrees of freedom, the p-value, the effect size with its interval and how the interval was computed, the follow-up test and its adjustment. For this example:

A Kruskal-Wallis test on all 72 plots found that insect counts were not distributed alike across the six sprays (H = 54.69, df = 5, p < 0.001), with a large effect (epsilon squared = 0.7703, 95% CI 0.7302 to 0.8336, two-sided percentile bootstrap, 1000 resamples, seed 20261101). In Dunn’s follow-up with Holm adjustment, sprays A, B and F formed one set and C, D and E another: all nine pairs across the two sets were separated (adjusted p from < 0.001 to 0.01276), with fewer insects under C, D and E, and none of the six pairs within a set was (adjusted p from 0.5857 to 1).

Naming the seed and the number of resamples matters here because the interval would move with either. Put the full table of 15 pairs, with each pair’s effect size and interval, in a table or supplement. If a reader asks why you did not use ANOVA, the answer is in the design: an ordered count outcome that Beall (1942) argued needs transforming before a classic ANOVA applies, and a question about which groups tend to be higher.

The research question here is about field plots in a published teaching dataset. The analysis says nothing about how any product performs outside these 72 plots. The APA reporting post has the reporting formats for the other common tests.

The same analysis as a Celeus run

The worked example above is a Celeus engine run: Kruskal-Wallis test, epsilon squared with a seeded two-sided bootstrap interval, and Dunn’s test with an adjustment you choose from Holm, Benjamini-Hochberg, Bonferroni or none. The run page shows the exact R script the engine wrote, before you download anything, and that script is the one below this article. It was rerun in the same software environment and gave the same Kruskal-Wallis statistic and p-value. The rank-biserial values per pair come from a separate effectsize call and are not in the script.

The guide to reproducible statistical analysis explains what a reviewer needs to rerun a result like this one, and the method reference for the Kruskal-Wallis test lists what the engine computes for it. For the two-group case, where the rank-based test is the Mann-Whitney test, start from the t-tests explainer.

# Results

- Test: Kruskal-Wallis test
- n: 72
- statistic: 54.69
- p-value: < 0.001
- effect size: epsilon-squared = 0.7703 (95% CI 0.7302 to 0.8336)
  (two-sided percentile-bootstrap CI at the run's confidence level)
- seed: 20261101

## Assumption checks
- normality by group (Shapiro-Wilk): context only, passed=FALSE (evidence of non-normality - normality is not assumed by this rank-based test)
- sample size: passed=TRUE (adequate n)

## Post-hoc (Dunn's test; correction: Holm)
- A vs B: adj p = 1 | rank-biserial r -0.1389 [-0.545, 0.32]
- A vs C: adj p = 2.347e-05 | rank-biserial r 0.9931 [0.9823, 0.9973]
- B vs C: adj p = 5.098e-06 | rank-biserial r 0.9931 [0.9823, 0.9973]
- A vs D: adj p = 0.01276 | rank-biserial r 0.9514 [0.8798, 0.9808]
- B vs D: adj p = 0.004823 | rank-biserial r 0.9583 [0.8964, 0.9836]
- C vs D: adj p = 0.5857 | rank-biserial r -0.7222 [-0.8818, -0.4143]
- A vs E: adj p = 0.001179 | rank-biserial r 1 [1, 1]
- B vs E: adj p = 0.0003451 | rank-biserial r 1 [1, 1]
- C vs E: adj p = 1 | rank-biserial r -0.4653 [-0.7511, -0.03252]
- D vs E: adj p = 1 | rank-biserial r 0.3264 [-0.1319, 0.6697]
- A vs F: adj p = 1 | rank-biserial r -0.1944 [-0.584, 0.2678]
- B vs F: adj p = 1 | rank-biserial r -0.03472 [-0.467, 0.411]
- C vs F: adj p = 3.334e-06 | rank-biserial r -1 [-1, -1]
- D vs F: adj p = 0.003838 | rank-biserial r -0.9583 [-0.9836, -0.8964]
- E vs F: adj p = 0.0002496 | rank-biserial r -1 [-1, -1]

groups: A (a), B (a), C (b), D (b), E (b), F (a) - groups sharing a letter are not significantly different at the adjusted threshold.

## Power
- Power is not calculated for this test. The matching parametric test gives a rough guide.

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Kruskal-Wallis test  |  engine 0.41.7  |  catalog c40-dc08
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261101)
.data_path <- "data/kruskal-wallis-test-in-r.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("count", "spray")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- kruskal.test(`count` ~ `spray`, data = data)
print(result)

# Effect size (epsilon_squared) and its confidence interval.
effect <- { set.seed(20261101); .es_d <- data; .es_d[["spray"]] <- as.factor(.es_d[["spray"]]); effectsize::rank_epsilon_squared(.es_d[["count"]], .es_d[["spray"]], ci = 0.95, alternative = "two.sided", iterations = 1000L) }
print(effect)

# Post-hoc pairwise comparisons (Dunn's test; correction: Holm).
data[["spray"]] <- as.factor(data[["spray"]])
posthoc <- { .ph_d <- data; .ph_d[["spray"]] <- as.factor(.ph_d[["spray"]]); PMCMRplus::kwAllPairsDunnTest(`count` ~ `spray`, data = .ph_d, p.adjust.method = "holm") }
print(posthoc)