r i = {{i}}

# Check design

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

# Fit a model for assumptions plots

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

# Estimate missing values

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

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

# Get anova table with estimated missing values

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

# CV

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

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

r if (lc$nmis == 0) {"There are no missing values for this trait; 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 genotypes

tapply(dfr[, trait.est], dfr[, geno], mean)

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

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

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

tapply(dfr[, trait.est], list(dfr[, geno], 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 genotypes and environments. The following plots help to evaluate this 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 genotypes 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 genotypes 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 genotypes 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) {"In the next two sections the least significance difference method and the multiple comparisons method of Tukey are used to evaluate differences among genotypes, both at the 5% level. However take into account that differences among genotypes can be obscured by the interaction effects, and that in the case of strong interaction, differences among genotypes must depend on the specific environments."} else {"Because the effect of genotypes was not significant in the ANOVA, multiple comparison tests are not presented."}

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

if (at[1, 5] < 0.05)
  agricolae::LSD.test(dfr[, trait.est], dfr[, geno], at[5, 1], at[5, 3])$groups

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

if (at[1, 5] < 0.05)
  agricolae::HSD.test(dfr[, trait.est], dfr[, geno], at[5, 1], at[5, 3])$groups

r if (at[1, 5] < 0.05) {paste0("### ", {{i+1}}, ".2.5. Variance components estimation")} else {paste0("### ", {{i+1}}, ".2.3. Variance components estimation")}

Under the assumption that all the factors (genotypes, environments, and blocks) have random effects, below it is shown the variance components estimation. Here the model is estimated by REML (Restricted Maximum Likelihood) and the original data without the estimation of missing values is used.

y <- dfr[, traits[i]]
fg <- dfr[, geno]
fe <- dfr[, env]
fr <- dfr[, rep]
ff <- as.formula(y ~ (1|fg) + (1|fg:fe) + (1|fe/fr))
model.reml <- lme4::lmer(ff)
vc <- data.frame(lme4::VarCorr(model.reml))
vg <- vc[vc[, 1] == "fg", 4]
vgxe <- vc[vc[, 1] == "fg:fe", 4]
vr <- vc[vc[, 1] == "Residual", 4]
vc[vc[, 1] == "fg", 1] <- geno
vc[vc[, 1] == "fe", 1] <- env
vc[vc[, 1] == "fg:fe", 1] <- paste0(geno, ":", env)
vc[vc[, 1] == "fr:fe", 1] <- paste0(rep, "(", env, ")")
rownames(vc) <- vc[, 1]
vc <- vc[, c(4, 5)]
colnames(vc) <- c("Variance", "Std.Dev.")
h2 <- vg / (vg + vgxe / lc$nl[2] + vr / lc$nl[2] / lc$nrep) * 100
vc

With these variance estimates, the broad sense heritability results r paste0(format(h2, digits = 4), "%").

{{i+1}}.3. Stability analysis

r if (at[4, 5] > 0.05 | lc$nl[2] <= 2) {"This analysis is not shown because:"}

r if (at[4, 5] > 0.05) {"- Interaction is non significant."} r if (lc$nl[2] <= 2) {"- There are only 2 environments. At least 3 are needed."}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"Because interaction is significant a stability analysis is presented."}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {paste0("### ", {{i+1}}, ".3.1. AMMI")}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### AMMI biplots"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2) {
  ammimodel <- suppressWarnings(ammi(trait.est, geno, env, rep, dfr))
  plot(ammimodel, bp.type = 1)
}
if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  plot(ammimodel, bp.type = 2)

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### Interaction effects"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ammimodel$Interaction_effects

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### PC-values for genotypes"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ammimodel$PC_values_genotypes

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### PC-values for environments"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ammimodel$PC_values_environments

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### PC contributions"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ammimodel$Contribution_PCs

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {paste0("### ", {{i+1}}, ".3.2. GGE")}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### GGE biplots"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2) {
  ggemodel <- suppressWarnings(ammi(trait.est, geno, env, rep, dfr, method = "gge"))
  plot(ggemodel, bp.type = 1)
}
if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  plot(ggemodel, bp.type = 2)

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### PC-values for genotypes"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ggemodel$PC_values_genotypes

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### PC-values for environments"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ggemodel$PC_values_environments

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### PC contributions"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  ggemodel$Contribution_PCs

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {paste0("### ", {{i+1}}, ".3.3. Regression Stability Analysis")}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  rsamodel <- suppressWarnings(rsa(traits[i], geno, env, rep, dfr))

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### ANOVA"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  rsamodel$ANOVA

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### Stability measures for genotypes"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  rsamodel$Stability_for_genotypes

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"Here: \n - a is the linear regression intercept, \n - b is the linear regression slope, \n - se is the slope standard error, \n - MSe is the error mean square, \n - MSentry is the variance of the means, and \n - MSinter is the variance of the interaction effects. \n The same is shown in the next section for each environment."}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### Stability measures for environments"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  rsamodel$Stability_for_environments

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {paste0("### ", {{i+1}}, ".3.4. Tai stability analysis")}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### Tai plot"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2) {
  taimodel <- suppressWarnings(tai(trait.est, geno, env, rep, dfr))
  plot(taimodel)
}

r if (at[4, 5] < 0.05 & lc$nl[2] > 2) {"#### Tai alpha and lambda values"}

if (at[4, 5] < 0.05 & lc$nl[2] > 2)
  taimodel$Tai_values


CIP-RIU/hidap documentation built on April 30, 2021, 9:21 p.m.