Kruskal-Wallis Test in R: kruskal.test()

Q: Do I need to install a package to run a Kruskal-Wallis test in R?

No. kruskal.test() is in the stats package, which loads with base R. You only need packages for the things base R omits: a post-hoc test ( rstatix , FSA or dunn.test ), an effect size ( rstatix ), or an exact/Monte Carlo null distribution ( coin ).

Q: Why do dunn.test and rstatix give me different p-values on the same data?

Because their defaults differ on two axes at once. dunn.test::dunn.test() defaults to no multiplicity adjustment and reports one-tailed p-values intended for comparison against α/2; rstatix::dunn_test() and FSA::dunnTest() both default to Holm adjustment and two-sided p-values. Neither is wrong, but only one of them is the number you probably meant to report. Always state the package and adjustment.

Most people arrive at kruskal.test() in R the same way: an ANOVA assumption check looked bad, and Kruskal-Wallis is the test everyone names as the non-parametric fallback. That reflex is where the mistake usually happens. Kruskal-Wallis fixes exactly one of the things that can go wrong with a one-way ANOVA, and it makes another one worse. Before you switch, you need to know which assumption failed — and if the answer is unequal variances, R already has a better one-line answer for you.

First: which ANOVA assumption actually failed?

Kruskal-Wallis is not a general-purpose repair for a misbehaving ANOVA. Kroeger and colleagues, writing in the American Journal of Clinical Nutrition, ran simulations under conditions typical of applied research and found that Kruskal-Wallis type I error rates “deviated from the expected significance level” when variances were unequal, with “greater deviation from the expected type I error rate…observed as the heterogeneity increased, especially in the presence of an imbalanced sample size.” They note the practice persists “sometimes because of the mistaken rationale that the test corrects for heteroscedasticity.” It does not.

So map the failure to the fix before you type anything:

Residuals are clearly non-normal, and a transform doesn’t help. This is the case Kruskal-Wallis is actually for. Use kruskal.test().
The outcome is ordinal (Likert-type items, pain scores, stage or grade). Means are not defined on that scale. Use kruskal.test().
Variances differ across groups, but the data are otherwise reasonable. Use Welch’s ANOVA, not Kruskal-Wallis. In R that is oneway.test() — and note its signature is oneway.test(formula, data, subset, na.action, var.equal = FALSE), so Welch is already the default. Writing oneway.test(score ~ group, data = df) gives you “an approximate method of Welch (1951)…which generalizes the commonly known 2-sample Welch test to the case of arbitrarily many samples.” Passing var.equal = TRUE is what turns it back into the classical F test.
Both problems at once, or badly unequal groups plus unequal spread. Neither classical test is safe; see the rank-based Behrens-Fisher options below.
Observations are repeated on the same subjects. Independence is violated, which Kruskal-Wallis also requires. Use the Friedman test or a mixed model.

If you haven’t formally checked the spread across groups yet, CASRAI’s guide to Levene’s test for equality of variances covers that check and how to read it. For the parametric side of this decision — aov(), car::Anova() and Type III sums of squares — see running ANOVA in R.

Don’t let a normality test pick the test for you

The natural-looking R workflow — run shapiro.test() on each group, then branch to aov() or kruskal.test() depending on the result — is a formally invalid procedure, not just an inelegant one. Rochon, Gondan and Kieser simulated exactly this two-stage strategy and found that conditioning the choice of test on a preliminary Shapiro-Wilk result pushes the realised type I error rate away from the nominal 5% in both directions, depending on the underlying distribution and the pre-test alpha. Their conclusion is blunt: “From a formal perspective, preliminary testing for normality is incorrect and should therefore be avoided.” They also point out the practical trap that at small sample sizes — the case where normality matters most — “the Shapiro-Wilk test lacks power to detect deviations from normality,” so it tends to wave through exactly the data you were worried about.

The defensible alternative is to decide from what you know about the measurement — its scale, its known skew, prior data on the same instrument — and from diagnostic plots, then commit to that analysis. Choose the test before you see the p-value it produces, and say in your methods section why.

Running kruskal.test()

The function ships with base R in the stats package; nothing to install. It takes two forms:

kruskal.test(x, ...)
## Default S3 method:
kruskal.test(x, g, ...)
## S3 method for class 'formula':
kruskal.test(formula, data, subset, na.action, ...)

In practice the formula interface is the one to use, because it keeps the response and the grouping variable in the same data frame:

kruskal.test(pain_score ~ protocol, data = trial)

The default method is the fallback for data that isn’t in a tidy frame — kruskal.test(x, g) with a numeric vector and a parallel grouping vector, or kruskal.test(list_of_vectors). R’s documentation states the formula argument is “a formula of the form response ~ group where response gives the data values and group a vector or factor of the corresponding groups.”

