#' Genotype plus genotype-by-environment model
#' @description
#' `r badge('stable')`
#'
#' Produces genotype plus genotype-by-environment model based on a multi-environment
#' trial dataset containing at least the columns for genotypes, environments and one
#' response variable or a two-way table.
#'
#'
#' @param .data The dataset containing the columns related to Environments, Genotypes
#' and the response variable(s).
#' @param env The name of the column that contains the levels of the environments.
#' @param gen The name of the column that contains the levels of the genotypes.
#' @param resp The response variable(s). To analyze multiple variables in a
#' single procedure a vector of variables may be used. For example `resp
#' = c(var1, var2, var3)`. Select helpers are also supported.
#' @param centering The centering method. Must be one of the `'none | 0'`, for no
#' centering; `'global | 1'`, for global centered (E+G+GE); `'environment | 2'` (default),
#' for environment-centered (G+GE); or `'double | 3'`, for double centered (GE).
#' A biplot cannot be produced with models produced without centering.
#' @param scaling The scaling method. Must be one of the `'none | 0'` (default), for no scaling;
#' or `'sd | 1'`, where each value is divided by the standard deviation of its corresponding
#' environment (column). This will put all environments roughly he same rang of values.
#'
#' @param svp The method for singular value partitioning. Must be one of the `'genotype | 1'`,
#' (The singular value is entirely partitioned into the genotype eigenvectors, also called row
#' metric preserving); `'environment | 2'`, default, (The singular value is entirely partitioned into the
#' environment eigenvectors, also called column metric preserving); or `'symmetrical | 3'`
#' (The singular value is symmetrically partitioned into the genotype and the environment eigenvectors
#' This SVP is most often used in AMMI analysis and other biplot analysis, but it is not ideal for
#' visualizing either the relationship among genotypes or that among the environments).
#' @param by One variable (factor) to compute the function by. It is a shortcut
#' to [dplyr::group_by()].This is especially useful, for example,
#' when the researcher want to produce GGE biplots for each level of a
#' categorical variable. In this case, an object of class gge_grouped is
#' returned.
#' @param ... Arguments passed to the function
#' [impute_missing_val()] for imputation of missing values in case
#' of unbalanced data.
#'
#' @return The function returns a list of class `gge` containing the following objects
#'
#' * **coordgen** The coordinates for genotypes for all components.
#'
#' * **coordenv** The coordinates for environments for all components.
#'
#' * **eigenvalues** The vector of eigenvalues.
#'
#' * **totalvar** The overall variance.
#'
#' * **labelgen** The name of the genotypes.
#'
#' * **labelenv** The names of the environments.
#'
#' * **labelaxes** The axes labels.
#'
#' * **ge_mat** The data used to produce the model (scaled and centered).
#'
#' * **centering** The centering method.
#'
#' * **scaling** The scaling method.
#'
#' * **svp** The singular value partitioning method.
#'
#' * **d** The factor used to generate in which the ranges of genotypes and environments
#' are comparable when singular value partitioning is set to 'genotype' or 'environment'.
#' * **grand_mean** The grand mean of the trial.
#' * **mean_gen** A vector with the means of the genotypes.
#' * **mean_env** A vector with the means of the environments.
#' * **scale_var** The scaling vector when the scaling method is `'sd'`.
#' @md
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @references Yan, W., and M.S. Kang. 2003. GGE biplot analysis: a graphical tool for breeders,
#' geneticists, and agronomists. CRC Press.
#'
#' @export
#' @examples
#' \donttest{
#' library(metan)
#' mod <- gge(data_ge, ENV, GEN, GY)
#' plot(mod)
#'
#' # GGE model for all numeric variables
#' mod2 <- gge(data_ge2, ENV, GEN, resp = everything())
#' plot(mod2, var = "ED")
#'
#' # If we have a two-way table with the mean values for
#' # genotypes and environments
#'
#' table <- make_mat(data_ge, GEN, ENV, GY) %>% round(2)
#' table
#' make_long(table) %>%
#' gge(ENV, GEN, Y) %>%
#' plot()
#'}
gge <- function(.data,
env,
gen,
resp,
centering = "environment",
scaling = "none",
svp = "environment",
by = NULL,
...) {
if (!missing(by)){
if(length(as.list(substitute(by))[-1L]) != 0){
stop("Only one grouping variable can be used in the argument 'by'.\nUse 'group_by()' to pass '.data' grouped by more than one variable.", call. = FALSE)
}
.data <- group_by(.data, {{by}})
}
if(is_grouped_df(.data)){
results <-
.data %>%
doo(gge,
env = {{env}},
gen = {{gen}},
resp = {{resp}},
centering = centering,
scaling = scaling,
svp = svp,
...)
return(set_class(results, c("tbl_df", "gge_group", "tbl", "data.frame")))
} else{
factors <-
.data %>%
select({{env}}, {{gen}}) %>%
mutate(across(everything(), as.factor))
vars <- .data %>% select({{resp}}, -names(factors))
vars %<>% select_numeric_cols()
factors %<>% set_names("ENV", "GEN")
listres <- list()
nvar <- ncol(vars)
for (var in 1:nvar) {
ge_mat <-
factors %>%
mutate(Y = vars[[var]]) %>%
make_mat(GEN, ENV, Y) %>%
as.matrix()
if(has_na(ge_mat)){
ge_mat <- impute_missing_val(ge_mat, verbose = FALSE, ...)$.data
warning("Data imputation used to fill the GxE matrix", call. = FALSE)
}
grand_mean <- mean(ge_mat)
mean_env <- colMeans(ge_mat)
mean_gen <- rowMeans(ge_mat)
scale_val <- apply(ge_mat, 2, sd)
labelgen <- rownames(ge_mat)
labelenv <- colnames(ge_mat)
if (any(is.na(ge_mat))) {
stop("missing data in input data frame")
}
if (any(apply(ge_mat, 2, is.numeric) == FALSE)) {
stop("not all columns are of class 'numeric'")
}
if (!(centering %in% c("none", "environment", "global", "double") |
centering %in% 0:3)) {
warning(paste("Centering method", centering, "not found; defaulting to environment centered"))
centering <- "environment"
}
if (!(svp %in% c("genotype", "environment", "symmetrical") |
svp %in% 1:3)) {
warning(paste("svp method", svp, "not found; defaulting to column metric preserving"))
svp <- "environment"
}
if (!(scaling %in% c("none", "sd") | scaling %in% 0:1)) {
warning(paste("scaling method", scaling, "not found; defaulting to no scaling"))
sd <- "none"
}
labelaxes <- paste("PC", 1:ncol(diag(svd(ge_mat)$d)), sep = "")
# Centering
if (centering == 1 | centering == "global") {
ge_mat <- ge_mat - mean(ge_mat)
}
if (centering == 2 | centering == "environment") {
ge_mat <- sweep(ge_mat, 2, colMeans(ge_mat))
}
if (centering == 3 | centering == "double") {
grand_mean <- mean(ge_mat)
mean_env <- colMeans(ge_mat)
mean_gen <- rowMeans(ge_mat)
for (i in 1:nrow(ge_mat)) {
for (j in 1:ncol(ge_mat)) {
ge_mat[i, j] <- ge_mat[i, j] + grand_mean - mean_env[j] -
mean_gen[i]
}
}
}
# Scaling
if (scaling == 1 | scaling == "sd") {
ge_mat <- sweep(ge_mat, 2, apply(ge_mat, 2, sd), FUN = "/")
}
# Singular value partitioning
if (svp == 1 | svp == "genotype") {
coordgen <- svd(ge_mat)$u %*% diag(svd(ge_mat)$d)
coordenv <- svd(ge_mat)$v
d1 <- (max(coordenv[, 1]) - min(coordenv[, 1]))/(max(coordgen[,
1]) - min(coordgen[, 1]))
d2 <- (max(coordenv[, 2]) - min(coordenv[, 2]))/(max(coordgen[,
2]) - min(coordgen[, 2]))
coordenv <- coordenv/max(d1, d2)
}
if (svp == 2 | svp == "environment") {
coordgen <- svd(ge_mat)$u
coordenv <- svd(ge_mat)$v %*% diag(svd(ge_mat)$d)
d1 <- (max(coordgen[, 1]) - min(coordgen[, 1]))/(max(coordenv[,
1]) - min(coordenv[, 1]))
d2 <- (max(coordgen[, 2]) - min(coordgen[, 2]))/(max(coordenv[,
2]) - min(coordenv[, 2]))
coordgen <- coordgen/max(d1, d2)
}
if (svp == 3 | svp == "symmetrical") {
coordgen <- svd(ge_mat)$u %*% diag(sqrt(svd(ge_mat)$d))
coordenv <- svd(ge_mat)$v %*% diag(sqrt(svd(ge_mat)$d))
}
eigenvalues <- svd(ge_mat)$d
totalvar <- round(as.numeric(sum(eigenvalues^2)), 2)
varexpl <- round(as.numeric((eigenvalues^2/totalvar) * 100),
2)
if (svp == "genotype" | svp == "environment") {
d <- max(d1, d2)
} else {
d <- NULL
}
tmp <- structure(
list(coordgen = coordgen, coordenv = coordenv, eigenvalues = eigenvalues,
totalvar = totalvar, varexpl = varexpl, labelgen = labelgen,
labelenv = labelenv, labelaxes = labelaxes, ge_mat = ge_mat,
centering = centering, scaling = scaling, svp = svp,
d = d, grand_mean = grand_mean, mean_gen = mean_gen,
mean_env = mean_env, scale_val = scale_val),
class = "gge")
if (nvar > 1) {
listres[[paste(names(vars[var]))]] <- tmp
} else {
listres[[paste(names(vars[var]))]] <- tmp
}
}
return(structure(listres, class = "gge"))
}
}
#' Create GGE, GT or GYT biplots
#'
#' Produces a ggplot2-based GGE-GT-GYT biplot based on a model fitted with the
#' functions [gge()], [gtb()], and [gytb()].
#'
#'@param x An object with classes `gge` `gtb`, or `gytb`.
#'@param var The variable to plot (useful for `gge` objects. Defaults to
#' `var = 1` the first variable of `x`.
#'@param type The type of biplot to produce.
#' 1. Basic biplot.
#' 2. Mean performance vs. stability (gge biplots) or the The Average Tester
#' Coordination view for genotype-trait and genotype-yield*trait biplots.
#'
#' 3. Which-won-where.
#' 4. Discriminativeness vs. representativeness.
#' 5. Examine an environment (or trait/yield*trait combination).
#'
#' 6. Ranking environments (or trait/yield*trait combination).
#' 7. Examine a genotype.
#' 8. Ranking genotypes.
#' 9. Compare two genotypes.
#' 10. Relationship among environments (or trait/yield*trait combination).
#' @param repel If `TRUE` (default), the text labels repel away from each
#' other and away from the data points.
#' @param repulsion Force of repulsion between overlapping text labels. Defaults
#' to `1`.
#' @param max_overlaps Exclude text labels that overlap too many things.
#' Defaults to 20.
#' @param sel_env,sel_gen The name of the environment (or trait/yield*trait
#' combination) and genotype to examine when `type = 5` and `type = 7`,
#' respectively. Must be a string which matches a environment or genotype
#' label.
#' @param sel_gen1,sel_gen2 The name of genotypes to compare between when
#' `type = 9`. Must be a string present in the genotype's name.
#' @param shape.gen,shape.env The shape for genotype and environment indication
#' in the biplot. Defaults to `shape.gen = 21` (circle) for genotypes and
#' `shape.env = 23` (rhombus) for environments. Values must be between
#' `21-25`: `21` (circle), `22` (square), `23` (rhombus),
#' `24` (up triangle), and `25` (low triangle).
#' @param line.type.gen The line type to highlith the genotype's vectors.
#' Defaults to `line.type.gen == "dotted`.
#' @param size.shape The size of the shape (both for genotypes and
#' environments). Defaults to `2.2`.
#' @param size.shape.win The size of the shape for winners genotypes when
#' `type = 3`. Defaults to `3.2`.
#' @param size.stroke,col.stroke The width and color of the border,
#' respectively. Default to `size.stroke = 0.3` and `col.stroke =
#' "black"`. The size of the shape will be `size.shape + size.stroke`
#' @param col.gen,col.env,col.line Color for genotype/environment labels and for
#' the line that passes through the biplot origin. Defaults to `col.gen =
#' 'blue'`, `col.env = 'forestgreen'`, and `col.line =
#' 'forestgreen'`.
#' @param col.alpha The alpha value for the color. Defaults to `1`. Values
#' must be between `0` (full transparency) to `1` (full color).
#' @param col.circle,col.alpha.circle The color and alpha values for the circle
#' lines. Defaults to `'gray'` and `0.4`, respectively.
#' @param leg.lab The labs of legend. Defaults to `NULL` is `c('Env', 'Gen')`.
#' @param size.text.gen,size.text.env,size.text.lab The size of the text for
#' genotypes, environments and labels, respectively.
#' @param size.text.win The text size to use for winner genotypes where
#' `type = 3` and for the two selected genotypes where `type = 9`.
#' Defaults to 4.5.
#' @param size.line The size of the line in biplots (Both for segments and circles).
#' @param axis_expand multiplication factor to expand the axis limits by to
#' enable fitting of labels. Defaults to 1.2
#' @param title Logical values (Defaults to `TRUE`) to include
#' automatically generated information in the plot such as singular value
#' partitioning, scaling and centering.
#' @param plot_theme The graphical theme of the plot. Default is
#' `plot_theme = theme_metan()`. For more details, see
#' [ggplot2::theme()].
#' @param ... Currently not used.
#' @return A ggplot2-based biplot.
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @references Yan, W., and M.S. Kang. 2003. GGE biplot analysis: a graphical
#' tool for breeders, geneticists, and agronomists. CRC Press.
#' @method plot gge
#' @importFrom ggforce geom_arc
#' @export
#' @return An object of class `gg, ggplot`.
#' @examples
#' \donttest{
#' library(metan)
#' mod <- gge(data_ge, ENV, GEN, GY)
#' plot(mod)
#' plot(mod,
#' type = 2,
#' col.gen = 'blue',
#' col.env = 'red',
#' size.text.gen = 2)
#' }
plot.gge <- function(x,
var = 1,
type = 1,
repel = TRUE,
repulsion = 1,
max_overlaps = 20,
sel_env = NA,
sel_gen = NA,
sel_gen1 = NA,
sel_gen2 = NA,
shape.gen = 21,
shape.env = 23,
line.type.gen = "dotted",
size.shape = 2.2,
size.shape.win = 3.2,
size.stroke = 0.3,
col.stroke = "black",
col.gen = "blue",
col.env = "forestgreen",
col.line = "forestgreen",
col.alpha = 1,
col.circle = "gray",
col.alpha.circle = 0.5,
leg.lab = NULL,
size.text.gen = 3.5,
size.text.env = 3.5,
size.text.lab = 12,
size.text.win = 4.5,
size.line = 0.5,
axis_expand = 1.2,
title = TRUE,
plot_theme = theme_metan(),
...) {
if(is.null(leg.lab)){
leg.lab <- case_when(
has_class(x, "gytb") ~ c("Y-Trait", "Gen"),
has_class(x, "gtb") ~ c("Trait", "Gen"),
TRUE ~ c("Env", "Gen")
)
} else{
leg.lab <- leg.lab
}
model <- x[[var]]
if (!has_class(model, c("gge", "gtb", "gytb"))) {
stop("The model must be of class 'gge'")
}
coord_gen <- model$coordgen[, c(1, 2)]
coord_env <- model$coordenv[, c(1, 2)]
varexpl <- model$varexpl[c(1, 2)]
labelgen <- model$labelgen
labelenv <- model$labelenv
labelaxes <- model$labelaxes[c(1, 2)]
Data <- model$ge_mat
centering <- model$centering
svp <- model$svp
scaling <- model$scaling
ngen <- nrow(coord_gen)
nenv <- nrow(coord_env)
if (centering == "none") {
stop("It is not possible to create a GGE biplot with a model produced without centering")
}
plotdata <- data.frame(rbind(data.frame(coord_gen,
type = "genotype",
label = labelgen),
data.frame(coord_env,
type = "environment",
label = labelenv)))
colnames(plotdata)[1:2] <- c("d1", "d2")
xlim <- c(min(plotdata$d1 * axis_expand),
max(plotdata$d1 * axis_expand))
ylim <- c(min(plotdata$d2 * axis_expand),
max(plotdata$d2 * axis_expand))
if (which(c(diff(xlim), diff(ylim)) == max(c(diff(xlim), diff(ylim)))) == 1) {
xlim1 <- xlim
ylim1 <- c(ylim[1] - (diff(xlim) - diff(ylim))/2,
ylim[2] + (diff(xlim) - diff(ylim))/2)
}
if (which(c(diff(xlim), diff(ylim)) == max(c(diff(xlim), diff(ylim)))) == 2) {
ylim1 <- ylim
xlim1 <- c(xlim[1] - (diff(ylim) - diff(xlim))/2,
xlim[2] + (diff(ylim) - diff(xlim))/2)
}
# Base plot
P1 <-
ggplot(data = plotdata, aes(x = d1, y = d2, group = "type")) +
scale_color_manual(values = c(col.env, col.gen)) +
scale_size_manual(values = c(size.text.env, size.text.gen)) +
xlab(paste(labelaxes[1], " (", round(varexpl[1], 2), "%)", sep = "")) +
ylab(paste(labelaxes[2]," (", round(varexpl[2], 2), "%)", sep = "")) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
scale_x_continuous(limits = xlim1, expand = c(0, 0)) +
scale_y_continuous(limits = ylim1, expand = c(0, 0)) +
coord_fixed() +
plot_theme %+replace%
theme(axis.text = element_text(size = size.text.lab, colour = "black"),
axis.title = element_text(size = size.text.lab, colour = "black"))
# Basic plot
if (type == 1) {
if(has_class(model, c("gtb", "gytb"))){
P2 <- P1 +
geom_segment(xend = 0,
yend = 0,
size = size.line,
aes(col = type),
data = plotdata,
show.legend = FALSE)
} else{
P2 <- P1 +
geom_segment(xend = 0,
yend = 0,
size = size.line,
col = alpha(col.line, col.alpha),
data = subset(plotdata, type == "environment"))
}
P2 <- P2 +
geom_segment(xend = 0,
yend = 0,
size = size.line,
col = alpha(col.line, col.alpha),
data = subset(plotdata, type == "environment")) +
geom_point(aes(d1, d2, fill = type, shape = type),
size = size.shape,
stroke = size.stroke,
color = col.stroke,
alpha = col.alpha) +
scale_shape_manual(labels = leg.lab,
values = c(shape.env, shape.gen)) +
scale_fill_manual(labels = leg.lab,
values = c(col.env, col.gen)) +
{if(repel)geom_text_repel(aes(col = type,
label = label,
size = type),
show.legend = FALSE,
col = c(rep(col.gen, ngen), rep(col.env, nenv)),
force = repulsion,
max.overlaps = max_overlaps)} +
{if(!repel)geom_text(aes(col = type,
label = label,
size = type),
show.legend = FALSE,
col = c(rep(col.gen, ngen), rep(col.env, nenv)),
hjust = "outward",
vjust = "outward")} +
plot_theme %+replace%
theme(axis.text = element_text(size = size.text.lab, colour = "black"),
axis.title = element_text(size = size.text.lab, colour = "black"))
if (title == TRUE) {
ggt <- case_when(
has_class(model, "gytb") ~ "GYT Biplot",
has_class(model, "gtb") ~ "GT Biplot",
TRUE ~ "GGE biplot"
) %>%
ggtitle()
}
}
# Mean vs. stability
if (type == 2) {
med1 <- mean(coord_env[, 1])
med2 <- mean(coord_env[, 2])
x1 <- NULL
for (i in 1:nrow(Data)) {
x <- solve(matrix(c(-med2, med1, med1, med2), nrow = 2),
matrix(c(0, med2 * coord_gen[i, 2] + med1 * coord_gen[i, 1]), ncol = 1))
x1 <- rbind(x1, t(x))
}
plotdata$x1_x <- NA
plotdata$x1_x[plotdata$type == "genotype"] <- x1[, 1]
plotdata$x1_y <- NA
plotdata$x1_y[plotdata$type == "genotype"] <- x1[, 2]
P2 <- P1 +
geom_segment(aes(xend = x1_x, yend = x1_y),
color = col.gen,
linetype = line.type.gen,
size = size.line,
data = subset(plotdata, type == "genotype")) +
geom_abline(intercept = 0,
slope = med2/med1,
color = col.line,
size = size.line) +
geom_abline(intercept = 0,
slope = -med1/med2,
color = col.line,
size = size.line) +
geom_segment(x = 0,
y = 0,
xend = med1,
yend = med2,
arrow = arrow(length = unit(0.3, "cm")),
size = size.line,
color = col.line) +
geom_text(aes(color = type,
label = label,
size = type),
show.legend = FALSE)
if (title == TRUE) {
ggt <- case_when(
has_class(model, "gytb") ~ "Average tester coordination view of the GYT biplot",
has_class(model, "gtb") ~ "Average tester coordination form of the GT biplot",
TRUE ~ "Mean vs. Stability"
) %>%
ggtitle()
}
}
# Which-won-where
if (type == 3) {
indice <- c(grDevices::chull(coord_gen[, 1], coord_gen[, 2]))
polign <- data.frame(coord_gen[indice, ])
indice <- c(indice, indice[1])
segs <- NULL
limx <- layer_scales(P1)$x$limits
limy <- layer_scales(P1)$y$limits
i <- 1
while (is.na(indice[i + 1]) == FALSE) {
m <- (coord_gen[indice[i], 2] - coord_gen[indice[i + 1], 2])/(coord_gen[indice[i], 1] - coord_gen[indice[i + 1], 1])
mperp <- -1/m
c2 <- coord_gen[indice[i + 1], 2] - m * coord_gen[indice[i + 1], 1]
xint <- -c2/(m - mperp)
xint <- ifelse(xint < 0, min(coord_env[, 1], coord_gen[, 1]), max(coord_env[, 1], coord_gen[, 1]))
yint <- mperp * xint
xprop <- ifelse(xint < 0, xint/limx[1], xint/limx[2])
yprop <- ifelse(yint < 0, yint/limy[1], yint/limy[2])
m3 <- which(c(xprop, yprop) == max(c(xprop, yprop)))
m2 <- abs(c(xint, yint)[m3])
if (m3 == 1 & xint < 0)
sl1 <- (c(xint, yint)/m2) * abs(limx[1])
if (m3 == 1 & xint > 0)
sl1 <- (c(xint, yint)/m2) * limx[2]
if (m3 == 2 & yint < 0)
sl1 <- (c(xint, yint)/m2) * abs(limy[1])
if (m3 == 2 & yint > 0)
sl1 <- (c(xint, yint)/m2) * limy[2]
segs <- rbind(segs, sl1)
i <- i + 1
}
rownames(segs) <- NULL
colnames(segs) <- NULL
segs <- data.frame(segs)
colnames(polign) <- c("X1", "X2")
winners <- plotdata[plotdata$type == "genotype", ][indice[-1], ] %>% add_cols(win = "yes")
others <- anti_join(plotdata, winners, by = "label") %>% add_cols(win = "no")
df_winners <- rbind(winners, others)
P2 <- P1 +
geom_polygon(data = polign,
aes(x = X1, y = X2),
fill = NA,
col = col.gen,
size = size.line) +
geom_segment(data = segs,
aes(x = X1, y = X2),
xend = 0,
yend = 0,
linetype = line.type.gen,
color = col.gen,
size = size.line) +
geom_point(data = subset(df_winners, win == "no"),
aes(d1, d2, fill = type, shape = type),
size = size.shape,
stroke = size.stroke,
alpha = col.alpha,
color = col.stroke,
show.legend = FALSE) +
{if(repel)geom_text_repel(data = subset(df_winners, win == "no"),
aes(col = type, label = label, size = type),
show.legend = FALSE,
force = repulsion,
max.overlaps = max_overlaps)} +
{if(!repel)geom_text(data = subset(df_winners, win == "no"),
aes(col = type, label = label, size = type),
show.legend = FALSE,
hjust = "outward",
vjust = "outward")} +
geom_point(data = subset(df_winners, win == "yes"),
aes(d1, d2, fill = type, shape = type),
size = size.shape.win,
stroke = size.stroke,
color = col.stroke,
alpha = col.alpha) +
{if(repel)geom_text_repel(data = subset(df_winners, win == "yes"),
aes(col = type, label = label),
show.legend = FALSE,
fontface = "bold",
size = size.text.win,
force = repulsion,
max.overlaps = max_overlaps)}
{if(!repel)geom_text(data = subset(df_winners, win == "yes"),
aes(col = type, label = label),
show.legend = FALSE,
fontface = "bold",
size = size.text.win,
hjust = "outward",
vjust = "outward")}
if("genotype" %in% others$type){
P2 <-
P2 +
scale_shape_manual(labels = leg.lab, values = c(shape.env, shape.gen)) +
scale_fill_manual(labels = leg.lab, values = c(col.env, col.gen))
} else{
P2 <-
suppressMessages(
P2 +
scale_shape_manual(labels = leg.lab[c(2,1)], values = c(shape.env, shape.gen)) +
scale_fill_manual(labels = leg.lab[c(2,1)], values = c(col.env, col.gen)) +
scale_color_manual(labels = leg.lab[c(2,1)], values = c(col.env, col.gen)) +
scale_size_manual(labels = leg.lab[c(2,1)], values = c(size.text.env, size.text.gen))
)
}
if (title == TRUE) {
ggt <- ggtitle("Which-won-where pattern")
ggt <- case_when(
has_class(x, "gytb") ~ "The which-won-where view of the GYT biplot",
has_class(x, "gtb") ~ "The which-won-where view of the GT biplot",
TRUE ~ "Which-won-where view of the GGE biplot"
) %>%
ggtitle()
}
}
# Discrimination vs. representativeness
if (type == 4) {
circles <- data.frame(x0 = 0,
y0 = 0,
start = 0,
end = pi * 2,
radio = 1:5 * max((max(coord_env[1, ]) - min(coord_env[1, ])),
(max(coord_env[2, ]) - min(coord_env[2, ])))/10)
P2 <- P1 +
geom_arc(data = circles,
aes(r = radio,
x0 = x0,
y0 = y0,
start = start,
end = end),
color = col.circle,
alpha = col.alpha.circle,
size = size.line,
inherit.aes = F,
na.rm = TRUE) +
geom_segment(xend = 0,
yend = 0,
x = mean(coord_env[, 1]),
y = mean(coord_env[, 2]),
arrow = arrow(ends = "first", length = unit(0.1, "inches")),
size = size.line,
color = col.line) +
geom_abline(intercept = 0,
slope = mean(coord_env[, 2])/mean(coord_env[, 1]),
color = col.line,
size = size.line) +
geom_point(aes(d1, d2, fill = type, shape = type),
size = size.shape,
stroke = size.stroke,
color = col.stroke,
alpha = col.alpha) +
scale_shape_manual(labels = leg.lab, values = c(shape.env, shape.gen)) +
scale_fill_manual(labels = leg.lab, values = c(col.env, col.gen)) +
geom_segment(data = subset(plotdata, type == "environment"),
xend = 0,
yend = 0,
color = alpha(col.line, col.alpha),
size = size.line) +
{if(repel)geom_text_repel(aes(col = type, label = label, size = type),
show.legend = FALSE,
color = c(rep(col.gen, ngen), rep(col.env, nenv)),
force = repulsion,
max.overlaps = max_overlaps)} +
{if(!repel)geom_text(aes(col = type, label = label, size = type),
show.legend = FALSE,
color = c(rep(col.gen, ngen), rep(col.env, nenv)),
hjust = "outward",
vjust = "outward")}
if (title == TRUE) {
ggt <- ggtitle("Discriminativeness vs. representativeness")
}
}
# Examine an environment
if (type == 5) {
if (!sel_env %in% labelenv) {
stop(paste("The environment", sel_env, "is not in the list of environment labels"))
}
venvironment <- labelenv == sel_env
x1 <- NULL
for (i in 1:nrow(Data)) {
x <- solve(matrix(c(-coord_env[venvironment, 2],
coord_env[venvironment, 1],
coord_env[venvironment, 1],
coord_env[venvironment, 2]),
nrow = 2),
matrix(c(0, coord_env[venvironment, 1] * coord_gen[i, 1] +
coord_env[venvironment, 2] * coord_gen[i, 2]), ncol = 1))
x1 <- rbind(x1, t(x))
}
plotdata$x1_x <- NA
plotdata$x1_x[plotdata$type == "genotype"] <- x1[, 1]
plotdata$x1_y <- NA
plotdata$x1_y[plotdata$type == "genotype"] <- x1[, 2]
P2 <- P1 + geom_segment(data = subset(plotdata, type == "genotype"),
aes(xend = x1_x, yend = x1_y),
col = col.gen,
linetype = line.type.gen) +
geom_abline(slope = coord_env[venvironment, 2]/coord_env[venvironment, 1],
intercept = 0,
color = alpha(col.line, col.alpha),
size = size.line) +
geom_abline(slope = -coord_env[venvironment, 1]/coord_env[venvironment, 2],
intercept = 0,
col = alpha(col.line, col.alpha),
size = size.line) +
geom_segment(data = subset(plotdata, type == "environment" & label == sel_env),xend = 0,
yend = 0,
col = alpha(col.line, col.alpha),
size = size.line,
arrow = arrow(ends = "first", length = unit(0.5, "cm"))) +
geom_text(data = subset(plotdata, type == "genotype"),
aes(label = label),
show.legend = FALSE,
color = col.gen,
size = size.text.gen) +
geom_text(data = subset(plotdata, type == "environment" & label %in% sel_env),
aes(label = label),
show.legend = FALSE,
col = col.env,
size = size.text.env,
fontface = "bold",
hjust = "outward",
vjust = "outward") +
geom_point(data = subset(plotdata, type == "environment" & label == sel_env),
aes(d1, d2),
shape = 1,
color = col.stroke,
size = 4)
if (title == TRUE) {
ggt <- ggtitle(paste("Environment:", sel_env))
}
}
# Ranking environments
if (type == 6) {
med1 <- mean(coord_env[, 1])
med2 <- mean(coord_env[, 2])
mod <- max((coord_env[, 1]^2 + coord_env[, 2]^2)^0.5)
xcoord <- sign(med1) * (mod^2/(1 + med2^2/med1^2))^0.5
ycoord <- (med2/med1) * xcoord
circles <- data.frame(x0 = xcoord,
y0 = ycoord,
start = 0,
end = pi * 2,
radio = 1:8 * ((xcoord - med1)^2 + (ycoord - med2)^2)^0.5/3)
P2 <- P1 + geom_abline(intercept = 0,
slope = med2/med1,
col = col.line,
size = size.line) +
geom_abline(intercept = 0,
slope = -med1/med2,
color = col.line,
size = size.line) +
geom_arc(data = circles,
aes(r = radio,
x0 = x0,
y0 = y0,
start = start,
end = end),
col = col.circle,
alpha = col.alpha.circle,
inherit.aes = F,
na.rm = TRUE) +
geom_segment(x = 0,
y = 0,
xend = xcoord,
yend = ycoord,
arrow = arrow(length = unit(0.15, "inches")),
size = size.line,
color = col.line) +
geom_point(fill = col.env,
shape = shape.env,
size = size.shape,
stroke = size.stroke,
color = col.stroke,
alpha = col.alpha,
data = subset(plotdata, type == "environment")) +
{if(repel)geom_text_repel(data = subset(plotdata, type == "environment"),
aes(label = label),
size = size.text.env,
show.legend = FALSE,
col = col.env,
force = repulsion,
max.overlaps = max_overlaps)} +
{if(!repel)geom_text(data = subset(plotdata, type == "environment"),
aes(label = label),
size = size.text.env,
show.legend = FALSE,
col = col.env,
hjust = "outward",
vjust = "outward")} +
geom_point(aes(xcoord, ycoord),
shape = 1,
size = 4,
color = col.stroke)
if (title == TRUE) {
ggt <- ggtitle("Ranking Environments")
}
}
# Examine a genotype
if (type == 7) {
if (!sel_gen %in% labelgen) {
stop(paste("The genotype", sel_gen, "is not in the list of genotype labels"))
}
vgenotype <- labelgen == sel_gen
x1 <- NULL
for (i in 1:ncol(Data)) {
x <- solve(matrix(c(-coord_gen[vgenotype, 2],
coord_gen[vgenotype, 1],
coord_gen[vgenotype, 1],
coord_gen[vgenotype, 2]),
nrow = 2),
matrix(c(0, coord_gen[vgenotype, 1] * coord_env[i, 1] +
coord_gen[vgenotype, 2] * coord_env[i, 2]),
ncol = 1))
x1 <- rbind(x1, t(x))
}
plotdata$x1_x <- NA
plotdata$x1_x[plotdata$type == "environment"] <- x1[, 1]
plotdata$x1_y <- NA
plotdata$x1_y[plotdata$type == "environment"] <- x1[, 2]
P2 <- P1 + geom_segment(data = subset(plotdata, type == "environment"),
aes(xend = x1_x, yend = x1_y),
col = col.line,
linetype = line.type.gen) +
geom_abline(slope = coord_gen[vgenotype, 2]/coord_gen[vgenotype, 1],
intercept = 0,
color = alpha(col.gen, col.alpha),
size = size.line) +
geom_abline(slope = -coord_gen[vgenotype, 1]/coord_gen[vgenotype, 2],
intercept = 0,
col = alpha(col.gen, col.alpha),
size = size.line) +
geom_segment(data = subset(plotdata, type == "genotype" & label == sel_gen),
xend = 0,
yend = 0,
col = alpha(col.gen, col.alpha),
arrow = arrow(ends = "first", length = unit(0.5, "cm")),
size = size.line) +
geom_text(data = subset(plotdata, type == "environment"),
aes(label = label),
show.legend = FALSE,
color = col.env,
size = size.text.gen) +
geom_text(data = subset(plotdata, type == "genotype" & label %in% sel_gen),
aes(label = label),
show.legend = FALSE,
color = col.gen,
size = size.text.gen,
fontface = "bold",
hjust = "outward",
vjust = "outward") +
geom_point(data = subset(plotdata, type == "genotype" & label == sel_gen),
aes(d1, d2),
shape = 1,
color = col.stroke,
size = 4)
if (title == TRUE) {
ggt <- ggtitle(paste("Genotype:", sel_gen))
}
}
# Ranking genotypes
if (type == 8) {
med1 <- mean(coord_env[, 1])
med2 <- mean(coord_env[, 2])
coordx <- 0
coordy <- 0
for (i in 1:nrow(Data)) {
x <- solve(matrix(c(-med2, med1, med1, med2), nrow = 2),
matrix(c(0, med2 * coord_gen[i, 2] + med1 * coord_gen[i, 1]), ncol = 1))
if (sign(x[1]) == sign(med1)) {
if (abs(x[1]) > abs(coordx)) {
coordx <- x[1]
coordy <- x[2]
}
}
}
circles <- data.frame(x0 = coordx,
y0 = coordy,
radio = 1:10 * ((coordx - med1)^2 + (coordy - med2)^2)^0.5/3)
P2 <- P1 +
geom_abline(intercept = 0,
slope = med2/med1,
col = col.gen,
size = size.line) +
geom_abline(intercept = 0,
slope = -med1/med2,
color = col.gen,
size = size.line) +
geom_arc(data = circles,
aes(r = radio,
x0 = x0,
y0 = y0,
start = 0,
end = pi * 2),
color = col.circle,
alpha = col.alpha.circle,
size = size.line,
inherit.aes = F) +
geom_segment(x = 0,
y = 0,
xend = coordx,
yend = coordy,
arrow = arrow(length = unit(0.15, "inches")),
size = size.line,
color = col.gen) +
geom_point(aes(d1, d2, fill = type, shape = type),
size = size.shape,
stroke = size.stroke,
color = col.stroke,
alpha = col.alpha) +
scale_shape_manual(labels = leg.lab, values = c(shape.env, shape.gen)) +
scale_fill_manual(labels = leg.lab, values = c(col.env, col.gen)) +
{if(repel)geom_text_repel(aes(col = type, label = label, size = type),
show.legend = FALSE,
color = c(rep(col.gen, ngen), rep(col.env, nenv)),
force = repulsion,
max.overlaps = max_overlaps)} +
{if(!repel)geom_text(aes(col = type, label = label, size = type),
show.legend = FALSE,
color = c(rep(col.gen, ngen), rep(col.env, nenv)),
hjust = "outward",
vjust = "outward")} +
geom_point(aes(coordx, coordy),
shape = 1,
color = col.stroke,
size = 3)
if (title == TRUE) {
ggt <- ggtitle("Ranking Genotypes")
}
}
# Compare two genotypes
if (type == 9) {
if (!sel_gen1 %in% labelgen) {
stop(paste("The genotype", sel_gen1, "is not in the list of genotype labels"))
}
if (!sel_gen2 %in% labelgen) {
stop(paste("The genotype", sel_gen2, "is not in the list of genotype labels"))
}
if (sel_gen1 == sel_gen2) {
stop(paste("It is not possible to compare a genotype to itself"))
}
vgenotype1 <- labelgen == sel_gen1
vgenotype2 <- labelgen == sel_gen2
P2 <- P1 +
geom_segment(x = plotdata$d1[plotdata$label == sel_gen1 & plotdata$type == "genotype"],
xend = plotdata$d1[plotdata$label == sel_gen2 & plotdata$type == "genotype"],
y = plotdata$d2[plotdata$label == sel_gen1 & plotdata$type == "genotype"],
yend = plotdata$d2[plotdata$label == sel_gen2 & plotdata$type == "genotype"],
col = col.gen,
size = size.line) +
geom_abline(intercept = 0,
slope = -(coord_gen[vgenotype1, 1] - coord_gen[vgenotype2, 1])/
(coord_gen[vgenotype1, 2] - coord_gen[vgenotype2, 2]),
color = col.gen,
size = size.line) +
geom_text(aes(label = label),
show.legend = FALSE,
data = subset(plotdata, type == "environment"),
col = col.env,
size = size.text.env) +
geom_text(data = subset(plotdata, type == "genotype" & !label %in% c(sel_gen1, sel_gen2)),
aes(label = label),
show.legend = FALSE,
col = alpha(col.gen, alpha = col.alpha),
size = size.text.gen) +
geom_text(data = subset(plotdata, type == "genotype" & label %in% c(sel_gen1, sel_gen2)),
aes(label = label),
show.legend = FALSE,
col = col.gen,
size = size.text.win,
hjust = "outward",
vjust = "outward") +
geom_point(data = subset(plotdata, type == "genotype" & label %in% c(sel_gen1, sel_gen2)),
shape = shape.gen,
fill = col.gen,
color = col.stroke,
size = size.shape) +
plot_theme
if (title == TRUE) {
ggt <- ggtitle(paste("Comparison of Genotype", sel_gen1, "with Genotype", sel_gen2))
}
}
# Relationship among environments
if (type == 10) {
P2 <- P1 +
geom_segment(data = subset(plotdata, type == "environment"),
xend = 0,
yend = 0,
col = alpha(col.line, col.alpha),
size = size.line) +
geom_point(data = subset(plotdata, type == "environment"),
shape = shape.env,
fill = col.env,
size = size.shape,
stroke = size.stroke,
color = col.stroke,
alpha = col.alpha) +
{if(repel)geom_text_repel(data = subset(plotdata, type == "environment"),
aes(label = label),
show.legend = FALSE,
col = col.env,
size = size.text.env,
force = repulsion,
max.overlaps = max_overlaps)} +
{if(!repel)geom_text(data = subset(plotdata, type == "environment"),
aes(label = label),
show.legend = FALSE,
col = col.env,
size = size.text.env,
hjust = "outward",
vjust = "outward")}
if (title == TRUE) {
if(any(inherits(x, "gtb"))){
ggt <- ggtitle("Relationship Among Traits")
} else{
ggt <- ggtitle("Relationship Among Environments")
}
}
}
if (title == T) {
scal_text <- ifelse(scaling == 1 | scaling == "sd", "Scaling = 1", "Scaling = 0")
cent_text <-
case_when(
centering == 1 | centering == "global" ~ "Centering = 1",
centering == 2 | centering == "environment" | centering == "trait" ~ "Centering = 2",
centering == 3 | centering == "double" ~ "Centering = 3",
FALSE ~ "No Centering"
)
svp_text <-
case_when(
svp == 1 | svp == "genotype" ~ "SVP = 1",
svp == 2 | svp == "environment" | svp == "trait" ~ "SVP = 2",
svp == 3 | svp == "symmetrical" ~ "SVP = 3"
)
annotationtxt <- paste(scal_text, ", ", cent_text, ", ", svp_text, sep = "")
P2 <- P2 + ggtitle(label = ggt, subtitle = annotationtxt)
}
return(P2)
}
#' Predict a two-way table based on GGE model
#'
#' Predict the means for a genotype-vs-environment trial based on a Genotype
#' plus Genotype-vs-Environment interaction (GGE) model.
#'
#' This function is used to predict the response variable of a two-way table
#' (for examples the yielding of g genotypes in e environments) based on GGE
#' model. This prediction is based on the number of principal components used.
#' For more details see Yan and Kang (2007).
#'
#' @param object An object of class `gge`.
#' @param naxis The the number of principal components to be used in the
#' prediction. Generally, two axis may be used. In this case, the estimated
#' values will be those shown in the biplot.
#' @param output The type of output. It must be one of the `'long'`
#' (default) returning a long-format table with the columns for environment
#' (ENV), genotypes (GEN) and response variable (Y); or `'wide'` to
#' return a two-way table with genotypes in the row, environments in the
#' columns, filled by the estimated values.
#' @param ... Currently not used.
#' @return A two-way table with genotypes in rows and environments in columns if
#' `output = "wide"` or a long format (columns ENV, GEN and Y) if
#' `output = "long"` with the predicted values by the GGE model.
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @references Yan, W., and M.S. Kang. 2003. GGE biplot analysis: a graphical
#' tool for breeders, geneticists, and agronomists. CRC Press.
#' @method predict gge
#' @export
#' @examples
#' \donttest{
#' library(metan)
#' mod <- gge(data_ge, GEN, ENV, c(GY, HM))
#' predict(mod)
#' }
#'
predict.gge <- function(object, naxis = 2, output = "wide", ...) {
if (has_class(object, "gtb")) {
stop("The object must be of class 'gge'.")
}
listres <- list()
varin <- 1
for (var in 1:length(object)) {
objectin <- object[[var]]
if (naxis > min(dim(objectin$coordenv))) {
stop("The number of principal components cannot be greater than min(g, e), in this case ",
min(dim(objectin$coordenv)))
}
# SVP
if (objectin$svp == "environment" | objectin$svp == 2) {
pred <- (objectin$coordgen[, 1:naxis] * (objectin$d)) %*%
t(objectin$coordenv[, 1:naxis])
}
if (objectin$svp == "genotype" | objectin$svp == 1) {
pred <- (objectin$coordgen[, 1:naxis] %*% t(objectin$coordenv[,
1:naxis] * (objectin$d)))
}
if (objectin$svp == "symmetrical" | objectin$svp == 3) {
pred <- (objectin$coordgen[, 1:naxis] %*% t(objectin$coordenv[,
1:naxis]))
}
# Scaling
if (objectin$scaling == "sd" | objectin$scaling == 1) {
pred <- sweep(pred, 2, objectin$scale_val, FUN = "*")
}
# Centering
if (objectin$centering == "global" | objectin$centering == 1) {
pred <- pred + objectin$grand_mean
}
if (objectin$centering == "environment" | objectin$centering ==
2) {
pred <- sweep(pred, 2, objectin$mean_env, FUN = "+")
}
if (objectin$centering == "double" | objectin$centering == 3) {
for (i in 1:nrow(pred)) {
for (j in 1:ncol(pred)) {
pred[i, j] <- pred[i, j] - objectin$grand_mean +
objectin$mean_env[j] + objectin$mean_gen[i]
}
}
}
rownames(pred) <- objectin$labelgen
colnames(pred) <- objectin$labelenv
if (output == "wide") {
temp <- as_tibble(pred, rownames = NA)
}
if (output == "long") {
temp <-
pred %>%
make_long() %>%
as_tibble()
}
listres[[paste(names(object[var]))]] <- temp
}
return(listres)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.