diff --git a/DESCRIPTION b/DESCRIPTION index 4bda4c1..d341408 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: statgenGxE Type: Package Title: Genotype by Environment (GxE) Analysis -Version: 1.0.5 -Date: 2022-08-11 +Version: 1.0.6 +Date: 2023-12-12 Authors@R: c(person(given = "Bart-Jan", family = "van Rossum", email = "bart-jan.vanrossum@wur.nl", @@ -58,7 +58,7 @@ Description: Analysis of multi environment data of plant breeding experiments License: GPL Encoding: UTF-8 LazyData: true -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 Roxygen: list(markdown = TRUE) Depends: R (>= 3.3) Imports: @@ -68,6 +68,7 @@ Imports: knitr, lme4, methods, + rlang, statgenSTA (>= 1.0.6), xtable Suggests: asreml(>= 3.0), diff --git a/NAMESPACE b/NAMESPACE index a02c49e..857c566 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -48,6 +48,7 @@ import(stats) importFrom(grDevices,topo.colors) importFrom(graphics,plot) importFrom(methods,getFunction) +importFrom(rlang,.data) importFrom(statgenSTA,createTD) importFrom(statgenSTA,report) importFrom(stats,predict) diff --git a/NEWS.md b/NEWS.md index e5aa40e..898c1d6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,14 @@ -# statgenGxE 1.0.4.1 +# statgenGxE 1.0.6 -* The predict function for gxeVarComp output is extended so all variables in the fitted model can now be used for making predictions. +* Functions no longer rely on soft-deprecated ggplot2 functions. +* A small bug is fixed that made plotting of `gxeVarComp` output impossible when using asreml for fitting the models. + +# statgenGxE 1.0.5 + +* The predict function for `gxeVarComp` output is extended so all variables in the fitted model can now be used for making predictions. * The plot functions for AMMI and GGE analysis now have an argument `rotatePC` allowing the specification of a trial that is aligned with the positive x-axis in the plot. -* The gxeVarComp function now has an argument `models` allowing a subset of the available models to be fitted. -* A small bug in gxeMegaEnv that sometimes caused NA for predicted values is fixed. +* The `gxeVarCov` function now has an argument `models` allowing a subset of the available models to be fitted. +* A small bug in `gxeMegaEnv` that sometimes caused NA for predicted values is fixed. * Some minor changes in order and capitalization of outputs. # statgenGxE 1.0.4 diff --git a/R/createAMMI.R b/R/createAMMI.R index c3c7ab2..d48f3a1 100644 --- a/R/createAMMI.R +++ b/R/createAMMI.R @@ -599,8 +599,8 @@ plotAMMI1 <- function(loadings, } lims <- c(min(xMin, yMin), max(xMax, yMax)) p <- ggplot2::ggplot(totDat, - ggplot2::aes_string(x = "x", y = "y", - color = ".group")) + + ggplot2::aes(x = .data[["x"]], y = .data[["y"]], + color = .data[[".group"]])) + ## Needed for a square plot output. ggplot2::coord_equal(xlim = lims, ylim = lims, clip = "off") + ggplot2::scale_color_manual(breaks = colGroups, values = colNamed, @@ -626,8 +626,8 @@ plotAMMI1 <- function(loadings, if (!is.null(textDat)) { ## Plot genotypes and environments as text. p <- p + ggplot2::geom_text(data = textDat, - ggplot2::aes_string(label = "rownames(textDat)", - size = ".size"), + ggplot2::aes(label = rownames(textDat), + size = .data[[".size"]]), show.legend = !is.null(colorEnvBy), vjust = 0) + ggplot2::scale_size(range = range(textDat[[".size"]]), guide = "none") } @@ -838,8 +838,9 @@ plotAMMI2 <- function(loadings, nrowGuide <- shapesGuide <- sizesGuide <- NULL } lims <- c(min(xMin, yMin), max(xMax, yMax)) - p <- ggplot2::ggplot(totDat, ggplot2::aes_string(x = primAxis, y = secAxis, - color = ".group")) + + p <- ggplot2::ggplot(totDat, ggplot2::aes(x = .data[[primAxis]], + y = .data[[secAxis]], + color = .data[[".group"]])) + ## Needed for a square plot output. ggplot2::coord_equal(xlim = lims, ylim = lims, clip = "off") + ggplot2::scale_color_manual(breaks = colGroups, values = colNamed, @@ -861,10 +862,10 @@ plotAMMI2 <- function(loadings, if (!is.null(textDat)) { ## Plot genotypes and environments as text. p <- p + ggplot2::geom_text(data = textDat, - ggplot2::aes_string(label = "rownames(textDat)", - size = ".size", - vjust = ".vjust", - hjust = ".hjust"), + ggplot2::aes(label = rownames(textDat), + size = .data[[".size"]], + vjust = .data[[".vjust"]], + hjust = .data[[".hjust"]]), show.legend = !is.null(colorEnvBy)) + ggplot2::scale_size(range = range(textDat[[".size"]]), guide = "none") } @@ -917,10 +918,11 @@ plotAMMI2 <- function(loadings, ## for ease of plotting. perpDat <- data.frame(xend = xNew, yend = yNew) ## Add perpendicular lines to plot as segments. - p <- p + ggplot2::geom_segment(ggplot2::aes_string(x = 0, y = 0, - xend = "xend", - yend = "yend"), - data = perpDat, col = "grey50", size = 0.6) + p <- p + ggplot2::geom_segment(ggplot2::aes(x = 0, y = 0, + xend = .data[["xend"]], + yend = .data[["yend"]]), + data = perpDat, col = "grey50", + linewidth = 0.6) } } if (plotEnv) { @@ -930,9 +932,9 @@ plotAMMI2 <- function(loadings, ## the plot otherwise. p <- p + ggplot2::geom_segment(data = envDat, - ggplot2::aes_string(x = 0, y = 0, - xend = primAxis, - yend = secAxis), + ggplot2::aes(x = 0, y = 0, + xend = .data[[primAxis]], + yend = .data[[secAxis]]), arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), show.legend = FALSE) } diff --git a/R/createFW.R b/R/createFW.R index 5b32f31..3bb49be 100644 --- a/R/createFW.R +++ b/R/createFW.R @@ -208,9 +208,9 @@ plot.FW <- function(x, scatterDat <- ggplot2::remove_missing(scatterDat, na.rm = TRUE) ## Create plot of mean x mse. No x axis because of position in grid. p1 <- ggplot2::ggplot(data = scatterDat, - ggplot2::aes_string(x = "GenMean", - y = "sqrt(MSdeviation)", - color = colorGenoBy)) + + ggplot2::aes(x = .data[["GenMean"]], + y = sqrt(.data[["MSdeviation"]]), + color = .data[[colorGenoBy]])) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = colGeno) + ggplot2::theme(axis.title.x = ggplot2::element_blank(), @@ -229,17 +229,18 @@ plot.FW <- function(x, p1 <- p1 + ggplot2::theme(legend.position = "none") ## Create plot of mean x sensitivity. p2 <- ggplot2::ggplot(data = scatterDat, - ggplot2::aes_string(x = "GenMean", y = "Sens", - color = colorGenoBy)) + + ggplot2::aes(x = .data[["GenMean"]], + y = .data[["Sens"]], + color = .data[[colorGenoBy]])) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = colGeno) + ggplot2::theme(legend.position = "none") + ggplot2::labs(x = "Mean", y = "Sensitivity") ## Create plot of mse x sensitivity. No y axis because of position in grid. p3 <- ggplot2::ggplot(data = scatterDat, - ggplot2::aes_string(x = "sqrt(MSdeviation)", - y = "Sens", - color = colorGenoBy)) + + ggplot2::aes(x = sqrt(.data[["MSdeviation"]]), + y = .data[["Sens"]], + color = .data[[colorGenoBy]])) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = colGeno) + ggplot2::theme(legend.position = "none", @@ -275,9 +276,6 @@ plot.FW <- function(x, lineDat <- ggplot2::remove_missing(lineDat, na.rm = TRUE) ## Set arguments for plot aesthetics. yVar <- ifelse(response == "observed", "genoMean", "fitted") - aesArgs <- list(x = "EnvMean", y = yVar, group = "genotype", - color = if (colorGenoBy == ".colorGenoBy") "genotype" else - enquote(colorGenoBy)) ## Order descending can be achieved by reversing the x-axis. if (order == "descending") { xTrans <- "reverse" @@ -285,12 +283,17 @@ plot.FW <- function(x, xTrans <- "identity" } plotLims <- range(c(lineDat[["EnvMean"]], lineDat[[yVar]])) + colorVar <- if (colorGenoBy == ".colorGenoBy") "genotype" else + enquote(colorGenoBy) ## Create plot. - p <- ggplot2::ggplot(data = lineDat, - do.call(ggplot2::aes_string, args = aesArgs)) + + p <- ggplot2::ggplot( + data = lineDat, + ggplot2::aes(x = .data[["EnvMean"]], y = .data[[yVar]], + group = .data[["genotype"]], + color = .data[[colorVar]])) + ggplot2::geom_point() + - ggplot2::geom_line(ggplot2::aes_string(y = "fitted"), - size = 0.5, alpha = 0.7) + + ggplot2::geom_line(ggplot2::aes(y = .data[["fitted"]]), + linewidth = 0.5, alpha = 0.7) + ggplot2::scale_x_continuous(trans = xTrans, sec.axis = ggplot2::dup_axis(name = "Environment", breaks = envEffs[["EnvMean"]], @@ -329,10 +332,12 @@ plot.FW <- function(x, trellisDat <- trellisDat[order(trellisDat[["genotype"]], trellisDat[["EnvMean"]]), ] p <- ggplot2::ggplot(data = trellisDat, - ggplot2::aes_string(x = "EnvMean", y = "genoMean")) + + ggplot2::aes(x = .data[["EnvMean"]], + y = .data[["genoMean"]])) + ggplot2::geom_point() + ggplot2::geom_line(data = trellisDat, - ggplot2::aes_string(x = "EnvMean", y = "fitted")) + + ggplot2::aes(x = .data[["EnvMean"]], + y = .data[["fitted"]])) + ggplot2::facet_wrap(facets = "genotype") + ggplot2::labs(x = "Environment", y = trait) + ggplot2::ggtitle(title) + @@ -356,8 +361,9 @@ plot.FW <- function(x, plotDat <- merge(plotDat, genoDat[, c("genotype", colorGenoBy)]) ## Create scatter plot of fitted values. p <- ggplot2::ggplot(data = plotDat, - ggplot2::aes_string(x = "trMin", y = "trMax", - color = colorGenoBy)) + + ggplot2::aes(x = .data[["trMin"]], + y = .data[["trMax"]], + color = .data[[colorGenoBy]])) + ggplot2::geom_point(na.rm = TRUE, show.legend = colorGenoBy != ".colorGenoBy") + ggplot2::scale_color_manual(values = colGeno) + @@ -394,7 +400,7 @@ plot.FW <- function(x, #' @export fitted.FW <- function(object, ...) { - return(object$fittedGeno) + return(object$fittedGeno) } #' Extract residuals. diff --git a/R/createStability.R b/R/createStability.R index cb06b5a..36b0e6c 100644 --- a/R/createStability.R +++ b/R/createStability.R @@ -133,9 +133,9 @@ plot.stability <- function(x, ## Create superiority plot. supDat <- merge(x$superiority, genoDat, by.x = "Genotype", by.y = "genotype") plots$p1 <- ggplot2::ggplot(data = supDat, - ggplot2::aes_string(x = "Mean", - y = "sqrt(Superiority)", - color = colorGenoBy)) + + ggplot2::aes(x = .data[["Mean"]], + y = sqrt(.data[["Superiority"]]), + color = .data[[colorGenoBy]])) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = colGeno) + ggplot2::labs(x = "Mean", y = "Square root of\n Cultivar superiority") @@ -144,9 +144,9 @@ plot.stability <- function(x, ## Create static plot. statDat <- merge(x$static, genoDat, by.x = "Genotype", by.y = "genotype") plots$p2 <- ggplot2::ggplot(data = statDat, - ggplot2::aes_string(x = "Mean", - y = "sqrt(Static)", - color = colorGenoBy)) + + ggplot2::aes(x = "Mean", + y = sqrt(.data[["Static"]]), + color = .data[[colorGenoBy]])) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = colGeno) + ggplot2::labs(x = "Mean", y = "Square root of\n Static stability") @@ -155,9 +155,9 @@ plot.stability <- function(x, ## Create Wricke plot. wrickeDat <- merge(x$wricke, genoDat, by.x = "Genotype", by.y = "genotype") plots$p3 <- ggplot2::ggplot(data = wrickeDat, - ggplot2::aes_string(x = "Mean", - y = "sqrt(Wricke)", - color = colorGenoBy)) + + ggplot2::aes(x = "Mean", + y = sqrt(.data[["Wricke"]]), + color = .data[[colorGenoBy]])) + ggplot2::geom_point() + ggplot2::scale_color_manual(values = colGeno) + ggplot2::labs(x = "Mean", y = "Square root of\n Wricke's ecovalence") diff --git a/R/createVarComp.R b/R/createVarComp.R index 4191787..4d6e35c 100644 --- a/R/createVarComp.R +++ b/R/createVarComp.R @@ -174,10 +174,12 @@ plot.varComp <- function(x, ## The idea is to annotate the y-axis just left of the x = 0. annoPosX <- -max(fullRandVC[[plotVar]]) / 5e5 p <- ggplot2::ggplot(fullRandVC, - ggplot2::aes_string(x = plotVar, y = "term")) + + ggplot2::aes(x = .data[[plotVar]], + y = .data[["term"]])) + ggplot2::geom_point(na.rm = TRUE, size = 2) + ## Add line from y-axis to points. - ggplot2::geom_segment(ggplot2::aes_string(xend = plotVar, yend = "term"), + ggplot2::geom_segment(ggplot2::aes(xend = .data[[plotVar]], + yend = .data[["term"]]), x = 0) + ggplot2::scale_x_continuous(expand = ggplot2::expansion(mult = c(0, 0.05))) + ## Set lower xlim to 0. This assures 0 is always displayed on the x-axis diff --git a/R/createVarCov.R b/R/createVarCov.R index 81e3799..fbebcbf 100644 --- a/R/createVarCov.R +++ b/R/createVarCov.R @@ -102,7 +102,9 @@ plot.varCov <- function(x, meltedCorMatLow <- meltedCorMat[as.numeric(meltedCorMat[["trial1"]]) > as.numeric(meltedCorMat[["trial2"]]), ] p <- ggplot2::ggplot(data = meltedCorMatLow, - ggplot2::aes_string("trial1", "trial2", fill = "cor")) + + ggplot2::aes(x = .data[["trial1"]], + y = .data[["trial2"]], + fill = .data[["cor"]])) + ggplot2::geom_tile(color = "white") + ggplot2::scale_y_discrete(position = "right") + ## Create a gradient scale. diff --git a/R/gxeVarComp.R b/R/gxeVarComp.R index a1e689d..da89e98 100644 --- a/R/gxeVarComp.R +++ b/R/gxeVarComp.R @@ -153,7 +153,7 @@ gxeVarComp <- function(TD, ## Random terms are genotype x fixedTerms. ## If there are no replicates or weights the final random term is the actual ## residual and therefore left out of the model. - if (hasReps || useWt) { + if (hasReps) { #} || useWt) { randTermIncl <- fixedTerms } else { randTermIncl <- fixedTerms[-length(fixedTerms)] @@ -208,6 +208,7 @@ gxeVarComp <- function(TD, randTermPos <- sapply(X = aovTermSets, FUN = setequal, randTermSet) ## Get MSS for current term. MSSRandTerm <- aovFullFixedMod[randTermPos, "Mean Sq"] + if (is.na(MSSRandTerm)) MSSRandTerm <- 0 ## For all other terms in the ANOVA table that have the current term ## as a subset the MSS cannot be higher. ## If it is the corresponding variance component is possibly zero. @@ -244,7 +245,7 @@ gxeVarComp <- function(TD, fullRandTxt <- paste("~", paste(c(fixedTerms, randTerms), collapse = "+")) fullRandMod <- tryCatchExt(asreml::asreml(fixed = formula(paste0("`", trait, "`~ 1")), random = formula(fullRandTxt), - family = asreml::asr_gaussian(dispersion = 1), + #family = asreml::asr_gaussian(dispersion = 1), data = TDTot, weights = "wt", maxiter = maxIter, trace = FALSE)) if (!is.null(fullRandMod$warning)) { @@ -306,16 +307,15 @@ gxeVarComp <- function(TD, } else if (engine == "asreml") { if (requireNamespace("asreml", quietly = TRUE)) { randTxt <- paste("~ ", paste(randTerms, collapse = "+")) - ## Put arguments for models in a list to make it easier to switch - ## between asreml3 and asreml4. Usually only one or two arguments differ. - ## Also some arguments are identical for all models - modArgs0 <- list(fixed = formula(fixedTxt), random = formula(randTxt), + ## Putting family argument in a list for some reason makes it impossible + ## to extract fitted values from asreml. + ## Since we only fit a single model here it can be called directly. + mr <- tryCatchExt( + asreml::asreml(fixed = formula(fixedTxt), + random = formula(randTxt), family = asreml::asr_gaussian(dispersion = 1), data = TDTot, weights = "wt", maxiter = maxIter, - trace = FALSE) - modArgs <- modArgs0 - ## Fit the actual model. - mr <- tryCatchExt(do.call(asreml::asreml, modArgs)) + trace = FALSE)) if (!is.null(mr$warning)) { ## Check if param 1% increase is significant. Remove warning if not. mr <- chkLastIter(mr) diff --git a/R/statgenGxE.R b/R/statgenGxE.R index b8d8828..8139140 100644 --- a/R/statgenGxE.R +++ b/R/statgenGxE.R @@ -6,6 +6,7 @@ #' @importFrom xtable xtable #' @importFrom utils hasName #' @importFrom statgenSTA createTD +#' @importFrom rlang .data NULL statgenDefaultOptions <- list( diff --git a/_pkgdown.yml b/_pkgdown.yml index fa7c59a..d2e22bc 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -30,9 +30,11 @@ navbar: - text: "statgenIBD" href: https://biometris.github.io/statgenIBD - text: "statgenMPP" - href: https://biometris.github.io/statgenIBD - - text: "LMMsolver" href: https://biometris.github.io/statgenMPP + - text: "statgenQTLxT" + href: https://biometris.github.io/statgenQTLxT + - text: "LMMsolver" + href: https://biometris.github.io/LMMsolver - icon: fa-github fa-lg text: "github" href: https://github.com/Biometris/statgenGxE diff --git a/vignettes/statgenGxE.Rmd b/vignettes/statgenGxE.Rmd index 6128671..0250215 100644 --- a/vignettes/statgenGxE.Rmd +++ b/vignettes/statgenGxE.Rmd @@ -55,9 +55,9 @@ The first step is loading the data into R. ```{r loadData} data(dropsPheno) ``` -dropsPheno contains data for the genotypic means (Best Linear Unbiased Estimators, BLUEs), with one value per genotype per experiment, for a selection of 10 experiments from the full data set and eight traits. These 10 experiments form a good representation of the full set of experiments covering the five scenarios described in @Millet2016. Throughout this vignette in all examples the trait grain.yield will be analyzed. For a more detailed description of the contents of the data see `help(dropsData)`. +dropsPheno contains data for the genotypic means (Best Linear Unbiased Estimators, BLUEs), with one value per genotype per experiment, for a selection of 10 experiments from the full data set and eight traits. These 10 experiments form a good representation of the full set of experiments covering the five scenarios described in @Millet2016. Throughout this vignette in all examples the trait grain.yield will be analyzed. For a more detailed description of the contents of the data see `help(dropsPheno)`. -The input for GxE analysis in the statgenGxE package is an object of class TD (Trial Data). For a detailed description on how to construct such an object see the vignette of the statgenSTA package (`vignette("Modeling field trials using statgenSTA")`). For GxE analysis it is enough to specify the genotype and trial option of the `createTD` function. +The input for GxE analysis in the statgenGxE package is an object of class TD (Trial Data). For a detailed description on how to construct such an object see the vignette of the statgenSTA package (`vignette("statgenSTA")`). For GxE analysis it is enough to specify the genotype and trial option of the `createTD` function. ```{r createTD} ## Create a TD object from dropsPheno. @@ -409,8 +409,8 @@ The models fitted are of the form $y_{ij} = \mu_j + \epsilon_{ij}$, where $y_{ij | Model | Description | var($g_{ij}$) | cov($g_{ij}$;$g_{ik}$) | Number of parameters| |:-----------|:-------------------|:------------|:------------|:-----------------| -| identity | identity | $\sigma_G^2$ | 0 | 1 | -| cs | compound symmetry | $\sigma_G^2 + \sigma_{GE}^2$ | $\sigma_{GE}^2$ | 2 | +| identity | identity | $\sigma_{GE}^2$ | 0 | 1 | +| cs | compound symmetry | $\sigma_G^2 + \sigma_{GE}^2$ | $\sigma_G^2$ | 2 | | diagonal | diagonal matrix (heteroscedastic) | $\sigma_{GE_j}^2$ | 0 | $J$ | | hcs | heterogeneous compound symmetry | $\sigma_G^2+\sigma_{GE_j}^2$ | $\sigma_G^2$ | $J+1$ | | outside | heterogeneity outside | $\sigma_{G_j}^2$ | $\theta$ | $J+1$ |