R/add.ps.R

In raw.means.plot: Raw-Means Plots for Experimental Designs

add.ps <- function(data, col.id, col.offset, col.x, col.value, fun.aggregate = "mean", ref.offset = 1, prefixes, test = "t", alternative = c("two.sided", "less", "greater"), distribution = approximate(100000), mu = 0, paired = FALSE, var.equal = FALSE, lty = 0, ... ) {
  round.ps <- function(x) {
		as.character(ifelse(x < 0.001, "<.001", ifelse(x < 0.01, paste("=", substr(round(x, 3),2,5), sep =""),ifelse((x < 0.05) & (round(x, 2) == 0.05), "<.05", ifelse(round(x,2) == 1, ">.99", paste("=", substr(round(x, 2), 2, 4), sep = ""))))))
	}
  if(!is.data.frame(data)) stop("data must be a data.frame")
  columns <- c(col.id, col.offset, col.x, col.value)
  if (any(!(columns %in% colnames(data)))) stop("column not matching the data")
  implemented.tests <- c("t", "stats.wilcoxon", "Wilcoxon", "permutation", "median")
  test = match.arg(test, implemented.tests)
  formula.agg <- as.formula(paste(col.value, "~", col.id, "+", col.offset, "+", col.x))
  d.new <- aggregate(formula.agg, data = data, FUN = fun.aggregate)
  #so far same as in rm.plot2
  l.offset <- levels(factor(d.new[,col.offset]))
  l.offset.test <- l.offset[-ref.offset]
  if (missing(prefixes)) prefixes <- paste("p(", l.offset.test, ")", sep = "")
  l.x <- levels(factor(d.new[,col.x]))

  for (c.offset in seq_along(l.offset.test)) {
    tmp.ps <- vector("numeric", length(l.x))
    for (c.x in seq_along(l.x)) {
      tmpx <- d.new[d.new[,col.offset] == l.offset[ref.offset] & d.new[,col.x] == l.x[c.x], ]
	  tmpx <- tmpx[order(tmpx[,col.id]),]
	  tmpy <- d.new[d.new[,col.offset] == l.offset.test[c.offset] & d.new[,col.x] == l.x[c.x], ]
	  tmpy <- tmpy[order(tmpy[,col.id]),]
      if (test == "t") {
        tmp.ps[c.x] <- t.test(tmpx[,col.value], tmpy[,col.value], alternative = alternative, mu = mu, paired = paired, var.equal = var.equal)$p.value
      } else {
        if (test == "stats.wilcoxon") tmp.ps[c.x] <- wilcox.test(tmpx[,col.value], tmpy[,col.value], alternative = alternative, mu = mu, paired = paired)$p.value
        else if (test %in% c("Wilcoxon", "permutation", "median")) tmp.ps[c.x] <- compare.2.vectors(tmpx[,col.value], tmpy[,col.value], alternative = alternative, paired = paired, tests = "nonparametric", coin = test, perm.distribution = distribution)[["nonparametric"]][2,"p"]
      }
    }
    tmp.labels <- paste(prefixes[c.offset], round.ps(tmp.ps), sep = "")
    axis(1,seq_along(l.x), labels = tmp.labels, line = c.offset, lty = lty,  ...)
  }      
}