[R] Multiple Comparisons-Kruskal-Wallis-Test: kruskal{agricolae} and kruskalmc{pgirmess} don't yield the same results although they should do (?)

[greatest.possible.newbie]

Hi there,

I am doing multiple comparisons for data that is not normally distributed.
For this purpose I tried both functions kruskal{agricolae} and
kruskalmc{pgirmess}. It confuses me that these functions do not yield the
same results although they are doing the same thing, don't they? Can anyone
tell my why this happens and which function I can trust?

kruskalmc() tells me that there are no differences between any of the groups
(i.e. the "difference" column of the results is filled only with FALSE).
kruskal() tells me that there are indeed differences (between group 4 and
11).

Here is my data and code:

x <- structure(c(-0.089, 0.093, -0.125, -0.253, 0.053, 0.029, 0.429,
0.139, 0.153, -0.035, 0.721, 0, -0.271, -0.014, 0.038, 0.107,
-0.064, -0.178, 0, -0.065, 0.232, 0, -0.036, 0.107, -0.018, -0.632,
0.189, -0.247, -0.982, 0.161, 0.307, 0.146, 0.21, -0.043, -0.029,
0.107, -0.043, -0.178, 0.036, 0.511, 0, 0.307, 0.389, -0.754,
0.152, -0.129, 0.65, 0.411, 0, 0.539, -0.122, 0.514, -0.15, -1.19,
-0.032, 0.103, 0, -0.522, -0.857, -0.036, 0.104, -0.357, 0.114,
-0.054, 0.04, 0.036, 0.086, 0.736, -0.097, 0.784, 0.053, 0.007,
-0.646, 0.185, 0.107, -0.115, 0, -0.036, -0.082, -0.113, -0.032,
0.354, -0.095, -0.328, -0.215, -0.036, 0, 0.357, 0, 0.108, -0.014,
0.307, -0.053, 0.318, -0.058, 0.268, -0.067, 0.071, 0.261, -0.018,
0.054, 0.086, 1.107, -0.617, 0.286, 0.072, 0.036, 0.179, -0.096,
0.143, 0.45, -0.21, 0.372, 0.061, -0.218, -0.214, 0.272, 0.108,
0.175, -0.017, 0.473, -0.575, 0.083, 0.025, 0.25, -0.843, -0.054,
-0.775, 0.036, -0.297, 0.8, 0.004, 0.189, 0.005, 0.103, 0.289,
-0.107, -0.096, 0, 0.015, -0.035, -0.125, 0.125, -0.071, -0.029,
-0.643, -0.008, 0.184, 0.303, 0, -0.164, 0.047, -0.062, -0.164,
-0.604, -0.178, 0.233, -0.154, -0.107, -0.14, -0.207, 0.211,
0.175, 0.714, 0, 0.286, -0.143, 0.018, 0.643, -0.036, 0.357,
0.071, 0.186, 0.104, -0.086, -0.611, -0.028, -0.025, 0.179, -0.032,
0.058, 0.04, -0.428, 0.447, 0.178, -0.061, 0.167, -0.071, 0.321,
0.082, 0.532, -0.22, 0.086, -0.107, 0.118, -0.139, 0.03, -0.228,
0.008, 0.178, -0.3, 0.018, -0.025, -0.329, 0.136, 0.304, 0.085,
-0.014, 0.07, 0.136, -0.218, 0.071, -0.178, 0.012, 0.229, 0.268,
-0.535, 0.164, -0.15, 0.097, 0.125, -0.536, 0.214, 0.222, -0.089,
-0.121, -0.155, -0.286, -0.282, -0.443, 0.071, -0.05, -0.04,
-0.075, -0.03, -0.357, -0.071, 0.641, 0.007, 0, 0.018, -0.573,
0.132, -0.33, -0.279, -0.639, -0.093, -0.5, -0.197, 0.303, 0.322,
-0.071, -0.071, 0.165, 0, 0, 0, 0, 0.054, -0.321, 0.093, 0.268,
-0.511, -0.3, 0.202, 0.328, -0.24, 0.871, -0.021, 0.211, 0.118,
-0.157, -0.357, 0.107, 0, 0.072, -0.357, 0.003, 0.147, 0.057,
-0.315, 0.053, 0.35, -0.107, 0.036, -0.143, 0.168, 0, 0.172,
0.321, 0.178, -0.526, -0.035, 0.247, 0.557, 0.168, 0.143, 0,
-0.432, 0.072, -0.065, 1, 0, 0.179, 0.218, 0.196, -0.122, -0.457,
0.072, 0, 0.247, -0.296, -0.118, 0.107, 0.136, 0.029, -0.058,
0.25, 0.139, -0.057, 0.15, 0.042, -0.703, -0.018, -0.318, -0.011,
-0.321, -0.6, -1.189, 0.225, -0.143, -0.112, 0.09, -0.071, -0.015,
0.828, 0.124, 0.582, 0.689, -0.107, 0.017, -0.15, 0.057, 0.143,
0.107, -0.204, -0.118, 0.021, 0.067, -0.035, -0.357, -1.015,
-0.039, 0.046, -0.036, 0.072, 0.204, 0.05, -0.038, 0, 0.057,
-0.05, 0.41, 0.143, -0.325, 0.332, -0.153, 0.157, -0.185, 0.206,
-0.086, -0.204, 0, 0.271, -0.143, 0, -0.357, 0.218, 0.036, 0.179,
-0.45, 0.072, 0.018, -0.259, 0, -0.53, -0.018, -0.054, 0.435,
0.378, 0.221, -0.921, -0.375, -0.54, 0.25, -0.16, 0, -0.007,
-0.357, 0.204, 0.129, 0, -0.45, 0, 0, 0.571, 0.392, 0, -0.465,
0.072, 0.072, -0.257, -0.007, 0.039, 0, 0.299, -0.526, 0.268,
-0.05, -0.45, 0.178, 0.025, 0.072, 0.107, 0.195, 0.089, 0.115,
0, 0.189, -0.036, -0.337, -0.45, -0.34, -0.065, 0.343, -0.071,
0.107, -0.111, 0, -0.411, -1.012, -0.108, -0.036, -0.036, 0.243,
-0.118, 0.045, 0.018, -0.607, 0.196, -0.425, -0.174, 0.068, -0.886,
-0.075, 0.143, 0.09, 0.482, -0.058, -0.2, 0.341, -0.014, 0.068,
0.107, 0.078, 0.107, -0.25, 0.382, -0.122, -0.466, -0.471, 0.046,
0.229, 0.329, 0.043, -0.207, 0.357, 0.029, 0.107, 0.286, 0, -0.15,
-0.389, -0.261, -0.135, 1.028, 0.384, 0.136, -0.147, -0.143,
0.486, -0.075, -0.036, 0.04, -0.214, 0, 0.286, 0.071, 0.25, 0.115,
1, 7, 5, 1, 1, 1, 11, 8, 3, 11, 9, 4, 9, 2, 1, 1, 4, 2, 2, 4,
2, 1, 1, 6, 4, 4, 8, 1, 4, 1, 9, 1, 11, 9, 3, 10, 1, 4, 1, 5,
2, 1, 1, 1, 1, 3, 7, 11, 1, 1, 2, 1, 1, 2, 1, 1, 8, 2, 1, 1,
1, 1, 1, 10, 3, 4, 1, 10, 2, 1, 1, 6, 5, 1, 1, 11, 1, 11, 3,
1, 1, 2, 11, 2, 1, 6, 7, 1, 9, 1, 3, 8, 1, 1, 4, 2, 5, 2, 1,
1, 3, 1, 1, 1, 10, 5, 3, 1, 1, 11, 4, 1, 1, 10, 1, 1, 6, 2, 1,
1, 5, 1, 2, 4, 9, 1, 1, 5, 1, 1, 1, 1, 5, 1, 1, 11, 5, 1, 1,
1, 1, 1, 1, 2, 1, 10, 1, 1, 9, 9, 1, 11, 1, 1, 1, 4, 2, 3, 1,
11, 1, 4, 1, 1, 1, 5, 1, 1, 4, 1, 1, 1, 1, 8, 2, 1, 1, 1, 1,
1, 7, 2, 2, 2, 2, 1, 11, 1, 8, 1, 1, 8, 2, 1, 1, 5, 1, 5, 1,
1, 2, 10, 7, 8, 1, 9, 1, 5, 6, 9, 11, 2, 1, 1, 1, 1, 3, 1, 1,
1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 11, 1, 1, 5, 5, 5, 11, 10,
3, 1, 8, 1, 1, 1, 6, 1, 1, 6, 1, 11, 1, 4, 1, 1, 8, 11, 8, 1,
3, 1, 1, 1, 1, 11, 11, 5, 9, 1, 1, 5, 5, 7, 11, 4, 3, 1, 10,
5, 8, 1, 1, 7, 1, 5, 1, 1, 1, 6, 5, 4, 1, 2, 1, 7, 1, 3, 7, 1,
8, 7, 1, 1, 2, 11, 1, 1, 1, 1, 1, 10, 1, 1, 1, 11, 3, 3, 1, 1,
9, 1, 1, 1, 11, 1, 11, 1, 1, 1, 11, 7, 1, 6, 1, 1, 1, 11, 10,
3, 4, 3, 4, 1, 1, 1, 11, 1, 8, 6, 1, 6, 5, 11, 1, 2, 1, 4, 1,
1, 1, 10, 1, 1, 6, 1, 1, 2, 1, 1, 1, 7, 2, 1, 8, 8, 1, 1, 1,
1, 5, 1, 11, 1, 1, 1, 9, 4, 1, 2, 2, 8, 1, 1, 6, 3, 1, 2, 1,
1, 1, 1, 1, 10, 9, 4, 1, 1, 6, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1,
1, 1, 2, 10, 1, 6, 2, 1, 1, 1, 1, 1, 1, 1, 10, 10, 6, 11, 1,
1, 1, 1, 11, 1, 1, 1, 1, 2, 1, 1, 1, 7, 6, 1, 1, 1, 9, 1, 7,
11, 4, 1, 1, 1, 1, 8, 1, 1, 6, 1, 2, 1, 3, 1, 11, 1, 5, 1, 6,
1, 1, 11, 6, 1, 1, 1, 2, 1, 1, 2, 6, 7, 7, 1, 8, 3, 1, 1, 1,
7, 1, 1), .Dim = c(500L, 2L))