What comes back is an htest object with, per the documentation, “statistic the Kruskal-Wallis rank sum statistic; parameter the degrees of freedom; p.value the p-value of the test.” You can pull them out directly — res$statistic, res$parameter, res$p.value — which is what you want when you are assembling a results table rather than reading console output.

Note carefully what R says the null hypothesis is. The help page states the function tests “the null that the location parameters of the distribution of x are the same in each group.” That phrasing is a location claim, and it is only justified when the group distributions have broadly similar shape and spread. When shapes genuinely differ, a significant result tells you the groups differ distributionally, not that their medians differ — a distinction covered in more depth in CASRAI’s guide to when to use the Kruskal-Wallis test and how to report it.

Five things kruskal.test() does without telling you

These are visible in the function’s source, and each one has bitten someone’s analysis.

It coerces your grouping variable with g <- factor(g). A numeric or character group column works fine. But this also re-levels, so if you subset a data frame and leave empty factor levels behind, they are dropped rather than raising an error — and your degrees of freedom quietly change. Check parameter against the number of groups you meant to compare.
It drops incomplete rows silently via OK <- complete.cases(x, g). Any row missing either the outcome or the group is removed with no message. The N used by the test is therefore not necessarily nrow(your_data) — which matters, because every Kruskal-Wallis effect size divides by N. Compute your own N from the complete cases, not from the data frame.
Tie correction is always applied, and cannot be switched off. The statistic is computed as ((12 * STATISTIC / (n * (n + 1)) - 3 * (n + 1)) / (1 - sum(TIES^3 - TIES) / (n^3 - n))), where the denominator is the standard tie correction. This is the right default — without it, ordinal data with many repeated values yields a downward-biased statistic — but it means an R result will not match software or a textbook that reports the uncorrected value.
It refuses only the degenerate case. The source contains if (k < 2L) stop("all observations are in the same group"). There is no minimum group size check. Kruskal-Wallis with two observations per group will run and return a p-value that the chi-square approximation does not really support.
There is no exact test. Base R gives you the chi-square approximation on k − 1 degrees of freedom and nothing else. For small samples, use the coin package, whose kruskal_test(formula, data, subset = NULL, weights = NULL, ...) accepts a distribution argument with "asymptotic" (the default), "approximate" for Monte Carlo resampling, and "exact", which is available for univariate two-sample problems. For a small three-group design, distribution = approximate(nresample = 10000) is the practical choice.

The post-hoc step: three R packages, three different p-value columns

A significant Kruskal-Wallis result is an omnibus finding — at least one group differs from at least one other. Locating the difference needs Dunn’s test, and this is where R will hand you three defensible-looking but non-comparable answers depending on which package you loaded.

dunn.test::dunn.test() has two defaults that will burn you. Its signature is dunn.test(x, g=NA, method=p.adjustment.methods, kw=TRUE, label=TRUE, ...), and p.adjustment.methods is c("none", "bonferroni", "sidak", "holm", "hs", "hochberg", "bh", "by"). Because "none" comes first, the default is no multiple-comparison adjustment at all — the documentation confirms “the default is no adjustment for multiple comparisons.” Second, it reports one-tailed p-values: the docs state “the default is to express p-value = P(Z ≥ |z|), and reject Ho if p ≤ α/2,” with the altp=TRUE option switching to “p-value = P(|Z| ≥ |z|), and reject Ho if p ≤ α.” These give “identical test results” only if you apply the matching threshold. Compare this package’s default output to 0.05 and you are running an unadjusted test at twice your stated alpha.

FSA::dunnTest() fixes both defaults. Its signature is dunnTest(x, g, method = dunn.test::p.adjustment.methods[c(4, 2:3, 5:8, 1)], two.sided = TRUE, altp = two.sided, ...) — the reordering puts "holm" first, and the documentation notes that in contrast to the underlying function, “the p-values are adjusted by default with the ‘holm’ method, and two-sided p-values are returned by default.”

rstatix::dunn_test() matches the commercial convention. Signature dunn_test(data, formula, p.adjust.method = "holm", ref.group = NULL, detailed = FALSE, effect.size = FALSE), and its documentation states the default “is to perform a two-sided Dunn test like the well known commercial softwares, such as SPSS and GraphPad.”

library(rstatix)

trial %>% kruskal_test(pain_score ~ protocol)
trial %>% dunn_test(pain_score ~ protocol, p.adjust.method = "holm")

Whichever you use, state the package and the adjustment method in your write-up, because the numbers are not interchangeable. Holm is the usual default rather than Bonferroni because it is uniformly less conservative at the same familywise error rate; CASRAI’s guide to the Bonferroni correction and multiple comparisons covers that trade-off.

Why not just loop pairwise wilcox.test()?

