createFitMod.R

#' @keywords internal
createFitMod <- function(models,
                         experimentName,
                         what,
                         useRepId,
                         useCheck,
                         spatial,
                         timePoints) {
  fitMod <- structure(models,
                      experimentName = experimentName,
                      timePoints = timePoints,
                      what = what,
                      useRepId = useRepId,
                      useCheck = useCheck,
                      spatial = spatial,
                      class = c("fitMod", "list"),
                      timestamp = Sys.time())
  return(fitMod)
}

#' Summary function for fitMod objects
#'
#' Function for creating a short summary of the contents of a TP object. The
#' summary consists of the name of the experiment, the number of time points,
#' the engine used to fit the models and, in case spatial models where fitted
#' using asreml, the selected spatial model.
#'
#' @param object An object of class fitMod.
#' @param ... Ignored.
#'
#' @return No return value, a summary is printed.
#'
#' @examples
#' \donttest{
#' ## Using the first example dataset (PhenovatorDat1):
#' ## Create an object of class TP.
#' phenoTP <- createTimePoints(dat = PhenovatorDat1,
#'                             experimentName = "Phenovator",
#'                             genotype = "Genotype",
#'                             timePoint = "timepoints",
#'                             repId = "Replicate",
#'                             plotId = "pos",
#'                             rowNum = "y", colNum = "x",
#'                             addCheck = TRUE,
#'                             checkGenotypes = c("check1", "check2",
#'                                                "check3", "check4"))
#'
#' ## Fit a SpATS model on few time points:
#' modPhenoSp <- fitModels(TP = phenoTP,
#'                         trait = "EffpsII",
#'                         timePoints = c(1, 6, 36))
#'
#' ## Create a summary.
#' summary(modPhenoSp)
#' }
#'
#' @family functions for spatial modeling
#'
#' @export
summary.fitMod <- function(object,
                           ...)  {
  experimentName <- attr(x = object, which = "experimentName")
  noTP <- length(object)
  engine <- class(object[[1]])
  if (engine == "asreml" && attr(x = object, which = "spatial")) {
    tpUsed <- min(max(noTP / 5, 10), noTP)
    bestSpat <- attr(x = object[[1]], which = "sumTab")[[1]][, "spatial"]
  }
  cat("Models in ", deparse(substitute(object)),
      " where fitted for experiment ", experimentName, ".\n\n", sep = "")
  cat("It contains", noTP, "time points.\n")
  cat("The models were fitted using ", engine, ".\n\n", sep = "")
  if (engine == "asreml" && attr(x = object, which = "spatial")) {
    cat("The selected spatial model is ", bestSpat, ".\n", sep = "")
    cat(tpUsed, "time points were used to select the best spatial model.\n")
  }
}

