[R-sig-eco] post hoc in Kruskal Wallis

[Richard Boyce]

Ken Aho was kind enough to alert me to the post hoc multiple comparison test KW.multi.comp that he wrote for his package, asbio. This works well for unbalanced designs, and with a minor modification, it also works well when you have tied ranks.

I've included the modified KW.multi.comp as a function. Please give it a whirl and see how it works.


KW.multi.comp <- function (Y, X, conf)
{
    ranks <- rank(Y)
    ni <- tapply(ranks, X, length)
    N <- length(ranks)
    r <- length(ni)
 
    y <- sort(match(ranks,ranks))
      t <- 0
      for (i in 1:(N - 1)) {
          tie.ranks <- sum(match(y, i),na.rm=T)
          if (tie.ranks > 1) t <- t + tie.ranks^3 - tie.ranks
          }
 
    mean.ranks <- tapply(ranks, X, mean)
    dif.mat <- outer(mean.ranks, mean.ranks, "-")
    diffs <- dif.mat[upper.tri(dif.mat)]
    ni.mat <- outer(1/ni, 1/ni, "+")
    ni.den <- ni.mat[upper.tri(ni.mat)]
    SE.diff <- sqrt(((N * (N + 1))/12 - t/(12 * (N - 1))) * ni.den)
    B <- qnorm(1 - ((1 - conf)/(r^2 - r)))
    p.val <- 2 * pnorm(abs(diffs)/SE.diff, lower = FALSE)
    p.adj <- ifelse(p.val * ((r^2 - r)/2) >= 1, 1, round(p.val *
        ((r^2 - r)/2), 6))
    hwidths <- B * SE.diff
    val <- round(cbind(diffs, diffs - hwidths, diffs + hwidths),
        5)
    Decision <- ifelse((val[, 2] > 0 & val[, 3] > 0) | val[,
        2] < 0 & val[, 3] < 0, "Reject H0", "FTR H0")
    val <- as.data.frame(cbind(val, Decision, p.adj))
    lvl <- outer(levels(X), levels(X), function(x1, x2) {
        paste(paste("Avg.rank", x1, sep = ""), paste("Avg.rank",
            x2, sep = ""), sep = "-")
    })
    dimnames(val) <- list(lvl[upper.tri(lvl)], c("Diff", "Lower",
        "Upper", "Decision", "Adj. p-value"))
    val
}

> Message: 1
> Date: Tue, 6 Dec 2011 16:33:52 -0500
> From: "Richard L. Boyce" <boycer at nku.edu>
> To: <r-sig-ecology at r-project.org>
> Subject: Re: [R-sig-eco] post hoc in Kruskal Wallis
> Message-ID: <a06240800cb042a042650@[192.168.200.253]>
> Content-Type: text/plain
> 
> As promised, here is a function for the post hoc test in Kruskal 
> Wallis when the design is unbalanced (group sizes are unequal). Note 
> that both my and Dave's routine won't work if you have ties--perhaps 
> I'll get versions made for those cases in the future.
> 
> npcompare.unequal <- function (var, class)
> {
>       k <- nlevels(class) #groups
>       n <- array(dim=k)
> 
>       for (i in 1:k) n[i] <- as.numeric(table(class)[i]) # no. in each group
> 
>       rvar <- rank(var) #ranks
> 
>       #prints 1st group, 2nd group, difference, mean 1st group rank, 
> mean 2nd group rank, standard error, Q, and p
>       cat("a b Ra Rb diff se Q p", '\n') #header column
> 
>       for (i in 1:(k-1)) {
>           left <- mean(rvar[as.numeric(class)==i])
>           for (j in (i+1):k) {
>               right <- mean(rvar[as.numeric(class)==j])
>               diff = abs(left-right)
>               se <- sqrt((length(var)*(length(var)+1))/12*(1/n[i]+1/n[j]))
>               p <- round(k*(k-1)*(1-pnorm(diff/se)),4)
>               if (p > 0.50) p <- c(">0.50")
>               cat(paste(i, j, round(left,2), round(right,2), 
> round(diff,2), round(se,4), round(diff/se,4), p, '\n'))
>           }
>       }
>       cat('\n')
> }
> 
> 
>> Message: 1
>> Date: Thu, 01 Dec 2011 10:27:20 -0700
>> From: Dave Roberts <dvrbts at ecology.msu.montana.edu>
>> To: r-sig-ecology at r-project.org
>> Subject: Re: [R-sig-eco] post hoc in Kruskal Wallis
>> Message-ID: <4ED7B8F8.1070801 at ecology.msu.montana.edu>
>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>> 
>> If the design is balanced this function will work
>> 
>> npcompare <- function (var,class)
>> {
>>      k <- length(table(class))
>>      n <- length(var)/k
>> 
>>      se <- sqrt((n*(n*k)*(n*k+1))/12)
>>      rvar <- rank(var)
>> 
>>      for (i in 1:(k-1)) {
>>          left <- sum(rvar[as.numeric(class)==i])
>>          for (j in (i+1):k) {
>>              right <- sum(rvar[as.numeric(class)==j])
>>              diff = abs(left-right)
>>              cat(paste(i,j,left,right,diff,round(1-pt(diff/se,k),4),'\n'))
>>          }
>>      }
>>      cat('\n')
>> }
>> 
>> It's not formatted very elegantly, but it works.
>> 
>> Dave R.
>> 
>> On 11/24/2011 08:01 AM, Richard L. Boyce wrote:
>>>  Zar, in Biostatistical Analysis, has a Tukey-like test that is highly
>>>  recommended. You may need to do it by hand, but it is quite easy.
>>> 
>>>>  Message: 5
>>>>  Date: Wed, 23 Nov 2011 18:21:54 +0100
>>>>  From: "Jakub Szymkowiak"<qbaszym at tlen.pl>
>>>>  To:<r-sig-ecology at r-project.org>
>>>>  Subject: [R-sig-eco]  post hoc in Kruskal Wallis
>>>>  Message-ID:<C8041F4B78EB487C8707E72114C3F4CF at komp>
>>>>  Content-Type: text/plain; format=flowed; charset="iso-8859-1";
>>>> 	reply-type=original
>>>> 
>>>>  Hi,
>>>>  does anyone know, how can I perform post-hoc tests (especially Least
>>>>  Significant Difference and Sheffe Test) for results from Kruskal-Wallis
>>>>  test? In KruskaI-Wallis test I found some significant differences between
>>>>  tested groups, but I want to know between which groups this difference is
>>>>  really signifficant.
>>>> 
>>>>  Cheers,
>>>>  Jakub
>>> 
>> 
>> --
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>> David W. Roberts                                     office 406-994-4548
>> Professor and Head                                      FAX 406-994-3190
>> Department of Ecology                         email droberts at montana.edu
>> Montana State University
>> Bozeman, MT 59717-3460
> 
> -- 
> ================================
> Richard L. Boyce
> Director, Environmental Science Program
> Associate Professor
> Department of Biological Sciences
> Northern Kentucky University
> Nunn Drive
> Highland Heights, KY  41099  USA
> 
> 859-572-1407 (tel.)
> 859-572-5639 (fax)
> boycer at nku.edu
> http://www.nku.edu/~boycer/
> =================================
> 
> "One of the advantages of being disorderly is that one is constantly 
> making exciting discoveries." - A.A. Milne