R/wmc.R

In Zachary-Wu/learn: learn something about development of R package

#' @description Nonparametric pairwise multiple comparisons using the Wilcoxon Signed Rank Test
#' @description Probability values are adjusted using the p.adjust function
#' @description The url of the source code is:
#' @description \url{https://www.statmethods.net/RiA/wmc.txt}
#' @title wmc
#'
#' @param formula
#' @param data
#' @param exact
#' @param sort
#' @param method
#' @importFrom stats wilcox.test aggregate
#' @return Nonparametric pairwise multiple comparisons using the Wilcoxon Signed Rank Test,
#' and probability values are adjusted using the p.adjust function.
#' @export
#'
#' @examples
wmc <- function(formula, data, exact=FALSE, sort=TRUE, method="holm"){

  # setup
  df <- model.frame(formula, data)
  y <- df[[1]]
  x <- as.factor(df[[2]])


  # reorder levels of x by median y
  if(sort){
    medians <- aggregate(y, by=list(x), FUN=median)[2]
    index <- order(medians)
    x <- factor(x, levels(x)[index])
  }

  groups <- levels(x)
  k <- length(groups)

  # summary statistics
  stats <- function(z)(c(N = length(z), Median = median(z), MAD = mad(z)))
  sumstats <- t(aggregate(y, by=list(x), FUN=stats)[2])
  rownames(sumstats) <- c("n", "median", "mad")
  colnames(sumstats) <- groups
  cat("Descriptive Statistics\n\n")
  print(sumstats)

  # multiple comparisons
  mc <- data.frame(Group.1=character(0),
                   Group.2=character(0),
                   W=numeric(0),
                   p.unadj=numeric(0),
                   p=numeric(0),
                   stars=character(0),
                   stringsAsFactors=FALSE)

  # perform Wilcoxon test
  row <- 0
  for(i in 1:k){
    for(j in 1:k){
      if (j > i){
        row <- row + 1
        y1 <- y[x==groups[i]]
        y2 <- y[x==groups[j]]
        test <- wilcox.test(y1, y2, exact=exact)
        mc[row,1] <- groups[i]
        mc[row,2] <- groups[j]
        mc[row,3] <- test$statistic
        mc[row,4] <- test$p.value
      }
    }
  }
  mc$p <- p.adjust(mc$p.unadj, method=method)

  # add stars
  mc$stars <- " "
  mc$stars[mc$p <   .1] <- "."
  mc$stars[mc$p <  .05] <- "*"
  mc$stars[mc$p <  .01] <- "**"
  mc$stars[mc$p < .001] <- "***"
  names(mc)[6] <- " "

  cat("\nMultiple Comparisons (Wilcoxon Rank Sum Tests)\n")
  cat(paste("Probability Adjustment = ", method, "\n\n", sep=""))
  print(mc[-4], right=TRUE)
  cat("---\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")
  return(invisible(NULL))

}