#Now perform the kruskalmc() test.

require(pgirmess)
kruskalmc(x[,1], x[,2])

#Multiple comparison test after Kruskal-Wallis
#p.value: 0.05
#Comparisons
# obs.dif critical.dif difference
#1-2 8.4929415 82.16065 FALSE
#1-3 2.6363784 104.15270 FALSE
#1-4 22.5503277 102.13541 FALSE
#1-5 13.2511610 100.24365 FALSE
#1-6 0.8574104 108.62307 FALSE
#[..........]
#8-9 17.8228070 165.54197 FALSE
#8-10 9.5230263 162.62538 FALSE
#8-11 26.2159443 137.28174 FALSE
#9-10 27.3458333 172.25281 FALSE
#9-11 8.3931373 148.56074 FALSE
#10-11 35.7389706 145.30369 FALSE

require(agricolae)
kruskal(x[,1], x[,2])

# [ ... ]
#Means with the same letter are not significantly different
#Groups, Treatments and mean of the ranks
#a 11 304.4
#ab 9 296
#ab 7 286.6
#ab 8 278.2
#ab 10 268.7
#ab 2 250.6
#ab 6 242.9
#ab 1 242.1
#ab 3 239.4
#ab 5 228.8
#b 4 219.5



--
View this message in context: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-tp4639004.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-h...@r-project.org mailing list
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.

