R/waasb.R

Defines functions waasb

Documented in waasb

#' Weighted Average of Absolute Scores
#' @description
#' `r badge('stable')`
#'
#' Compute the Weighted Average of Absolute Scores (Olivoto et al., 2019) for
#' quantifying the stability of *g* genotypes conducted in *e*
#' environments using linear mixed-effect models.
#'
#' The weighted average of absolute scores is computed considering all
#' Interaction Principal Component Axis (IPCA) from the Singular Value Decomposition (SVD) of the
#' matrix of genotype-environment interaction (GEI) effects generated by a linear
#' mixed-effect model, as follows:
#' \loadmathjax
#' \mjsdeqn{WAASB_i = \sum_{k = 1}^{p} |IPCA_{ik} \times EP_k|/ \sum_{k = 1}^{p}EP_k}
#'
#' where \mjseqn{WAASB_i} is the weighted average of absolute scores of the
#' *i*th genotype; \mjseqn{IPCA_{ik}} is the score of the *i*th genotype
#' in the *k*th Interaction Principal Component Axis (IPCA); and \mjseqn{EP_k} is the explained variance of the *k*th
#' IPCA for *k = 1,2,..,p*, considering \mjseqn{p = min(g - 1; e - 1)}.
#'
#' The nature of the effects in the model is
#' chosen with the argument `random`. By default, the experimental design
#' considered in each environment is a randomized complete block design. If
#' `block` is informed, a resolvable alpha-lattice design (Patterson and
#' Williams, 1976) is implemented. The following six models can be fitted
#' depending on the values of `random` and `block` arguments.
#'   *  **Model 1:** `block = NULL` and `random = "gen"` (The
#'   default option). This model considers a Randomized Complete Block Design in
#'   each environment assuming genotype and genotype-environment interaction as
#'   random effects. Environments and blocks nested within environments are
#'   assumed to fixed factors.
#'
#'   *  **Model 2:** `block = NULL` and `random = "env"`. This
#'   model considers a Randomized Complete Block Design in each environment
#'   treating environment, genotype-environment interaction, and blocks nested
#'   within environments as random factors. Genotypes are assumed to be fixed
#'   factors.
#'
#'   *  **Model 3:** `block = NULL` and `random = "all"`. This
#'   model considers a Randomized Complete Block Design in each environment
#'   assuming a random-effect model, i.e., all effects (genotypes, environments,
#'   genotype-vs-environment interaction and blocks nested within environments)
#'   are assumed to be random factors.
#'
#'   *  **Model 4:** `block` is not `NULL` and `random =
#'   "gen"`. This model considers an alpha-lattice design in each environment
#'   assuming genotype, genotype-environment interaction, and incomplete blocks
#'   nested within complete replicates as random to make use of inter-block
#'   information (Mohring et al., 2015). Complete replicates nested within
#'   environments and environments are assumed to be fixed factors.
#'
#'   *  **Model 5:** `block` is not `NULL` and `random =
#'   "env"`. This model considers an alpha-lattice design in each environment
#'   assuming genotype as fixed. All other sources of variation (environment,
#'   genotype-environment interaction, complete replicates nested within
#'   environments, and incomplete blocks nested within replicates) are assumed
#'   to be random factors.
#'
#'   *  **Model 6:** `block` is not `NULL` and `random =
#'   "all"`. This model considers an alpha-lattice design in each environment
#'   assuming all effects, except the intercept, as random factors.
#'
#' @param .data The dataset containing the columns related to Environments,
#'   Genotypes, replication/block and 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 rep The name of the column that contains the levels of the
#'   replications/blocks.
#' @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)`.
#' @param block Defaults to `NULL`. In this case, a randomized complete
#'   block design is considered. If block is informed, then an alpha-lattice
#'   design is employed considering block as random to make use of inter-block
#'   information, whereas the complete replicate effect is always taken as
#'   fixed, as no inter-replicate information was to be recovered (Mohring et
#'   al., 2015).
#'@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 compute the indexes by mega-environments. In
#'  this case, an object of class waasb_grouped is returned.
#'  [mtsi()] can then be used to compute the mtsi index within each
#'  mega-environment.
#' @param mresp  The new maximum value after rescaling the response variable. By
#'   default, all variables in `resp` are rescaled so that de maximum value
#'   is 100 and the minimum value is 0 (i.e., `mresp = NULL`). It must be a
#'   character vector of the same length of `resp` if rescaling is assumed
#'   to be different across variables, e.g., if for the first variable smaller
#'   values are better and for the second one, higher values are better, then
#'   `mresp = c("l, h")` must be used. Character value of length 1 will be
#'   recycled with a warning message.
#' @param wresp The weight for the response variable(s) for computing the WAASBY
#'   index. By default, all variables in `resp` have equal weights for mean
#'   performance and stability (i.e., `wresp = 50`). It must be a numeric
#'   vector of the same length of `resp` to assign different weights across
#'   variables, e.g., if for the first variable equal weights for mean
#'   performance and stability are assumed and for the second one, a higher
#'   weight for mean performance (e.g. 65) is assumed, then `wresp = c(50,
#'   65)` must be used. Numeric value of length 1 will be recycled with a
#'   warning message.
#' @param random The effects of the model assumed to be random. Defaults to
#'   `random = "gen"`. See **Details** to see the random effects
#'   assumed depending on the experimental design of the trials.
#' @param prob The probability for estimating confidence interval for BLUP's
#'   prediction.
#' @param ind_anova Logical argument set to `FALSE`. If `TRUE` an
#'   within-environment ANOVA is performed.
#' @param verbose Logical argument. If `verbose = FALSE` the code will run
#'   silently.
#' @param ... Arguments passed to the function
#'   [impute_missing_val()] for imputation of missing values in the
#'   matrix of BLUPs for genotype-environment interaction, thus allowing the
#'   computation of the WAASB index.
#' @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.
#' 111:2949-2960.
#' \doi{10.2134/agronj2019.03.0220}
#'
#' Mohring, J., E. Williams, and H.-P. Piepho. 2015. Inter-block information: to
#' recover or not to recover it? TAG. Theor. Appl. Genet. 128:1541-54.
#' \doi{10.1007/s00122-015-2530-0}
#'
#' Patterson, H.D., and E.R. Williams. 1976. A new class of resolvable
#' incomplete block designs. Biometrika 63:83-92.
#'
#'
#' @return An object of class `waasb` with the following items for each
#'   variable:
#'
#' * **individual** A within-environments ANOVA considering a
#'   fixed-effect model.
#'
#' * **fixed** Test for fixed effects.
#'
#' * **random** Variance components for random effects.
#'
#' * **LRT** The Likelihood Ratio Test for the random effects.
#'
#' * **model** A tibble with the response variable, the scores of all
#' IPCAs, the estimates of Weighted Average of Absolute Scores, and WAASBY (the
#' index that considers the weights for stability and mean performance in the
#' genotype ranking), and their respective ranks.
#'
#' * **BLUPgen** The random effects and estimated BLUPS for genotypes (If
#' `random = "gen"` or `random = "all"`)
#'
#' * **BLUPenv** The random effects and estimated BLUPS for environments,
#' (If `random = "env"` or `random = "all"`).
#'
#' * **BLUPint** The random effects and estimated BLUPS of all genotypes in
#' all environments.
#'
#' * **PCA** The results of Principal Component Analysis with the
#' eigenvalues and explained variance of the matrix of genotype-environment
#' effects estimated by the linear fixed-effect model.
#'
#' * **MeansGxE** The phenotypic means of genotypes in the environments.
#'
#' * **Details** A list summarizing the results. The following information
#' are shown: `Nenv`, the number of environments in the analysis;
#' `Ngen` the number of genotypes in the analysis; `mresp` The value
#' attributed to the highest value of the response variable after rescaling it;
#' `wresp` The weight of the response variable for estimating the WAASBY
#' index. `Mean` the grand mean; `SE` the standard error of the mean;
#' `SD` the standard deviation. `CV` the coefficient of variation of
#' the phenotypic means, estimating WAASB, `Min` the minimum value observed
#' (returning the genotype and environment), `Max` the maximum value
#' observed (returning the genotype and environment); `MinENV` the
#' environment with the lower mean, `MaxENV` the environment with the
#' larger mean observed, `MinGEN` the genotype with the lower mean,
#' `MaxGEN` the genotype with the larger.
#'
#' * **ESTIMATES** A tibble with the genetic parameters (if `random =
#' "gen"` or `random = "all"`) with the following columns: `Phenotypic
#' variance` the phenotypic variance; `Heritability` the broad-sense
#' heritability; `GEr2` the coefficient of determination of the interaction
#' effects; `h2mg` the heritability on the mean basis;
#' `Accuracy` the selective accuracy; `rge` the genotype-environment
#' correlation; `CVg` the genotypic coefficient of variation; `CVr`
#' the residual coefficient of variation; `CV ratio` the ratio between
#' genotypic and residual coefficient of variation.
#'
#'  * **residuals** The residuals of the model.
#'
#'  * **formula** The formula used to fit the model.
#' @md
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @seealso [mtsi()] [waas()]
#'   [get_model_data()] [plot_scores()]
#' @export
#' @examples
#' \donttest{
#' library(metan)
#' #===============================================================#
#' # Example 1: Analyzing all numeric variables assuming genotypes #
#' # as random effects with equal weights for mean performance and #
#' # stability                                                     #
#' #===============================================================#
#'model <- waasb(data_ge,
#'               env = ENV,
#'               gen = GEN,
#'               rep = REP,
#'               resp = everything())
#' # Distribution of random effects (first variable)
#' plot(model, type = "re")
#'
#' # Genetic parameters
#' get_model_data(model, "genpar")
#'
#'
#'
#' #===============================================================#
#' # Example 2: Analyzing variables that starts with "N"           #
#' # assuming environment as random effects with higher weight for #
#' # response variable (65) for the three traits.                  #
#' #===============================================================#
#'
#'model2 <- waasb(data_ge2,
#'                env = ENV,
#'                gen = GEN,
#'                rep = REP,
#'                random = "env",
#'                resp = starts_with("N"),
#'                wresp = 65)
#'
#'
#' # Get the index WAASBY
#' get_model_data(model2, what = "WAASBY")
#'
#' # Plot the scores (response x WAASB)
#' plot_scores(model2, type = 3)
#'
#' #===============================================================#
#' # Example 3: Analyzing GY and HM assuming a random-effect model.#
#' # Smaller values for HM and higher values for GY are better.    #
#' # To estimate WAASBY, higher weight for the GY (60%) and lower  #
#' # weight for HM (40%) are considered for mean performance.      #
#' #===============================================================#
#'
#' model3 <- waasb(data_ge,
#'                 env = ENV,
#'                 gen = GEN,
#'                 rep = REP,
#'                 resp = c(GY, HM),
#'                 random = "all",
#'                 mresp = c("h, l"),
#'                 wresp = c(60, 40))
#'
#' # Get Likelihood-ratio test
#' get_model_data(model3, "lrt")
#'
#' # Get the random effects
#' get_model_data(model3, what = "ranef")
#' }
#'
waasb <- function(.data,
                  env,
                  gen,
                  rep,
                  resp,
                  block = NULL,
                  by = NULL,
                  mresp = NULL,
                  wresp = NULL,
                  random = "gen",
                  prob = 0.05,
                  ind_anova = FALSE,
                  verbose = TRUE,
                  ...) {
    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)){
        if(!missing(block)){
            results <-
                .data %>%
                doo(waasb,
                    env = {{env}},
                    gen = {{gen}},
                    rep = {{rep}},
                    resp = {{resp}},
                    block = {{block}},
                    mresp = mresp,
                    wresp = wresp,
                    random = random,
                    prob = prob,
                    ind_anova = ind_anova,
                    verbose = verbose,
                    ...)
        } else{
            results <-
                .data %>%
                doo(waasb,
                    env = {{env}},
                    gen = {{gen}},
                    rep = {{rep}},
                    resp = {{resp}},
                    mresp = mresp,
                    wresp = wresp,
                    random = random,
                    prob = prob,
                    ind_anova = ind_anova,
                    verbose = verbose,
                    ...)
        }
        return(set_class(results, c("tbl_df",  "waasb_group", "tbl",  "data.frame")))
    }
    if (!random %in% c("env", "gen", "all")) {
        stop("The argument 'random' must be one of the 'gen', 'env', or 'all'.", call. = FALSE)
    }
    if(is.numeric(mresp)){
        stop("Using a numeric vector in 'mresp' is deprecated as of metan 1.9.0. use 'h' or 'l' instead.\nOld code: 'mresp = c(100, 100, 0)'.\nNew code: 'mresp = c(\"h, h, l\")", call. = FALSE)
    }
    block_test <- missing(block)
    if(!missing(block)){
        factors  <- .data %>%
            select({{env}},
                   {{gen}},
                   {{rep}},
                   {{block}}) %>%
            mutate(across(everything(), as.factor))
    } else{
        factors  <- .data %>%
            select({{env}},
                   {{gen}},
                   {{rep}}) %>%
            mutate(across(everything(), as.factor))
    }
    vars <- .data %>% select({{resp}}, -names(factors))
    vars %<>% select_numeric_cols()
    if(!missing(block)){
        factors %<>% set_names("ENV", "GEN", "REP", "BLOCK")
    } else{
        factors %<>% set_names("ENV", "GEN", "REP")
    }
    model_formula <-
        case_when(
            random == "gen" & block_test ~ paste("Y ~ ENV/REP + (1 | GEN) + (1 | GEN:ENV)"),
            random == "env" & block_test ~ paste("Y ~ GEN + (1 | ENV/REP) + (1 | GEN:ENV)"),
            random == "all" & block_test ~ paste("Y ~ (1 | GEN) + (1 | ENV/REP) + (1 | GEN:ENV)"),
            random == "gen" & !block_test ~ paste("Y ~  (1 | GEN) + ENV / REP + (1|BLOCK:(REP:ENV))  + (1 | GEN:ENV)"),
            random == "env" & !block_test ~ paste("Y ~ 0 + GEN + (1| ENV/REP/BLOCK)  + (1 | GEN:ENV)"),
            random == "all" & !block_test ~ paste("Y ~  (1 | GEN) + (1|ENV/REP/BLOCK) + (1 | GEN:ENV)")
        )
    lrt_groups <-
        strsplit(
            case_when(
                random == "gen" & block_test ~ c("COMPLETE GEN GEN:ENV"),
                random == "env" & block_test ~ c("COMPLETE REP(ENV) ENV GEN:ENV"),
                random == "all" & block_test ~ c("COMPLETE GEN REP(ENV) ENV GEN:ENV"),
                random == "gen" & !block_test ~ c("COMPLETE GEN BLOCK(ENV:REP) GEN:ENV"),
                random == "env" & !block_test ~ c("COMPLETE BLOCK(ENV:REP) REP(ENV) ENV GEN:ENV"),
                random == "all" & !block_test ~ c("COMPLETE GEN BLOCK(ENV:REP) REP(ENV) ENV GEN:ENV")
            ), " ")[[1]]
    mod1 <- random == "gen" & block_test
    mod2 <- random == "gen" & !block_test
    mod3 <- random == "env" & block_test
    mod4 <- random == "env" & !block_test
    mod5 <- random == "all" & block_test
    mod6 <- random == "all" & !block_test
    nvar <- ncol(vars)
    if (is.null(mresp)) {
        mresp <- replicate(nvar, 100)
        minresp <- 100 - mresp
    } else {
        mresp <- unlist(strsplit(mresp, split="\\s*(\\s|,)\\s*")) %>% all_lower_case()
        if(!any(mresp %in% c("h", "l", "H", "L"))){
            if(!mresp[[1]] %in% c("h", "l")){
                stop("Argument 'mresp' must have only h or l.", call. = FALSE)
            } else{
            warning("Argument 'mresp' must have only h or l. Setting mresp = ", mresp[[1]],
                    " to all the ", nvar, " variables.", call. = FALSE)
                mresp <- replicate(nvar, mresp[[1]])
            }
        }
        if (length(mresp) != nvar) {
            warning("Invalid length in 'mresp'. Setting mresp = ", mresp[[1]],
                    " to all the ", nvar, " variables.", call. = FALSE)
            mresp <- replicate(nvar, mresp[[1]])
        }

        mresp <- ifelse(mresp == "h", 100, 0)
        minresp <- 100 - mresp
    }
    if (is.null(wresp)) {
        PesoResp <- replicate(nvar, 50)
        PesoWAASB <- 100 - PesoResp
    } else {
        PesoResp <- wresp
        PesoWAASB <- 100 - PesoResp
        if (length(wresp) != nvar) {

            warning("Invalid length in 'wresp'. Setting wresp = ", wresp[[1]],
                    " to all the ", nvar, " variables.", call. = FALSE)
            PesoResp <- replicate(nvar, wresp[[1]])
            PesoWAASB <- 100 - PesoResp
        }
        if (min(wresp) < 0 | max(wresp) > 100) {
            stop("The range of the numeric vector 'wresp' must be equal between 0 and 100.")
        }
    }
    listres <- list()
    vin <- 0
    if (verbose == TRUE) {
        pb <- progress(max = nvar, style = 4)
    }
    for (var in 1:nvar) {
        data <- factors %>%
            mutate(Y = vars[[var]])
        check_labels(data)
        if(has_na(data)){
            data <- remove_rows_na(data)
            has_text_in_num(data)
        }
        if(!is_balanced_trial(data, ENV, GEN, Y) && random == "env"){
            warning("Fitting a model with unbalanced data considering genotype as fixed effect is not suggested.", call. = FALSE)
        }
        Nenv <- nlevels(data$ENV)
        Ngen <- nlevels(data$GEN)
        Nrep <- nlevels(data$REP)
        minimo <- min(Nenv, Ngen) - 1
        vin <- vin + 1
        ovmean <- mean(data$Y)
        if (minimo < 2) {
            cat("\nWarning. The analysis is not possible.")
            cat("\nThe number of environments and number of genotypes must be greater than 2\n")
        }
        if(ind_anova == TRUE){
            if(missing(block)){
                individual <- data %>% anova_ind(ENV, GEN, REP, Y)
            } else{
                individual <- data %>% anova_ind(ENV, GEN, REP, Y, block = BLOCK)
            }
        } else{
            individual = NULL
        }
        Complete <- suppressWarnings(suppressMessages(lmerTest::lmer(model_formula, data = data)))
        LRT <- suppressWarnings(suppressMessages(lmerTest::ranova(Complete, reduce.terms = FALSE) %>%
                                                     mutate(model = lrt_groups) %>%
                                                     column_to_first(model)))
        fixed <- anova(Complete)
        var_eff <-
            lme4::VarCorr(Complete) %>%
            as.data.frame() %>%
            select_cols(1, 4) %>%
            arrange(grp) %>%
            rename(Group = grp, Variance = vcov) %>%
            add_cols(Percent = (Variance / sum(Variance)) * 100)
        if(random %in% c("gen", "all")){
            GV <- as.numeric(var_eff[which(var_eff[1] == "GEN"), 2])
            IV <- as.numeric(var_eff[which(var_eff[1] == "GEN:ENV"), 2])
            RV <- as.numeric(var_eff[which(var_eff[1] == "Residual"), 2])
            FV <- sum(var_eff$Variance)
            h2g <- GV/FV
            h2mg <- GV/(GV + IV/Nenv + RV/(Nenv * Nrep))
            GEr2 <- IV/(GV + IV + RV)
            AccuGen <- sqrt(h2mg)
            rge <- IV/(IV + RV)
            CVg <- (sqrt(GV)/ovmean) * 100
            CVr <- (sqrt(RV)/ovmean) * 100
            CVratio <- CVg/CVr
            PROB <- ((1 - (1 - prob))/2) + (1 - prob)
            t <- qt(PROB, Nrep)
            Limits <- t * sqrt(((1 - AccuGen) * GV))
            genpar <- tibble(Parameters = c("Phenotypic variance", "Heritability", "GEIr2", "h2mg",
                                            "Accuracy", "rge", "CVg", "CVr", "CV ratio"),
                             Values = c(FV, h2g, GEr2, h2mg, AccuGen, rge, CVg, CVr, CVratio))
        } else{
            genpar <- NULL
        }
        bups <- lme4::ranef(Complete)
        bINT <-
            data.frame(Names = rownames(bups[["GEN:ENV"]])) %>%
            separate(Names, into = c("GEN", "ENV"), sep = ":") %>%
            add_cols(BLUPge = bups[["GEN:ENV"]][[1]]) %>%
            as_factor(1:2)
        intmatrix <- as.matrix(make_mat(bINT, GEN, ENV, BLUPge))
        if(has_na(intmatrix)){
            intmatrix <- impute_missing_val(intmatrix, verbose = verbose, ...)$.data
            warning("Data imputation used to fill the GxE matrix", call. = FALSE)
        }
        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 <- means_by(data, ENV, GEN)
        MGEN <- MEDIAS %>% means_by(GEN) %>% add_cols(type = "GEN")
        MGEN <- cbind(MGEN, SCOREG)
        MENV <- MEDIAS %>% means_by(ENV) %>% add_cols(type = "ENV")
        MENV <- cbind(MENV, SCOREE)
        MEDIAS <- suppressMessages(dplyr::mutate(MEDIAS,
                                                 envPC1 = left_join(MEDIAS, MENV %>% select(ENV, PC1))$PC1,
                                                 genPC1 = left_join(MEDIAS, MGEN %>% select(GEN, PC1))$PC1,
                                                 nominal = left_join(MEDIAS, MGEN %>% select(GEN, Y))$Y + genPC1 * envPC1))
        MGEN %<>% rename(Code = GEN)
        MENV %<>% rename(Code = ENV)
        Escores <- rbind(MGEN, MENV) %>%
            column_to_first(type)
        Pesos <- data.frame(Percent = Eigenvalue$Proportion)
        WAASB <- Escores %>%
            select(contains("PC")) %>%
            abs() %>%
            t() %>%
            as.data.frame() %>%
            mutate(Percent = Pesos$Percent)
        WAASAbs <-
            mutate(Escores, WAASB = sapply(WAASB[, -ncol(WAASB)], weighted.mean, w = WAASB$Percent)) %>%
            group_by(type) %>%
            mutate(PctResp = (mresp[vin] - minresp[vin])/(max(Y) - min(Y)) * (Y - max(Y)) + mresp[vin],
                   PctWAASB = (0 - 100)/(max(WAASB) - min(WAASB)) * (WAASB - max(WAASB)) + 0,
                   wRes = PesoResp[vin],
                   wWAASB = PesoWAASB[vin],
                   OrResp = rank(-Y),
                   OrWAASB = rank(WAASB),
                   OrPC1 = rank(abs(PC1)),
                   WAASBY = ifelse(is.na(((PctResp * wRes) + (PctWAASB * wWAASB))/(wRes + wWAASB)), PctResp, ((PctResp * wRes) + (PctWAASB * wWAASB))/(wRes + wWAASB)),
                   OrWAASBY = rank(-WAASBY)) %>%
            ungroup()
        Details <-
            rbind(ge_details(data, ENV, GEN, Y),
                  tribble(~Parameters,  ~Y,
                          "wresp", PesoResp[vin],
                          "mresp", mresp[vin],
                          "Ngen", Ngen,
                          "Nenv", Nenv)) %>%
            rename(Values = Y)
        if(mod1){
            ran_ef <- c("GEN, GEN:ENV")
            fix_ef <- c("ENV, REP(ENV)")
            data_factors <- data %>% select_non_numeric_cols()
            BLUPgen <-
                data.frame(GEN = MGEN$Code,
                           BLUPg = bups$GEN$`(Intercept)`) %>%
                add_cols(Predicted = BLUPg + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted),
                         LL = Predicted - Limits,
                         UL = Predicted + Limits) %>%
                column_to_first(Rank)
            BLUPint <-
                suppressWarnings(
                    left_join(data_factors, bINT, by = c("ENV", "GEN")) %>%
                        left_join(BLUPgen, by = "GEN") %>%
                        select(ENV, GEN, REP, BLUPg, BLUPge) %>%
                        add_cols(`BLUPg+ge` = BLUPge + BLUPg,
                                 Predicted = predict(Complete))
                )
            BLUPenv <- NULL
        } else if(mod2){
            ran_ef <- c("GEN, BLOCK(ENV:REP), GEN:ENV")
            fix_ef <- c("ENV, REP(ENV)")
            data_factors <- data %>% select_non_numeric_cols()
            BLUPgen <-
                data.frame(GEN = MGEN$Code,
                           BLUPg = bups$GEN$`(Intercept)`) %>%
                add_cols(Predicted = BLUPg + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted),
                         LL = Predicted - Limits,
                         UL = Predicted + Limits) %>%
                column_to_first(Rank)
            blupBRE <-
                data.frame(Names = rownames(bups$`BLOCK:(REP:ENV)`)) %>%
                separate(Names, into = c("BLOCK", "REP", "ENV"), sep = ":") %>%
                add_cols(BLUPbre = bups$`BLOCK:(REP:ENV)`[[1]]) %>%
                as_factor(1:3)
            BLUPint <-
                suppressWarnings(
                    left_join(data_factors, bINT, by = c("ENV", "GEN")) %>%
                        left_join(BLUPgen, by = "GEN") %>%
                        left_join(blupBRE, by = c("ENV", "REP", "BLOCK")) %>%
                        select(ENV, REP, BLOCK, GEN, BLUPg, BLUPge, BLUPbre) %>%
                        add_cols(`BLUPg+ge+bre` = BLUPge + BLUPg + BLUPbre,
                                 Predicted = `BLUPg+ge+bre` + left_join(data_factors, data %>% means_by(ENV, REP), by = c("ENV", "REP"))$Y)
                )
            BLUPenv <- NULL
        } else if (mod3){
            ran_ef <- c("REP(ENV), ENV, GEN:ENV")
            fix_ef <- c("GEN")
            data_factors <- data %>% select_non_numeric_cols()
            BLUPgen <- NULL
            BLUPenv <- data.frame(ENV = MENV$Code,
                                  BLUPe = bups$ENV$`(Intercept)`) %>%
                add_cols(Predicted = BLUPe + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted)) %>%
                column_to_first(Rank)
            blupRWE <-
                data.frame(Names = rownames(bups$`REP:ENV`)) %>%
                separate(Names, into = c("REP", "ENV"), sep = ":") %>%
                add_cols(BLUPre = bups$`REP:ENV`[[1]]) %>%
                as_factor(1:2)
            BLUPint <-
                suppressWarnings(
                    left_join(data_factors, bINT, by = c("ENV", "GEN")) %>%
                        left_join(BLUPenv, by = "ENV") %>%
                        left_join(blupRWE, by = c("ENV", "REP")) %>%
                        select(ENV, GEN, REP, BLUPe, BLUPge, BLUPre) %>%
                        add_cols(`BLUPge+e+re` = BLUPge + BLUPe + BLUPre,
                                 Predicted = `BLUPge+e+re` + left_join(data_factors, MGEN %>% select(Code, Y), by = c("GEN" = "Code"))$Y)
                )
        } else if (mod4){
            ran_ef <- c("BLOCK(ENV:REP), REP(ENV), ENV, GEN:ENV")
            fix_ef <- c("GEN")
            data_factors <- data %>% select_non_numeric_cols()
            BLUPgen <- NULL
            BLUPenv <-
                data.frame(ENV = MENV$Code,
                           BLUPe = bups$ENV$`(Intercept)`) %>%
                add_cols(Predicted = BLUPe + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted)) %>%
                column_to_first(Rank)
            blupRWE <-
                data.frame(Names = rownames(bups$`REP:ENV`)) %>%
                separate(Names, into = c("REP", "ENV"), sep = ":") %>%
                add_cols(BLUPre = bups$`REP:ENV`[[1]]) %>%
                as_factor(1:2)
            blupBRE <-
                data.frame(Names = rownames(bups$`BLOCK:(REP:ENV)`)) %>%
                separate(Names, into = c("BLOCK", "REP", "ENV")) %>%
                add_cols(BLUPbre = bups$`BLOCK:(REP:ENV)`[[1]]) %>%
                as_factor(1:3)
            genCOEF <- summary(Complete)[["coefficients"]] %>%
                as_tibble(rownames = NA) %>%
                rownames_to_column("GEN") %>%
                replace_string(GEN, pattern = "GEN") %>%
                rename(Y = Estimate) %>%
                as_factor(1)
            BLUPint <-
                suppressWarnings(
                    left_join(data_factors, bINT, by = c("ENV", "GEN")) %>%
                        left_join(BLUPenv, by = "ENV") %>%
                        left_join(blupRWE, by = c("ENV", "REP")) %>%
                        left_join(blupBRE, by = c("ENV", "REP", "BLOCK")) %>%
                        select(ENV, REP, BLOCK, GEN, BLUPe, BLUPge, BLUPre, BLUPbre) %>%
                        add_cols(`BLUPe+ge+re+bre` = BLUPge + BLUPe + BLUPre + BLUPbre,
                                 Predicted = `BLUPe+ge+re+bre` + left_join(data_factors, genCOEF, by = "GEN")$Y)
                )
        } else if (mod5){
            ran_ef <- c("GEN, REP(ENV), ENV, GEN:ENV")
            fix_ef <- c("-")
            data_factors <- data %>% select_non_numeric_cols()
            BLUPgen <-
                data.frame(GEN = MGEN$Code,
                           BLUPg = bups$GEN$`(Intercept)`) %>%
                add_cols(Predicted = BLUPg + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted),
                         LL = Predicted - Limits,
                         UL = Predicted + Limits) %>%
                column_to_first(Rank)
            BLUPenv <- data.frame(ENV = MENV$Code,
                                  BLUPe = bups$ENV$`(Intercept)`) %>%
                add_cols(Predicted = BLUPe + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted)) %>%
                column_to_first(Rank)
            blupRWE <- data.frame(Names = rownames(bups$`REP:ENV`)) %>%
                separate(Names, into = c("REP", "ENV"), sep = ":") %>%
                add_cols(BLUPre = bups$`REP:ENV`[[1]]) %>%
                arrange(ENV) %>%
                as_factor(1:2)
            BLUPint <-
                suppressWarnings(
                    left_join(data_factors, bINT, by = c("ENV", "GEN")) %>%
                        left_join(BLUPgen, by = "GEN") %>%
                        left_join(BLUPenv, by = "ENV") %>%
                        left_join(blupRWE, by = c("ENV", "REP")) %>%
                        select(GEN, ENV, REP, BLUPe, BLUPg, BLUPge, BLUPre) %>%
                        add_cols(`BLUPg+e+ge+re` = BLUPge + BLUPe + BLUPg + BLUPre,
                                 Predicted = `BLUPg+e+ge+re` + ovmean)
                )
        } else if (mod6){
            ran_ef <- c("GEN, BLOCK(ENV:REP), REP(ENV), ENV, GEN:ENV")
            fix_ef <- c("-")
            data_factors <- data %>% select_non_numeric_cols()
            BLUPgen <-
                data.frame(GEN = MGEN$Code,
                           BLUPg = bups$GEN$`(Intercept)`) %>%
                add_cols(Predicted = BLUPg + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted),
                         LL = Predicted - Limits,
                         UL = Predicted + Limits) %>%
                column_to_first(Rank)
            BLUPenv <- data.frame(ENV = MENV$Code,
                                  BLUPe = bups$ENV$`(Intercept)`) %>%
                add_cols(Predicted = BLUPe + ovmean) %>%
                arrange(-Predicted) %>%
                add_cols(Rank = rank(-Predicted)) %>%
                column_to_first(Rank)
            blupRWE <- data.frame(Names = rownames(bups$`REP:ENV`)) %>%
                separate(Names, into = c("REP", "ENV"), sep = ":") %>%
                add_cols(BLUPre = bups$`REP:ENV`[[1]]) %>%
                arrange(ENV) %>%
                as_factor(1:2)
            blupBRE <-
                data.frame(Names = rownames(bups$`BLOCK:(REP:ENV)`)) %>%
                separate(Names, into = c("BLOCK", "REP", "ENV"), sep = ":") %>%
                add_cols(BLUPbre = bups$`BLOCK:(REP:ENV)`[[1]]) %>%
                as_factor(1:3)
            BLUPint <-
                suppressWarnings(
                    left_join(data_factors, bINT, by = c("ENV", "GEN")) %>%
                        left_join(BLUPgen, by = "GEN") %>%
                        left_join(BLUPenv, by = "ENV") %>%
                        left_join(blupRWE, by = c("ENV", "REP")) %>%
                        left_join(blupBRE, by = c("ENV", "REP", "BLOCK")) %>%
                        select(GEN, ENV, REP, BLOCK, BLUPg, BLUPe, BLUPge, BLUPre, BLUPbre) %>%
                        add_cols(`BLUPg+e+ge+re+bre` = BLUPg + BLUPge + BLUPe + BLUPre + BLUPbre,
                                 Predicted = `BLUPg+e+ge+re+bre` + ovmean)
                )
        }
        residuals <- data.frame(fortify.merMod(Complete))
        residuals$reff <- BLUPint$BLUPge
        temp <- structure(list(individual = individual[[1]],
                               fixed = fixed %>% rownames_to_column("SOURCE") %>% as_tibble(),
                               random = var_eff,
                               LRT = LRT,
                               model = as_tibble(WAASAbs),
                               BLUPgen = BLUPgen,
                               BLUPenv = BLUPenv,
                               BLUPint = BLUPint,
                               PCA = as_tibble(Eigenvalue),
                               modellme = Complete,
                               MeansGxE = as_tibble(MEDIAS),
                               Details = as_tibble(Details),
                               ESTIMATES = genpar,
                               residuals = as_tibble(residuals),
                               formula = model_formula), class = "waasb")
        if (verbose == TRUE) {
            run_progress(pb,
                         actual = var,
                         text = paste("Evaluating trait", names(vars[var])))
        }
        listres[[paste(names(vars[var]))]] <- temp
    }
    if (verbose == TRUE) {
        message("Method: REML/BLUP\n", appendLF = FALSE)
        message("Random effects: ", ran_ef, "\n", appendLF = FALSE)
        message("Fixed effects: ", fix_ef, "\n", appendLF = FALSE)
        message("Denominador DF: Satterthwaite's method\n", appendLF = FALSE)
        cat("---------------------------------------------------------------------------\n")
        cat("P-values for Likelihood Ratio Test of the analyzed traits\n")
        cat("---------------------------------------------------------------------------\n")
        print.data.frame(sapply(listres, function(x){
            x$LRT[["Pr(>Chisq)"]]
        }) %>%
            as.data.frame() %>%
            add_cols(model = listres[[1]][["LRT"]][["model"]]) %>%
            column_to_first(model), row.names = FALSE, digits = 3)
        cat("---------------------------------------------------------------------------\n")
        if (length(which(unlist(lapply(listres, function(x) {
            x[["LRT"]] %>% dplyr::filter(model == "GEN:ENV") %>% pull(`Pr(>Chisq)`)
        })) > prob)) > 0) {
            cat("Variables with nonsignificant GxE interaction\n")
            cat(names(which(unlist(lapply(listres, function(x) {
                x[["LRT"]][which(x[["LRT"]][[1]] == "GEN:ENV"), 7]
            })) > prob)), "\n")
            cat("---------------------------------------------------------------------------\n")
            chek_na_waasb <- listres %>% map(~.x[["model"]] %>% has_na)
            if(any(chek_na_waasb == TRUE)){
            cat("The following traits had p-value for GE interaction = 1\n")
            cat(names(which(chek_na_waasb == TRUE)), "\n")
            cat("WAASBY value for these traits is based on mean performance only (PctResp)\n")
            cat("---------------------------------------------------------------------------\n")
            }
        } else {
            cat("All variables with significant (p < 0.05) genotype-vs-environment interaction\n")
        }
    }
    invisible(structure(listres, class = "waasb"))
}