#' Plot function for class fitMod
#'
#' Plotting function for objects of class \code{fitMod}. Seven different types
#' of plots can be made for an object of class \code{fitMod}. A detailed
#' description and optional extra parameters for the different plots is given
#' in the sections below.
#'
#' @section rawPred plot:
#' Plots the raw data (colored dots) overlayed with the predicted values from
#' the fitted model (black dots). For each genotype a plot is made per
#' plot/plant over time. These plots are put together in a 5x5 grid. By using
#' the parameter \code{genotypes} a selection of genotypes can be plotted.
#' Extra parameter options:
#' \describe{
#' \item{genotypes}{A character vector indicating the genotypes to be plotted.}
#' \item{plotChecks}{Should the check genotypes be included in the plot?}
#' \item{plotLine}{Should the data be displayed as lines? Default is
#' \code{FALSE}.}
#' }
#'
#' @section corrPred plot:
#' Plots the spatially corrected data (colored dots) overlayed with the
#' predicted values from the fitted model (black dors). For each genotype a plot
#' is made per plot/plant over time. These plots are put together in a 5x5 grid.
#' By using the parameter \code{genotypes} a selection of genotypes can be
#' plotted. Extra parameter options:
#' \describe{
#' \item{genotypes}{A character vector indicating the genotypes to be plotted.}
#' \item{plotChecks}{Should the check genotypes be included in the plot?}
#' \item{plotLine}{Should the data be displayed as lines? Default is
#' \code{FALSE}.}
#' }
#'
#' @section herit plot:
#' Plots the heritability over time. This plot is only available when genotype
#' is fitted as random factor in the model. If \code{geno.decomp} is used when
#' fitting the model, heritabilities are plotted for each level of geno.decomp
#' in a single plot. Extra parameter options:
#' \describe{
#' \item{yLim}{A numerical vector of length two, used for setting the limits of
#' the y-axis of the plot. If values outside of the plotting range are given,
#' then these are ignored.}
#' }
#'
#' @section effDim plot:
#' Plots the effective dimension over time for models fitted using SpATS.
#' Extra parameter options:
#' \describe{
#' \item{whichED}{A character vector indicating which effective dimensions
#' should be plotted. This should be a subset of "colId", "rowId", "fCol",
#' "fRow", "fColRow", "colfRow", "fColfRow" and "surface". When
#' \code{useRepId = TRUE}, the effective dimensions of "colId" and "rowId"
#' become "RepId:colId" and "RepId:rowId". Default all effective dimensions
#' are plotted.}
#' \item{EDType}{A character string specifying if the effective dimension
#' ("dimension") or the ratio of effective dimensions ("ratio") should be
#' plotted. Default the dimensions are plotted.}
#' \item{yLim}{A numerical vector of length two, used for setting the limits of
#' the y-axis of the plot. If values outside of the plotting range are given,
#' then these are ignored.}
#' }
#'
#' @section variance plot:
#' Plots the residual, column and row variances over time for the fitted models.
#' Extra parameter options:
#' \describe{
#' \item{yLim}{A numerical vector of length two, used for setting the limits of
#' the y-axis of the plot. If values outside of the plotting range are given,
#' then these are ignored.}
#' }
#'
#' @section timeLapse plot:
#' Creates a time lapse of the spatial trends of models fitted using SpATS over
#' time.
#'
#' @section spatial plot:
#' Creates five plots per time point, spatial plots of the raw data,
#' fitted values, residuals and either BLUEs or BLUPs, and a histogram of the
#' BLUEs or BLUPs. When SpATS was used for modeling an extra plot with the
#' fitted spatial trend is included Extra parameter options:
#' \describe{
#' \item{spaTrend}{A character string indicating how the spatial trend should
#' be displayed. Either "raw" for raw values, or "percentage" for displaying
#' as a percentage of the original phenotypic values.}
#' }
#'
#' @inheritParams plot.TP
#'
#' @param x An object of class fitMod.
#' @param outFile A character string indicating the .pdf file or .gif file
#' (for \code{plotType} = "timeLapse") to which the plots should be written.
#'
#' @return Depending on the plot type either a ggplot object or a list of
#' ggplot objects is invisibly returned.
#'
#' @examples
#' \donttest{
#' ## Using the first example dataset (PhenovatorDat1):
#' ## Create an object of class TP.
#' phenoTP <- createTimePoints(dat = PhenovatorDat1,
#'                             experimentName = "Phenovator",
#'                             genotype = "Genotype",
#'                             timePoint = "timepoints",
#'                             repId = "Replicate",
#'                             plotId = "pos",
#'                             rowNum = "y", colNum = "x",
#'                             addCheck = TRUE,
#'                             checkGenotypes = c("check1", "check2",
#'                                                "check3", "check4"))
#'
#' ## Fit a SpATS model on three points:
#' modPhenoSp <- fitModels(TP = phenoTP,
#'                         trait = "EffpsII",
#'                         timePoints = c(1, 6, 36))
#'
#' ## Plot the spatial trends for one time point:
#' plot(modPhenoSp,
#'      timePoints = 36,
#'      plotType = "spatial",
#'      spaTrend = "percentage")
#' }
#'
#' \dontrun{
#' ## Create a time lapse of all available time points:
#' plot(modPhenoSp,
#'      plotType = "timeLapse",
#'      outFile = "TimeLapse_modPhenoSp.gif")
#' }
#'
#' \donttest{
#' ## Plot the corrected values for a subset of four genotypes:
#' plot(modPhenoSp,
#'      plotType = "corrPred",
#'      genotypes = c("check1", "check2", "G007", "G058") )
#'
#' ## Plot the effective dimensions of all available time points in the model
#' ## for a subset of effective dimensions:
#' plot(modPhenoSp,
#'      plotType = "effDim",
#'      whichED = c("colId", "rowId", "fColRow","colfRow"),
#'      EDType = "ratio")
#' }
#'
#' @family functions for spatial modeling
#'
#' @export
plot.fitMod <- function(x,
                        ...,
                        plotType = c("rawPred", "corrPred", "herit", "effDim",
                                     "variance", "timeLapse", "spatial"),
                        timePoints = names(x),
                        title = NULL,
                        output = TRUE,
                        outFile = NULL,
                        outFileOpts = NULL) {
  ## Checks.
  timePoints <- chkTimePoints(x, timePoints)
  plotType <- match.arg(plotType)
  experimentName <- attr(x = x, which = "experimentName")
  dotArgs <- list(...)
  if (!is.null(title) && (!is.character(title) || length(title) > 1)) {
    stop("title should be NULL or a character string.\n")
  }
  ## Restrict x to selected time points.
  fitMods <- x[timePoints]
  ## Get engine from fitted models.
  engine <- class(fitMods[[1]])
  if (engine == "asreml" && plotType == "effDim") {
    stop("Effective dimensions can only be plotted for models fitted ",
         "with SpATS.\n")
  }
  if (engine == "asreml" && plotType == "spatial" &&
      !attr(x = fitMods, which = "spatial")) {
    stop("spatial plots can only be made when setting spatial = TRUE when ",
         "fitting the asreml models.\n")
  }
  ## Get geno.decomp from fitted models.
  if (engine == "SpATS") {
    geno.decomp <- fitMods[[1]]$model$geno$geno.decomp
  } else if (engine == "asreml") {
    if ("geno.decomp" %in% all.vars(fitMods[[1]]$formulae$random)) {
      geno.decomp <- "geno.decomp"
    } else {
      geno.decomp <- NULL
    }
  }
  ## Get trait from fitted models.
  if (engine == "SpATS") {
    trait <- fitMods[[1]]$model$response
  } else if (engine == "asreml") {
    ## Trait is always the only lhs variable in the fixed part.
    trait <- all.vars(update(fitMods[[1]]$formulae$fixed, .~0))
  }
  ## Get check from fitted models.
  if (engine == "SpATS") {
    useCheck <- grepl(pattern = "check", x = deparse(fitMods[[1]]$model$fixed))
  } else if (engine == "asreml") {
    useCheck <- "check" %in% all.vars(update(fitMods[[1]]$formulae$fixed, 0~.))
  }
  ## Get what from fitted models.
  what <- attr(x = fitMods, which = "what")
  if (!is.null(outFile) && plotType != "timeLapse") {
    chkFile(outFile, fileType = "pdf")
    output <- TRUE
    outFileOpts <- c(list(file = outFile), outFileOpts)
    on.exit(dev.off(), add = TRUE)
    do.call(pdf, args = outFileOpts)
  }
  if (plotType == "rawPred") {
    genotypes <- dotArgs$genotypes
    plotLine <- isTRUE(dotArgs$plotLine)
    if (!is.null(genotypes) && !is.character(genotypes)) {
      stop("genotypes should be NULL or a character vector.\n")
    }
    if (is.null(title)) title <-
      paste(experimentName, "- genotypic prediction + raw data")
    if (isTRUE(dotArgs$plotChecks) && useCheck) {
      totPred <- getGenoPred(fitMods, predictChecks = TRUE)
      ## Get check predictions.
      genoPred <- totPred$genoPred
      checkPred <- totPred$checkPred
      ## Rename check column to genotype so rbinding is possible.
      colnames(checkPred)[colnames(checkPred) == "check"] <- "genotype"
      preds <- rbind(genoPred, checkPred)
    } else {
      ## Get genotypic predictions.
      preds <- getGenoPred(fitMods)$genoPred
    }
    ## Construct full raw data from models.
    if (engine == "SpATS") {
      raw <- Reduce(f = rbind, x = lapply(X = fitMods, FUN = `[[`, "data"))
    } else if (engine == "asreml") {
      raw <- Reduce(f = rbind, x = lapply(X = fitMods, FUN = function(fitMod) {
        fitMod$call$data
      }))
    }
    ## Remove observations from raw where genotype is missing.
    ## These where included when fitting spatial asreml models to create
    ## a full grid.
    raw <- raw[!is.na(raw[["genotype"]]), ]
    if (!is.null(geno.decomp) && engine == "SpATS" && !useCheck) {
      ## Genotype was converted to an interaction term of genotype and
      ## geno.decomp in the proces of fitting the model. That needs to be
      ## undone to get the genotype back in the output again.
      genoStart <- nchar(as.character(raw[["geno.decomp"]])) + 2
      raw[["genotype"]] <- as.factor(substring(raw[["genotype"]],
                                               first = genoStart))
    }
    ## Remove check genotypes from raw data.
    if (!isTRUE(dotArgs$plotChecks) && useCheck) {
      raw <- droplevels(raw[!is.na(raw[["genoCheck"]]), ])
    }
    ## Restrict genotypes.
    if (!is.null(genotypes)) {
      if (!all(genotypes %in% preds[["genotype"]])) {
        stop("All genotypes should be in ", deparse(substitute(x)), ".\n")
      }
      preds <- preds[preds[["genotype"]] %in% genotypes, ]
      preds <- droplevels(preds)
      raw <- raw[raw[["genotype"]] %in% genotypes, ]
      raw <- droplevels(raw)
    }
    ## Restrict raw to neccessary columns so adding missing combinations works
    ## properly. With extra columns unneeded combinations might be added.
    raw <- raw[colnames(raw) %in% c("timePoint", "genotype", "plotId", trait,
                                    if (!is.null(geno.decomp)) "geno.decomp",
                                    if (useCheck) c("check", "genoCheck"))]
    ## Add combinations missing in data to raw.
    raw <- addMissVals(dat = raw, trait = trait)
    p <- xyFacetPlot(baseDat = raw, overlayDat = preds, yVal = trait,
                     yValOverlay = "predicted.values",
                     facetVal = c("genotype",
                                  if (!is.null(geno.decomp)) "geno.decomp"),
                     title = title, yLab = trait, output = output,
                     plotLine = plotLine)
  } else if (plotType == "corrPred") {
    genotypes <- dotArgs$genotypes
    plotLine <- isTRUE(dotArgs$plotLine)
    if (!is.null(genotypes) && !is.character(genotypes)) {
      stop("genotypes should be NULL or a character vector.\n")
    }
    if (is.null(title))
      title <- paste(experimentName, "- genotypic prediction + corrected data")
    if (isTRUE(dotArgs$plotChecks) && useCheck) {
      totPred <- getGenoPred(fitMods, predictChecks = TRUE)
      ## Get check predictions.
      genoPred <- totPred$genoPred
      checkPred <- totPred$checkPred
      ## Rename check column to genotype so rbinding is possible.
      colnames(checkPred)[colnames(checkPred) == "check"] <- "genotype"
      preds <- rbind(genoPred, checkPred)
    } else {
      ## Get genotypic predictions.
      preds <- getGenoPred(fitMods)$genoPred
    }
    ## Get spatial corrected values.
    corrected <- suppressWarnings(getCorrected(fitMods))
    ## Remove check genotypes from corrected data.
    if (!isTRUE(dotArgs$plotChecks) && useCheck) {
      corrected <- droplevels(corrected[corrected[["check"]] == "noCheck", ])
    }
    ## Restrict genotypes.,
    if (!is.null(genotypes)) {
      if (!all(genotypes %in% preds[["genotype"]])) {
        stop("All genotypes should be in ", deparse(substitute(x)), ".\n")
      }
      preds <- preds[preds[["genotype"]] %in% genotypes, ]
      preds <- droplevels(preds)
      corrected <- corrected[corrected[["genotype"]] %in% genotypes, ]
      corrected <- droplevels(corrected)
    }
    newTrait <- paste0(trait, "_corr")
    corrected <- corrected[c("timeNumber", "timePoint", "genotype",
                             newTrait, "plotId",
                             if (!is.null(geno.decomp)) "geno.decomp")]
    ## Add combinations missing in data to corrected.
    corrected <- addMissVals(dat = corrected, trait = newTrait)
    p <- xyFacetPlot(baseDat = corrected, overlayDat = preds, yVal = newTrait,
                     yValOverlay = "predicted.values",
                     facetVal = c("genotype",
                                  if (!is.null(geno.decomp)) "geno.decomp"),
                     title = title,
                     yLab = trait, output = output, plotLine = plotLine)
  } else if (plotType == "herit") {
    if (is.null(title)) title <- paste(experimentName, "- Heritability")
    ## Get heritabilities.
    herit <- getHerit(fitMods)
    ## Convert to long format needed by ggplot.
    herit <- reshape2::melt(herit, measure.vars = setdiff(colnames(herit),
                                                          c("timeNumber",
                                                            "timePoint")),
                            variable.name = "herit", value.name = "h2")
    ## Manually modify limit of y-axis.
    yLim <- c(min(dotArgs$yLim[1], herit[["h2"]]),
              max(dotArgs$yLim[2], herit[["h2"]]))
    p <- ggplot2::ggplot(herit,
                         ggplot2::aes_string(x = "timePoint", y = "h2",
                                             group = "herit",
                                             color = "herit")) +
      ggplot2::geom_point(size = 3, na.rm = TRUE) +
      plotTheme() +
      ggplot2::theme(axis.title.y = ggplot2::element_text(angle = 0,
                                                          vjust = 0.5)) +
      ggplot2::ylim(yLim) +
      ggplot2::labs(title = title)
    ## Compute the number of breaks for the time scale.
    ## If there are less than 4 time points use positions of the time points.
    ## Otherwise use 3.
    nTp <- length(unique(herit[["timePoint"]]))
    if (nTp < 5) {
      p <- p + ggplot2::scale_x_datetime(breaks = unique(herit[["timePoint"]]),
                                         labels = scales::date_format("%B %d"))
    } else {
      ## Format the time scale to Month + day.
      p <- p + ggplot2::scale_x_datetime(breaks = prettier(n = 3),
                                         labels = scales::date_format("%B %d"))
    }
    if (nTp > 1) {
      p <- p + ggplot2::geom_line(size = 0.5, na.rm = TRUE)
    }
    if (output) {
      plot(p)
    }
  } else if (plotType == "effDim") {
    useRepId <- attr(x = fitMods, which = "useRepId")
    colVarId <- ifelse(useRepId, "repId:colId", "colId")
    rowVarId <- ifelse(useRepId, "repId:rowId", "rowId")
    whichEDopts <- c(colVarId, rowVarId, "fCol", "fRow", "fColRow", "colfRow",
                     "fColfRow", "surface")
    if (is.null(dotArgs$which)) {
      whichED <- whichEDopts
    } else {
      whichED <- match.arg(dotArgs$whichED, choices = whichEDopts,
                           several.ok = TRUE)
    }
    EDType <- match.arg(dotArgs$EDType, choices = c("dimension", "ratio"))
    if (is.null(title)) title <- paste(experimentName, "- Effective dimension")
    ## Get effective dimensions.
    effDim <- getEffDims(fitMods, EDType = EDType)
    ## Convert to long format needed by ggplot.
    effDim <- reshape2::melt(effDim, measure.vars = whichED,
                             variable.name = "effDim", value.name = "ED")
    ## Manually modify limit of y-axis.
    yLim <- c(min(dotArgs$yLim[1], effDim[["ED"]]),
              max(dotArgs$yLim[2], effDim[["ED"]]))
    p <- ggplot2::ggplot(effDim,
                         ggplot2::aes_string(x = "timePoint", y = "ED",
                                             group = "effDim",
                                             color = "effDim")) +
      ggplot2::geom_point(size = 3, na.rm = TRUE) +
      plotTheme() +
      ggplot2::theme(axis.title.y = ggplot2::element_text(angle = 0,
                                                          vjust = 0.5)) +
      ggplot2::ylim(yLim) +
      ggplot2::labs(title = title, color = "Effective dimension")
    ## Compute the number of breaks for the time scale.
    ## If there are less than 4 time points use positions of the time points.
    ## Otherwise use 3.
    nTp <- length(unique(effDim[["timePoint"]]))
    if (nTp < 5) {
      p <- p + ggplot2::scale_x_datetime(breaks = unique(effDim[["timePoint"]]),
                                         labels = scales::date_format("%B %d"))
    } else {
      ## Format the time scale to Month + day.
      p <- p + ggplot2::scale_x_datetime(breaks = prettier(n = 3),
                                         labels = scales::date_format("%B %d"))
    }
    if (nTp > 1) {
      p <- p + ggplot2::geom_line(size = 0.5, na.rm = TRUE)
    }
    if (output) {
      plot(p)
    }
  } else if (plotType == "variance") {
    if (is.null(title)) title <- paste(experimentName, "- Variances")
    ## Get variances.
    variance <- getVar(fitMods)
    ## Get variance columns from variance, i.e. all columns starting with var.
    varCols <- colnames(variance)[grepl(pattern = "^var",
                                        x = colnames(variance))]
    ## Construct labels for variances.
    varLabs <- c(if ("varGen" %in% varCols) "Genotypic" else
      substring(varCols[1:(length(varCols) - 3)], first = 17),
      "Residual", "Columns", "Rows")
    ## Convert to long format needed by ggplot.
    variance <- reshape2::melt(variance, measure.vars = varCols,
                               variable.name = "var")
    ## Manually modify limit of y-axis.
    yLim <- c(min(dotArgs$yLim[1], variance[["value"]]),
              max(dotArgs$yLim[2], variance[["value"]]))
    p <- ggplot2::ggplot(variance,
                         ggplot2::aes_string(x = "timePoint", y = "value",
                                             group = "var", color = "var")) +
      ggplot2::geom_point(size = 3, na.rm = TRUE) +
      ggplot2::scale_color_discrete(labels = varLabs) +
      plotTheme() +
      ggplot2::theme(axis.title.y = ggplot2::element_text(angle = 0,
                                                          vjust = 0.5)) +
      ggplot2::ylim(yLim) +
      ggplot2::labs(title = title, color = "variance",
                    y = expression(sigma ^ 2))
    ## Compute the number of breaks for the time scale.
    ## If there are less than 4 time points use positions of the time points.
    ## Otherwise use 3.
    nTp <- length(unique(variance[["timePoint"]]))
    if (nTp < 5) {
      p <- p + ggplot2::scale_x_datetime(breaks = unique(variance[["timePoint"]]),
                                         labels = scales::date_format("%B %d"))
    } else {
      ## Format the time scale to Month + day.
      p <- p + ggplot2::scale_x_datetime(breaks = prettier(n = 3),
                                         labels = scales::date_format("%B %d"))
    }
    if (nTp > 1) {
      p <- p + ggplot2::geom_line(size = 0.5, na.rm = TRUE)
    }
    if (output) {
      plot(p)
    }
  } else if (plotType == "spatial") {
    p <- lapply(X = fitMods, FUN = spatPlot, trait = trait, what = what,
                geno.decomp = geno.decomp, useCheck = useCheck, engine = engine,
                output = output, ... = ..., experimentName = experimentName,
                title = title)
  } else if (plotType == "timeLapse") {
    chkFile(outFile, fileType = "gif")
    timeLapsePlot(fitMods, outFile = outFile, ...)
  }
  if (!plotType == "timeLapse") {
    invisible(p)
  }
}

