R/npmc.test.R

In helophilus/ColsTools: A variety of convenience tools and short-cuts

#' multiple comparisons for non-parametric data after, e.g., Kruskal-Wallis test
#'
#' \code{npmc.test(formula,data,display=1)}
#'
#' @param \code{formula} of the type \code{resp~explan} where <i>resp</i> is the response variable and explan the explanatory variable.
#' @param \code{data} is the name of the data frame containing the variables
#' @param \code{display} setting this to zero produces no output table
#'
#' @details Code for carrying out Dunn/Nemenyi test after significant Kruskal Wallis test using formula
#' If group sizes are equal Nemenyi is used, otherwise Dunn.
#'
#' See Zar's Biostatistical Analysis section 11.6 for details
#' The test used is printed unless display=0


npmc.test <- function(formula,data=NULL,display=1){
  options(warn=-1)
  attach(data)
  category=get(all.vars(formula)[2])
  resp=get(all.vars(formula)[1])
  category <- factor(category)
  detach(data)
  L <- nlevels(category)
  lvl <- levels(category)
  ns=numeric(L)
  rk <- rank(resp)		# rank the values
  lengths <- tapply(rk,category,length)[]
  means <- tapply(rk,category,mean)[]
  rsums <- tapply(rk,category,sum)[]
  medians <- tapply(resp,category,median)
  Nem=ifelse(length(unique(lengths))==1,1,0)# distinguishes which test is needed
  if(display!=0){
    if(Nem==1){(cat("\n","Nemenyi Test Results  ","\n","\n"))#
    } else {(cat("\n","Dunn Test Results  ","\n","\n"))}	# Correct test title
  }
  Colbl<-c("median1","median2","     Q","     SE"," p-value"," sig")# column headers for output table
  Tl <- length(resp)			# find N
  cons <- Tl*(Tl+1)/12			# pre-calc part of later calcs
  X=L-1;Y=L-1				# set start for permutation loop
  table<-NULL; RLABS <- NULL	# set up vectors for output values
  Ties=ifelse(length(rk)==length(unique(rk)),0,1) # establish if ranks are tied
  if(Ties==1 & Nem==0 & display==1)(cat("SE calculation based upon Tied Ranks formula","\n","\n"))
  ng=unique(lengths)		#
  SE=sqrt(ng^2*L*(ng*L+1)/12)	# calculate SE for Nemenyi Test overwritten if Dunn
  for (a in 1:X) {		# start permutation outer loop
    for (b in 1:Y) {	# start permutation inner loop
      b2 <- a + b	# b2 is used as second permute
      if(Nem==1){			#
        diff=abs(rsums[a]-rsums[b2])	#
        Q=diff/SE			# calc Nemenyi Q
        Pr=ptukey(Q,L,Inf,lower=F)	# calc Nemenyi Pr
      } else {			# Dunn Test calcs begin
        if(Ties==1){	# calc additional bit for tied ranks
          U=unique(sort(resp[duplicated(resp)]))
          m=length(U)
          ti=numeric(m)
          for(h in 1:m){
            ti[h]=length(resp[resp==U[h]])
          }
          sumt=sum(ti^3-ti)	# result of the tied ranks calculation bit
        }
        first <- lengths[a]	# length of the first compared group
        second <- lengths[b2]	# length of the second compared group
        if (Ties==0){
          SE<- sqrt(cons*(1/first + 1/second)) # calc SE not tied
        } else {
          SE<- sqrt((cons-(sumt/(12*Tl-1)))*(1/first + 1/second))
        }					# calc SE for ties
        rankdiff <- abs(means[a]-means[b2])
        Q<-rankdiff/SE					# calc Dunn Q
        Pr<-pnorm(-Q)*(L*(L-1))				# calc Dunn Pr
      }		# Finish Dunn bit and go to common stuff
      sig="";
      if(Pr<0.05)(sig="*");if(Pr<0.01)(sig="**");if(Pr<0.001)(sig="***")
      if(Pr>1)(Pr=1)
      RLAB <- paste(lvl[a]," vs ",lvl[b2]," ")
      ROW<-(c(medians[a],medians[b2],round(Q,3),round(SE,3),round(Pr,4),sig))
      table <- rbind(table,ROW)
      RLABS <- rbind(RLABS,RLAB)
    }
    Y<-Y-1 }
  options(warn=0)
  OUT <- matrix(table,nc=6,dimnames=list(RLABS,Colbl))
  OUTDF=as.data.frame(OUT)
  if(display!=0){
    print(OUTDF)
  }else{
    return(invisible(OUTDF))
  }
}