post hoc in Kruskal Wallis

Tue, Dec 13, 2011 7:25 AM

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

post hoc in Kruskal Wallis

Thread (5 messages)