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

r i = {{i}}

# Check design

lc <- ck.f(dfr, vars[i], c(trt, env), rep)

# Fit a model for assumptions plots

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

# Estimate missing values

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

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

# Get anova table with estimated missing values

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

# CV

cv <- (at[5, 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, ANOVA, regression stability analysis and Tai sections.")

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

{{i+1}}.1.1. Means by treatments

tapply(dfr[, y.est], dfr[, trt], mean)

{{i+1}}.1.2. Means by environments

tapply(dfr[, y.est], dfr[, env], mean)

{{i+1}}.1.3. Means by treatments and environments

tapply(dfr[, y.est], list(dfr[, trt], dfr[, env]), 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 treatments and environments. The following plots help to evaluate these assumptions:

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

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

For this analysis it is assumed that treatments and environments have fixed effects and that the blocks are random.

at

The coefficient of variation for this experiment is r format(cv, digits = 4)%. The p-values for the model are: r format(at[1, 5], digits = 4) for treatments r if (at[1, 5] < 0.05) {" which is significant at the 5% level, "} else {" which is not significant at the 5% level, "} r format(at[2, 5], digits = 4) for environments r if (at[2, 5] < 0.05) {" which is significant at the 5% level,"} else {" which is not significant at the 5% level,"} and r format(at[4, 5], digits = 4) for the treatments by environments interaction r if (at[4, 5] < 0.05) {" which is significant at the 5% level."} else {" which is not significant at the 5% level."}

r if (at[1, 5] < 0.05 | at[4, 5] < 0.05) {paste0("### ", {{i+1}}, ".2.3. LSD test for treatments")}

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

if (at[4, 5] < 0.05) {
  for (j in 1:lc$nl[2]) {
    tmp <- dfr[dfr[, env] == unique(dfr[, env])[j], ]
    y <- tmp[, y.est]
    print(paste("LSD test for environment", unique(dfr[, env])[j]))
    print(agricolae::LSD.test(y, tmp[, trt], at[5, 1], at[5, 3])$groups)
  }
}

r if (at[1, 5] < 0.05 | at[4, 5] < 0.05) {paste0("### ", {{i+1}}, ".2.4. Tukey test for treatments")}

if (at[1, 5] < 0.05 & at[4, 5] > 0.05) {
  y <- dfr[, y.est]
  agricolae::HSD.test(y, dfr[, trt], at[5, 1], at[5, 3])$groups
}

if (at[4, 5] < 0.05) {
  for (j in 1:lc$nl[2]) {
    tmp <- dfr[dfr[, env] == unique(dfr[, env])[j], ]
    y <- tmp[, y.est]
    print(paste("HSD Tukey test for environment", unique(dfr[, env])[j]))
    print(agricolae::HSD.test(y, tmp[, trt], at[5, 1], at[5, 3])$groups)
  }
}

r if (at[2, 5] < 0.05 & at[1, 5] > 0.05 & at[4, 5] > 0.05) {paste0("### ", {{i+1}}, ".2.3. LSD test for environments")} r if (at[4, 5] < 0.05) {paste0("### ", {{i+1}}, ".2.5. LSD test for environments")} r if (at[1, 5] < 0.05 & at[2, 5] < 0.05 & at[4, 5] > 0.05) {paste0("### ", {{i+1}}, ".2.5. LSD test for environments")}

if (at[2, 5] < 0.05 & at[4, 5] > 0.05) {
  y <- dfr[, y.est]
  agricolae::LSD.test(y, dfr[, env], at[3, 1], at[3, 3])$groups
}

if (at[4, 5] < 0.05) {
  for (j in 1:lc$nl[1]) {
    tmp <- dfr[dfr[, trt] == unique(dfr[, trt])[j], ]
    y <- tmp[, y.est]
    print(paste("LSD test for treatments", unique(dfr[, trt])[j]))
    print(agricolae::LSD.test(y, tmp[, env], at[3, 1], at[3, 3])$groups)
  }
}

r if (at[2, 5] < 0.05 & at[1, 5] > 0.05 & at[4, 5] > 0.05) {paste0("### ", {{i+1}}, ".2.4. Tukey test for environments")} r if (at[4, 5] < 0.05) {paste0("### ", {{i+1}}, ".2.6. Tukey test for environments")} r if (at[1, 5] < 0.05 & at[2, 5] < 0.05 & at[4, 5] > 0.05) {paste0("### ", {{i+1}}, ".2.6. Tukey test for environments")}

if (at[2, 5] < 0.05 & at[4, 5] > 0.05) {
  y <- dfr[, y.est]
  agricolae::HSD.test(y, dfr[, env], at[3, 1], at[3, 3])$groups
}

if (at[4, 5] < 0.05) {
  for (j in 1:lc$nl[1]) {
    tmp <- dfr[dfr[, trt] == unique(dfr[, trt])[j], ]
    y <- tmp[, y.est]
    print(paste("HSD Tukey test for treatments", unique(dfr[, trt])[j]))
    print(agricolae::HSD.test(y, tmp[, env], at[3, 1], at[3, 3])$groups)
  }
}