# 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))
}
