![[banner]](images/banner.png "Batsto Lake, Wharton State Forest, New Jersey")
Examples in Summary and Analysis of Extension Program Evaluation
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
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
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")
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.