Kruskal-Wallis Test

In plain language

The Kruskal-Wallis test answers the question: "Do three or more groups have meaningfully different outcome distributions?" — without assuming the data follow a bell-curve shape. It works by replacing every observed value with its rank (1st smallest, 2nd smallest, and so on) across all groups combined, then checking whether the groups ended up with surprisingly unequal shares of high and low ranks. If the answer is yes, a follow-up procedure called Dunn's test identifies which specific pairs of groups differ. It is the standard nonparametric substitute for one-way ANOVA and is commonly used in health outcomes research to compare length of stay, cost, or patient-reported scores across three or more treatment groups when the numbers are skewed or the sample is small.

What the Kruskal-Wallis test does and why it exists

One-way ANOVA tests whether the means of three or more groups are equal by partitioning variance — but it assumes the outcome is approximately normally distributed within each group and that the variances are roughly equal. In real-world evidence and health outcomes research, the outcomes most worth comparing across groups — length of stay, total cost per episode, patient-reported pain scores, time to treatment escalation — are routinely right-skewed, bounded, or ordinal. For these outcomes at small-to-moderate sample sizes, the ANOVA normality assumption is not tenable, and the Kruskal-Wallis H test provides a valid alternative.

The Kruskal-Wallis test works by collapsing all N observations from k groups into a single pooled ranking: every value is replaced by its rank from 1 (smallest) to N (largest), ignoring group membership. The test statistic H measures how unequally the rank mass is distributed across groups relative to what would be expected if group labels were assigned at random. Under the null hypothesis of no difference among groups, H follows an approximate chi-square distribution with k − 1 degrees of freedom, provided each group has at least five observations. For very small groups (n < 5 per group), exact tables or permutation-based p-values should be used instead.

The H statistic: mechanics and tie correction

The H statistic is:

H = [12 / (N(N+1))] * sum_j[R_j^2 / n_j] - 3(N+1)

where N is the total number of observations, k is the number of groups, n_j is the sample size of group j, and R_j is the sum of ranks assigned to group j. Under the null, all groups share the same population distribution and the expected rank sum for group j is n_j * (N+1)/2, so H measures the weighted squared deviation of observed rank sums from their expected values.

When ties are present — multiple observations with the same value — they are assigned the average of the ranks they would have received if they were distinct. Ties shrink the variance of the rank distribution and can inflate H if uncorrected. The standard tie correction divides H by:

C = 1 - [sum_g(t_g^3 - t_g)] / (N^3 - N)

where t_g is the number of observations in tie group g. The corrected statistic is H/C. In practice the correction has little effect unless there are many large tie groups, but it should always be applied and is the default in all major software implementations.

The chi-square approximation and degrees of freedom

The chi-square approximation for H is valid when each group has at least five observations. The degrees of freedom are k − 1: for three groups, df = 2; for four groups, df = 3. A p-value below the chosen alpha (typically 0.05) leads to rejection of the null hypothesis that all k groups share the same distribution. The p-value alone does not indicate which pairs differ — only that at least one pair does.

Stochastic dominance, not medians

The Kruskal-Wallis test shares the same subtle interpretation as the Mann-Whitney U test (its two-group special case): it tests a form of stochastic dominance, not equality of medians. The null hypothesis is that a randomly chosen observation from any one group has an equal probability of being larger than a randomly chosen observation from any other group. This is equivalent to testing median equality only under the additional assumption that the distributions across groups differ solely in location — that is, they have the same shape and spread and merely shift left or right. When the groups have different shapes (one is bimodal, another is heavily skewed) or different spreads, a significant Kruskal- Wallis result signals that the distributions differ in some way, but it does not cleanly identify a median shift. Analysis reports should acknowledge this: describe the test as examining rank distributions rather than medians, and accompany the p-value with group medians (and interquartile ranges) for interpretive context.

From omnibus rejection to post-hoc pairwise comparisons

