r i = {{i}}

{{i+1}}. Analysis for trait r traits[i]

# Check design

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

# Fit a model for assumptions plots

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

# Estimate missing values

trait.est <- paste0(traits[i], ".est")

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

# Get anova table with estimated missing values

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

# CV

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

# ANOVA with subsampling

show.anova.sub <- FALSE

if (!is.null(eu))
  if (sum(is.na(dfr2[, traits[i]])) == 0 & dim(dfr)[1] < dim(dfr2)[1]) {
    model2 <- aov(dfr2[, traits[i]] ~ dfr2[, trt] + dfr2[, rep] + dfr2[, eu] %in% dfr2[, trt])
    model2$terms[[2]] <- traits[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. 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}}.2. 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 must 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}}.3. r trt.lab.c means

r if(at[1, 5] < 0.05) {paste("Below are the sorted means for each", trt.lab, "with letters indicating if there are significant differences using the least significance difference method and the multiple comparisons method of Tukey, both at the 5% level.")} 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) {paste0("### ", {{i+1}}, ".3.1. LSD test")}

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

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

if (at[1, 5] < 0.05)
  agricolae::HSD.test(means, dfr[, trt], at[3, 1], at[3, 3])$groups
if (at[1, 5] > 0.05)
    tapply(dfr[, trait.est], dfr[, trt], mean)

r if (lc$ng < 10 & at[1, 5] < 0.05) {paste0("### ", {{i+1}}, ".3.3. Plot of means")}

r if(lc$ng < 10) {paste0("It is always good to have some visualization of the data. Because the number of" , trt.lab.s, " in your experiment is not so big, we can plot the data for each ", trt.lab, ":")}

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

{{i+1}}.4. 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[, traits[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


AGROFIMS/hagrofims documentation built on May 6, 2020, 7:43 p.m.