Merge branch 'develop'

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 = "[email protected]", 
@@ -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),

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)

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

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)
   }

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,22 +276,24 @@ 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"
     } else {
       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.

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")

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

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.

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)

R/statgenGxE.R

@@ -6,6 +6,7 @@
 #' @importFrom xtable xtable
 #' @importFrom utils hasName
 #' @importFrom statgenSTA createTD
+#' @importFrom rlang .data
 NULL
 
 statgenDefaultOptions <- list(

_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