#' Helper function for creating spatial plots.
#'
#' @noRd
#' @keywords internal
spatPlot <- function(fitMod,
                     trait,
                     what,
                     geno.decomp,
                     useCheck,
                     engine,
                     output = TRUE,
                     ...) {
  dotArgs <- list(...)
  ## Get plot type for spatial trend from args.
  spaTrend <- match.arg(dotArgs$spaTrend, choices = c("raw", "percentage"))
  ## Extract data from model.
  if (engine == "SpATS") {
    modDat <- fitMod$data
  } else if (engine == "asreml") {
    modDat <- fitMod$call$data
  }
  ## Check spatial information in modDat.
  if (!chkRowCol(modDat)) {
    return(NULL)
  }
  ## Extract time point from model data.
  timePoint <- modDat[["timePoint"]][1]
  ## Extract raw data.
  raw <- modDat[c("genotype", trait, "rowNum", "colNum", geno.decomp)]
  ## Extract fitted values from model.
  fitted <- fitted(fitMod)
  ## Extract predictions (BLUEs or BLUPs) from model.
  if (useCheck) {
    totPred <- predictGeno(fitMod, predictChecks = TRUE)
    ## Get check predictions.
    genoPred <- totPred$predGeno
    checkPred <- totPred$predCheck
    ## Rename check column to genotype so rbinding is possible.
    colnames(checkPred)[colnames(checkPred) == "check"] <- "genotype"
    pred <- rbind(genoPred, checkPred)
  } else {
    ## Get genotypic predictions.
    pred <- predictGeno(fitMod)$predGeno
  }
  pred <- pred[c("genotype", "predicted.values", geno.decomp)]
  if (!is.null(geno.decomp) && engine == "SpATS") {
    ## Genotype was converted to an interaction term of genotype and
    ## geno.decomp in the proces of fitting the model. That needs to be
    ## undone to get the genotype back in the output again.
    genoStart <- nchar(as.character(raw[["geno.decomp"]])) + 2
    raw[["genotype"]] <- as.factor(substring(raw[["genotype"]],
                                             first = genoStart))
  }
  ## Create plot data by merging extracted data together and renaming some
  ## columns.
  plotDat <- cbind(raw, fitted)
  plotDat <- merge(plotDat, pred, by = c("genotype", geno.decomp))
  plotDat[["predicted.values"]][is.na(plotDat[["fitted"]])] <- NA
  plotDat[["resid"]] <- plotDat[[trait]] - plotDat[["fitted"]]
  ## Get limits for row and columns.
  yMin <- min(plotDat[["rowNum"]])
  yMax <- max(plotDat[["rowNum"]])
  xMin <- min(plotDat[["colNum"]])
  xMax <- max(plotDat[["colNum"]])
  ## Execute this part first since it needs plotData without missings
  ## removed.
  ## Code mimickes code from SpATS package but is adapted to create a
  ## data.frame useable by ggplot.
  if (engine == "SpATS") {
    plotDat <- plotDat[order(plotDat[["colNum"]], plotDat[["rowNum"]]), ]
    nCol <- xMax - xMin + 1
    nRow <- yMax - yMin + 1
    p1 <- 100 %/% nCol + 1
    p2 <- 100 %/% nRow + 1
    ## Get spatial trend from SpATS object.
    spatTr <- SpATS::obtain.spatialtrend(fitMod, grid = c(nCol * p1, nRow * p2))
    ## spatial trend contains values for all data points, so NA in original
    ## data need to be removed. The kronecker multiplication is needed to
    ## convert the normal row col pattern to the smaller grid extending the
    ## missing values.
    ## First a matrix M is created containing information for all
    ## columns/rows in the field even if they are completely empty.
    M <- matrix(nrow = nRow, ncol = nCol,
                dimnames = list(yMin:yMax, xMin:xMax))
    for (i in seq_len(nrow(plotDat))) {
      M[as.character(plotDat[i, "rowNum"]),
        as.character(plotDat[i, "colNum"])] <-
        ifelse(is.na(plotDat[i, trait]), NA, 1)
    }
    spatTrDat <- kronecker(M, matrix(data = 1, ncol = p1, nrow = p2)) *
      spatTr$fit
    ## Melt to get the data in ggplot shape. Rows and columns in the
    ## spatial trend coming from SpATS are swapped so therefore use t()
    plotDatSpat <- reshape2::melt(t(spatTrDat),
                                  varnames = c("colNum", "rowNum"))
    ## Add true values for columns and rows for plotting.
    plotDatSpat[["colNum"]] <- spatTr$col.p
    plotDatSpat[["rowNum"]] <- rep(x = spatTr$row.p, each = p1 * nCol)
    ## Remove missings from data.
    plotDatSpat <- ggplot2::remove_missing(plotDatSpat, na.rm = TRUE)
  }
  ## Now missing values can be removed from plotDat.
  plotDat <- ggplot2::remove_missing(plotDat, na.rm = TRUE)
  ## Code taken from plot.SpATS and simplified.
  ## Set colors and legends.
  colors <- topo.colors(100)
  legends <- c("Raw data", "Fitted data", "Residuals",
               "Fitted Spatial Trend",
               ifelse(what == "fixed", "Genotypic BLUEs",
                      "Genotypic BLUPs"), "Histogram")
  ## Compute range of values in response + fitted data so same scale
  ## can be used over plots.
  zlim <- range(c(plotDat[[trait]], plotDat[["fitted"]]), na.rm = TRUE)
  ## Create empty list for storing plots
  plots <- vector(mode = "list")
  ## Create main plot title.
  plotTitle <- ifelse(!is.null(dotArgs$title), dotArgs$title,
                      paste(dotArgs$experimentName, "-", trait, "-", timePoint))
  ## Create separate plots.
  plots$p1 <- fieldPlot(plotDat = plotDat, fillVar = trait,
                        title = legends[1], colors = colors, zlim = zlim)
  plots$p2 <- fieldPlot(plotDat = plotDat, fillVar = "fitted",
                        title = legends[2], colors = colors, zlim = zlim)
  plots$p3 <- fieldPlot(plotDat = plotDat, fillVar = "resid",
                        title = legends[3], colors = colors)
  ## Spatial plot only for SpATS.
  if (engine == "SpATS") {
    ## Get tickmarks from first plot to be used as ticks.
    ## Spatial plot tends to use different tickmarks by default.
    xTicks <- ggplot2::ggplot_build(plots[[1]])$layout$panel_params[[1]]$x$breaks
    if (spaTrend == "raw") {
      plots$p4 <- fieldPlot(plotDat = plotDatSpat, fillVar = "value",
                            title = legends[4], colors = colors,
                            xTicks = xTicks)
    } else {
      phenoMean <- mean(modDat[[trait]], na.rm = TRUE)
      plotDatSpat[["value"]] <- plotDatSpat[["value"]] / phenoMean
      zlim <- c(-1, 1) * max(c(abs(plotDatSpat[["value"]]), 0.1))
      plots$p4 <- fieldPlotPcts(plotDat = plotDatSpat, fillVar = "value",
                                title = legends[4], zlim = zlim,
                                colors = colorRampPalette(c("red", "yellow",
                                                            "blue"),
                                                          space = "Lab")(100),
                                xTicks = xTicks)
    }
  }
  plots$p5 <- fieldPlot(plotDat = plotDat, fillVar = "predicted.values",
                        title = legends[5], colors = colors)
  plots$p6 <- ggplot2::ggplot(data = plotDat) +
    ggplot2::geom_histogram(ggplot2::aes_string(x = "predicted.values"),
                            fill = "white", col = "black", bins = 10,
                            boundary = 0) +
    ## Remove empty space between ticks and actual plot.
    ggplot2::scale_x_continuous(expand = c(0, 0)) +
    ggplot2::scale_y_continuous(expand = c(0, 0)) +
    ## No background. Center and resize title. Resize axis labels.
    ggplot2::theme(panel.background = ggplot2::element_blank(),
                   plot.title = ggplot2::element_text(hjust = 0.5, size = 10),
                   axis.title = ggplot2::element_text(size = 9)) +
    ggplot2::labs(y = "Frequency", x = legends[5], title = legends[6])
  if (output) {
    ## do.call is needed since grid.arrange doesn't accept lists as input.
    do.call(gridExtra::grid.arrange,
            args = c(Filter(f = Negate(f = is.null), x = plots),
                     list(ncol = 3, top = plotTitle)))
  }
  return(plots)
}

