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

r i = {{i}}

# Check design

lc <- ck.f(dfr, vars[i], factors, rep)

# Fit a model for assumptions plots

expr <- paste(vars[i], '~', factors[1])
for (j in 2:lc$nf)
  expr <- paste(expr, '*', factors[j])

if (is.null(rep))
  ff <- as.formula(expr)

if (!is.null(rep)) {
  expr <- paste(expr, '+', rep)
  ff <- as.formula(expr)
}

model <- aov(ff, dfr)

# Estimate missing values

y.est <- paste0(vars[i], ".est")

if (lc$nmis > 0) {
  dfr[, y.est] <- mve.f(dfr, vars[i], factors, rep, maxp)[, y.est]
} else {
  dfr[, y.est] <- dfr[, vars[i]]
}

# Get anova table with estimated missing values

at <- suppressWarnings(aov.f(dfr, vars[i], factors, rep, maxp))

# CV

rr <- dim(at)[1]
cv <- (at[rr, 3])^0.5 / mean(dfr[, y.est]) * 100

{{i+1}}. Analysis for variable `r vars[i]`

r if (lc$nmis == 0) {"There are no missing values for this variable; the design is balanced."}

r if (lc$nmis > 0) paste0("There are some missing values (", format(lc$pmis * 100, digits = 3), "%) and they have been estimated for the descriptive statistics and ANOVA.")

{{i+1}}.1. Descriptive statistics

{{i+1}}.1.1. Means by individual factor levels

for (j in 1:lc$nf) {
  print(tapply(dfr[, y.est], dfr[, factors[j]], mean))
  cat("\n")
}

{{i+1}}.1.2. Means by factor levels combinations

# Create expression for list of factors

lf.expr <- 'list(dfr[, factors[1]]'

for (j in 2:lc$nf)
  lf.expr <- paste0(lf.expr, ', dfr[, factors[', j, ']]')

lf.expr <- paste0(lf.expr, ')')

# Compute means over replications

tapply(dfr[, y.est], eval(parse(text = lf.expr)), mean)

{{i+1}}.2. ANOVA

{{i+1}}.2.1. Checking assumptions

As it was stated in section 1, it is supposed that the error has a normal distribution with the same variance for all the combinations among the levels of the factors. The following plots help to evaluate these assumptions:

par(mfrow = c(1, 2))
tryCatch(
  suppressWarnings(plot(model, which = 1)),
    error = function(e) {}
)
tryCatch(
  suppressWarnings(plot(model, which = 2)),
    error = function(e) {}
)

Funnel shapes for the first plot may suggest heterogeneity of variances while departures from the theoretical normal line are symptoms of lack of normality.

{{i+1}}.2.2. ANOVA table

at

The coefficient of variation for this experiment is r format(cv, digits = 4)%.

nf <- lc$nf

if (is.null(rep)) {
  df.res <- at[4, 1]
  ms.res <- at[4, 3]
  pv.gen <- at[1, 5]
  pv.loc <- at[2, 5]
  pv.int <- at[3, 5]
  } else {
  df.res <- at[5, 1]
  ms.res <- at[5, 3]
  pv.gen <- at[1, 5]
  pv.loc <- at[2, 5]
  pv.int <- at[4, 5]
}

r if (!is.nan(pv.gen) & !is.nan(pv.int)) {if (nf == 2 & (pv.gen < 0.05 | pv.int < 0.05 | pe | se)) {paste0("### ", {{i+1}}, ".2.3. LSD test for ", factors[1])}}

if (!is.nan(pv.gen) & !is.nan(pv.int)) {

  if (nf == 2 & ((pv.gen < 0.05 & pv.int > 0.05) | pe)) {
    y <- dfr[, y.est]
    print(agricolae::LSD.test(y, dfr[, factors[1]], df.res, ms.res)$groups)
    cat("\n")
  }

  if (nf == 2 & pv.int < 0.05 | se) {
    for (j in 1:lc$nl[2]) {
      tmp <- dfr[dfr[, factors[2]] == unique(dfr[, factors[2]])[j], ]
      y <- tmp[, y.est]
      print(paste("LSD test at", factors[2], unique(dfr[, factors[2]])[j]))
      print(agricolae::LSD.test(y, tmp[, factors[1]], df.res, ms.res)$groups)
      cat("\n")
    }
  }
}

