In reyzaguirre/pepa: Package for the Execution of Pre Cooked Analysis

library(knitr)
opts_chunk$set(echo = FALSE, comment = NA)

1. ANOVA for variable `r response`

This is a simple analysis of the data, you may want to pursue this further using more sophisticated models.

## Keep copy of full data

dfr.full <- params$dfr

## Data

dfr <- params$dfr
response <- params$response
treat <- params$treat
rep <- params$rep
eu <- params$eu
dap <- params$dap

dfr[, treat] <- as.character(dfr[, treat])
dfr[, rep] <- as.character(dfr[, rep])
dfr[, eu] <- as.character(dfr[, eu])
dfr[, dap] <- as.character(dfr[, dap])

subsamp <- FALSE
if (max(table(dfr[, treat], dfr[, dap], dfr[, eu])) > 1)
  subsamp <- TRUE

## Model crd without subsampling

if (is.null(rep) & !subsamp) {
  ff <- as.formula(paste(response, "~", treat, "*", dap))
  model <- aov(ff, dfr)
  at <- anova(model)
  rownames(at)[4] <- "Exp. Error"
  at
}

## Model rcbd without subsampling

if (!is.null(rep) & !subsamp) {
  ff <- as.formula(paste(response, "~", treat, "*", dap, "+", rep))
  model <- aov(ff, dfr)
  at <- anova(model)
  rownames(at)[5] <- "Exp. Error"
  at
}

## Model crd with subsampling

if (is.null(rep) & subsamp) {
  ff <- as.formula(paste(response, "~", treat, "*", dap, "+",
                         eu, ":", treat, ":", dap))
  model <- aov(ff, dfr)
  at <- anova(model)
  rownames(at)[4:5] <- c("Exp. Error", "Sampling Error")
  at[1:3, 4] <- at[1:3, 3] / at[4, 3]
  at[1:3, 5] <- pf(at[1:3, 4], at[1:3, 1], at[4, 1], lower.tail = FALSE)
  at
}

## Model rcbd with subsampling

if (!is.null(rep) & subsamp) {
  ff <- as.formula(paste(response, "~", treat, "*", dap, "+", rep, "+",
                         eu, ":", treat, ":", dap))
  model <- aov(ff, dfr)
  at <- anova(model)
  rownames(at)[5:6] <- c("Exp. Error", "Sampling Error")
  at[1:4, 4] <- at[1:4, 3] / at[5, 3]
  at[1:4, 5] <- pf(at[1:4, 4], at[1:4, 1], at[5, 1], lower.tail = FALSE)
  at
}

2. Dispersion plot with a fitted smoothing function

## Data

dfr <- dfr.full[, c(response, treat, eu, dap)]

colnames(dfr) <- c("y", "Treatment", "eu", "dap")

if (average)
  dfr <- docomp("mean", "y", c("Treatment", "eu", "dap"), dfr = dfr)

## ggplot

if (se) {
  ggplot(dfr, aes(dap, y, colour = Treatment)) +
    geom_point() + 
    labs(x = dap, y = response) +
    geom_smooth(se = TRUE)
} else {
  ggplot(dfr, aes(dap, y, colour = Treatment)) +
    geom_point() + 
    labs(x = dap, y = response) +
    geom_smooth(se = FALSE)
}