#' Helper function for creating field plots.
#'
#' @noRd
#' @keywords internal
fieldPlot <- function(plotDat,
                      fillVar,
                      title,
                      colors,
                      zlim = range(plotDat[fillVar]),
                      xTicks = ggplot2::waiver(),
                      ...) {
  p <- ggplot2::ggplot(data = plotDat,
                       ggplot2::aes_string(x = "colNum", y = "rowNum",
                                           fill = fillVar)) +
    ggplot2::geom_tile() +
    ## Remove empty space between ticks and actual plot.
    ggplot2::scale_x_continuous(expand = c(0, 0), breaks = xTicks) +
    ggplot2::scale_y_continuous(expand = c(0, 0)) +
    ## Adjust plot colors.
    ggplot2::scale_fill_gradientn(limits = zlim, colors = colors) +
    ## No background. Center and resize title. Resize axis labels.
    ## Remove legend title and resize legend entries.
    ggplot2::theme(panel.background = ggplot2::element_blank(),
                   plot.title = ggplot2::element_text(hjust = 0.5, size = 10),
                   axis.title = ggplot2::element_text(size = 9),
                   legend.title = ggplot2::element_blank(),
                   legend.text = ggplot2::element_text(size = 8,
                                                       margin = ggplot2::margin(l = 5))) +
    ggplot2::ggtitle(title)
  return(p)
}


