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

r i = {{i}}

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

# Check design

lc <- ck.rcbd(dfr, vars[i], trt, rep)

# Fit a model for assumptions plots

model <- aov(dfr[, vars[i]] ~ dfr[, trt] + dfr[, rep])

# Estimate missing values

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

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

# Get anova table with estimated missing values

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

# CV

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

# ANOVA with subsampling

show.anova.sub <- FALSE

if (!is.null(eu))
  if (sum(is.na(dfr2[, vars[i]])) == 0 & dim(dfr)[1] < dim(dfr2)[1]) {
    model2 <- aov(dfr2[, vars[i]] ~ dfr2[, trt] + dfr2[, rep] + dfr2[, eu] %in% dfr2[, trt])
    model2$terms[[2]] <- vars[i]
    at2 <- anova(model2)
    at2[1, 4] <- at2[1, 3] / at2[3, 3]
    at2[2, 4] <- at2[2, 3] / at2[3, 3]
    at2[1, 5] <- pf(at2[1, 4], at2[1, 1], at2[3, 1], lower.tail = FALSE)
    at2[2, 5] <- pf(at2[2, 4], at2[2, 1], at2[3, 1], lower.tail = FALSE)
    rownames(at2) <- c(trt, rep, "Exp. Error", "Sampling Error")
    show.anova.sub <- TRUE
  }

{{i+1}}.1. Exploratory analysis

It is always good to have some visualization of your data. Below a histogram and a boxplot are shown.

par(mfrow = c(1, 2))
hist(dfr[, vars[i]], main = paste("Histogram of", vars[i]), xlab = vars[i])
boxplot(dfr[, vars[i]], ylab = vars[i])

r if(lc$ng < 10) {paste0(" Since the number of " , trt.lab.s, " in your experiment is not so large, we can plot the data for each ", trt.lab, ":")}

if (lc$ng < 10) msdplot(dfr, vars[i], trt, conf = 1, pch = 4)

{{i+1}}.2. ANOVA

You have fitted a linear model for a RCBD. The ANOVA table for your model is:

at

r if(lc$nmis > 0) paste0("You have some missing values (", format(lc$pmis * 100, digits = 3), "%) and they have been estimated before running ANOVA.")

The coefficient of variation for this experiment is r format(cv, digits = 4)%. The p-value for r trt.lab.s is r format(at[1, 5], digits = 4) r if(at[1, 5] < 0.05) {"which is significant at the 5% level."} else {"which is not significant at the 5% level."}

r if(show.anova.sub) {"At the subsample level, the ANOVA table is:"}

if (show.anova.sub)
  at2

r if(show.anova.sub) {if (at2[3, 5] > 0.05) {"The p-value for the experimental error is non-significant which implies that plot to plot variation is low. Thus, fewer plots with more subsamples could be used."}}

r if(show.anova.sub) {if (at2[3, 5] < 0.05) {"The p-value for the experimental error is significant which implies that plot to plot variation is high. Thus, more plots could be used."}}

{{i+1}}.3. Assumptions

Don't forget the assumptions of the model. It is supposed that the errors are independent with a normal distribution and with the same variance for all the r trt.lab.s. The following plots can help you evaluate this:

par(mfrow = c(1, 2))
plot(model, which = 1)
plot(model, which = 2)

Any trend in the residuals in the left plot would violate the assumption of independence while a trend in the variability of the residuals --for instance a funnel shape-- suggests heterogeneity of variances. Departures from the theoretical normal line on the right plot are symptoms of lack of normality.

{{i+1}}.4. `r trt.lab.c` means

r if(at[1, 5] < 0.05 | mc) {paste("Below are the sorted means for each", trt.lab, "using the Fisher's Least Significant Difference method and the multiple comparisons method of Tukey, both at the 5% level. Letters indicate if there are significant differences.")} else {paste("Because the effect of", trt.lab.s, "was not significant in the ANOVA, multiple comparison tests are not presented. The means of your", trt.lab.s, "are:")}

r if (at[1, 5] < 0.05 | mc) {paste0("### ", {{i+1}}, ".4.1. LSD test")}

if (at[1, 5] < 0.05 | mc) {
  means <- dfr[, y.est]
  agricolae::LSD.test(means, dfr[, trt], at[3, 1], at[3, 3])$groups
}

r if (at[1, 5] < 0.05 | mc) {paste0("### ", {{i+1}}, ".4.2. Tukey test")}

if (at[1, 5] < 0.05 | mc)
  agricolae::HSD.test(means, dfr[, trt], at[3, 1], at[3, 3])$groups

if (at[1, 5] > 0.05 & !mc)
    tapply(dfr[, y.est], dfr[, trt], mean)

{{i+1}}.5. Variance components

Below are the variance components for this model, under the assumption that r trt.lab.s and blocks are random. Here the model is fitted using REML and missing values are not estimated.

y <- dfr[, vars[i]]
fg <- dfr[, trt]
fr <- dfr[, rep]
ff <- as.formula(y ~ (1|fg) + (1|fr))
model <- lme4::lmer(ff)
vc <- data.frame(lme4::VarCorr(model))
vc[vc[, 1] == "fg", 1] <- trt
vc[vc[, 1] == "fr", 1] <- rep
rownames(vc) <- vc[, 1]
vc <- vc[, c(4, 5)]
colnames(vc) <- c("Variance", "Std.Dev.")
vc