[banner]

An R Companion for the Handbook of Biological Statistics

Salvatore S. Mangiafico

Kruskal–Wallis Test

Examples in Summary and Analysis of Extension Program Evaluation


SAEEPER: Kruskal–Wallis Test

 

Packages used in this chapter

The following commands will install these packages if they are not already installed:


if(!require(dplyr)){install.packages("dplyr")}
if(!require(FSA)){install.packages("FSA")}
if(!require(DescTools)){install.packages("DescTools")}
if(!require(rcompanion)){install.packages("rcompanion")}
if(!require(multcompView)){install.packages("multcompView")}

When to use it

See the Handbook for information on this topic.

 

Null hypothesis

This example shows just summary statistics, histograms by group, and the Kruskal–Wallis test.  An example with plots, post-hoc tests, and alternative tests is shown in the “Example” section below.

 

Kruskal–Wallis test example


### --------------------------------------------------------------
### Kruskal–Wallis test, hypothetical example, p. 159
### --------------------------------------------------------------

Input =("
Group      Value
Group.1      1
Group.1      2
Group.1      3
Group.1      4
Group.1      5
Group.1      6
Group.1      7
Group.1      8
Group.1      9
Group.1     46
Group.1     47
Group.1     48
Group.1     49
Group.1     50
Group.1     51
Group.1     52
Group.1     53
Group.1    342
Group.2     10
Group.2     11
Group.2     12
Group.2     13
Group.2     14
Group.2     15
Group.2     16
Group.2     17
Group.2     18
Group.2     37
Group.2     58
Group.2     59
Group.2     60
Group.2     61
Group.2     62
Group.2     63
Group.2     64
Group.2    193
Group.3     19
Group.3     20
Group.3     21
Group.3     22
Group.3     23
Group.3     24
Group.3     25
Group.3     26
Group.3     27
Group.3     28
Group.3     65
Group.3     66
Group.3     67
Group.3     68
Group.3     69
Group.3     70
Group.3     71
Group.3     72
")

Data = read.table(textConnection(Input),header=TRUE)


### Specify the order of factor levels

library(dplyr)

Data =
mutate(Data,
       Group = factor(Group, levels=unique(Group)))


Medians and descriptive statistics

As noted in the Handbook, each group has identical medians and means.


library(FSA)

Summarize(Value ~ Group,
          data = Data)


    Group  n mean       sd min    Q1 median    Q3 max
1 Group.1 18 43.5 77.77513   1  5.25   27.5 49.75 342
2 Group.2 18 43.5 43.69446  10 14.25   27.5 60.75 193
3 Group.3 18 43.5 23.16755  19 23.25   27.5 67.75  72


Histograms for each group


library(lattice)

histogram(~ Value | Group,
          data=Data,
          layout=c(1,3))      #  columns and rows of individual plots
         

Rplot

Kruskal–Wallis test

In this case, there is a significant difference in the distributions of values among groups, as is evident both from the histograms and from the significant Kruskal–Wallis test.  Only in cases where the distributions in each group are similar can a significant Kruskal–Wallis test be interpreted as a difference in medians.


kruskal.test(Value ~ Group,
             data = Data)


Kruskal-Wallis chi-squared = 7.3553, df = 2, p-value = 0.02528


#     #     #


How the test works

Assumptions

See the Handbook for information on these topics.

 

Example

The Kruskal–Wallis test is performed on a data frame with the kruskal.test function in the native stats package.  Shown first is a complete example with plots, post-hoc tests, and alternative methods, for the example used in R help.  It is data measuring if the mucociliary efficiency in the rate of dust removal is different among normal subjects, subjects with obstructive airway disease, and subjects with asbestosis.  For the original citation, use the ?kruskal.test command.  For both the submissive dog example and the oyster DNA example from the Handbook, a Kruskal–Wallis test is shown later in this chapter.

 

Kruskal–Wallis test example


### --------------------------------------------------------------
### Kruskal–Wallis test, asbestosis example from R help for
###   kruskal.test
### --------------------------------------------------------------

Input =("
Obs Health     Efficiency
1   Normal     2.9
2   Normal     3.0
3   Normal     2.5
4   Normal     2.6
5   Normal     3.2
6   OAD        3.8
7   OAD        2.7
8   OAD        4.0
9   OAD        2.4
10  Asbestosis 2.8
11  Asbestosis 3.4
12  Asbestosis 3.7
13  Asbestosis 2.2
14  Asbestosis 2.0
")

Data = read.table(textConnection(Input),header=TRUE)


### Specify the order of factor levels


library(dplyr)
         
Data =
mutate(Data,
       Health = factor(Health, levels=unique(Health)))


Medians and descriptive statistics


library(FSA)

Summarize(Efficiency ~ Health,
          data = Data)


      Health n  mean        sd min    Q1 median   Q3 max
1     Normal 5 2.840 0.2880972 2.5 2.600   2.90 3.00 3.2
2        OAD 4 3.225 0.7932003 2.4 2.625   3.25 3.85 4.0
3 Asbestosis 5 2.820 0.7362065 2.0 2.200   2.80 3.40 3.7


Graphing the results

 

Stacked histograms of values across groups

 

library(lattice)

histogram(~ Efficiency | Health,
          data=Data,
          layout=c(1,3))      #  columns and rows of individual plots
         

Rplot

Stacked histograms for each group in a Kruskal–Wallis test.  If the distributions are similar, then the Kruskal–Wallis test will test for a difference in medians.

 

 

Simple boxplots of values across groups

 

boxplot(Efficiency ~ Health,
        data = Data,
        ylab="Efficiency",
        xlab="Health")

 

Rplot

 

Kruskal–Wallis test

 

kruskal.test(Efficiency ~ Health,
             data = Data)

 

Kruskal-Wallis chi-squared = 0.7714, df = 2, p-value = 0.68

Dunn test for multiple comparisons

If the Kruskal–Wallis test is significant, a post-hoc analysis can be performed to determine which levels of the independent variable differ from each other level.  Probably the most popular test for this is the Dunn test, which is performed with the dunnTest function in the FSA package.  Adjustments to the p-values could be made using the method option to control the familywise error rate or to control the false discovery rate.  See ?p.adjust for details.

Zar (2010) states that the Dunn test is appropriate for groups with unequal numbers of observations.

If there are several values to compare, it can be beneficial to have R convert this table to a compact letter display for you.  The cldList function in the rcompanion package can do this.


### Order groups by median

Data$Health = factor(Data$Health,
                      levels=c("OAD", "Normal", "Asbestosis"))

### Dunn test

library(FSA)

PT = dunnTest(Efficiency ~ Health,
              data=Data,
              method="bh")    # Can adjust p-values;
                              # See ?p.adjust for options

PT


Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the False Discovery Rate method.

           Comparison         Z   P.unadj     P.adj
1        OAD - Normal 0.6414270 0.5212453 0.7818680
2    OAD - Asbestosis 0.8552360 0.3924205 1.0000000
3 Normal - Asbestosis 0.2267787 0.8205958 0.8205958



PT = PT$res

PT

library(rcompanion)

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


Error: No significant differences.


Nemenyi test for multiple comparisons

Zar (2010) suggests that the Nemenyi test is not appropriate for groups with unequal numbers of observations.


library(DescTools)

PT = NemenyiTest(x = Data$Efficiency,
                 g = Data$Health,
                 dist="tukey")

PT

Nemenyi's test of multiple comparisons for independent samples (tukey) 

                  mean.rank.diff   pval   
OAD-Normal                   1.8 0.7972   
Asbestosis-Normal           -0.6 0.9720   
Asbestosis-OAD              -2.4 0.6686   


library(rcompanion)

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


Error: No significant differences.


Pairwise Mann–Whitney U-tests

Another post-hoc approach is to use pairwise Mann–Whitney U-tests.  To prevent the inflation of type I error rates, adjustments to the p-values can be made using the p.adjust.method option to control the familywise error rate or to control the false discovery rate. See ?p.adjust for details.

If there are several values to compare, it can be beneficial to have R convert this table to a compact letter display for you.  The multcompLetters function in the multcompView package can do this, but first the table of p-values must be converted to a full table.


PT = pairwise.wilcox.test(Data$Efficiency,
                          Data$Health,
                          p.adjust.method="none")
                              # Can adjust p-values;
                              # See ?p.adjust for options


PT



Pairwise comparisons using Wilcoxon rank sum test

           Normal OAD
OAD        0.73   -  
Asbestosis 1.00   0.41



PT = PT$p.value

library(rcompanion)

PT1 = fullPTable(PT)

PT1


              Normal       OAD Asbestosis
Normal     1.0000000 0.7301587  1.0000000
OAD        0.7301587 1.0000000  0.4126984
Asbestosis 1.0000000 0.4126984  1.0000000


library(multcompView)

multcompLetters(PT1,
                compare="<",
                threshold=0.05,
                Letters=letters,
                reversed = FALSE)


    Normal        OAD Asbestosis
       "a"        "a"        "a"

    ### Values sharing the same letter are not significantly different

#     #     #


Kruskal–Wallis test example


### --------------------------------------------------------------
### Kruskal–Wallis test, submissive dog example, pp. 161–162
### --------------------------------------------------------------

Input =("
Dog          Sex      Rank
 Merlino      Male     1
 Gastone      Male     2
 Pippo        Male     3
 Leon         Male     4
 Golia        Male     5
 Lancillotto  Male     6
 Mamy         Female   7
 Nanà         Female   8
 Isotta       Female   9
 Diana        Female  10
 Simba        Male    11
 Pongo        Male    12
 Semola       Male    13
 Kimba        Male    14
 Morgana      Female  15
 Stella       Female  16
 Hansel       Male    17
 Cucciola     Male    18
 Mammolo      Male    19
 Dotto        Male    20
 Gongolo      Male    21
 Gretel       Female  22
 Brontolo     Female  23
 Eolo         Female  24
 Mag          Female  25
 Emy          Female  26
 Pisola       Female  27
 ")

Data = read.table(textConnection(Input),header=TRUE)

kruskal.test(Rank ~ Sex,
             data = Data)


Kruskal-Wallis chi-squared = 4.6095, df = 1, p-value = 0.03179

#     #     #


Graphing the results

Graphing of the results is shown above in the “Example” section.

 

Similar tests

One-way anova is presented elsewhere in this book.

 

How to do the test

Kruskal–Wallis test example


### --------------------------------------------------------------
### Kruskal–Wallis test, oyster DNA example, pp. 163–164
### --------------------------------------------------------------

Input =("
 Markername  Markertype  fst
 CVB1        DNA        -0.005
 CVB2m       DNA         0.116
 CVJ5        DNA        -0.006
 CVJ6        DNA         0.095
 CVL1        DNA         0.053
 CVL3        DNA         0.003
 6Pgd        protein    -0.005
 Aat-2       protein     0.016
 Acp-3       protein     0.041
 Adk-1       protein     0.016
 Ap-1        protein     0.066
 Est-1       protein     0.163
 Est-3       protein     0.004
 Lap-1       protein     0.049
 Lap-2       protein     0.006
 Mpi-2       protein     0.058
 Pgi         protein    -0.002
 Pgm-1       protein     0.015
 Pgm-2       protein     0.044
 Sdh         protein     0.024
")

Data = read.table(textConnection(Input),header=TRUE)


kruskal.test(fst ~ Markertype,
             data = Data)


Kruskal-Wallis chi-squared = 0.0426, df = 1, p-value = 0.8365

#     #     #


Power Analysis

See the Handbook for information on this topic.

 

References

Zar, J.H. 2010. Biostatistical Analysis, 5th ed.  Pearson Prentice Hall: Upper Saddle River, NJ.