[peter dalgaard]

On Aug 3, 2012, at 07:42 , greatest.possible.newbie wrote:

> I am doing multiple comparisons for data that is not normally distributed.
> For this purpose I tried both functions kruskal{agricolae} and
> kruskalmc{pgirmess}. It confuses me that these functions do not yield the
> same results although they are doing the same thing, don't they? Can anyone
> tell my why this happens and which function I can trust?
>
> kruskalmc() tells me that there are no differences between any of the groups
> (i.e. the "difference" column of the results is filled only with FALSE).
> kruskal() tells me that there are indeed differences (between group 4 and
> 11).

Trust nothing if you don't understand the issues involved. There's a reason that special code is required for multiple comparisons, and a rather complicated and inexact theory for dealing with it. There are special complications with rank tests because the standard theory assumes a global null (no group differences at all).

One of the approaches is to perform the testing pairwise and adjust the p-values for multiple comparisons. Notice that the kruskal() function has a p.adj argument which defaults to "none"!

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd....@cbs.dk Priv: PDa...@gmail.com

[greatest.possible.newbie]

Thank you for your answer.
The p.adj argument in the kruskal()-function doesn't seem to change
anything... Not even the "bonferroni"-method although it is described as the
most conservative one (multiplying all p-values with the number of
comparisons). I suppose the kruskal()-function is not working properly...
On the other hand I doubt the method behind the kruskalmc()-function as this
function doesn't even turn out to detect significant differences between the
grouping variable (which is obviously a severe error).
Do you think it is justifiable to use the kruskal()-function without
p-adjustment, i.e. doing only pairwise tests like you can do with the
kruskal.test()-function although I obviously want to do multiple
comparisons?

