FDB1: A brew of utility and convenience functions

##' @title Wrapper around \code{\link[stats]{aov}} with better output.
##'
##' @description Wrapper around \code{\link[stats]{aov}} for repeated measures anova.
##'
##' @details
##' No details yet.
##'
##' @param formula RHS of formula should be in partheses, () between any predictor and unit: y~(x1+x2)
##' @param unit Factor for stratification units, typicall subject or item, used to build error term
##' @param cv Say something here
##' @param data Data.frame passed to \code{\link[stats]{aov}}
##' @param dbg Enable verbose output for debugging
##' @param ... Additional arguments passed to \code{\link[stats]{aov}}
##' @return An anova table (data.frame)
##' @author David Braze \email{davebraze@@gmail.com}
##' @seealso \code{\link[stats]{aov}}
##' @family aovtable
##' @export
aovtable <-
function(formula=NULL, unit=NULL, cv=NULL, data=NULL, dbg=FALSE, ...) {

  if(is.null(formula)) stop ("formula must be specified\n")
  if(!inherits(formula, "formula")) stop("Bad formula.")
  if(is.null(data)) stop ("data must be specified\n")

  if (!is.null(unit)){
    form <- formula(paste(deparse(formula(formula)), "*", unit))
    mod <- aov(formula=form, data=data)

    if (dbg) {
      print(form)
      print (summary(mod))
    }

    mss <- sum(fitted(mod)^2)             # model SS
    rss <- sum(residuals(mod)^2)          # residual SS
    r.squared <- mss/(mss+rss)
    tab <- summary(mod)
    tab <- tab[[1]]

    rxp <- paste("(", unit, "|Residuals)", sep="")
    tab.f <- tab[-grep(rxp, rownames(tab)),]
    ## print(tab)
    result <- matrix(nrow=dim(tab.f)[1], ncol=6)
    colnames(result) <- c("Df1", "SS1", "MS1", "Df2", "SS2", "MS2")
    rownames(result) <- rownames(tab.f)

    nl <- regexpr(" ", rownames(result))
    rownames(result) <- strtrim(rownames(result), max(nl))

    result <- data.frame(result)
    result[,"Df1"] <- tab.f[,"Df"]
    result[,"SS1"] <- tab.f[,"Sum Sq"]
    result[,"MS1"] <- tab.f[,"Mean Sq"]

    for (rr in rownames(result)){
      e <- grep(paste("^", gsub(" +", "", rr), ":",
                      unit, sep=""), rownames(tab))
      if(length(e)>0) {
        result[rr,4:6] <- tab[e,1:3]
      } else {
        result[rr,4:6] <- tab[grep("Residuals", rownames(tab)),1:3]
      }
      result[rr,"F value"] <- result[rr,"MS1"]/result[rr,"MS2"]
      result[rr,"Pr(>F)"] <- pf(result[rr, "F value"], result[rr,"Df1"], result[rr, "Df2"],
                                lower.tail=FALSE)
      result[rr,"eta sq"] <- result[rr, "SS1"]/(result[rr,"SS1"]+result[rr, "SS2"])
    } # end for()
    attr(result, "units") <- unit
  } else {
    ## print(noquote("  Units: None"))
    mod <- aov(formula=formula, data=data)

    mss <- sum(fitted(mod)^2)             # model SS
    rss <- sum(residuals(mod)^2)          # residual SS
    r.squared <- mss/(mss+rss)
    tab <- summary(mod)
    tab <- tab[[1]]
    ## print(tab)

    tab.f <- tab[-grep("Residuals", rownames(tab)),]
    result <- matrix(nrow=dim(tab.f)[1], ncol=6)
    colnames(result) <- c("Df1", "SS1", "MS1", "Df2", "SS2", "MS2")
    rownames(result) <- rownames(tab.f)

    nl <- regexpr(" ", rownames(result))
    rownames(result) <- strtrim(rownames(result), max(nl))
    result <- data.frame(result)

    result[,"Df1"] <- tab.f[,"Df"]
    result[,"SS1"] <- tab.f[,"Sum Sq"]
    result[,"MS1"] <- tab.f[,"Mean Sq"]

    for (rr in rownames(result)){
      result[rr,4:6] <- tab[grep("Residuals", rownames(tab)),1:3]
      result[rr,"F value"] <- result[rr,"MS1"]/result[rr,"MS2"]
      result[rr,"Pr(>F)"] <- pf(result[rr, "F value"], result[rr,"Df1"], result[rr, "Df2"],
                                lower.tail=FALSE)
      result[rr,"eta sq"] <- result[rr, "SS1"]/(result[rr,"SS1"]+result[rr, "SS2"])
    } # end for()
    attr(result, "units") <- ""

  } # end if ... else
  attr(result, "R2") <- sprintf("%0.5f", r.squared)
  attr(result, "formula") <- deparse(formula)
  result <- list(table=result, model=mod)
  class(result) <- c("aovtable", class(result))
  return(result)
}