r if (!is.nan(pv.gen) & !is.nan(pv.int)) {if (nf == 2 & (pv.gen < 0.05 | pv.int < 0.05 | pe | se)) {paste0("### ", {{i+1}}, ".2.4. Tukey test for ", factors[1])}}

if (!is.nan(pv.gen) & !is.nan(pv.int)) {

  if (nf == 2 & ((pv.gen < 0.05 & pv.int > 0.05) | pe)) {
    y <- dfr[, y.est]
    print(agricolae::HSD.test(y, dfr[, factors[1]], df.res, ms.res)$groups)
    cat("\n")
  }

  if (nf == 2 & pv.int < 0.05 | se) {
    for (j in 1:lc$nl[2]) {
      tmp <- dfr[dfr[, factors[2]] == unique(dfr[, factors[2]])[j], ]
      y <- tmp[, y.est]
      print(paste("HSD Tukey test at", factors[2], unique(dfr[, factors[2]])[j]))
      print(agricolae::HSD.test(y, tmp[, factors[1]], df.res, ms.res)$groups)
      cat("\n")
    }
  }
}

r if (!is.nan(pv.gen) & !is.nan(pv.loc) & !is.nan(pv.int)) {if (nf == 2 & (pv.loc < 0.05 & pv.gen > 0.05 & pv.int > 0.05 & !pe & !se)) {paste0("### ", {{i+1}}, ".2.3. LSD test for ", factors[2])}} r if (!is.nan(pv.gen) & !is.nan(pv.loc) & !is.nan(pv.int)) {if (nf == 2 & (pv.int < 0.05 | (pv.gen < 0.05 & pv.loc < 0.05) | pe | se)) {paste0("### ", {{i+1}}, ".2.5. LSD test for ", factors[2])}}

if (!is.nan(pv.gen) & !is.nan(pv.loc) & !is.nan(pv.int)) {

  if (nf == 2 & ((pv.loc < 0.05 & pv.int > 0.05) | pe)) {
    y <- dfr[, y.est]
    print(agricolae::LSD.test(y, dfr[, factors[2]], df.res, ms.res)$groups)
    cat("\n")
  }

  if (nf == 2 & (pv.int < 0.05 | se)) {
    for (j in 1:lc$nl[1]) {
      tmp <- dfr[dfr[, factors[1]] == unique(dfr[, factors[1]])[j], ]
      y <- tmp[, y.est]
      print(paste("LSD test at", factors[1], unique(dfr[, factors[1]])[j]))
      print(agricolae::LSD.test(y, tmp[, factors[2]], df.res, ms.res)$groups)
      cat("\n")
    }
  }
}

r if (!is.nan(pv.gen) & !is.nan(pv.loc) & !is.nan(pv.int)) {if (nf == 2 & (pv.loc < 0.05 & pv.gen > 0.05 & pv.int > 0.05 & !pe & !se)) {paste0("### ", {{i+1}}, ".2.4. Tukey test for ", factors[2])}} r if (!is.nan(pv.gen) & !is.nan(pv.loc) & !is.nan(pv.int)) {if (nf == 2 & (pv.int < 0.05 | (pv.gen < 0.05 & pv.loc < 0.05) | pe | se)) {paste0("### ", {{i+1}}, ".2.6. Tukey test for ", factors[2])}}

if (!is.nan(pv.gen) & !is.nan(pv.loc) & !is.nan(pv.int)) {

  if (nf == 2 & ((pv.loc < 0.05 & pv.int > 0.05) | pe)) {
    y <- dfr[, y.est]
    print(agricolae::HSD.test(y, dfr[, factors[2]], df.res, ms.res)$groups)
    cat("\n")
  }

  if (nf == 2 & (pv.int < 0.05 | se)) {
    for (j in 1:lc$nl[1]) {
      tmp <- dfr[dfr[, factors[1]] == unique(dfr[, factors[1]])[j], ]
      y <- tmp[, y.est]
      print(paste("HSD Tukey test at", factors[1], unique(dfr[, factors[1]])[j]))
      print(agricolae::HSD.test(y, tmp[, factors[2]], df.res, ms.res)$groups)
      cat("\n")
    }
  }
}