Quant gen vignette

.Rbuildignore

@@ -5,4 +5,5 @@
 ^pkgdown$
 ^man/figures/.*\.jpg
 ^man/figures/.*\.png
+^vignettes/F2*.Rmd
 ^README.Rmd

R/Class-SimParamBee.R

@@ -599,12 +599,12 @@ isSimParamBee <- function(x) {
 #' @param n integer, number of samples
 #' @param average numeric, average number of workers
 #' @param lowerLimit numeric, returned numbers will be above this value
-#' @param queenTrait numeric, trait that represents queen's effect on the colony
-#'   phenotype (defined in \code{\link{SimParamBee}} - see examples); if
-#'   \code{0} then this effect is 0
-#' @param workersTrait numeric, trait that represents workers's effect on the
-#'   colony phenotype (defined in \code{\link{SimParamBee}} - see examples); if
-#'   \code{0} then this effect is 0
+#' @param queenTrait numeric (column position) or character (column name), trait
+#'   that represents queen's effect on the colony phenotype (defined in
+#'   \code{\link{SimParamBee}} - see examples); if \code{0} then this effect is 0
+#' @param workersTrait numeric (column position) or character (column name), trait
+#'   that represents workers's effect on the colony phenotype (defined in
+#'   \code{\link{SimParamBee}} - see examples); if \code{0} then this effect is 0
 #' @param checkProduction logical, does the phenotype depend on the production
 #'   status of colony; if yes and production is not \code{TRUE}, the result is
 #'   above \code{lowerLimit}
@@ -620,7 +620,8 @@ isSimParamBee <- function(x) {
 #'   traits influencing the colony phenotype and their parameters (mean and
 #'   variances) via \code{\link{SimParamBee}} (see examples).
 #'
-#' @seealso \code{\link{SimParamBee}} field \code{nWorkers}
+#' @seealso \code{\link{SimParamBee}} field \code{nWorkers} and
+#'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return numeric, number of workers
 #'
@@ -698,12 +699,12 @@ nWorkersColonyPhenotype <- function(colony, queenTrait = 1, workersTrait = NULL,
 #' @param n integer, number of samples
 #' @param average numeric, average number of drones
 #' @param lowerLimit numeric, returned numbers will be above this value
-#' @param queenTrait numeric, trait that represents queen's effect on the colony
-#'   phenotype (defined in \code{\link{SimParamBee}} - see examples); if
-#'   \code{0} then this effect is 0
-#' @param workersTrait numeric, trait that represents workers's effect on the
-#'   colony phenotype (defined in \code{\link{SimParamBee}} - see examples); if
-#'   \code{0} then this effect is 0
+#' @param queenTrait numeric (column position) or character (column name), trait
+#'   that represents queen's effect on the colony phenotype (defined in
+#'   \code{\link{SimParamBee}} - see examples); if \code{0} then this effect is 0
+#' @param workersTrait numeric (column position) or character (column name), trait
+#'   that represents workers's effect on the colony phenotype (defined in
+#'   \code{\link{SimParamBee}} - see examples); if \code{0} then this effect is 0
 #' @param checkProduction logical, does the phenotype depend on the production
 #'   status of colony; if yes and production is not \code{TRUE}, the result is
 #'   above \code{lowerLimit}
@@ -724,7 +725,8 @@ nWorkersColonyPhenotype <- function(colony, queenTrait = 1, workersTrait = NULL,
 #'   When \code{x} is \code{\link{Pop-class}}, only \code{workersTrait} is not
 #'   used, that is, only \code{queenTrait} is used.
 #'
-#' @seealso \code{\link{SimParamBee}} field \code{nDrones}
+#' @seealso \code{\link{SimParamBee}} field \code{nDrones} and
+#'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return numeric, number of drones
 #'
@@ -808,12 +810,12 @@ nDronesColonyPhenotype <- function(x, queenTrait = 1, workersTrait = NULL,
 #' @param n integer, number of samples
 #' @param average numeric, average number of virgin queens
 #' @param lowerLimit numeric, returned numbers will be above this value
-#' @param queenTrait numeric, trait that represents queen's effect on the colony
-#'   phenotype (defined in \code{\link{SimParamBee}} - see examples); if
-#'   \code{NULL} then this effect is 0
-#' @param workersTrait numeric, trait that represents workers's effect on the
-#'   colony phenotype (defined in \code{\link{SimParamBee}} - see examples); if
-#'   \code{NULL} then this effect is 0
+#' @param queenTrait numeric (column position) or character (column name), trait
+#'   that represents queen's effect on the colony phenotype (defined in
+#'   \code{\link{SimParamBee}} - see examples); if \code{NULL} then this effect is 0
+#' @param workersTrait numeric (column position) or character (column name), trait
+#'   that represents workers's effect on the colony phenotype (defined in
+#'   \code{\link{SimParamBee}} - see examples); if \code{NULL} then this effect is 0
 #' @param checkProduction logical, does the phenotype depend on the production
 #'   status of colony; if yes and production is not \code{TRUE}, the result is
 #'   above \code{lowerLimit}
@@ -833,7 +835,8 @@ nDronesColonyPhenotype <- function(x, queenTrait = 1, workersTrait = NULL,
 #'   influencing the colony phenotype and their parameters (mean and variances)
 #'   via \code{\link{SimParamBee}} (see examples).
 #'
-#' @seealso \code{\link{SimParamBee}} field \code{nVirginQueens}
+#' @seealso \code{\link{SimParamBee}} field \code{nVirginQueens} and
+#'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return numeric, number of virgin queens
 #'
@@ -853,9 +856,8 @@ nDronesColonyPhenotype <- function(x, queenTrait = 1, workersTrait = NULL,
 #' # Example for nVirginQueensColonyPhenotype()
 #' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100)
 #' SP <- SimParamBee$new(founderGenomes)
-#' # TODO: Find some means and variances (Vg and Ve) for swarming in honeybees
-#' #       https://github.com/HighlanderLab/SIMplyBee/issues/259
-#' meanP <- c(10, 0)
+#' # Setting trait scale such that mean is 10 split into queen and workers effects
+#' meanP <- c(5, 5 / SP$nWorkers)
 #' # setup variances such that the total phenotype variance will match the mean
 #' varA <- c(3 / 2, 3 / 2 / SP$nWorkers)
 #' corA <- matrix(data = c(
@@ -1189,38 +1191,51 @@ downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) {
 #'   own "caste functions" that satisfy your needs within this mapping function
 #'   (see \code{queenFUN}, \code{workersFUN}, and \code{dronesFUN} below)
 #'   or provide a complete replacement of this mapping function! For example,
-#'   this mappign function does not cater for indirect (social) genetic effects
+#'   this mapping function does not cater for indirect (social) genetic effects
 #'   where colony individuals value impacts value of other colony individuals.
-#'   Note though that you can achieve this also via multiple correlated traits,
-#'   such as a queen and a workers trait.
+#'   Note though that you can achieve this impact also via multiple correlated
+#'   traits, such as a queen and a workers trait.
 #'
 #' @param colony \code{\link{Colony-class}}
 #' @param value character, one of \code{pheno} or \code{gv}
-#' @param queenTrait numeric, trait that represents queen's contribution
-#'   to the colony value; if \code{NULL} then this contribution is 0
-#' @param queenFUN function, function that will be applied to the queen's value
-#' @param workersTrait numeric, trait that represents workers' contribution
-#'   to the colony value; if \code{NULL} then this contribution is 0
-#' @param workersFUN function, function that will be applied to the worker values
-#' @param dronesTrait numeric, trait that represents drones' contribution
-#'   to the colony value; if \code{NULL} then this contribution is 0
-#' @param dronesFUN function, function that will be applied to the drone values
+#' @param queenTrait numeric (column position) or character (column name),
+#'   trait(s) that represents queen's contribution to colony value(s); if
+#'   \code{NULL} then this contribution is 0; you can pass more than one trait
+#'   here, but make sure that \code{combineFUN} works with these trait dimensions
+#' @param queenFUN function, function that will be applied to queen's value
+#' @param workersTrait numeric (column position) or character (column name),
+#'   trait(s) that represents workers' contribution to colony value(s); if
+#'   \code{NULL} then this contribution is 0; you can pass more than one trait
+#'   here, but make sure that \code{combineFUN} works with these trait dimensions
+#' @param workersFUN function, function that will be applied to workers values
+#' @param dronesTrait numeric (column position) or character (column name),
+#'   trait(s) that represents drones' contribution to colony value(s); if
+#'   \code{NULL} then this contribution is 0; you can pass more than one trait
+#'   here, but make sure that \code{combineFUN} works with these trait dimensions
+#' @param dronesFUN function, function that will be applied to drone values
+#' @param traitName, the name of the colony trait(s), say, honeyYield; you can pass
+#'   more than one trait name here, but make sure to match them with
+#'   \code{combineFUN} trait dimensions
 #' @param combineFUN, function that will combine the queen, worker, and drone
 #'   contributions - this function should be defined as \code{function(q, w, d)}
-#'   where \code{q} represents queen's, \code{q} represents workers',
+#'   where \code{q} represents queen's, \code{q} represents workers', and
 #'   \code{d} represents drones' contribution.
 #' @param checkProduction logical, does the value depend on the production
-#'   status of colony; if yes and production is not \code{TRUE}, the return
+#'   status of colony; if yes and production is \code{FALSE}, the return
 #'   is \code{notProductiveValue} - this will often make sense for colony
-#'   phenotype value only
-#' @param notProductiveValue scalar, returned value when colony is not productive
+#'   phenotype value only; you can pass more than one logical value here (one
+#'   per trait coming out of \code{combineFUN})
+#' @param notProductiveValue numeric, returned value when colony is not productive;
+#'   you can pass more than one logical value here (one per trait coming out of
+#'   \code{combineFUN})
 #' @param simParamBee \code{\link{SimParamBee}}, global simulation parameters
 #' @param ... other arguments of \code{mapCasteToColonyValue} (for its aliases)
 #'
 #' @seealso \code{\link{SimParamBee}} field \code{colonyValueFUN} and functions
 #'   \code{\link{calcColonyValue}}, \code{\link{calcColonyPheno}},
 #'   \code{\link{calcColonyGv}}, \code{\link{getEvents}},
-#'   \code{\link{pheno}}, and \code{\link{gv}}
+#'   \code{\link{pheno}}, and \code{\link{gv}}, as well as
+#'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @details This is a utility/mapping function meant to be called by
 #'   \code{\link{calcColonyValue}}. It only works on a single colony - use
@@ -1237,8 +1252,8 @@ downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) {
 #' # 2) workers' effect on colony honey yield, say via foraging ability phenotype
 #' # The traits will have a negative genetic correlation of -0.5 and heritability
 #' # of 0.25 (on an individual level)
-#' mean <- c(20, 0)
 #' nWorkers <- 10
+#' mean <- c(10, 10 / nWorkers)
 #' varA <- c(1, 1 / nWorkers)
 #' corA <- matrix(data = c(
 #'   1.0, -0.5,
@@ -1278,17 +1293,16 @@ downsizePUnif <- function(colony, n = 1, min = 0.8, max = 0.9) {
 #' getWorkersGv(colony)
 #'
 #' @export
-# TODO: Calculate inheritance, selection and production criteria in the Colony #23
-#       https://github.com/HighlanderLab/SIMplyBee/issues/23
 # TODO: Do we need to do anything to add GxE to colony values? #353
 #       https://github.com/HighlanderLab/SIMplyBee/issues/353
 # TODO: Develop theory for colony genetic values under non-linearity/non-additivity #403
 #       https://github.com/HighlanderLab/SIMplyBee/issues/403
 mapCasteToColonyValue <- function(colony,
                                   value = "pheno",
                                   queenTrait = 1, queenFUN = function(x) x,
-                                  workersTrait = 2, workersFUN = sum,
+                                  workersTrait = 2, workersFUN = colSums,
                                   dronesTrait = NULL, dronesFUN = NULL,
+                                  traitName = NULL,
                                   combineFUN = function(q, w, d) q + w,
                                   checkProduction = TRUE, notProductiveValue = 0,
                                   simParamBee = NULL) {
@@ -1336,8 +1350,22 @@ mapCasteToColonyValue <- function(colony,
     dronesEff <- dronesFUN(tmp)
   }
   colonyValue <- combineFUN(q = queenEff, w = workersEff, d = dronesEff)
-  if (checkProduction && !colony@production) {
-    colonyValue <- notProductiveValue
+  nColTrt <- length(colonyValue)
+  colnames(colonyValue) <- traitName
+  if (any(checkProduction) && !isProductive(colony)) {
+    if (length(checkProduction) == 1 && nColTrt != 1) {
+      checkProduction <- rep(checkProduction, times = nColTrt)
+    }
+    if (length(notProductiveValue) == 1 && nColTrt != 1) {
+      notProductiveValue <- rep(notProductiveValue, times = nColTrt)
+    }
+    if (length(checkProduction) != nColTrt) {
+      stop("Dimension of checkProduction does not match the number of traits from combineFUN()!")
+    }
+    if (length(checkProduction) != length(notProductiveValue)) {
+      stop("Dimensions of checkProduction and notProductiveValue must match!")
+    }
+    colonyValue[checkProduction] <- notProductiveValue[checkProduction]
   }
   return(colonyValue)
 }