#' Helper function for creating field plots with percentages.
#'
#' @noRd
#' @keywords internal
fieldPlotPcts <- function(plotDat,
                          fillVar,
                          title,
                          colors,
                          zlim = range(plotDat[fillVar]),
                          xTicks = ggplot2::waiver(),
                          scaleLim = Inf,
                          ...) {
  p <- ggplot2::ggplot(
    data = plotDat,
    ggplot2::aes_string(x = "colNum", y = "rowNum",
                        fill = fillVar,
                        color = if (is.infinite(scaleLim)) NULL else "''")) +
    ggplot2::geom_tile(na.rm = TRUE) +
    ## Remove empty space between ticks and actual plot.
    ggplot2::scale_x_continuous(expand = c(0, 0), breaks = xTicks) +
    ggplot2::scale_y_continuous(expand = c(0, 0)) +
    ## Adjust plot colors.
    ggplot2::scale_fill_gradientn(limits = zlim, colors = colors, name = NULL,
                                  labels = scales::percent,
                                  breaks = seq(zlim[1], zlim[2], length.out = 5)) +
    ggplot2::scale_color_manual(values = NA) +
    ggplot2::guides(fill = ggplot2::guide_colorbar(order = 1),
                    color = ggplot2::guide_legend("Larger than scale limit",
                                                  override.aes = list(fill = "grey50",
                                                                      color = "grey50"))) +
    ## No background. Center and resize title. Resize axis labels.
    ## Remove legend title and resize legend entries.
    ggplot2::theme(panel.background = ggplot2::element_blank(),
                   plot.title = ggplot2::element_text(hjust = 0.5, size = 10),
                   axis.title = ggplot2::element_text(size = 9)) +
    ggplot2::ggtitle(title)
  return(p)
}