kruskal(x[,1],x[,2],p.adj="bonferroni")
#Yields exactely the same results.

#Groups, Treatments and mean of the ranks
#a 11 304.4
#ab 9 296
#ab 7 286.6
#ab 8 278.2
#ab 10 268.7
#ab 2 250.6
#ab 6 242.9
#ab 1 242.1
#ab 3 239.4
#ab 5 228.8
#b 4 219.5

kruskalmc(x[,2],x[,2])

#Multiple comparison test after Kruskal-Wallis
#p.value: 0.05
#Comparisons
# obs.dif critical.dif difference

#[......]
#6-9 54.0 162.02688 FALSE
#6-10 69.5 159.04584 FALSE
#6-11 94.5 133.02196 FALSE
#7-8 18.0 160.00778 FALSE
#7-9 35.0 169.78370 FALSE
#7-10 50.5 166.94123 FALSE
#7-11 75.5 142.36796 FALSE
#8-9 17.0 165.54197 FALSE
#8-10 32.5 162.62538 FALSE
#8-11 57.5 137.28174 FALSE
#9-10 15.5 172.25281 FALSE
#9-11 40.5 148.56074 FALSE
#10-11 25.0 145.30369 FALSE



--
View this message in context: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-tp4639004p4639027.html

Sent from the R help mailing list archive at Nabble.com.

[peter dalgaard]

On Aug 3, 2012, at 11:33 , greatest.possible.newbie wrote:

> Thank you for your answer.
> The p.adj argument in the kruskal()-function doesn't seem to change
> anything... Not even the "bonferroni"-method although it is described as the
> most conservative one (multiplying all p-values with the number of
> comparisons). I suppose the kruskal()-function is not working properly...

Apparently, the grouping logic doesn't care about p.adj. Not the most fortunate design in my opinion, but try looking at the output with group=FALSE.

> On the other hand I doubt the method behind the kruskalmc()-function as this
> function doesn't even turn out to detect significant differences between the
> grouping variable (which is obviously a severe error).

That's not obvious! Did you check all group comparisons? How big are the groups?

Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd....@cbs.dk Priv: PDa...@gmail.com

[David L Carlson]

Generally multiple comparisons are conducted after a test for a significant
difference among any of the groups. For your data

> kruskal.test(x[,1]~x[,2])

Kruskal-Wallis rank sum test

data: x[, 1] by x[, 2]
Kruskal-Wallis chi-squared = 11.0098, df = 10, p-value = 0.3568

There are no significant differences between the groups, so there is no
reason to use a multiple comparison test.

----------------------------------------------
David L Carlson
Associate Professor of Anthropology
Texas A&M University
College Station, TX 77843-4352

[peter dalgaard]

On Aug 3, 2012, at 18:49 , David L Carlson wrote:

> Generally multiple comparisons are conducted after a test for a significant
> difference among any of the groups. For your data
>
>> kruskal.test(x[,1]~x[,2])
>
> Kruskal-Wallis rank sum test
>
> data: x[, 1] by x[, 2]
> Kruskal-Wallis chi-squared = 11.0098, df = 10, p-value = 0.3568
>
> There are no significant differences between the groups, so there is no
> reason to use a multiple comparison test.

There's a point to that, but on the other hand, if multiple comparison methods control the familywise error rate by themselves, further "guarding" by a global test really just complicates things even further.

I forgot that we actually had the data... A few points can be made.

First, this is a really unbalanced design:

> table(x[,2])

1 2 3 4 5 6 7 8 9 10 11
267 39 23 24 25 21 17 19 15 16 34

This makes both the grouping method inherently suspect since it assumes at least roughly similar group sizes. However, since it doesn't even attempt to correct for multiple tests, the point is a bit moot.

If we forget the grouping technique, kruskal() gives
> kruskal(x[,1], x[,2], group=F)
....
Comparison between treatments mean of the ranks

Difference pvalue sig LCL UCL
1 - 10 -26.6163390 0.473942 -99.5913928 46.35871
1 - 11 -62.3553095 0.018026 * -113.9832721 -10.72735
...
11 - 3 64.9916880 0.095914 . -11.5557888 141.53916
11 - 4 84.9056373 0.027780 * 9.3153935 160.49588
11 - 5 75.6064706 0.047290 * 0.9077143 150.30523
11 - 6 61.4978992 0.125302 -17.1938261 140.18962
....

