how to get significance codes after Kruskal Wallis post hoc test

Question

Is there a way to get significance codes after a pairwise comparisons to a Kruskall wallis test? With significance codes I mean letter codes that are assigned to populations to indicate where differences are significant.

With a traditional anova such a test can be performed using HSD.test from the agricolae library but for non parametric counterparts of anova I have not been able to find anything.

A small toy example:

dv  <-  c(runif(100, 5.0, 10))
iv  <-  as.factor( c(rep("I", 10),  rep("II", 10),  rep("III", 10),  rep("IV", 10), rep("V", 10),
                    rep("VI", 10), rep("VII", 10), rep("VIII", 10), rep("IX", 10), rep("X", 10)))

df  <-  data.frame(dv, iv)

# with anova
library(agricolae)
aov.000  <-  aov(dv ~ iv,  data=df)
HSD.test(aov.000, "iv")

# after KW test: 
(kt  <-  kruskal.test(dv ~ iv,  data=df))

library(coin)
library(multcomp)
NDWD <- oneway_test(dv ~ iv, data = df,
        ytrafo = function(data) trafo(data, numeric_trafo = rank),
        xtrafo = function(data) trafo(data, factor_trafo = function(x)
            model.matrix(~x - 1) %*% t(contrMat(table(x), "Tukey"))),
        teststat = "max", distribution = approximate(B=1000))

### global p-value
print(pvalue(NDWD))

### sites (I = II) != (III = IV) at alpha = 0.01 (page 244)
print(pvalue(NDWD, method = "single-step"))

possible duplicate http://stackoverflow.com/questions/2478272/kruskal-wallis-test-with-details-on-pairwise-comparisons — Jonas Tundo, Apr 22 '14 at 11:07

thijs van den bergh · Accepted Answer · 2013-06-06T17:31:20.110

6

because it might be of use to others, the following seems to work (using the multcompView library):

library(multcompView)
mat  <-  data.frame(print(pvalue(NDWD, method = "single-step")))
(a   <-  c(mat[, 1]));  names(a)  <-  rownames(mat)
multcompLetters(a)

Alternatively the following will work:

test  <-  pairwise.wilcox.test(dv, iv, p.adj="bonferroni", exact=FALSE)
# test  <-  pairwise.wilcox.test(et.ef, s.t, p.adj="holm", exact=FALSE)

library(multcompView)
test$p.value
library(reshape)
(a <- melt(test$p.value))
a.cc  <-  na.omit(a)
a.pvals  <-  a.cc[, 3]
names(a.pvals)  <-  paste(a.cc[, 1], a.cc[, 2], sep="-")
a.pvals
multcompLetters(a.pvals)

edited Jun 06 '13 at 17:31

answered Jun 05 '13 at 09:43

thijs van den bergh

608
1
9
17

Great! This will be handy for me someday. – Lucas Fortini Jun 06 '13 at 17:06
1

I get the following error `Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) : cannot coerce class ""ci"" to a data.frame` after running `mat <- data.frame(print(pvalue(NDWD, method = "single-step")))` – toto_tico Nov 13 '15 at 04:18

score 3 · Answer 2 · answered Dec 17 '18 at 18:56

You can also use the cldList function from the rcompanion package (see https://rcompanion.org/rcompanion/d_06.html). Example:

k_test <- k_test$res

library(rcompanion)

cldList(comparison = k_test$Comparison,
    p.value    = PT$P.adj,
    threshold  = 0.05)


Error: No significant differences.

I used it in combination with the Dunn post-hoc and it worked perfectly.

score 2 · Answer 3 · answered Jun 04 '13 at 17:42

You can at least do it graphically using the multicomp package:

dv  <-  c(runif(100, 5.0, 10))
iv  <-  as.factor( c(rep("I", 10),  rep("II", 10),  rep("III", 10),  rep("IV", 10), rep("V", 10),
                rep("VI", 10), rep("VII", 10), rep("VIII", 10), rep("IX", 10), rep("X", 10)))
df  <-  data.frame(dv, iv)
anova_results  <-  aov(dv ~ iv,  data=df)
library(multcomp)
tuk <- glht(anova_results, linfct = mcp(iv = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))
plot(tuk.cld)
par(opar)

Of course given your randomly generated data, the resulting plot is not very interesting, but will give you the groupings-

enter image description here

This is one of my plots, using the same approach:

enter image description here Lastly, if you do not want the graphics, you can dig into the package and easily find the string that stores the grouping information to be used elsewhere.

I dont want to use an anova though. The example can indeed be analysed with anaova (everything is normally distributed and variances are equal) but for the real data I am working on I need to use a non parametric test. — thijs van den bergh, Jun 05 '13 at 09:34
Yes, I am sorry to have not caught the fact you wanted a non parametric test (it has been a while since I did KW tests). As of now, I do not know a way to do a similar plot on non-parametric contrast tests. However, as you may be aware, especially for simple designs such as a one-way comparison, you may try to perform the anova and contrast tests on ranks. If anything, you can compare your initial KW tests with anova on ranks to see how different the results are, and if not, apply the anova on ranks to be able to perform the post-hoc contrast tests. — Lucas Fortini, Jun 05 '13 at 18:05

score 1 · Answer 4 · answered Sep 22 '16 at 10:11

If you want the compact letter display for the Kruskal test, the same library agricolae seems to allow it with the function kruskal. Using your own data:

library(agricolae)
kruskal(df$dv, df$iv, group=TRUE, p.adj="bonferroni")$groups
####    trt  means M
#### 1  VI    59.2 a
#### 2  VII   57.0 a
#### 3  IX    56.4 a
#### 4  II    55.0 a
#### ...

(well, in this example the groups are not considered different...)

how to get significance codes after Kruskal Wallis post hoc test

4 Answers4