#' Helper function for creating time lapse plots.
#'
#' @importFrom grDevices colorRampPalette
#' @noRd
#' @keywords internal
timeLapsePlot <- function(fitMods,
                          outFile = "spatialTrends.gif",
                          ...) {
  dotArgs <- list(...)
  if (!is.null(dotArgs$scaleLim)) {
    scaleLim <- dotArgs$scaleLim / 100
  } else {
    scaleLim <- Inf
  }
  animation::saveGIF({
    ## Get trait from fitted models.
    trait <- fitMods[[1]]$model$response
    ## First get plot data for all fields so a single zLim can be extracted
    ## for all plots.
    plotSpatDats <- lapply(X = fitMods, FUN = function(fitMod) {
      ## Get data from fitted model
      plotDat <- fitMod$data
      ## Get min and max values for rows and columns.
      yMin <- min(plotDat$rowNum)
      yMax <- max(plotDat$rowNum)
      xMin <- min(plotDat$colNum)
      xMax <- max(plotDat$colNum)
      ## Execute this part first since it needs plotData without missings
      ## removed.
      ## Code mimickes code from SpATS package but is adapted to create a
      ## data.frame useable by ggplot.
      plotDat <- plotDat[order(plotDat$colNum, plotDat$rowNum), ]
      nCol <- xMax - xMin + 1
      nRow <- yMax - yMin + 1
      p1 <- 100 %/% nCol + 1
      p2 <- 100 %/% nRow + 1
      ## Get spatial trend from SpATS object.
      spatTr <- SpATS::obtain.spatialtrend(fitMod,
                                           grid = c(nCol * p1, nRow * p2))
      ## spatial trend contains values for all data points, so NA in original
      ## data need to be removed. The kronecker multiplication is needed to
      ## convert the normal row col pattern to the smaller grid extending the
      ## missing values.
      ## First a matrix M is created containing information for all
      ## columns/rows in the field even if they are completely empty.
      M <- matrix(nrow = nRow, ncol = nCol,
                  dimnames = list(yMin:yMax, xMin:xMax))
      for (i in seq_len(nrow(plotDat))) {
        M[as.character(plotDat[i, "rowNum"]),
          as.character(plotDat[i, "colNum"])] <-
          ifelse(is.na(plotDat[i, trait]), NA, 1)
      }
      spatTrDat <- kronecker(M, matrix(data = 1, ncol = p1, nrow = p2)) *
        spatTr[["fit"]]
      ## Melt to get the data in ggplot shape. Rows and columns in the
      ## spatial trend coming from SpATS are swapped so therefore use t()
      plotDatSpat <- reshape2::melt(t(spatTrDat),
                                    varnames = c("colNum", "rowNum"))
      ## Add true values for columns and rows for plotting.
      plotDatSpat[["colNum"]] <- spatTr[["col.p"]]
      plotDatSpat[["rowNum"]] <- rep(x = spatTr[["row.p"]], each = p1 * nCol)
      ## Remove missings from data.
      plotDatSpat <- ggplot2::remove_missing(plotDatSpat, na.rm = TRUE)
      ## Divide by mean value of trait to get trend as percentage.
      plotDatSpat[["mean"]] <- mean(plotDat[[trait]], na.rm = TRUE)
      plotDatSpat[["value"]] <- plotDatSpat[["value"]] / plotDatSpat[["mean"]]
      ## Set values outside scale limits to NA.
      plotDatSpat[["value"]][plotDatSpat[["value"]] > scaleLim |
                               plotDatSpat[["value"]] < -scaleLim] <- NA
      return(plotDatSpat)
    })
    ## Extract all zVals to use identical limits for spatial pattern
    ## This enables proper comparison of plots over timePoints.
    zVals <- unlist(sapply(X = plotSpatDats, `[[`, "value"))
    if (is.infinite(scaleLim)) {
      zLim <- c(-1, 1) * max(c(abs(zVals), 0.1), na.rm = TRUE)
    } else {
      zLim <- c(-scaleLim, scaleLim)
    }
    ## Create a plot of the spatial trend per time point.
    for (i in seq_along(plotSpatDats)) {
      p <- fieldPlotPcts(plotDat = plotSpatDats[[i]], fillVar = "value",
                         title = names(plotSpatDats)[i],
                         colors = colorRampPalette(c("red", "yellow", "blue"),
                                                   space = "Lab")(100),
                         zlim = zLim, scaleLim = scaleLim)
      plot(p)
    }
  }, movie.name = outFile, autobrowse = FALSE, loop = 1)
}

#' Function for extracting objects of class fitMod that keeps class.
#'
#' @param x An object of class fitMod.
#' @param i An index specifying the element to extract of replace.
#' @param ... Ignored.
#'
#' @noRd
#' @export
`[.fitMod` <- function(x,
                       i,
                       ...) {
  timePoints <- chkTimePoints(x, i)
  timePointsX <- attr(x, which = "timePoints")
  timePointsR <- timePointsX[timePointsX[["timePoint"]] %in% timePoints, ]
  if (nrow(timePointsR) > 0) {
    class(x) <- "list"
    r <- x[timePointsR[["timePoint"]]]
    attr(r, "timePoints") <- timePointsR
    attr(r, "experimentName") <- attr(x, "experimentName")
    attr(r, "what") <- attr(x, "what")
    attr(r, "useRepId") <- attr(x, "useRepId")
    attr(r, "useCheck") <- attr(x, "useCheck")
    attr(r, "spatial") <- attr(x, "spatial")
    attr(r, "class") <- c("fitMod", "list")
    attr(r, "timestamp") <- attr(x, "timestamp")
  } else {
    r <- NULL
  }
  return(r)
}