I'm about to apply Kruskal-Wallis test (non-parametric ANOVA) on three groups of unequal length. I was taught/advised to apply Krukal-Wallis only if:

  • dependent variable is at least at ordinal level of measurement
  • group's $ n > 5 $ (otherwise H statistic is not $ \chi ^2 $ distributed, so exact p-value cannot be calculated, etc.)
  • Wikipedia page says that it requires an identically-shaped and scaled distribution for each group

Apart from that, what are general considerations when applying Kruskal-Wallis test? Should I check homoscedascity with Fligner-Killeen's test, and what should I do in case of small $ n $, or unequal group sizes?

There are exactRankTests::wilcox.exact and PASWR::wilcoxE.test functions in R, but I can't seem to find an analogous one for Kruskal-Wallis...

Advice, anyone?

  • 1
    $\begingroup$ A permutation test might be appropriate here. $\endgroup$ – chl Mar 27 '11 at 19:52
  • $\begingroup$ I added a little something: "unequal groups". $\endgroup$ – aL3xa Mar 27 '11 at 21:12
  • $\begingroup$ Can you tell us what your three $n$s are? $\endgroup$ – onestop Mar 27 '11 at 22:09
  • $\begingroup$ Well, if I do that, then this question can easily be marked as an assignment/homework. =) But, OK, say: 10, 50, 100. $\endgroup$ – aL3xa Mar 27 '11 at 22:37

With small, and possibly unequal group sizes, I'd go with chl's and onestop's suggestion and do a Monte-Carlo permutation test. For the permutation test to be valid, you need exchangeability under $H_{0}$. If all distributions have the same shape (and are therefore identical under $H_{0}$), this is true.

Here's a first try at looking at the case of 3 groups and no ties. First, let's compare the asymptotic $\chi^{2}$ distribution function against a MC-permutation one for given group sizes (this implementation will break for larger group sizes).

P  <- 3                     # number of groups
Nj <- c(4, 8, 6)            # group sizes
N  <- sum(Nj)               # total number of subjects
IV <- factor(rep(1:P, Nj))  # grouping factor
alpha <- 0.05               # alpha-level

# there are N! permutations of ranks within the total sample, but we only want 5000
nPerms <- min(factorial(N), 5000)

# random sample of all N! permutations
# sample(1:factorial(N), nPerms) doesn't work for N! >= .Machine$integer.max
permIdx <- unique(round(runif(nPerms) * (factorial(N)-1)))
nPerms  <- length(permIdx)
H       <- numeric(nPerms)  # vector to later contain the test statistics

# function to calculate test statistic from a given rank permutation
getH <- function(ranks) {
    Rj <- tapply(ranks, IV, sum)
    (12 / (N*(N+1))) * sum((1/Nj) * (Rj-(Nj*(N+1) / 2))^2)

# all test statistics for the random sample of rank permutations (breaks for larger N)
# numperm() internally orders all N! permutations and returns the one with a desired index
library(sna)                # for numperm()
for(i in seq(along=permIdx)) { H[i] <- getH(numperm(N, permIdx[i]-1)) }

# cumulative relative frequencies of test statistic from random permutations
pKWH   <- cumsum(table(round(H, 4)) / nPerms)
qPerm  <- quantile(H, probs=1-alpha)  # critical value for level alpha from permutations
qAsymp <- qchisq(1-alpha, P-1)        # critical value for level alpha from chi^2

# illustration of cumRelFreq vs. chi^2 distribution function and resp. critical values
plot(names(pKWH), pKWH, main="Kruskal-Wallis: permutation vs. asymptotic",
     type="n", xlab="h", ylab="P(H <= h)", cex.lab=1.4)
points(names(pKWH), pKWH, pch=16, col="red")
curve(pchisq(x, P-1), lwd=2, n=200, add=TRUE)
abline(h=0.95, col="blue")                         # level alpha
abline(v=c(qPerm, qAsymp), col=c("red", "black"))  # critical values
legend(x="bottomright", legend=c("permutation", "asymptotic"),
       pch=c(16, NA), col=c("red", "black"), lty=c(NA, 1), lwd=c(NA, 2))

enter image description here

Now for an actual MC-permutation test. This compares the asymptotic $\chi^{2}$-derived p-value with the result from coin's oneway_test() and the cumulative relative frequency distribution from the MC-permutation sample above.

> DV1 <- round(rnorm(Nj[1], 100, 15), 2)  # data group 1
> DV2 <- round(rnorm(Nj[2], 110, 15), 2)  # data group 2
> DV3 <- round(rnorm(Nj[3], 120, 15), 2)  # data group 3
> DV  <- c(DV1, DV2, DV3)                 # all data
> kruskal.test(DV ~ IV)                   # asymptotic p-value
Kruskal-Wallis rank sum test
data:  DV by IV 
Kruskal-Wallis chi-squared = 7.6506, df = 2, p-value = 0.02181

> library(coin)                           # for oneway_test()
> oneway_test(DV ~ IV, distribution=approximate(B=9999))
Approximative K-Sample Permutation Test
data:  DV by IV (1, 2, 3) 
maxT = 2.5463, p-value = 0.0191

> Hobs <- getH(rank(DV))                  # observed test statistic

# proportion of test statistics at least as extreme as observed one (+1)
> (pPerm <- (sum(H >= Hobs) + 1) / (length(H) + 1))
[1] 0.0139972
  • $\begingroup$ Thank you so much for this thorough post (and accent on R). $\endgroup$ – aL3xa Mar 28 '11 at 12:38
  • You need not check homoscedasticity. Kruskal and Wallis stated in their original paper that the “test may be fairly insensitive to differences in variability”.
  • If there is no exact test available, you can use bootstrap.
  • $\begingroup$ Yeah, right... bootstrap is always an option. Thanks! $\endgroup$ – aL3xa Mar 27 '11 at 19:41
  • 1
    $\begingroup$ @GaBorgulya About your first point, it seems to me that non-parametric tests make the assumption that there's a shift in central tendency, but that overall the distributions that are compared have the same shape. $\endgroup$ – chl Mar 27 '11 at 19:54
  • $\begingroup$ So what should I use to test basic assumptions? Does that mean that I need to apply a homoscedascity test? $\endgroup$ – aL3xa Mar 27 '11 at 20:05
  • 2
    $\begingroup$ Unfortunatey boostrapping doesn't work well with very small samples either. @chl's suggestion of a permutation test is a much better idea. Monte Carlo permutation tests are a good option when exact tests haven't been programmed or would take unfeasibly long. $\endgroup$ – onestop Mar 27 '11 at 22:07
  • $\begingroup$ @chl The $H_0$ of KW is that the multiple samples come from the same distribution of any shape. Looking at the test statistic I agree that the test is more powerful in case of shifts only. But although the original paper(permalink) also states that the $H_A$ "is that the samples come from populations of approximately the same form, but shifted", I don't see where the test statistic relies on that assumption. $\endgroup$ – GaBorgulya Mar 27 '11 at 22:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.