R/anova.plots.R

In fionahevans/agric: Functions for accessing & handling agricultural data.

stderr <- function(x) {
  n <- sqrt(length(x[!is.na(x)]))
  sd(x, na.rm=T)/n
}

#' Dyamite plot
#' 
#' Plot an ANOVA style barchart with standard error bars
#'
#' @param x string; name of covariate in data frame.
#' @param y string; name of predictor in data frame.
#' @param data data frame; ninput data frame.
#'
#' @author Fiona Evans
#' 
#' @return NULL.
#' @export
dynamite.plot <- function(x, y, data) {
  
  # Dynamite plot
  means <- tapply(data[, y], as.factor(data[, x]), mean, na.rm=T)
  se <- tapply(data[, y], as.factor(data[, x]), stderr)
  
  n <- nlevels(as.factor(data[, x]))
  
  a <- aov(data[, y] ~ as.factor(data[, x]))
  
  p <- TukeyHSD(a)[[1]][1:(n-1), 4]
  txt <- c("")
  for (i in 1:length(p)) {
    if (p[i] < 0.1) this <- "*"
    if (p[i] < 0.01) this <- "**"
    if (p[i] < 0.001) this <-  "***"
    txt <- c(txt, this)
  }
  (bp <- barplot(means, xlab=y,ylab=x, ylim=c(0, 1.1*max(means+se, na.rm=T))))
  box()
  
  # Add SE
  ii <- diff(range(bp))/(n-1)/20
  for (i in 1:n) {
    lines(x=c(bp[i], bp[i]), y=c(means[i] - 2.96*se[i], means[i] + 2.96*se[i]))
    lines(x=c(bp[i]-ii, bp[i]+ii), y=c(means[i] - 2.96*se[i], means[i] - 2.96*se[i]))
    lines(x=c(bp[i]-ii, bp[i]+ii), y=c(means[i] + 2.96*se[i], means[i] + 2.96*se[i]))
  }
  
  # Add significance
  text(txt,x=bp,y=1.05*(means+2.96*abs(se)))
}

#' Anova plot
#' 
#' Plot an ANOVA style barchart with standard error bars
#'
#' @param x string; name of covariate in data frame.
#' @param y string; name of predictor in data frame.
#' @param data data frame; ninput data frame.
#'
#' @author Fiona Evans
#' 
#' @return NULL.
anova.plot <- function(x, y, data) {
  
  # Dynamite plot
  means <- tapply(data[, y], as.factor(data[, x]), mean, na.rm=T)
  se <- tapply(data[, y], as.factor(data[, x]), stderr)
  
  n <- nlevels(as.factor(data[, x]))
  
  a <- aov(data[, y] ~ as.factor(data[, x]))
  
  p <- TukeyHSD(a)[[1]][1:(n-1), 4]
  txt <- c("")
  for (i in 1:length(p)) {
    if (p[i] < 0.1) this <- "*"
    if (p[i] < 0.01) this <- "**"
    if (p[i] < 0.001) this <-  "***"
    txt <- c(txt, this)
  }
  (bp <- barplot(means, xlab=y,ylab=x, ylim=c(0, 1.1*max(means+se, na.rm=T))))
  box()
  
  # Add SE
  ii <- diff(range(bp))/(n-1)/20
  for (i in 1:n) {
    lines(x=c(bp[i], bp[i]), y=c(means[i] - 2.96*se[i], means[i] + 2.96*se[i]))
    lines(x=c(bp[i]-ii, bp[i]+ii), y=c(means[i] - 2.96*se[i], means[i] - 2.96*se[i]))
    lines(x=c(bp[i]-ii, bp[i]+ii), y=c(means[i] + 2.96*se[i], means[i] + 2.96*se[i]))
  }
  
  # Add significance
  text(txt,x=bp,y=1.05*(means+2.96*abs(se)))
}