# Nonparametric pairwise multiple comparisons using the Wilcoxon Signed Rank Test # Probability values are adjusted using the p.adjust function 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)) }