Because it discards the information the omnibus test was built on. Dunn’s test, as rstatix’s documentation puts it, “incorporates the pooled variance estimate implied by the null hypothesis of the Kruskal-Wallis test” and “retains the dependent ranking that produced the Kruskal-Wallis test statistic” — whereas the Wilcoxon rank sum test fails on both counts as a post-hoc procedure. Pairwise wilcox.test() re-ranks each pair in isolation, which is why it can disagree with both the omnibus result and Dunn’s test, most visibly when group sizes are unequal.

Effect size

Report one. An omnibus H and a p-value say nothing about magnitude, and with a large N a trivial rank difference clears any threshold. rstatix::kruskal_effsize() takes a method argument of "eta2" or "epsilon2"; for the default eta-squared it documents the formula as “eta2[H] = (H – k + 1)/(n – k); where H is the value obtained in the Kruskal-Wallis test; k is the number of groups; n is the total number of observations,” and offers bootstrap confidence intervals via ci = TRUE. Its stated interpretation bands are “0.01- < 0.06 (small effect), 0.06 – < 0.14 (moderate effect) and >= 0.14 (large effect).” Treat those as a calibration convention, not a verdict.

trial %>% kruskal_effsize(pain_score ~ protocol, ci = TRUE)

When ranks still aren’t the right answer

If you have unequal variances and a genuinely non-normal or ordinal outcome, you are in the nonparametric Behrens-Fisher problem, and neither kruskal.test() nor oneway.test() is a clean fit. The rankFD package targets exactly this. Its main function is rankFD(formula, data, alpha = 0.05, CI.method = c("logit", "normal"), effect = c("unweighted", "weighted"), hypothesis = c("H0F", "H0p"), ...), and the documentation explains that testing hypotheses in terms of relative effects — hypothesis = "H0p" — “allows for variance heteroscedasticity even under the null hypothesis of no treatment effect and thus covers the Nonparametric Behrens-Fisher problem.” It provides a Wald-type statistic, an ANOVA-type statistic, and multiple contrast tests, and reports relative effects with confidence intervals rather than only a p-value.

Reporting a Kruskal-Wallis analysis from R

A complete report names the software path, not just the numbers, because as shown above the numbers depend on it. Include:

The test statistic with degrees of freedom and the exact p-value, taken from statistic, parameter and p.value.
The N actually analysed, after complete.cases() dropped incomplete rows — and how many rows that removed.
Medians and interquartile ranges per group as descriptives, not means and standard deviations.
An effect size with its method (eta2 or epsilon2) named.
If significant: the post-hoc test, the R package that produced it, and the adjustment method.
Why Kruskal-Wallis rather than ANOVA — stated as a property of the measurement, not as the outcome of a normality test.

Frequently asked questions

Do I need to install a package to run a Kruskal-Wallis test in R?

No. kruskal.test() is in the stats package, which loads with base R. You only need packages for the things base R omits: a post-hoc test (rstatix, FSA or dunn.test), an effect size (rstatix), or an exact/Monte Carlo null distribution (coin).

Why does my kruskal.test() result differ from another program’s?

The two usual causes are ties and missing data. R always applies the tie correction, with no option to disable it, so a result computed without that correction will differ. R also silently drops rows where either the outcome or the group is NA, so if the other program handled missingness differently you are comparing tests run on different N.

Why do dunn.test and rstatix give me different p-values on the same data?

Because their defaults differ on two axes at once. dunn.test::dunn.test() defaults to no multiplicity adjustment and reports one-tailed p-values intended for comparison against α/2; rstatix::dunn_test() and FSA::dunnTest() both default to Holm adjustment and two-sided p-values. Neither is wrong, but only one of them is the number you probably meant to report. Always state the package and adjustment.

My ANOVA failed Levene’s test. Should I use Kruskal-Wallis?

Usually not. Unequal variance is the one problem Kruskal-Wallis does not solve — Kroeger et al. found its type I error rate deviates further from nominal as heterogeneity increases, especially with imbalanced group sizes. If normality is broadly acceptable, oneway.test(y ~ g, data = df) gives you Welch’s ANOVA by default. If both assumptions fail, look at rankFD.

Should I run shapiro.test() first and switch to kruskal.test() if it’s significant?

No. That two-stage procedure changes the type I error rate of whichever test you end up running, and Rochon, Gondan and Kieser concluded that “preliminary testing for normality is incorrect and should therefore be avoided.” Decide from the measurement scale, prior knowledge and diagnostic plots, then commit.

Can I use kruskal.test() on repeated measures?

No. It requires independent observations, and repeated measurements on the same subject are not independent. Use friedman.test() for a complete balanced block design, or a mixed-effects model if the design is unbalanced or has covariates.

What if kruskal.test() is significant but no Dunn comparison survives adjustment?

This is a normal and reportable outcome, not an error. The omnibus test pools evidence across all groups and can detect a spread of ranks that no single adjusted pairwise contrast is powered to isolate. Report both honestly rather than dropping the adjustment to manufacture a pairwise result.