[R] How to get significance codes after Kruskal Wallis test
David L Carlson
dcarlson at tamu.edu
Mon Sep 28 18:28:32 CEST 2015
It may be, but the documentation for kruskal() seems to have fallen behind updates to the function. The kruskal() manual page indicates values returned by the function that have no connection to what is actually returned and M is not described at all.
There is another version of the test, DunnTest() in the package DescTools. It appears to be less conservative than kruskalmc(). DunnTest() finds 6 significant differences between the 10 pairs at .05 to the 4 identified by kruskalmc().
> library(DescTools)
> DunnTest(Heliocide~Treatment, dta)
Dunn's test of multiple comparisons using rank sums : holm
mean.rank.diff pval
1_7d-1_2d 6.045455 0.5618
3_2d-1_2d 3.833333 0.5648
9_2d-1_2d 21.500000 0.0083 **
C-1_2d -17.045455 0.0342 *
3_2d-1_7d -2.212121 0.5648
9_2d-1_7d 15.454545 0.0602 .
C-1_7d -23.090909 0.0040 **
9_2d-3_2d 17.666667 0.0342 *
C-3_2d -20.878788 0.0083 **
C-9_2d -38.545455 3.2e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
David C
From: Michael Eisenring [mailto:Michael.Eisenring at gmx.ch]
Sent: Monday, September 28, 2015 10:40 AM
To: David L Carlson
Cc: David Winsemius; 'r-help'
Subject: Aw: RE: [R] How to get significance codes after Kruskal Wallis test
Dear David,
Thanks for your answer.
Another member of the R list pointed out that one can actually use the function
kk<-kruskal(dta$Heliocide,dta$Treatment,group=TRUE,p.adj="bonferroni")
R gives then in the Results a section$groups
$groups
trt means M
1 9_2d 47.500000 a
2 1_7d 32.045455 b
3 3_2d 29.833333 b
4 1_2d 26.000000 b
5 C 8.954545 c
I guess the codes under $groups in Column M are the significance codes I am looking for.
Thanks,
Mike
Gesendet: Montag, 28. September 2015 um 07:28 Uhr
Von: "David L Carlson" <dcarlson at tamu.edu>
An: "David Winsemius" <dwinsemius at comcast.net>, "Michael Eisenring" <michael.eisenring at gmx.ch>
Cc: 'r-help' <r-help at r-project.org>
Betreff: RE: [R] How to get significance codes after Kruskal Wallis test
The kruskalmc() function in package pgirmess performs a multiple comparisons analysis using the kruskal-wallis test. It indicates which pairs are significantly different, but it does not summarize the results in a compact letter display.
> library(pgirmess)
> kruskalmc(Heliocide~Treatment, dta)
Multiple comparison test after Kruskal-Wallis
p.value: 0.05
Comparisons
obs.dif critical.dif difference
1_2d-1_7d 6.045455 19.11021 FALSE
1_2d-3_2d 3.833333 18.69015 FALSE
1_2d-9_2d 21.500000 19.60240 TRUE
1_2d-C 17.045455 19.11021 FALSE
1_7d-3_2d 2.212121 19.11021 FALSE
1_7d-9_2d 15.454545 20.00331 FALSE
1_7d-C 23.090909 19.52123 TRUE
3_2d-9_2d 17.666667 19.60240 FALSE
3_2d-C 20.878788 19.11021 TRUE
9_2d-C 38.545455 20.00331 TRUE
-------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77840-4352
-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of David Winsemius
Sent: Saturday, September 26, 2015 11:07 AM
To: Michael Eisenring
Cc: 'r-help'
Subject: Re: [R] How to get significance codes after Kruskal Wallis test
On Sep 26, 2015, at 6:48 AM, Michael Eisenring wrote:
> Thank you very much Kristina,
>
> Unfortunately that's not what I am looking for.
>
> I am just very surprised if there would be no possibility to get the significance codes for Kruskal Wallis (I would have suggested that this is a pretty common test.)
Well, it is modestly common test but its really a global test, and not a pairwise one.
> I found another option called kruskal() which does pairwise comparison, but without significance codes.
>
> Maybe another R-list member knows more.
>
I'm not sure I "know more", and it's very possible I "know less". In particular I don't really know what the term "significance code" actually means. (I'm hoping it's not a request for "significance stars", a feature which is roundly deprecated by more knowledgeable R users.)
However, looking at the agricolae::kruskal function's help page and submitting the data in the non-formulaic manner it expects, I get this output very similar in form the the HSD.test output that it appeared you considered satisfied satisfactory:
(an.dta3<-kruskal(dta$Heliocide, dta$Treatment))
#--------------------
$statistics
Chisq p.chisq
30.25246 4.348055e-06
$parameters
Df ntr t.value
4 5 2.007584
$means
dta$Heliocide std r Min Max
1_2d 1992.7707 1747.1879 12 334.53973 4929.372
1_7d 2368.8057 1187.9285 11 767.22881 4624.945
3_2d 2640.1286 2659.5800 12 615.91181 8559.142
9_2d 5338.6711 1579.4428 10 3328.89713 8014.897
C 397.9086 443.6019 11 75.73956 1588.431
$rankMeans
dta$Treatment dta$Heliocide r
1 1_2d 26.000000 12
2 1_7d 32.045455 11
3 3_2d 29.833333 12
4 9_2d 47.500000 10
5 C 8.954545 11
$comparison
NULL
$groups
trt means M
1 9_2d 47.500000 a
2 1_7d 32.045455 b
3 3_2d 29.833333 b
4 1_2d 26.000000 b
5 C 8.954545 c
From context I am guessing that the "significance codes" you ask for are the items in the M column of the "groups" element of the list output.
--
David.
>
>
> Thank you,
>
> Mike
>
>
>
> Von: Kristina Wolf [mailto:kmwolf at ucdavis.edu]
> Gesendet: Freitag, 25. September 2015 23:26
> An: Michael Eisenring <michael.eisenring at gmx.ch>
> Cc: r-help <r-help at r-project.org>
> Betreff: Re: [R] How to get significance codes after Kruskal Wallis test
>
>
>
> Perhaps look into the function friedman.test.with.post.hoc()
>
> There is more information here: http://www.r-statistics.com/wp-content/uploads/2010/02/Friedman-Test-with-Post-Hoc.r.txt
>
>
>
> Note, this does not handle NA's though, and technically it is for blocked designs, but maybe it will lead you somewhere useful or could be adapted?
>
>
> ~ Kristina
>
> Kristina Wolf
> Ph.D. Candidate, Graduate Group in Ecology
> M.S. Soil Science
>
>
>
> On Fri, Sep 25, 2015 at 10:01 PM, Michael Eisenring <michael.eisenring at gmx.ch <mailto:michael.eisenring at gmx.ch> > wrote:
>
> Is there a way to get significance codes after a pairwise comparisons to a
> Kruskall wallis test? With significance codes I mean letter codes (a, b,c)
> that are assigned to treatments 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.
>
> Can anyone help me?
>
> Thanks mike
>
>
>
> I added two example codes.
>
> First code represents an ANOVA and a HSD.test() giving me significant codes
>
> #FIRST CODE USING ANOVA
>
> library(agricolae)
> an.dta<-aov(Gossypol~Treatment,data=dta)
> summary(an.dta)
>
> HSD.test(an.dta,"Treatment")
> # The level by alpha default is 0.05.
> outT<-HSD.test(an.dta,"Treatment", group=T)
> outT
>
> #I receive significant codes.
>
>
> #SECOND CODE USING KRUSKAL WALLIs
>
> library(agricolae)
> an.dta2<-kruskal.test(Heliocide~Treatment,dta)
> summary(an.dta2)
>
> HSD.test(an.dta2,"Treatment")
>
> #ERROR MESSAGE no significance codes, why??
>
>
>
> #DATA FOR CODES
>
>
> structure(list(Treatment = structure(c(1L, 3L, 4L, 2L, 1L, 3L,
> 4L, 2L, 5L, 1L, 3L, 2L, 5L, 1L, 3L, 4L, 2L, 5L, 1L, 3L, 4L, 2L,
> 5L, 1L, 3L, 4L, 2L, 5L, 1L, 3L, 4L, 2L, 5L, 1L, 3L, 4L, 2L, 5L,
> 1L, 3L, 4L, 2L, 5L, 1L, 3L, 4L, 2L, 5L, 1L, 3L, 4L, 2L, 5L, 1L,
> 3L, 5L), .Label = c("1_2d", "1_7d", "3_2d", "9_2d", "C"), class = "factor"),
>
> Code = structure(c(1L, 2L, 3L, 4L, 18L, 19L, 20L, 21L, 22L,
> 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L,
> 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L,
> 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 5L, 6L,
> 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L), .Label =
> c("1_2d_1c",
> "1_2d_3c", "1_2d_9c", "1_7d_1c", "10_2d_1c", "10_2d_3c",
> "10_2d_9c", "10_7d_1c", "10_C", "11_2d_1c", "11_2d_3c", "11_2d_9c",
> "11_7d_1c", "11_C", "12_2d_1c", "12_2d_3c", "12_C", "2_2d_1c",
> "2_2d_3c", "2_2d_9c", "2_7d_1c", "2_C", "3_2d_1c", "3_2d_3c",
> "3_7d_1c", "3_C", "4_2d_1c", "4_2d_3c", "4_2d_9c", "4_7d_1c",
> "4_C", "5_2d_1c", "5_2d_3c", "5_2d_9c", "5_7d_1c", "5_C",
> "6_2d_1c", "6_2d_3c", "6_2d_9c", "6_7d_1c", "6_C", "7_2d_1c",
> "7_2d_3c", "7_2d_9c", "7_7d_1c", "7_C", "8_2d_1c", "8_2d_3c",
> "8_2d_9c", "8_7d_1c", "8_C", "9_2d_1c", "9_2d_3c", "9_2d_9c",
> "9_7d_1c", "9_C"), class = "factor"), Glands = c(165, 289.3333333,
> 319.3333333, 472, 334.6666667, 259, 373.3333333, 525.6666667,
> 275.3333333, 230.6666667, 346.3333333, 377.6666667, 255.3333333,
> 217.6666667, 266, 300.3333333, 354.3333333, 225.3333333,
> 294, 359, 359, 222.6666667, 103, 246.6666667, 324.6666667,
> 277, 460, 163.6666667, 226.3333333, 228, 357.6666667, 505,
> 142.6666667, 324, 278.6666667, 317.3333333, 335.6666667,
> 193.6666667, 188, 255, 252, 393.3333333, 248.3333333, 353,
> 320.6666667, 228.3333333, 497, 165.6666667, 209.3333333,
> 162.3333333, 280, 337, 169.6666667, 231.6666667, 257.6666667,
> 218.6666667), Tannin = c(0.334252451, 1.376077586, 0.896849593,
> 0.888621795, 0.464285714, 0.830236486, 0.870881783, 0.768489583,
> 0.647727273, 0.81372549, 0.51380814, 0.859923246, 0.495265152,
> 0.699932796, 1.09375, 0.785037879, 0.892650463, 0.518963675,
> 1.05859375, 0.447916667, 1.269097222, 1.147522523, 0.391276042,
> 0.883400538, 1.523989899, 0.907930108, 0.749155405, 0.450126263,
> 0.562239583, 0.911151961, 0.611111111, 1.610677083, 0.446428571,
> 0.601151316, 1.073635057, 1.359923246, 1.00154321, 0.90933642,
> 0.012054398, 1.102083333, 1.017361111, 1.052372685, 0.958607456,
> 1.224702381, 0.982291667, 1.045138889, 1.611607143, 0.662574405,
> 1.385416667, 0.464518229, 0.994444444, 1.239583333, 0.877514368,
> 0.74453125, 0.804315476, 1.024066092), H.polone = c(6754.067177,
> 22380.26652, 23622.79158, 23733.77678, 13099.20833, 23564.74907,
> 2725.016387, 18751.03986, 4283.098494, 23008.35336, 10205.56354,
> 19787.63361, 4302.050374, 7400.640798, 22442.86044, 34315.09631,
> 16498.66728, 14170.13252, 9509.1073, 6265.29637, 20671.56905,
> 14517.15648, 2643.950729, 4974.607571, 14782.87029, 13918.82361,
> 12526.27863, 1236.908141, 4854.469195, 4076.396504, 9603.950212,
> 13762.57476, 2298.727719, 3514.186757, 5705.140289, 14178.21668,
> 14277.39878, 2656.552509, 8184.633961, 9931.163373, 21474.90732,
> 18522.74376, 9884.406532, 17242.54114, 8431.506608, 14601.11606,
> 15748.4912, 2849.90903, 16747.27644, 9396.645481, 21996.95822,
> 5767.358748, 5767.358748, 14207.1734, 10353.21833, 2859.51171
> ), Gossypol = c(1036.331811, 4171.427741, 6039.995102, 5909.068158,
> 4140.242559, 4854.985845, 6982.035521, 6132.876396, 948.2418407,
> 3618.448997, 3130.376482, 5113.942098, 1180.171957, 1500.863038,
> 4576.787021, 5629.979049, 3378.151945, 3589.187889, 2508.417927,
> 1989.576826, 5972.926124, 2867.610671, 450.7205451, 1120.955,
> 3470.09352, 3575.043632, 2952.931863, 349.0864019, 1013.807628,
> 910.8879471, 3743.331903, 3350.203452, 592.3403778, 1517.045807,
> 1504.491931, 3736.144027, 2818.419785, 723.885643, 1782.864308,
> 1414.161257, 3723.629772, 3747.076592, 2005.919344, 4198.569251,
> 2228.522959, 3322.115942, 4274.324792, 720.9785449, 2874.651764,
> 2287.228752, 5654.858696, 1247.806111, 1247.806111, 2547.326207,
> 2608.716056, 1079.846532), Heliocide = c(711.1776124, 8559.141828,
> 8014.897387, 3972.305107, 3227.467943, 5778.242027, 3628.427557,
> 3177.426984, 325.1764586, 3774.732152, 3111.880146, 4624.945228,
> 160.8912744, 336.4018128, 5207.091788, 6360.856306, 1740.091298,
> 1588.430761, 3509.141442, 685.6917982, 4664.118976, 1477.26149,
> 75.73956465, 402.1570283, 3703.317553, 4235.211434, 1730.465296,
> 91.53557346, 334.5397274, 698.1713846, 3328.897126, 1742.69355,
> 231.9097243, 513.7933372, 774.6461158, 4687.003829, 1692.296924,
> 179.1968506, 1022.628651, 1199.898583, 6132.303567, 1971.798098,
> 413.3375988, 4072.908467, 615.911814, 4906.642605, 3160.349616,
> 117.642134, 4929.371855, 616.8755006, 7428.352411, 767.2288107,
> 767.2288107, 1078.928494, 730.6740868, 425.9053258), Damage..cm. =
> c(0.4955,
> 1.516, 4.409, 3.2665, 0.491, 2.3035, 3.51, 1.8115, 0, 0.4435,
> 1.573, 1.8595, 0, 0.142, 2.171, 4.023, 4.9835, 0, 0.6925,
> 1.989, 5.683, 3.547, 0, 0.756, 2.129, 9.437, 3.211, 0, 0.578,
> 2.966, 4.7245, 1.8185, 0, 1.0475, 1.62, 5.568, 9.7455, 0,
> 0.8295, 2.411, 7.272, 4.516, 0, 0.4035, 2.974, 8.043, 4.809,
> 0, 0.6965, 1.313, 5.681, 3.474, 0, 0.5895, 2.559, 0)), .Names =
> c("Treatment",
> "Code", "Glands", "Tannin", "H.polone", "Gossypol", "Heliocide",
> "Damage..cm."), class = "data.frame", row.names = c(NA, -56L))
>
>
> [[alternative HTML version deleted]]
>
David Winsemius
Alameda, CA, USA
______________________________________________
R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list