Modified fitModels to better handle combinations of geno.decomp and c…

…ovariates

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method("[",TP)
 S3method(plot,TP)
 S3method(plot,fitMod)
 export(createTimePoints)

R/createFitMod.R

@@ -56,7 +56,7 @@ plot.fitMod <- function(x,
       corrected <- droplevels(corrected)
     }
     ## Add combinations missing in data to corrected.
-    corrected <- addMissVals(dat = corrected, trait = trait)
+    corrected <- addMissVals(dat = corrected, trait = "newTrait")
     xyFacetPlot(baseDat = corrected, overlayDat = preds, yVal = "newTrait",
                 yValOverlay = "predicted.values", title = title, yLab = trait)
   } else if (plotType == "herit") {

R/createTP.R

@@ -527,3 +527,16 @@ plot.TP <- function(x,
   }
   invisible(p)
 }
+
+#' Function for extracting for objects of class TP that keeps class.
+#'
+#' @export
+`[.TP` <- function(x, i, ...) {
+  r <- NextMethod("[")
+  attr(r, "class") <- attr(x, "class")
+  return(r)
+}
+
+
+
+

R/fitModels.R

@@ -4,39 +4,66 @@ fitModels <- function(TP,
                       covariates = NULL,
                       geno.decomp = NULL,
                       useCheck = FALSE) {
-  # In these analyses, the genotype is always included as random
-  # I have to check with Coté if we can also include it as fixed
-  genotype.as.random <- TRUE
-  ## All covariates should be factors.
-  for (covar in covariates) {
-    TP <- lapply(X = TP, FUN = function(timePoint) {
-      if (!is.factor(timePoint[[covar]])) {
-        timePoint[[covar]] <- as.factor(timePoint[[covar]])
+  ## Checks
+  if (!inherits(TP, "TP")) {
+    stop("TP should be an object of class TP.\n")
+  }
+  if (!all(sapply(X = TP, FUN = hasName, name = trait))) {
+    stop(trait, " should be a column for all timePoints.\n")
+  }
+  if (!is.null(covariates)) {
+    for (covar in covariates) {
+      if (!all(sapply(X = TP, FUN = hasName, name = covar))) {
+        stop(covar, " should be a column for all timePoints.\n")
       }
+    }
+  }
+  if (!is.null(geno.decomp)) {
+    for (gd in geno.decomp) {
+      if (!all(sapply(X = TP, FUN = hasName, name = gd))) {
+        stop(gd, " should be a column for all timePoints.\n")
+      }
+    }
+  }
+  ## If geno.decomp is used genotype and covariates have to be replaced by
+  ## an interaction of genotype and covariates with the geno.decomp variables.
+  ## Construct an interaction of all variables in geno.decomp.
+  if (length(geno.decomp) > 1) {
+    TP <- lapply(X = TP, FUN = function(timePoint) {
+      timePoint[["geno.decomp"]] <- interaction(timePoint[geno.decomp],
+                                                sep = "_")
       return(timePoint)
     })
+    ## Set geno.decomp to newly constructed variable.
+    geno.decomp <- "geno.decomp"
   }
-  ## If geno.decomp is used genotype and covariates have to be replaced by
-  ## an interaction of genotype and covariates with the geno.decomp variable.
-  ## Add these new variables to the data.
+  ## Replace genotype and covariates by their interaction with geno.decomp.
   if (!is.null(geno.decomp)) {
     TP <- lapply(X = TP, FUN = function(timePoint) {
       timePoint[["genotype"]] <- interaction(timePoint[[geno.decomp]],
-                                         timePoint[["genotype"]], sep = "_")
+                                             timePoint[["genotype"]], sep = "_")
       for (covar in covariates) {
         timePoint[[covar]] <- interaction(timePoint[[geno.decomp]],
                                           timePoint[[covar]], sep = "_")
       }
       return(timePoint)
     })
   }
+  ## All covariates should be factors. If not convert them to factor.
+  for (covar in covariates) {
+    TP <- lapply(X = TP, FUN = function(timePoint) {
+      if (!is.factor(timePoint[[covar]])) {
+        timePoint[[covar]] <- as.factor(timePoint[[covar]])
+      }
+      return(timePoint)
+    })
+  }
   ## Construct formula for fixed part.
-  if (!is.null(covariates)) {
-    fixedForm <- formula(paste("~", paste(covariates, collapse = "+"),
-                         if (useCheck) "+ check"))
-    if (!is.null(geno.decomp)) {
-      geno.decomp <- covariates
-    }
+  ## Fixed part consists of covariates, geno.decomp and check.
+  if (!is.null(c(covariates, geno.decomp))) {
+    fixedForm <- formula(paste("~", paste(c(covariates, geno.decomp),
+                                          collapse = "+"),
+                               if (useCheck) "+ check"))
   } else {
     fixedForm <- if (useCheck) formula("~ check") else NULL
   }
@@ -45,19 +72,20 @@ fitModels <- function(TP,
     message(timePoint[["timePoint"]][1])
     ## Only keep columns needed for analysis.
     modCols <- c("timePoint", "plotId", "genotype", "genoCheck", "check",
-                 "colId", "rowId", "colNum", "rowNum", covariates, trait)
+                 "colId", "rowId", "colNum", "rowNum", covariates, geno.decomp,
+                 trait)
     modDat <- timePoint[colnames(timePoint) %in% modCols]
     modDat <- droplevels(modDat)
     ## number of segments for SpATS.
     nseg = c(nlevels(modDat[["colId"]]), nlevels(modDat[["rowId"]]))
-    ## Fit the model.
+    ## Fit and return the model.
     SpATS::SpATS(response = trait, fixed = fixedForm,
                  random = ~ colId + rowId,
                  spatial = ~ SpATS::PSANOVA(colNum, rowNum, nseg = nseg,
                                             nest.div = c(2, 2)),
                  genotype = if (useCheck) "genoCheck" else "genotype",
-                 genotype.as.random = genotype.as.random,
-                 geno.decomp = geno.decomp, data = modDat,
+                 genotype.as.random = TRUE, geno.decomp = geno.decomp,
+                 data = modDat,
                  control = list(maxit = 50, tolerance = 1e-03, monitoring = 0))
   })
   return(createFitMod(fitMods))

testScript.R

@@ -1,19 +1,17 @@
-folder <- "."
-setwd(paste0(folder, "output/"))
+setwd(paste0("./output/"))
 
 inDat <- data.table::fread("../data-raw/Original_PAM_reshape.csv",
-                               data.table = FALSE)
+                           data.table = FALSE)
 
 # creating a unique ID per plant using the row and col coordinate
 # (in principle the column "Sowing_Position" was also a unique ID but I wanted to see the position)
 inDat$ID <- interaction(inDat[["x"]], inDat[["y"]], sep = "_")
 inDat <- inDat[!is.na(inDat$pheno), ]
 # Create an indicator for each plot (according to the row and column position)
-inDat$pos <- paste0("c", inDat$x, "r", inDat$y)
+inDat$pos <- paste0("c", inDat[["x"]], "r", inDat[["y"]])
 
 # I removed a plant that has very few measurements
 inDat <- inDat[inDat$pos != "c1r54",]
-inDat <- droplevels(inDat)
 
 inTP <- createTimePoints(dat = inDat, genotype = "Genotype",
                          timePoint = "timepoints", repId = "Sowing_Block",
@@ -36,7 +34,7 @@ dev.off()
 plot(inTP, plotType = "raw", traits = "pheno",
      genotypes = c("col", "ely", "evo1", "ler"))
 
-fitMods <- fitModels(TP = inTP, trait = "pheno",
+fitMods <- fitModels(TP = inTP[1:3], trait = "pheno",
                      covariates = c("repId", "Image_pos"))
 
 genoPreds <- getGenoPred(fitMods, outFile = "BLUPs_PAM_modRep.csv")
@@ -53,13 +51,12 @@ plot(fitMods, plotType = "timeLapse", outFile = "spatialtrends.gif")
 
 ## Second example
 inDat2 <- data.table::fread("../data-raw/Data_modif_ZA17_anonymous.csv",
-                           data.table = FALSE)
-
+                            data.table = FALSE)
 # Create an indicator for each plot (according to the row and column position)
 inDat2$pos <- paste0("c", inDat2$Line, "r", inDat2$Position)
 
 inTP2 <- createTimePoints(dat = inDat2, genotype = "geno", timePoint = "Date",
-                  plotId = "pos", rowNum = "Position", colNum = "Line")
+                          plotId = "pos", rowNum = "Position", colNum = "Line")
 
 plot(inTP2, plotType = "layout", timePoints = "2017-04-13")
 plot(inTP2, plotType = "cor", traits = "LA_Estimated")
@@ -69,8 +66,11 @@ pdf("Phenovator_ZA17_raw_data_na.pdf", height = 8, width = 12)
 plot(inTP2, plotType = "raw", traits = "LA_Estimated", geno.decomp = "Scenario")
 dev.off()
 
-fitMods2 <- fitModels(TP = inTP2, trait = "LA_Estimated",
-                      geno.decomp = "Scenario", covariates = "population")
+fitMods2 <- fitModels(TP = inTP2[1:3], trait = "LA_Estimated",
+                      geno.decomp = c("Scenario", "population"))
+## This crashes:
+fitMods2b <- fitModels(TP = inTP2[1:3], trait = "LA_Estimated",
+                       geno.decomp = "Scenario", covariates = "population")
 
 genoPreds2 <- getGenoPred(fitMods2, outFile = "BLUPs_ZA17_LeafArea.csv")
 spatCorr2 <- getCorrected(fitMods2, outFile = "Corrected_ZA17_LeafArea.csv")