#' Weighting between stability and mean performance
#' @description
#' `r badge('stable')`
#'
#' This function computes the WAASY or WAASBY indexes (Olivoto et al., 2019)
#' considering different scenarios of weights for stability and mean
#' performance.
#'
#' After fitting a model with the functions [waas()] or
#' [waasb()] it is possible to compute the superiority indexes WAASY
#' or WAASBY in different scenarios of weights for stability and mean
#' performance. The number of scenarios is defined by the arguments
#' `increment`. By default, twenty-one different scenarios are computed. In
#' this case, the the superiority index is computed considering the following
#' weights: stability (waasb or waas) = 100; mean performance = 0. In other
#' words, only stability is considered for genotype ranking. In the next
#' iteration, the weights becomes 95/5 (since increment = 5). In the third
#' scenario, the weights become 90/10, and so on up to these weights become
#' 0/100. In the last iteration, the genotype ranking for WAASY or WAASBY
#' matches perfectly with the ranks of the response variable.
#'
#' @param model An object computed with [waas()], [waasb()], or [mps()].
#' @param mresp A numeric value that will be the new maximum value after
#' rescaling. By default, the variable in `resp` is rescaled so that the
#' original maximum and minimum values are 100 and 0, respectively. Let us
#' consider that for a specific trait, say, lodging incidence, lower values
#' are better. In this case, you should use `mresp = 0` to rescale the
#' response variable so that the lowest values will become 100 and the highest
#' values 0.
#' @param increment The increment in the weight ratio for stability and mean
#' performance. Se the **Details** section for more information.
#' @param saveWAASY Automatically save the WAASY values when the weight for
#' stability is `saveWAASY`.
#' @param prob The p-value for considering an interaction principal component
#' axis significant. must be multiple of `increment`. If this assumption
#' is not valid, an error will be occur.
#' @param progbar A logical argument to define if a progress bar is shown.
#' Default is `TRUE`.
#' @references
#' Olivoto, T., A.D.C. L{\'{u}}cio, J.A.G. da silva, V.S. Marchioro, V.Q. de
#' Souza, and E. Jost. 2019. Mean performance and stability in multi-environment
#' trials I: Combining features of AMMI and BLUP techniques. Agron. J.
#' \doi{10.2134/agronj2019.03.0220}
#'
#' @return An object of class `wsmp` with the following items for each
#' variable
#' * When computed with [waas()] or [waasb()].
#' * **scenarios** A list with the model for all computed scenarios.
#' * **WAASY** The values of the WAASY estimated when the weight for the
#' stability in the loop match with argument `saveWAASY`.
#' * **hetdata, hetcomb** The data used to produce the heatmaps.
#' * **Ranks** All the values of WAASY estimated in the different
#' scenarios of WAAS/GY weighting ratio.
#' * When computed with [mps()]
#' * **hetcomb** showing the rank for mean performance and stability in the different weights.
#' @md
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @seealso [resca()], [mps()], [mtmps()]
#' @export
#' @importFrom dplyr bind_cols
#' @examples
#' \donttest{
#' library(metan)
#' # using the WAASB as statistic and BLUP as mean performance
#' # the same as using waasb()
#'
#' model <- mps(data_ge2,
#' env = ENV,
#' gen = GEN,
#' rep = REP,
#' resp = PH)
#' scenarios <- wsmp(model)
#'}
wsmp <- function(model,
mresp = 100,
increment = 5,
saveWAASY = 50,
prob = 0.05,
progbar = TRUE) {
if(!class(model) %in% c("waas", "waasb", "mps")){
stop("The model must be an object of class 'waas' or 'waasb'")
}
if(class(model) %in% c("waas", "waasb")){
test <- 100%%increment == 0
test2 <- saveWAASY%%increment == 0
if (!mresp %in% c(100, 0)) {
stop("The value 'mresp' must be 0 or 100.")
}
if (test == FALSE) {
stop("The argument 'increment = ", increment, "' is invalid. Please, note that this value must result in an integer in the expression '100 / increment'. Please, consider changing the values.")
}
if (test2 == FALSE) {
stop("The argument 'saveWAASY = ", saveWAASY, "' must be divisible by 'increment' (",
increment, "). Please, consider changing the values.")
}
datain <- model
if (inherits(model, "waasb")) {
dfs <- list()
for (k in 1:length(model)) {
PesoWAAS <- 100
PesoResp <- 0
minresp <- 100 - mresp
data <- datain[[k]][["residuals"]] %>%
select(ENV, GEN, REP, Y)
nam <- names(datain[k])
Nenv <- length(unique(data$ENV))
Ngen <- length(unique(data$GEN))
minimo <- min(Nenv, Ngen) - 1
ncomb <- (100/increment) + 1
totalcomb <- ncomb * minimo
CombWAASY <- data.frame(type = matrix(".", (Ngen + Nenv), 1))
ovmean <- mean(data$Y)
WAASY.Values <- list()
initial <- 0
model <- suppressWarnings(suppressMessages(lme4::lmer(Y ~ REP %in% ENV + ENV + (1 | GEN) + (1 | GEN:ENV),
data = data)))
summ <- summary(model)
bups <- lme4::ranef(model)
blups <- data.frame(Names = rownames(bups$`GEN:ENV`))
blups <- blups %>% data.frame(do.call("rbind", strsplit(as.character(blups$Names),
":", fixed = TRUE))) %>% dplyr::select(-Names) %>%
dplyr::select(-X1, everything()) %>% dplyr::mutate(BLUPge = bups[[1]]$`(Intercept)`) %>%
dplyr::rename(Code = X2, GEN = X1) %>% dplyr::arrange(Code)
intmatrix <- by(blups[, 3], blups[, c(2, 1)], function(x) sum(x, na.rm = TRUE))
s <- svd(intmatrix)
U <- s$u[, 1:minimo]
LL <- diag(s$d[1:minimo])
V <- s$v[, 1:minimo]
Eigenvalue <-
data.frame(Eigenvalue = s$d[1:minimo]^2) %>%
add_cols(Proportion = s$d[1:minimo]^2/sum(s$d[1:minimo]^2) * 100,
Accumulated = cumsum(Proportion),
PC = paste("PC", 1:minimo, sep = "")) %>%
column_to_first(PC)
SCOREG <- U %*% LL^0.5
SCOREE <- V %*% LL^0.5
colnames(SCOREG) <- colnames(SCOREE) <- paste("PC", 1:minimo, sep = "")
MEDIAS <- mean_by(data, ENV, GEN)
MGEN <-
mean_by(data, GEN) %>%
add_cols(type = "GEN")
MGEN <- cbind(MGEN, SCOREG) %>%
rename(Code = GEN)
MENV <- mean_by(data, ENV) %>%
add_cols(type = "ENV")
MENV <- cbind(MENV, SCOREE) %>%
rename(Code = ENV)
Escores <- rbind(MGEN, MENV) %>%
column_to_first(type)
Pesos <- data.frame(Percent = Eigenvalue$Proportion)
if (progbar == TRUE) {
pb <- progress(max = totalcomb, style = 4)
}
for (k in 1:ncomb) {
WAASB <-
Escores %>%
select_cols(contains("PC")) %>%
abs() %>%
t() %>%
as.data.frame() %>%
add_cols(Percent = Pesos$Percent)
WAASAbsInicial <- Escores %>%
mutate(WAASB = sapply(WAASB[, -ncol(WAASB)], weighted.mean, w = WAASB$Percent)) %>%
group_by(type) %>%
mutate(PctResp = (mresp - minresp)/(max(Y) - min(Y)) * (Y - max(Y)) + mresp,
PctWAASB = (minresp - mresp)/(max(WAASB) - min(WAASB)) * (WAASB - max(WAASB)) + minresp,
wRes = PesoResp, wWAASB = PesoWAAS, OrResp = rank(-Y),
OrWAASB = rank(WAASB), OrPC1 = rank(abs(PC1)),
WAASBY = ((PctResp * wRes) + (PctWAASB * wWAASB))/(wRes + wWAASB),
OrWAASBY = rank(-WAASBY)) %>%
dplyr::ungroup()
inicial <- as.data.frame(WAASAbsInicial$OrWAASB)
colnames(inicial) <- paste0(minimo, "PC")
SigPC2 <- 1
for (j in 1:minimo) {
WAASB <- Escores %>%
select(contains("PC")) %>%
abs() %>%
t() %>%
as.data.frame() %>%
slice(1:SigPC2) %>%
mutate(Percent = Pesos$Percent[1:SigPC2])
WAASAbs <- Escores %>%
mutate(WAASB = sapply(WAASB[, -ncol(WAASB)], weighted.mean, w = WAASB$Percent)) %>%
group_by(type) %>%
mutate(PctResp = (mresp - minresp)/(max(Y) - min(Y)) * (Y - max(Y)) + mresp,
PctWAASB = (minresp - mresp)/(max(WAASB) - min(WAASB)) * (WAASB - max(WAASB)) + minresp,
wRes = PesoResp, wWAASB = PesoWAAS, OrResp = rank(-Y),
OrWAASB = rank(WAASB), OrPC1 = rank(abs(PC1)),
WAASBY = ((PctResp * wRes) + (PctWAASB * wWAASB))/(wRes + wWAASB),
OrWAASBY = rank(-WAASBY)) %>%
dplyr::ungroup()
results <- as.data.frame(WAASAbs$OrWAASB)
names(results) <- paste0(SigPC2, "PCA")
final <- cbind(results, inicial)
inicial <- final
SigPC2 <- SigPC2 + 1
ProcdAtua <- j
initial <- initial + 1
if (progbar == TRUE) {
run_progress(pb,
actual = initial,
text = paste("Ranks considering ", PesoResp, " for Y and ", PesoWAAS, " for WAASB", sep = ""))
}
}
initial <- initial
WAAS <- WAASAbsInicial
WAAS$type <- ifelse(WAAS$type == "GEN", "Genotype", "Environment")
CombWAASY[[sprintf("%.0f/%.0f", PesoWAAS, PesoResp)]] <- WAASAbsInicial$OrWAASBY
WAASY.Values[[paste("WAAS/GY", PesoWAAS, "/", PesoResp)]] <- data.frame(WAAS)
PesoResp <- PesoResp + increment
PesoWAAS <- PesoWAAS - increment
if (PesoWAAS + increment == saveWAASY) {
genotypes <- WAAS %>%
dplyr::filter(type == "Genotype") %>%
dplyr::select(Code, wRes, wWAASB, WAASBY) %>%
dplyr::arrange(WAASBY) %>%
mutate(Mean = ifelse(WAASBY < mean(WAASBY), "below", "above"))
}
}
Rank <- final[, -(SigPC2)]
Names <- WAASAbsInicial %>% select(type, Code, OrResp, OrPC1, OrWAASB)
Rank <- cbind(Names, Rank) %>% as_tibble()
hetdata <-
Rank %>%
dplyr::filter(type == "GEN") %>%
column_to_rownames("Code") %>%
select(contains("PCA")) %>%
as_tibble(rownames = NA)
CombWAASY %<>%
select(-type) %>%
mutate(type = Names$type, Code = Names$Code) %>%
select(type, Code, everything()) %>%
dplyr::filter(type == "GEN") %>%
column_to_rownames("Code") %>%
select(contains("/")) %>%
as_tibble(rownames = NA)
PC1 <- Pesos[1, 1]
PC2 <- Pesos[2, 1]
mean <- mean(WAAS$Y)
dfs[[paste(nam)]] <- structure(list(scenarios = WAASY.Values,
WAASY = genotypes,
hetcomb = CombWAASY,
hetdata = hetdata,
Ranks = Rank),
class = "wsmp")
}
}
if (inherits(model, "waas")) {
dfs <- list()
for (k in 1:length(model)) {
PesoWAAS <- 100
PesoResp <- 0
minresp <- 100 - mresp
data <- datain[[k]][["augment"]] %>%
select(ENV, GEN, REP, Y)
nam <- names(datain[k])
Nenv <- length(unique(data$ENV))
Ngen <- length(unique(data$GEN))
minimo <- min(Nenv, Ngen) - 1
ncomb <- (100/increment) + 1
CombWAASY <- data.frame(type = matrix(".", (Ngen + Nenv), 1))
WAASY.Values <- list()
model <- performs_ammi(data, ENV, GEN, REP, Y, verbose = FALSE)[[1]]
anova <- model$anova
PC <- model$PCA
MeansGxE <- model$MeansGxE
totalcomb <- ncomb * nrow(PC)
initial <- 0
if (progbar == TRUE) {
pb <- progress(max = totalcomb, style = 4)
}
for (k in 1:ncomb) {
Escores <- model$model
SigPC1 <- nrow(PC[which(PC[, 6] < prob), ])
Pesos <- as.data.frame(model$PCA[7][c(1:SigPC1), ])
colnames(Pesos) <- "Percent"
WAAS <-
Escores %>%
select(contains("PC")) %>%
abs() %>%
t() %>%
as.data.frame() %>%
slice(1:SigPC1) %>%
mutate(Percent = Pesos$Percent)
WAASAbsInicial <-
Escores %>%
mutate(WAAS = sapply(WAAS[, -ncol(WAAS)], weighted.mean, w = WAAS$Percent)) %>%
group_by(type) %>%
mutate(PctResp = (mresp - minresp)/(max(Y) - min(Y)) * (Y - max(Y)) + mresp,
PctWAAS = (minresp - mresp)/(max(WAAS) - min(WAAS)) * (WAAS - max(WAAS)) + minresp,
wRes = PesoResp, wWAAS = PesoWAAS, OrResp = rank(-Y),
OrWAAS = rank(WAAS), OrPC1 = rank(abs(PC1)),
WAASY = ((PctResp * wRes) + (PctWAAS * wWAAS))/(wRes + wWAAS),
OrWAASY = rank(-WAASY)) %>%
dplyr::ungroup()
inicial <- as.data.frame(WAASAbsInicial$OrWAAS)
colnames(inicial) <- paste0(SigPC1, "PCA")
SigPC2 <- 1
for (j in 1:nrow(PC)) {
WAAS <-
Escores %>%
select(contains("PC")) %>%
abs() %>%
t() %>%
as.data.frame() %>%
slice(1:SigPC2) %>%
mutate(Percent = Pesos$Percent[1:SigPC2])
WAASAbs <- Escores %>% mutate(WAAS = sapply(WAAS[, -ncol(WAAS)], weighted.mean, w = WAAS$Percent)) %>%
group_by(type) %>%
mutate(PctResp = (mresp - minresp)/(max(Y) - min(Y)) * (Y - max(Y)) + mresp,
PctWAAS = (minresp - mresp)/(max(WAAS) - min(WAAS)) * (WAAS - max(WAAS)) + minresp,
wRes = PesoResp, wWAAS = PesoWAAS, OrResp = rank(-Y),
OrWAAS = rank(WAAS), OrPC1 = rank(abs(PC1)),
WAASY = ((PctResp * wRes) + (PctWAAS * wWAAS))/(wRes + wWAAS),
OrWAASY = rank(-WAASY)) %>%
dplyr::ungroup()
results <- as.data.frame(WAASAbs$OrWAAS)
names(results) <- paste0(SigPC2, "PCA")
final <- cbind(results, inicial)
inicial <- final
SigPC2 %<>% +1
ProcdAtua <- j
initial <- initial + 1
if (progbar == TRUE) {
run_progress(pb,
actual = initial,
text = paste("Ranks considering ", PesoResp, " for GY and ", PesoWAAS, " for WAASB", sep = ""))
# pb$tick(tokens = list(what = paste("Ranks considering ", PesoResp, "% for GY and ", PesoWAAS, "% for WAAS", sep = "")))
}
}
initial <- initial
WAAS <- WAASAbsInicial
WAAS$type <- ifelse(WAAS$type == "GEN", "Genotype", "Environment")
CombWAASY[[sprintf("%.0f/%.0f", PesoWAAS, PesoResp)]] <- WAASAbsInicial$OrWAASY
WAASY.Values[[paste("WAAS/GY", PesoWAAS, "/", PesoResp)]] <- as_tibble(WAAS)
PesoResp <- PesoResp + increment
PesoWAAS <- PesoWAAS - increment
if (PesoWAAS + increment == saveWAASY) {
genotypes <- WAAS %>%
dplyr::filter(type == "Genotype") %>%
dplyr::select(Code, wRes, wWAAS, WAASY) %>%
dplyr::arrange(WAASY) %>%
mutate(Mean = ifelse(WAASY < mean(WAASY), "below", "above"))
}
}
Rank <- final[, -(SigPC2)]
Names <- WAASAbsInicial %>%
select(type, Code, OrResp, OrPC1, OrWAAS)
Rank <- cbind(Names, Rank) %>% as_tibble()
hetdata <-
Rank %>%
dplyr::filter(type == "GEN") %>%
column_to_rownames("Code") %>%
select(contains("PCA")) %>%
as_tibble(rownames = NA)
CombWAASY %<>%
select(-type) %>%
mutate(type = Names$type, Code = Names$Code) %>%
select(type, Code, everything()) %>%
dplyr::filter(type == "GEN") %>%
column_to_rownames("Code") %>%
select(contains("/")) %>%
as_tibble(rownames = NA)
PC1 <- Pesos[1, 1]
PC2 <- Pesos[2, 1]
mean <- mean(WAAS$Y)
dfs[[paste(nam)]] <- structure(list(scenarios = WAASY.Values,
WAASY = genotypes,
hetcomb = CombWAASY,
hetdata = hetdata,
Ranks = Rank),
class = "wsmp")
}
}
return(structure(dfs, class = "wsmp"))
} else{
nvar <- ncol(model$performance_res)
ncomb <- (100/increment) + 1
wstab <- 100
wmp <- 0
vars <- list()
for (i in 1:nvar){
stab <- model$stability_res[[i]]
mper <- model$performance_res[[i]]
tmp <- list()
for(j in 1:ncomb){
tmp[[paste0(wstab, "/", wmp)]] <- as.numeric(((stab * wstab) + (mper * wmp)) / (wstab + wmp))
wmp <- wmp + increment
wstab <- wstab - increment
}
bind <- as.data.frame(bind_cols(tmp))
rownames(bind) <- rownames(model$performance_res)
vars[[paste(colnames(model$stability_res[[i]]))]][["hetcomb"]] <- bind
wstab <- 100
wmp <- 0
}
return(structure(vars, class = "wsmp"))
}
}
#' Plot heat maps with genotype ranking
#'
#' Plot heat maps with genotype ranking in two ways.
#'
#' The first type of heatmap shows the genotype ranking depending on the number
#' of principal component axis used for estimating the WAASB index. The second type of heatmap shows the
#' genotype ranking depending on the WAASB/GY ratio. The ranks obtained with a
#' ratio of 100/0 considers exclusively the stability for the genotype ranking.
#' On the other hand, a ratio of 0/100 considers exclusively the productivity
#' for the genotype ranking. Four clusters of genotypes are shown by label colors (red) unproductive and
#' unstable genotypes; (blue) productive, but unstable genotypes; (black) stable, but
#' unproductive genotypes; and (green), productive and stable genotypes.
#'
#' @param x An object returned by the function `wsmp`.
#' @param var The variable to plot. Defaults to `var = 1` the first
#' variable of `x`.
#' @param type `1 = Heat map Ranks`: this graphic shows the genotype
#' ranking considering the WAASB index estimated with different numbers of
#' Principal Components; `2 = Heat map WAASY-GY ratio`: this graphic
#' shows the genotype ranking considering the different combinations in the
#' WAASB/GY ratio.
#' @param y.lab The label of y axis. Default is 'Genotypes'.
#' @param x.lab The label of x axis. Default is 'Number of axes'.
#' @param size.lab The size of the labels.
#' @param ... Currently not used.
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @method plot wsmp
#' @export
#' @return An object of class `gg`.
#' @examples
#' \donttest{
#' library(metan)
#' model <- waasb(data_ge2,
#' env = ENV,
#' gen = GEN,
#' rep = REP,
#' resp = PH) %>%
#' wsmp()
#'
#' p1 <- plot(model)
#' p2 <- plot(model, type = 2)
#' arrange_ggplot(p1, p2, ncol = 1)
#' }
#'
plot.wsmp <- function(x,
var = 1,
type = 1,
y.lab = NULL,
x.lab = NULL,
size.lab = 12,
...) {
nam_dat <- ifelse(type == 1, "hetdata", "hetcomb")
dat <- x[[var]][[nam_dat]]
if(is.null(dat)){
warning("object `x` seems to be computed with `mps()`. Switching to `type = 2`.", call. = FALSE)
dat <- x[[var]][["hetcomb"]]
}
grps <-
dist(dat) %>%
hclust() %>%
cutree(4) %>%
as.data.frame() %>%
rownames_to_column("GEN") %>%
set_names("GEN", "GROUP") %>%
arrange(GROUP) %>%
add_cols(color =
case_when(GROUP == 1 ~ "green",
GROUP == 2 ~ "red",
GROUP == 3 ~ "blue",
GROUP == 4 ~ "black"))
data <-
dat %>%
rownames_to_column("GEN") %>%
metan::as_factor(1) %>%
pivot_longer(-GEN) %>%
add_row_id()
if(type == 1){
x.lab = ifelse(is.null(x.lab), paste0("Number of IPCAs"), x.lab)
y.lab = ifelse(is.null(y.lab), paste0("Genotypes"), y.lab)
}
if(type == 2){
x.lab = ifelse(is.null(x.lab), paste0("WAASB/GY ratio"), x.lab)
y.lab = ifelse(is.null(y.lab), paste0("Genotypes"), y.lab)
}
if(type == 1){
p <-
ggplot(data, aes(reorder(name, desc(row_id)), factor(GEN, levels = grps$GEN), fill= value))
}
if(type == 2){
p <-
ggplot(data, aes(reorder(name, row_id), factor(GEN, levels = grps$GEN), fill= value))
}
p <- p +
geom_raster() +
theme_metan() +
theme(legend.position = "right")+
scale_y_discrete(expand = expansion(mult = c(0,0)))+
scale_x_discrete(expand = expansion(mult = c(0,0)))+
scale_fill_viridis_c() +
guides(fill = guide_colourbar(label = TRUE,
draw.ulim = TRUE,
draw.llim = TRUE,
frame.colour = "black",
ticks = TRUE,
nbin = 10,
label.position = "right",
barwidth = 1.3,
barheight = 10,
direction = 'vertical'))+
theme(axis.text.x = element_text(size = size.lab, angle = 90, hjust = 1, vjust = 0.5))+
theme(axis.text.y = element_text(size = size.lab, colour = grps$color),
axis.title = element_text(size = size.lab)) +
labs(x = x.lab, y = y.lab)
return(p)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.