The Kruskal-Wallis test is an omnibus test: it asks whether any group differs from any other, not which specific pairs differ. When H is significant, post-hoc pairwise comparisons are needed to locate the source of the difference.

The correct post-hoc procedure is Dunn's test (Dunn 1964), which compares each pair of groups using the difference in their mean ranks from the pooled ranking — not separate pairwise Mann-Whitney tests. Applying repeated pairwise Mann-Whitney tests without multiplicity adjustment inflates the family-wise type-I error rate: with k = 4 groups and six pairwise comparisons at alpha = 0.05, the probability of at least one false positive is approximately 1 − (0.95)^6 ≈ 0.26. Dunn's test controls this by using the pooled variance from the full ranking, and the resulting p-values are then adjusted for multiple comparisons using the Bonferroni correction (conservative), the Holm step-down procedure (uniformly more powerful than Bonferroni), or the Benjamini-Hochberg false discovery rate correction when many comparisons are expected. The choice among these corrections should be pre-specified; defaulting to Bonferroni is the safest option when confirmatory inference is required and the number of comparisons is small (k < 6).

RWE applications and common patterns

The Kruskal-Wallis test appears in HEOR and RWE in several recurring patterns:

• Length of stay or cost across three or more lines of therapy: comparing index

• Patient-reported outcome scores across site or payer type: ordinal PRO scores (pain,

• Unadjusted descriptive comparison in Table 2 of an observational study: presenting

• Multi-site quality comparisons: comparing median time-to-diagnosis, door-to-balloon

Clustering violation: a critical warning for EHR and registry data

The Kruskal-Wallis test assumes that all observations are independent. In EHR, registry, and claims data, patients are routinely nested within providers, practices, or facilities. Applying Kruskal-Wallis to individual patient rows without accounting for this clustering produces artificially narrow confidence intervals and inflated test statistics, because patients within the same provider share unmeasured factors that make their outcomes more similar than observations from different providers would be. When clustering is present, either aggregate to the provider level before testing, use a mixed-effects model, or apply a nonparametric test that accounts for the clustered design (such as a permutation test with cluster-level exchangeability). This is not a theoretical concern — provider-level intraclass correlation (ICC) for cost and utilization outcomes in claims data commonly exceeds 0.05 and can exceed 0.20 for specialty care, making the independence assumption materially wrong.

Pros, cons, and trade-offs

Pros: - Valid under broad distributional assumptions; no requirement for normality, equal variances, or any specific distributional family within groups. - Robust to outliers: an extreme value in one group is capped at rank N and cannot distort H the way it distorts an ANOVA F statistic. - Well-suited to ordinal outcomes (Likert scales, PRO scores, severity grades) where the interval assumption required by ANOVA is not justified. - Exact in small samples when combined with permutation-based critical values; the chi-square approximation is adequate for n_j ≥ 5. - Directly extensible to the Jonckheere-Terpstra test when the k groups have a natural order (dose levels, line of therapy), which is more powerful than Kruskal-Wallis for detecting monotone dose-response relationships.

Cons: - Lower statistical power than one-way ANOVA when the normality assumption actually holds, because the rank transformation discards information about the magnitude of differences. - The omnibus H statistic requires post-hoc testing (Dunn's test) to identify which pairs differ; each step adds complexity and reduces power relative to a planned contrast. - Cannot adjust for covariates directly: the rank-based framework does not extend naturally to multivariable regression. Quantile regression or a ranked-outcome GLM can partly address this limitation but are less familiar to most audiences. - The test statistic does not produce an inherently interpretable effect size on the original measurement scale; the epsilon-squared or eta-squared rank-based effect size must be computed separately and reported alongside H and the p-value. - Tests stochastic dominance rather than median equality under heterogeneous shapes — requires careful communication to avoid the common mischaracterization as a median test.

When to use

• Three or more independent groups where the outcome is continuous but non-normal,
• Sensitivity analysis alongside a one-way ANOVA primary analysis, to verify that
• Descriptive comparison in Table 2 or Table 3 of an observational RWE study,
• Ordinal outcomes: PRO scales, severity classifications, adherence categories —
• Any multi-group comparison where one ANOVA assumption is clearly violated: if

When NOT to use

• *Repeated-measures or pre-post-post designs where the same patients appear in multiple

• When adjusted estimates are needed: Kruskal-Wallis cannot include covariates. If the

• When the groups have a meaningful natural order (dose-response, line of therapy):

• When the inference target is the mean cost difference: for budget-impact analyses,

• When clustering is present and unaccounted for: applying Kruskal-Wallis row-by-row to

• *As a substitute for the Jonckheere-Terpstra test when detecting monotone trends in

Interpreting the output

In the worked example, nine patients across three hospitals (three per hospital) had inpatient lengths of stay ranging from 1 to 9 days. After pooled ranking, Hospital A received rank sum R_A = 6, Hospital B received R_B = 15, and Hospital C received R_C = 24. The rank-sum contribution is 6²/3 + 15²/3 + 24²/3 = 12 + 75 + 192 = 279, yielding H = (12/90) × 279 − 30 = 37.2 − 30 = 7.2. On df = 2, the chi-squared critical value at alpha = 0.05 is approximately 5.99, and p ≈ 0.027.

(1) Formal interpretation. H = 7.2 exceeds the critical value of approximately 5.99 on df = 2, placing p ≈ 0.027 below alpha = 0.05. Under the null hypothesis that all three groups share the same distribution, a rank imbalance this extreme arises in approximately 2.7% of random rank assignments. The omnibus test rejects the null, indicating that at least one hospital's LOS rank distribution is stochastically different from the others. This result pertains to stochastic dominance — whether one group tends to produce larger values — not specifically to differences in medians unless the three distributions are assumed to differ only in location. Pairwise identification of which hospitals differ requires Dunn's test with an appropriate multiplicity correction.

(2) Practical interpretation. The rank distribution of length of stay differs across the three hospitals in a way unlikely to arise by chance. Hospital C consistently occupies the highest ranks (longest stays) and Hospital A the lowest, which is visible directly from the rank sums (6 vs 15 vs 24). Whether this reflects differences in patient case mix, clinical practice, or discharge protocols cannot be determined from an unadjusted Kruskal- Wallis test — the comparison is descriptive. Dunn's test with Bonferroni or Holm adjustment for the three pairwise comparisons is the appropriate next step to identify which specific hospital pairs drive the overall signal.

Worked example

Scenario

A health systems analyst is comparing the length of inpatient stay (days) for patients admitted with a primary diagnosis of heart failure across three hospitals — Community Hospital (Group A), Regional Medical Center (Group B), and Academic Medical Center (Group C) — using a small convenience sample of three patients per site. The analyst wants to test whether length of stay differs across the three facilities using a method that does not assume normality, and then manually verify the Kruskal-Wallis H statistic from first principles.

Dataset

Length of stay (days) for nine heart failure admissions across three hospitals. All nine values are pooled and ranked together for the Kruskal-Wallis test. There are no tied values in this example.

Steps

• Pool all 9 observations and assign ranks 1 through 9 in ascending order of LOS. Since all values are distinct, there are no ties. Value 1 gets rank 1, value 2 gets rank 2, value 3 gets rank 3, value 4 gets rank 4, value 5 gets rank 5, value 6 gets rank 6, value 7 gets rank 7, value 8 gets rank 8, value 9 gets rank 9.

• Compute the rank sum for each group. Group A (Community): patients A1, A2, A3 received ranks 1, 2, 3 respectively. R_A = 1 + 2 + 3 = 6. Group B (Regional): patients B1, B2, B3 received ranks 4, 5, 6. R_B = 4 + 5 + 6 = 15. Group C (Academic): patients C1, C2, C3 received ranks 7, 8, 9. R_C = 7 + 8 + 9 = 24.

• Verify that rank sums total to N(N+1)/2 = 910/2 = 45. Check: 6 + 15 + 24 = 45. Correct — all ranks are accounted for.

• Compute the bracketed term in the H formula. Each group has n_j = 3 observations. For Group A: R_A squared equals 36; divide by n_A = 3 to get 36/3 = 12. For Group B: R_B squared equals 225; divide by n_B = 3 to get 225/3 = 75. For Group C: R_C squared equals 576; divide by n_C = 3 to get 576/3 = 192. Bracketed sum = 12 + 75 + 192 = 279.

• Apply the H formula with N = 9: the leading constant is 12 divided by N(N+1), which is 12/(910) = 12/90. Multiply by the bracketed sum: (12/90)279 = 3348/90 = 37.2. Subtract 3(N+1) = 3*10 = 30. So H = 37.2 - 30 = 7.2.

• Compare H = 7.2 to the chi-square critical value with df = k - 1 = 3 - 1 = 2. The chi-square critical value at alpha = 0.05 with df = 2 is 5.991. Since 7.2 > 5.991, we reject the null hypothesis (p approximately 0.027). At least one hospital's length-of-stay distribution differs from the others.

• Because the omnibus test is significant, Dunn's test with Bonferroni correction would be applied to identify which pairs differ (A vs B, A vs C, B vs C). With three pairs, the Bonferroni-adjusted alpha threshold is 0.05/3 = 0.0167. Given the clean rank separation in this example, the A-vs-C comparison yields the smallest adjusted p-value. In practice, Dunn's test is run in software rather than computed by hand.

Result

N = 9, k = 3, n per group = 3. Rank sums: R_A = 6, R_B = 15, R_C = 24; total = 6 + 15 + 24 = 45. Bracketed sum = 12 + 75 + 192 = 279. H = (12/90)*279 - 30 = 3348/90 - 30 = 37.2 - 30 = 7.2. With df = 2, chi-square critical value at 0.05 is 5.991; H = 7.2 exceeds 5.991, so p is approximately 0.027. The test rejects the null; post-hoc Dunn's test with multiplicity adjustment identifies which hospital pairs differ.

Runnable example

import math
from scipy import stats

# ── Motivating dataset: LOS (days) across three hospitals ──
group_a = [1, 2, 3]   # Community Hospital
group_b = [4, 5, 6]   # Regional Medical Center
group_c = [7, 8, 9]   # Academic Medical Center

# ── 1. Kruskal-Wallis omnibus test ──
H, p_kw = stats.kruskal(group_a, group_b, group_c)
print(f"Kruskal-Wallis H = {H:.4f}, p = {p_kw:.4f}, df = {3 - 1}")
# Expected: H = 7.2, p ≈ 0.0273

# ── 2. Epsilon-squared effect size (eta^2 for Kruskal-Wallis) ──
N = len(group_a) + len(group_b) + len(group_c)
k = 3
eps_sq = (H - k + 1) / (N - k)
print(f"Epsilon-squared (eta^2) = {eps_sq:.4f}")
# Effect size interpretation: 0.01 small, 0.06 medium, 0.14 large

# ── 3. Descriptive summary: medians and IQRs ──
import statistics
for name, grp in [("Community", group_a), ("Regional", group_b), ("Academic", group_c)]:
    med = statistics.median(grp)
    q1 = sorted(grp)[len(grp)//4]
    q3 = sorted(grp)[3*len(grp)//4]
    print(f"{name}: median = {med}, IQR = [{q1}, {q3}]")

# ── 4. Post-hoc Dunn's test (only run after significant omnibus result) ──
# scikit-posthocs returns a symmetric p-value matrix
try:
    import scikit_posthocs as sp
    import pandas as pd
    data = group_a + group_b + group_c
    labels = ["A"] * 3 + ["B"] * 3 + ["C"] * 3
    df = pd.DataFrame({"los": data, "hospital": labels})
    # Bonferroni adjustment (p_adjust = "bonferroni" or "holm")
    dunn_bonf = sp.posthoc_dunn(df, val_col="los", group_col="hospital",
                                p_adjust="bonferroni")
    print("\nDunn post-hoc p-values (Bonferroni-adjusted):")
    print(dunn_bonf.to_string())
    dunn_holm = sp.posthoc_dunn(df, val_col="los", group_col="hospital",
                                p_adjust="holm")
    print("\nDunn post-hoc p-values (Holm-adjusted):")
    print(dunn_holm.to_string())
    print("\nNote: Apply post-hoc tests only after a significant omnibus result.")
    print("Never run unadjusted pairwise Mann-Whitney tests as a post-hoc substitute.")
except ImportError:
    print("\nscikit-posthocs not installed — run: pip install scikit-posthocs")
    print("Alternative: use scipy.stats.mannwhitneyu with Bonferroni correction,")
    print("BUT use Dunn's test (not Mann-Whitney) for correct family-wise error control.")

# ── 5. Manual H verification (matches the worked example exactly) ──
N_all = 9
n_j = 3
R = [6, 15, 24]   # rank sums computed in the worked example
H_manual = (12 / (N_all * (N_all + 1))) * sum(r**2 / n_j for r in R) - 3 * (N_all + 1)
print(f"\nManual H = (12/90) * {sum(r**2/n_j for r in R):.1f} - 30 = {H_manual:.4f}")
# Expected: H_manual = 7.2

# ── Motivating dataset: LOS (days) across three hospitals ──
los    <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
hospital <- factor(c("Community", "Community", "Community",
                     "Regional",  "Regional",  "Regional",
                     "Academic",  "Academic",  "Academic"))

# ── 1. Kruskal-Wallis omnibus test (base R) ──
kw_res <- kruskal.test(los ~ hospital)
print(kw_res)
# Expected: H = 7.2, df = 2, p ≈ 0.0273

# ── 2. Epsilon-squared effect size ──
N <- length(los)
k <- nlevels(hospital)
H <- kw_res$statistic
eps_sq <- (H - k + 1) / (N - k)
cat(sprintf("Epsilon-squared (eta^2) = %.4f\n", eps_sq))

# ── 3. Group descriptive statistics ──
tapply(los, hospital, function(x) {
  cat(sprintf("%s: median = %.1f, IQR = [%.1f, %.1f]\n",
              deparse(substitute(x)), median(x), quantile(x, 0.25), quantile(x, 0.75)))
})
by(los, hospital, function(x) c(median = median(x), IQR = IQR(x)))

# ── 4a. Dunn post-hoc via dunn.test package ──
if (requireNamespace("dunn.test", quietly = TRUE)) {
  cat("\nDunn post-hoc (Bonferroni adjustment):\n")
  dunn.test::dunn.test(los, hospital, method = "bonferroni")
  cat("\nDunn post-hoc (Holm adjustment):\n")
  dunn.test::dunn.test(los, hospital, method = "holm")
} else {
  message("Install: install.packages('dunn.test')")
}

# ── 4b. Alternative: FSA::dunnTest (returns a tidy data frame) ──
if (requireNamespace("FSA", quietly = TRUE)) {
  cat("\nDunn post-hoc via FSA::dunnTest (Bonferroni):\n")
  print(FSA::dunnTest(los ~ hospital, data = data.frame(los, hospital),
                      method = "bonferroni"))
} else {
  message("Install: install.packages('FSA')")
}

# ── Note on incorrect practice ──
cat("\nIMPORTANT: Do NOT use repeated pairwise wilcox.test calls as a post-hoc.\n")
cat("Use dunn.test or FSA::dunnTest which use the pooled variance from the full ranking.\n")