so three comparisons are formally significant at level 0.05, but this is without correction for the multiple comparisons. This roughly amounts to multiplying all p values by 55 (actually 55, 54, 53, ... in the Holm method), which of course doesn't leave anything significant:

> kruskal(x[,1], x[,2], group=F, p.adj="holm")

P value adjustment method: holm
Comparison between treatments mean of the ranks

Difference pvalue sig LCL UCL
1 - 10 -26.6163390 1.00000 -150.58158 97.34890
1 - 11 -62.3553095 0.99143 -150.05751 25.34689
...
11 - 3 64.9916880 1.00000 -65.04215 195.02552
11 - 4 84.9056373 1.00000 -43.50211 213.31339
11 - 5 75.6064706 1.00000 -51.28688 202.49982
11 - 6 61.4978992 1.00000 -72.17844 195.17424
...

which is in no way at variance with the global test or the fact that kruskalmc shows no differences to be significant. It also roughly fits results of the stock pairwise.wilcox.test:

> pairwise.wilcox.test(x[,1], x[,2])

Pairwise comparisons using Wilcoxon rank sum test

data: x[, 1] and x[, 2]

1 2 3 4 5 6 7 8 9 10
2 1.00 - - - - - - - - -
3 1.00 1.00 - - - - - - - -
4 1.00 1.00 1.00 - - - - - - -
5 1.00 1.00 1.00 1.00 - - - - - -
6 1.00 1.00 1.00 1.00 1.00 - - - - -
7 1.00 1.00 1.00 1.00 1.00 1.00 - - - -
8 1.00 1.00 1.00 1.00 1.00 1.00 1.00 - - -
9 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 - -
10 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 -
11 0.91 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

P value adjustment method: holm

> pairwise.wilcox.test(x[,1], x[,2], p.adj="none")

Pairwise comparisons using Wilcoxon rank sum test

data: x[, 1] and x[, 2]

1 2 3 4 5 6 7 8 9 10
2 0.705 - - - - - - - - -
3 0.916 0.748 - - - - - - - -
4 0.469 0.557 0.343 - - - - - - -
5 0.675 0.587 0.733 0.920 - - - - - -
6 0.985 0.727 0.805 0.608 0.869 - - - - -
7 0.226 0.349 0.311 0.153 0.187 0.362 - - - -
8 0.282 0.562 0.288 0.135 0.325 0.524 0.874 - - -
9 0.173 0.246 0.296 0.105 0.162 0.351 0.820 0.728 - -
10 0.519 0.650 0.539 0.172 0.407 0.759 0.857 0.855 0.737 -
11 0.016 0.105 0.091 0.041 0.055 0.121 0.727 0.420 0.922 0.499

P value adjustment method: none

(They don't get exactly the same p-values because pairwise.wilcox.test is based on ranks recomputed separately for each pair of groups.)

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd....@cbs.dk Priv: PDa...@gmail.com

______________________________________________
R-h...@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help

[greatest.possible.newbie]

I see.. So apparently the different functions are doing the same thing! :-)
Besides I didn't know the groups should have about the same size.
Thank you four your time Mr. Dalgaard.



--
View this message in context: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-tp4639004p4639431.html

Sent from the R help mailing list archive at Nabble.com.

[Tham Tran]

Dear Mr. Dalgaard,

Could you help me know the name of post-hoc multi-comparaison test mentioned
in kruskal function of agricolae package?

Thank you in advance.

Tham Tran



--
View this message in context: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-tp4639004p4690739.html

Sent from the R help mailing list archive at Nabble.com.

[David Winsemius]

On May 16, 2014, at 2:07 PM, Tham Tran wrote:

> Dear Mr. Dalgaard,

You do realize that was a posting from 2012, right?

>
> Could you help me know the name of post-hoc multi-comparaison test mentioned
> in kruskal function of agricolae package?

There are multiple such tests mentioned on that function help page.

Did you review the material above the Nabble posting? (This is why including the context in every posting is the recommended way to pose a question on r-help.) See also the rest of the Posting Guide.

PLEASE do read the posting guide http://www.R-project.org/posting-guide.html

>
> Thank you in advance.
>
> Tham Tran
>
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Multiple-Comparisons-Kruskal-Wallis-Test-kruskal-agricolae-and-kruskalmc-pgirmess-don-t-yield-the-sa-tp4639004p4690739.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-h...@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
>
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA