10
$\begingroup$

I have box plots of 13 groups that I show in one plot. The groups have unbalanced populations and are not normally distributed. I want to show which pairs are statistically similar (i.e., have kruskal.test p-value < 0.05) by putting a,b,c, etc. on top of boxes that match. Here is a pseudo code to show what I have:

A = c(1, 5, 8, 17, 16, 3, 24, 19, 6) 
B = c(2, 16, 5, 7, 4, 7, 3) 
C = c(1, 1, 3, 7, 9, 6, 10, 13) 
D = c(2, 15, 2, 9, 7) 
junk = list(g1=A, g2=B, g3=C, g4=D) 
boxplot(junk)

Here is an plot I found that does what I want (except I have 13 groups in one row):

Improve this question
$\endgroup$

1 Answer 1

Reset to default
8
$\begingroup$

The simplest code that comes to my mind is shown below. I'm pretty certain there's some already existing function(s) to do that on CRAN but I'm too lazy to search for them, even on R-seek.

dd <- data.frame(y=as.vector(unlist(junk)), 
                 g=rep(paste("g", 1:4, sep=""), unlist(lapply(junk, length))))

aov.res <- kruskal.test(y ~ g, data=dd)
alpha.level <- .05/nlevels(dd$g)  # Bonferroni correction, but use 
                                  # whatever you want using p.adjust()

# generate all pairwise comparisons
idx <- combn(nlevels(dd$g), 2)

# compute p-values from Wilcoxon test for all comparisons
pval.res <- numeric(ncol(idx))
for (i in 1:ncol(idx))
  # test all group, pairwise
  pval.res[i] <- with(dd, wilcox.test(y[as.numeric(g)==idx[1,i]], 
                                      y[as.numeric(g)==idx[2,i]]))$p.value

# which groups are significantly different (arranged by column)
signif.pairs <- idx[,which(pval.res<alpha.level)]

boxplot(y ~ g, data=dd, ylim=c(min(dd$y)-1, max(dd$y)+1))
# use offset= to increment space between labels, thanks to vectorization
for (i in 1:ncol(signif.pairs))
    text(signif.pairs[,i], max(dd$y)+1, letters[i], pos=4, offset=i*.8-1)

Here is an example of what the above code would produce (with significant differences between the four groups):

enter image description here

Instead of the Wilcoxon test, one could rely on the procedure implemented in the kruskalmc() function from the pgirmess package (see a description of the procedure used here).

Also, be sure to check Rudolf Cardinal's R tips about R: basic graphs 2 (see in particular, Another bar graph, with annotations).

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks, chl. I'm not a statistician so I'll need to learn your answer more, but it looks like it does what I need. $\endgroup$
    – Ilik
    Commented Dec 21, 2011 at 22:19
  • $\begingroup$ If there's something unclear, don't hesitate to ask. There is not much statistical stuff here: I used Wilcoxon tests for pairwise group comparisons and corrected individual p-values to limit the overall risk of false positive to 5%. The npmc package includes additional facilities for handling non-parametric multiple comparisons, but there's also the coin framework. The rest is purely R code using base R graphics; it could be done with lattice or ggplot2 as well. $\endgroup$
    – chl
    Commented Dec 22, 2011 at 0:10
  • $\begingroup$ I apologies for my ignorance in statistics. So I tried your code and I first noticed you calculate the kruskal.test but don't use the result (aov.res). I now understand the wilcox.test is a special case for kruskal for two samples. But then I tried to change the values in my groups to make them (intuitively) different and see what comes up. A = c(10, 50, 18, 17, 16, 31, 24, 19, 6) B = c(10, 50, 18, 17, 16, 30, 25, 18, 7) C = c(1, 1, 3, 7, 9, 6, 10, 13) D = c(200, 158, 206, 119, 77). g1 &g2 are now significantly different?? and g3 & g4 are not? not sure I understand why. $\endgroup$
    – Ilik
    Commented Dec 22, 2011 at 20:33
  • $\begingroup$ @Ilik Something went wrong in my code (the with(dd, instruction was in the wrong place resulting in a weird test!). Thanks for catching that! Yes, I don't use results from the KW test but it's always a good idea to check it first, otherwise multiple post-hoc tests would be meaningless (or at least, suggestive of data snooping). Note I've corrected the code and nothing turned out to be significant, but I left the original image for the sake of clarity with significant results. $\endgroup$
    – chl
    Commented Dec 22, 2011 at 21:35
  • $\begingroup$ For anyone interested, I've changed a little the code to write the results to the boxplot: MyText = rep('',nlevels(dd$g)) for (i in 1:ncol(signif.pairs)) { MyText[signif.pairs[,i]] = paste(MyText[signif.pairs[,i]],letters[i],sep='') } text(c(1:nlevels(dd$g)), rep(max(dd$y)+1,nlevels(dd$g)),MyText) $\endgroup$
    – Ilik
    Commented Dec 23, 2011 at 23:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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