#' Several types of residual plots
#'
#' Residual plots for a output model of class `waas` and `waasb`. Six types
#' of plots are produced: (1) Residuals vs fitted, (2) normal Q-Q plot for the
#' residuals, (3) scale-location plot (standardized residuals vs Fitted
#' Values), (4) standardized residuals vs Factor-levels, (5) Histogram of raw
#' residuals and (6) standardized residuals vs observation order. For a `waasb`
#' object, normal Q-Q plot for random effects may also be obtained declaring
#' `type = 're'`
#'
#'
#' @param x An object of class `waasb`.
#' @param var The variable to plot. Defaults to `var = 1` the first
#'   variable of `x`.
#' @param type One of the `"res"` to plot the model residuals (default),
#'   `type = 're'` to plot normal Q-Q plots for the random effects, or
#'   `"vcomp"` to create a bar plot with the variance components.
#' @param position The position adjustment when `type = "vcomp"`. Defaults
#'   to `"fill"`, which shows relative proportions at each trait by
#'   stacking the bars and then standardizing each bar to have the same height.
#'   Use `position = "stack"` to plot the phenotypic variance for each
#'   trait.
#' @param rotate Logical argument. If `rotate = TRUE` the plot is rotated,
#'   i.e., traits in y axis and value in the x axis.
#' @param conf Level of confidence interval to use in the Q-Q plot (0.95 by
#' default).
#' @param out How the output is returned. Must be one of the 'print' (default)
#' or 'return'.
#' @param n.dodge The number of rows that should be used to render the x labels.
#'   This is useful for displaying labels that would otherwise overlap.
#' @param check.overlap Silently remove overlapping labels, (recursively)
#'   prioritizing the first, last, and middle labels.
#' @param labels Logical argument. If `TRUE` labels the points outside
#' confidence interval limits.
#' @param plot_theme The graphical theme of the plot. Default is
#'   `plot_theme = theme_metan()`. For more details, see
#'   [ggplot2::theme()].
#' @param alpha The transparency of confidence band in the Q-Q plot. Must be a
#' number between 0 (opaque) and 1 (full transparency).
#' @param fill.hist The color to fill the histogram. Default is 'gray'.
#' @param col.hist The color of the border of the the histogram. Default is
#' 'black'.
#' @param col.point The color of the points in the graphic. Default is 'black'.
#' @param col.line The color of the lines in the graphic. Default is 'red'.
#' @param col.lab.out The color of the labels for the 'outlying' points.
#' @param size.line The size of the line in graphic. Defaults to 0.7.
#' @param size.text The size for the text in the plot. Defaults to 10.
#' @param width.bar The width of the bars if `type = "contribution"`.
#' @param size.lab.out The size of the labels for the 'outlying' points.
#' @param size.tex.lab The size of the text in axis text and labels.
#' @param size.shape The size of the shape in the plots.
#' @param bins The number of bins to use in the histogram. Default is 30.
#' @param which Which graphics should be plotted. Default is `which =
#' c(1:4)` that means that the first four graphics will be plotted.
#' @param ncol,nrow The number of columns and rows of the plot pannel. Defaults
#'   to `NULL`
#' @param ... Additional arguments passed on to the function
#'  [patchwork::wrap_plots()].
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @importFrom dplyr distinct_all
#' @importFrom tibble tribble
#' @method plot waasb
#' @export
#' @examples
#'\donttest{
#' library(metan)
#' model2 <- waasb(data_ge,
#'                 resp = GY,
#'                 gen = GEN,
#'                 env = ENV,
#'                 rep = REP)
#' plot(model2)
#'}

