1
votes

All,

I would like to perform the equivalent of TukeyHSD on the rank ordering median shift test that such as kruskal wallis

X=matrix(c(1,1,1,1,2,2,2,4,4,4,4,4,1,3,6,9,4,6,8,10,1,2,1,3),ncol=2)
anova=aov(X[,2]~factor(X[,1]))
TukeyHSD(anova)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = X[, 2] ~ factor(X[, 1]))
##
## $`factor(X[, 1])`
## diff lwr upr p adj
## 2-1 1.25 -5.927068 8.427068 0.8794664
## 4-1 -1.35 -7.653691 4.953691 0.8246844
## 4-2 -2.60 -9.462589 4.262589 0.5617125
kruskal.test(X[,2]~factor(X[,1]))
##
## Kruskal-Wallis rank sum test
##
## data: X[, 2] by factor(X[, 1])
## Kruskal-Wallis chi-squared = 1.7325, df = 2, p-value = 0.4205

I would like now to analyze the contrasts. Please help. Thanks. Rik

2
You can do post-hoc analysis with the kruskalmc function from the pgirmess package.Jaap
Hi @Rik, if any answer solves your problem can you click on "accept it" so that other people can see it? thanksagenis

2 Answers

0
votes

If you want to do multiple comparisons after a Kruskal-Wallis test, you need the kruskalmc function from the pgirmess package. Before you can implement this function, you will need to transform your matrix to a dataframe. In your example:

# convert matrix to dataframe
dfx <- as.data.frame(X)

# the Kruskal-Wallis test & output
kruskal.test(dfx$V2~factor(dfx$V1))


    Kruskal-Wallis rank sum test

data:  dfx$V2 by factor(dfx$V1)
Kruskal-Wallis chi-squared = 1.7325, df = 2, p-value = 0.4205

# the post-hoc tests & output
kruskalmc(V2~factor(V1), data = dfx)

Multiple comparison test after Kruskal-Wallis 
p.value: 0.05 
Comparisons
    obs.dif critical.dif difference
1-2    1.75     6.592506      FALSE
1-4    1.65     5.790265      FALSE
2-4    3.40     6.303642      FALSE
0
votes

If you want the compact letter display similar to what is outputed from TukeyHSD, for the Kruskal test, the library agricolae allows it with the function kruskal. Using your own data:

  library(agricolae)
print( kruskal(X[, 2], factor(X[, 1]), group=TRUE, p.adj="bonferroni") )
#### ...
#### $groups
####   trt means M
#### 1   2  8.50 a
#### 2   1  6.75 a
#### 3   4  5.10 a

(well, in this example the groups are not considered different, same result than the other answer..)