plot.waasb <- function(x,
                       var = 1,
                       type = "res",
                       position = "fill",
                       rotate = FALSE,
                       conf = 0.95,
                       out = "print",
                       n.dodge = 1,
                       check.overlap = FALSE,
                       labels = FALSE,
                       plot_theme = theme_metan(),
                       alpha = 0.2,
                       fill.hist = "gray",
                       col.hist = "black",
                       col.point = "black",
                       col.line = "red",
                       col.lab.out = "red",
                       size.line = 0.7,
                       size.text = 10,
                       width.bar = 0.75,
                       size.lab.out = 2.5,
                       size.tex.lab = 10,
                       size.shape = 1.5,
                       bins = 30,
                       which = c(1:4),
                       ncol = NULL,
                       nrow = NULL,
                       ...) {
    if(!type  %in% c("res", 're', "vcomp")){
        stop("Argument type = '", match.call()[["type"]], "' invalid. Use one of 'res', 're', or 'vcomp'", call. = FALSE)
    }
    if(type %in% c("vcomp", "re") && !class(x)  %in% c("waasb", "gamem")){
        stop("Arguments 're' and 'vcomp' valid for objects of class 'waasb' or 'gamem'. ")
    }
    if(is.numeric(var)){
        var_name <- names(x)[var]
    } else{
        var_name <- var
    }
    if(!var_name %in% names(x)){
        stop("Variable not found in ", match.call()[["x"]] , call. = FALSE)
    }
    if (type == "re" & max(which) >= 5) {
        stop("When type =\"re\", 'which' must be a value between 1 and 4")
    }
    if(type == "vcomp"){
        list <- lapply(x, function(x){
            x[["random"]] %>% select(Group, Variance)
        })
        vcomp <- suppressWarnings(
            lapply(seq_along(list),
                   function(i){
                       set_names(list[[i]], "Group", names(list)[i])
                   }) %>%
                reduce(full_join, by = "Group") %>%
                pivot_longer(-Group))
        p1 <-
            ggplot(vcomp, aes(x = name, y = value, fill = Group)) +
            geom_bar(stat = "identity",
                     position = position,
                     color = "black",
                     size = size.line,
                     width = width.bar) +
            scale_y_continuous(expand = expansion(c(0, ifelse(position == "fill", 0, 0.05))))+
            theme_bw()+
            theme(legend.position = "bottom",
                  axis.ticks = element_line(size = size.line),
                  axis.ticks.length = unit(0.2, "cm"),
                  panel.grid = element_blank(),
                  legend.title = element_blank(),
                  strip.background = element_rect(fill = NA),
                  text = element_text(size = size.text, colour = "black"),
                  axis.text = element_text(size = size.text, colour = "black")) +
            scale_x_discrete(guide = guide_axis(n.dodge = n.dodge, check.overlap = check.overlap))+
            labs(x = "Traits",
                 y = ifelse(position == "fill", "Proportion of phenotypic variance", "Phenotypic variance"))
        if(rotate == TRUE){
            p1 <- p1 + coord_flip()
        }
        return(p1)
    }

    if (type == "res") {
        x <- x[[var]]
        df <- data.frame(x$residuals)
        df$id <- rownames(df)
        df <- data.frame(df[order(df$.scresid), ])
        P <- ppoints(nrow(df))
        df$z <- qnorm(P)
        n <- nrow(df)
        Q.x <- quantile(df$.scresid, c(0.25, 0.75), na.rm = TRUE)
        Q.z <- qnorm(c(0.25, 0.75))
        b <- diff(Q.x)/diff(Q.z)
        coef <- c(Q.x[1] - b * Q.z[1], b)
        zz <- qnorm(1 - (1 - conf)/2)
        SE <- (coef[2]/dnorm(df$z)) * sqrt(P * (1 - P)/n)
        fit.value <- coef[1] + coef[2] * df$z
        df$upper <- fit.value + zz * SE
        df$lower <- fit.value - zz * SE
        df$label <- ifelse(df$.scresid > df$.scresid | df$.scresid <
                               df$lower, rownames(df), "")
        df$factors <- paste(df$ENV, df$GEN)
        # Residuals vs .fitted
        p1 <- ggplot(df, aes(.fitted, .resid)) +
            geom_point(col = col.point, size = size.shape) +
            geom_smooth(se = F, method = "loess", col = col.line) +
            geom_hline(yintercept = 0, linetype = 2, col = "gray") +
            labs(x = "Fitted values", y = "Residual") +
            ggtitle("Residual vs fitted") + plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1))
        if (labels != FALSE) {
            p1 <- p1 +
                ggrepel::geom_text_repel(aes(.fitted, .resid, label = (label)),
                                         color = col.lab.out,
                                         size = size.lab.out)
        } else {
            p1 <- p1
        }
        # normal qq
        p2 <- ggplot(df, aes(z, .scresid)) +
            geom_point(col = col.point, size = size.shape) +
            geom_abline(intercept = coef[1],
                        slope = coef[2],
                        size = 1,
                        col = col.line) +
            geom_ribbon(aes_(ymin = ~lower, ymax = ~upper),
                        alpha = 0.2) +
            labs(x = "Theoretical quantiles", y = "Sample quantiles") +
            ggtitle("Normal Q-Q") +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1))
        if (labels != FALSE) {
            p2 <- p2 + ggrepel::geom_text_repel(aes(z, .scresid, label = (label)),
                                                color = col.lab.out,
                                                size = size.lab.out)
        } else {
            p2 <- p2
        }
        # scale-location
        p3 <- ggplot(df, aes(.fitted, sqrt(abs(.resid)))) +
            geom_point(col = col.point, size = size.shape) +
            geom_smooth(se = F, method = "loess", col = col.line) +
            labs(x = "Fitted Values", y = expression(sqrt("|Standardized residuals|"))) +
            ggtitle("Scale-location") +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1))
        if (labels != FALSE) {
            p3 <- p3 + ggrepel::geom_text_repel(aes(.fitted, sqrt(abs(.resid)),
                                                    label = (label)),
                                                color = col.lab.out,
                                                size = size.lab.out)
        } else {
            p3 <- p3
        }
        # Residuals vs Factor-levels
        p4 <- ggplot(df, aes(factors, .scresid)) +
            geom_point(col = col.point, size = size.shape) +
            geom_hline(yintercept = 0, linetype = 2, col = "gray") +
            labs(x = "Factor levels", y = "Standardized residuals") +
            ggtitle("Residuals vs factor-levels") +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  panel.grid.major.x = element_blank(),
                  axis.text.x = element_text(color = "white"),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1))
        if (labels != FALSE) {
            p4 <- p4 + ggrepel::geom_text_repel(aes(factors,
                                                    .scresid, label = (label)),
                                                color = col.lab.out,
                                                size = size.lab.out)
        } else {
            p4 <- p4
        }
        # Histogram of residuals
        p5 <- ggplot(df, aes(x = .resid)) +
            geom_histogram(bins = bins,
                           colour = col.hist,
                           fill = fill.hist,
                           aes(y = ..density..)) +
            stat_function(fun = dnorm,
                          color = col.line,
                          size = 1,
                          args = list(mean = mean(df$.resid),
                                      sd = sd(df$.resid))) +
            labs(x = "Raw residuals", y = "Density") +
            ggtitle("Histogram of residuals") +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1))
        # Residuals vs order
        p6 <- ggplot(df, aes(as.numeric(id), .scresid, group = 1)) +
            geom_point(col = col.point, size = size.shape) +
            geom_line(col = col.line) +
            geom_hline(yintercept = 0,
                       linetype = 2,
                       col = col.line) +
            labs(x = "Observation order", y = "Standardized residuals") +
            ggtitle("Residuals vs observation order") +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1))
        p7 <- ggplot(df, aes(.fitted, Y)) +
            geom_point(col = col.point, size = size.shape) +
            facet_wrap(~GEN) +
            geom_abline(intercept = 0, slope = 1, col = col.line) +
            labs(x = "Fitted values", y = "Observed values") +
            ggtitle("1:1 line plot") +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  panel.grid.major.x = element_blank(),
                  panel.grid.major.y = element_blank(),
                  plot.title = element_text(size = size.tex.lab, hjust = 0, vjust = 1),
                  panel.spacing = unit(0, "cm"))
        plots <- list(p1, p2, p3, p4, p5, p6, p7)
        p1 <- wrap_plots(plots[c(which)],
                         ncol = ncol,
                         nrow = nrow,
                         ...) +
            plot_annotation(title = var_name)
        return(p1)
    }
    if (type == "re") {
        x <- x[[var]]
        blups <-
            x$BLUPint %>%
            select_cols(contains("BLUP"))
        fact <-x$BLUPint %>% select_non_numeric_cols()
        qlist <- list()
        for (i in 1:ncol(blups)) {
            df <-
                data.frame(blups[i]) %>%
                distinct_all() %>%
                add_row_id(var = "id") %>%
                arrange(across(2))
            P <- ppoints(nrow(df))
            df$z <- qnorm(P)
            n <- nrow(df)
            Q.x <- quantile(df[[2]], c(0.25, 0.75), na.rm = TRUE)
            Q.z <- qnorm(c(0.25, 0.75))
            b <- diff(Q.x)/diff(Q.z)
            coef <- c(Q.x[1] - b * Q.z[1], b)
            zz <- qnorm(1 - (1 - conf)/2)
            SE <- (coef[2]/dnorm(df$z)) * sqrt(P * (1 - P)/n)
            fit.value <- coef[1] + coef[2] * df$z
            df %<>% add_cols(upper = fit.value + zz * SE,
                             lower = fit.value - zz * SE,
                             label = ifelse(df[[2]] > upper | df[[2]] < lower, id, ""),
                             intercept = coef[1],
                             slope = coef[2],
                             var = paste(names(blups[i]))
            ) %>%
                set_names("id",    "blup", "z",     "upper", "lower", "label", "intercept", "slope", "var")
            qlist[[paste(names(blups[i]))]] <- df
        }

        df <- do.call(rbind, qlist)
        # normal qq GEI effects
        p1 <- ggplot(df, aes(z, blup)) +
            geom_point(col = col.point, size = size.shape) +
            geom_abline(aes(intercept = intercept,
                            slope = slope),
                        size = 1, col = col.line) +
            geom_ribbon(aes_(ymin = ~lower, ymax = ~upper),
                        alpha = 0.2) +
            labs(x = "Theoretical quantiles", y = "Sample quantiles")+
            facet_wrap( ~var,
                        scales = "free",
                        ncol = ncol,
                        nrow = nrow) +
            plot_theme %+replace%
            theme(axis.text = element_text(size = size.tex.lab, colour = "black"),
                  axis.title = element_text(size = size.tex.lab, colour = "black"),
                  plot.title.position = "plot",
                  plot.title = element_text(size = size.tex.lab + 1, hjust = 0, vjust = 1, face = "bold"))
        if (labels != FALSE) {
            p1 <- p1 + ggrepel::geom_text_repel(aes(z, blup, label = (label)),
                                                color = col.lab.out,
                                                size = size.lab.out)+
                labs(title = var_name)
        } else {
            p1 <- p1 + labs(title = var_name)
        }
        return(p1)
    }
}




#' Print an object of class waasb
#'
#' Print a `waasb` object in two ways. By default, the results are shown in
#' the R console. The results can also be exported to the directory.
#'
#'
#' @param x An object of class `waasb`.
#' @param export A logical argument. If `TRUE|T`, a *.txt file is exported
#'   to the working directory
#' @param blup A logical argument. If `TRUE|T`, the blups are shown.
#' @param file.name The name of the file if `export = TRUE`
#' @param digits The significant digits to be shown.
#' @param ... Options used by the tibble package to format the output. See
#'   [`tibble::print()`][tibble::formatting] for more details.
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @method print waasb
#' @export
#' @examples
#'\donttest{
#' library(metan)
#' model <- waasb(data_ge,
#'   resp = c(GY, HM),
#'   gen = GEN,
#'   env = ENV,
#'   rep = REP
#' )
#' print(model)
#' }
print.waasb <- function(x, export = FALSE, blup = FALSE, file.name = NULL, digits = 4, ...) {
    if (!class(x) == "waasb") {
        stop("The object must be of class 'waasb'")
    }
    if (export == TRUE) {
        file.name <- ifelse(is.null(file.name) == TRUE, "waasb print", file.name)
        sink(paste0(file.name, ".txt"))
    }
    opar <- options(pillar.sigfig = digits)
    on.exit(options(opar))
    for (i in 1:length(x)) {
        var <- x[[i]]
        cat("Variable", names(x)[i], "\n")
        cat("---------------------------------------------------------------------------\n")
        cat("Individual fixed-model analysis of variance\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$individual$individual)
        cat("---------------------------------------------------------------------------\n")
        cat("Fixed effects\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$fixed)
        cat("---------------------------------------------------------------------------\n")
        cat("Random effects\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$random)
        cat("---------------------------------------------------------------------------\n")
        cat("Likelihood ratio test\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$LRT)
        cat("---------------------------------------------------------------------------\n")
        cat("Variance components and genetic parameters\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$ESTIMATES)
        cat("---------------------------------------------------------------------------\n")
        cat(" Principal component analysis of the G x E interaction matrix\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$PCA)
        cat("---------------------------------------------------------------------------\n")
        if (blup == TRUE) {
            cat("BLUPs for genotypes\n")
            print(var$BLUPgen)
            cat("---------------------------------------------------------------------------\n")
            cat("BLUPs for genotypes-vs-environments\n")
            cat("---------------------------------------------------------------------------\n")
            print(var$BLUPgge)
            cat("---------------------------------------------------------------------------\n")
        }
        cat("Some information regarding the analysis\n")
        cat("---------------------------------------------------------------------------\n")
        print(var$Details)
        cat("\n\n\n")
    }
    if (export == TRUE) {
        sink()
    }
}










#' Predict method for waasb fits
#'
#' Obtains predictions from an object fitted with [waasb()].
#'
#'
#' @param object An object of class `waasb`
#' @param ... Currently not used
#' @return A tibble with the predicted values for each variable in the model
#' @author Tiago Olivoto \email{tiagoolivoto@@gmail.com}
#' @method predict waasb
#' @export
#' @examples
#'\donttest{
#' library(metan)
#' model <- waasb(data_ge,
#'                env = ENV,
#'                gen = GEN,
#'                rep = REP,
#'                resp = c(GY, HM))
#' predict(model)
#' }
#'
predict.waasb <- function(object, ...) {
    if (class(object) != "waasb") {
        stop("The objectin must be an objectin of the class 'waasb'")
    }
    factors <- object[[1]][["BLUPint"]] %>% select_non_numeric_cols()
    numeric <- sapply(object, function(x){
        x[["BLUPint"]][["Predicted"]]
    })
    return(cbind(factors, numeric) %>% as_tibble())
}

Try the metan package in your browser

Any scripts or data that you put into this service are public.

metan documentation built on Nov. 10, 2021, 9:11 a.m.