F_Quantitative_Genetics.Rmd

---
title: "Quantitative genetics"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quantitative genetics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
editor_options:
  markdown:
    wrap: 80
    canonical: true
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  include = TRUE
)
```

# Introduction

This vignette describes and demonstrates how SIMplyBee implements quantitative
genetics principles for honeybees. Specifically, it describes three different
examples where we simulate:

1.  Honey yield - a single colony trait,

2.  Honey yield and Calmness - two colony traits, and

3.  Colony strength and Honey yield - two colony traits where one trait impacts
    the other one via the number of workers.

We start by loading SIMplyBee and quickly simulating genomes for some founder
honeybees. Specifically, we will simulate genomes for 20 individuals with 16
chromosomes and 1000 segregating sites per chromosome.

```{r founderGenomes}
library(package = "SIMplyBee")
library(package = "ggplot2")
founderGenomes <- quickHaplo(nInd = 20, nChr = 16, segSites = 1000)
```

# Honey yield

This section shows how to simulate one colony trait, honey yield, that is
influenced by the queen and workers as well as the environment. We will achieve
this by:

a.  setting base population quantitative genetic parameters,
b.  inspecting individual values in the base population,
c.  inspecting individual values in a colony,
d.  calculating colony value,
e.  calculating multi-colony values, and
f.  selecting on colony values.

## Base population quantitative genetic parameters

AlphaSimR, and hence SIMplyBee, simulates each individual with its corresponding
genome, and quantitative genetic and phenotypic values. To enable this
simulation, we must set base population quantitative genetic parameters for the
traits of interest in the global simulation parameters via `SimParamBee`. We
must set:

1)  the number of traits,
2)  the number of quantitative trait loci (QTL) that affect the traits,
3)  the distribution of QTL effects,
4)  trait means, and
5)  trait genetic and environmental variances - if we simulate multiple traits,
    we must also specify genetic and environmental covariances between the
    traits.

In honeybees, the majority of traits are influenced by the queen and workers.
There are many biological mechanisms for these queen and workers effects.
Depending on which caste is the main driver of the trait (the queen or workers),
we also talk about direct and indirect effects. For example, for honey yield,
workers directly affect honey yield by foraging, while the queen indirectly
affects honey yield by stimulating workers via pheromone production. The queen
and workers effects for a trait can be genetically and environmentally
independent or correlated (usually negatively).

Here, we will simulate two traits to represent the queen and workers effects on
honey yield. From this point onward we will use the terms the queen effect and
queen trait interchangeably. The same applies to workers effect and workers
trait. These two effects (=traits) will give rise to honey yield trait. We will
assume that colony honey yield is approximately normally distributed with the
mean of 20 kg and variance of 4 $kg^2$, which implies that most colonies will
have honey yield between 14 kg and 26 kg (see
`hist(rnorm(n = 1000, mean = 20, sd = sqrt(4)))`). Traits like honey yield have
a complex polygenic genetic architecture, so we will assume that this trait is
influenced by 100 QTL per chromosome (with 16 chromosomes, this gives us 1600
QTL in total).

We will first initiate global simulation parameters and set the mean of queen
effects to 10 kg with genetic variance of 1 $kg^2$, while we will set the mean
of workers effects to 10 kg with genetic variance of 1 $kg^2$. The mean and the
variance for the worker effect are proportionally scaled by the expected number
of workers in a colony. The mean and variance for the queen effect is assumed
larger than for the workers effect, because there is one queen and many workers
in colony and we assume that workers effects "accumulate". Deciding how to split
the colony mean between queen and workers effects will depend on the individual
to colony mapping function, which we will describe in the Colony value
sub-section.

```{r SimParamBee_mean_and_varA}
# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

nQtlPerChr <- 100

# Genetic parameters for queen and workers effects - each represented by a trait
mean <- c(10, 10 / SP$nWorkers)
varA <- c(1, 1 / SP$nWorkers)
```

We next set genetic correlation between the queen and workers effects to -0.5 to
reflect the commonly observed antagonistic relationship between these effects.
With all the quantitative genetic parameters defined, we now add two additive
traits to global simulation parameters and name them `queenTrait` and
`workerTrait`. These parameters drive the simulation of QTL effects. Read about
all the other trait simulation options in AlphaSimR via:
`vignette(topic = "traits", package="AlphaSimR")`.

```{r SimParamBee_corA_and_addTrait}
corA <- matrix(data = c( 1.0, -0.5, 
                        -0.5,  1.0), nrow = 2, byrow = TRUE)
SP$addTraitA(nQtlPerChr = nQtlPerChr, mean = mean, var = varA, corA = corA,
             name = c("queenTrait", "workersTrait"))
```

Finally, we set the environmental variance of the queen and workers effects to 3
$kg^2$ and we again scale the worker variance by the expected number of workers.
Contrary to the negative genetic correlation, we here assume that environmental
correlation between the queen and workers effects is slightly positive, 0.3.
This is just an example! These parameters should be based on literature or
simulation scenarios of interest.

```{r SimParamBee_varE_and_corR}
varE <- c(3, 3 / SP$nWorkers)
corE <- matrix(data = c(1.0, 0.3, 
                        0.3, 1.0), nrow = 2, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)
```

## Individual values in the base population

Now we create a base population of virgin queens. Since we defined two traits,
all honeybees in the simulation will have genetic and phenotypic values for both
traits. The genetic values are stored in the `gv` slot of each `Pop` object,
while phenotypic values are stored in the `pheno` slot.

```{r basePop_virgin_queens, echo = FALSE, fig.height = 5, fig.width = 6}
# Base population virgin queens
basePop <- createVirginQueens(founderGenomes, n = 20)
head(basePop@gv)
head(basePop@pheno)
oldpar <- par(mfrow=c(2,2))
limQ <- range(c(basePop@gv[, "queenTrait"], basePop@pheno[, "queenTrait"]))
brkQ <- seq(from = limQ[1], to = limQ[2], length.out = 10)
limW <- range(c(basePop@gv[, "workersTrait"], basePop@pheno[, "workersTrait"]))
brkW <- seq(from = limW[1], to = limW[2], length.out = 10)
hist(basePop@gv[, "queenTrait"], xlab = "Genetic value", main = "Queen effect", xlim = limQ, breaks = brkQ)
hist(basePop@gv[, "workersTrait"], xlab = "Genetic value", main = "Workers effect", xlim = limW, breaks = brkW)
hist(basePop@pheno[, "queenTrait"], xlab = "Phenotypic value", main = "Queen effect", xlim = limQ, breaks = brkQ)
hist(basePop@pheno[, "workersTrait"], xlab = "Phenotypic value", main = "Workers effect", xlim = limW, breaks = brkW)
par(oldpar)
```

Note that these are virgin queens, yet we obtained queen and workers effect
values for them! Is this wrong? No! Virgin queens carry DNA with genes that are
differentially expressed in different castes, which would be only showed in
their phenotype. Hence, virgin queens have genetic values for the queen and
worker effects, but they might never actually express these effects. In this
simulation virgin queens also obtained phenotypic values for both of the
effects. This is technically incorrect because virgin queens don't express genes
for the worker effect at all, and they also do not express the queen effect, not
until they become the queen of a colony. We can treat these phenotypic values
for virgin queens as values that we could see if these virgin queens would
express these traits. We will show later in the Colony value sub-section how we
use these traits from different castes. If existence of these phenotypic values
for certain castes is a hindrance, we can always remove them for population or
colony objects by modifying the corresponding slots as required.

As with the virgin queens, drones also carry DNA with genes that are expressed
in different castes. Therefore, drones will also have the queen and workers
effect genetic (and phenotypic values) for honey yield even though they do not
contribute to this trait in a colony.

```{r basePop_drones, echo = FALSE, fig.height = 2.7, fig.width = 5}
# Base population drones
drones <- createDrones(x = basePop[1:5], nInd = 3)
head(drones@gv)
oldpar <- par(mfrow=c(1,2))
hist(drones@gv[, "queenTrait"], xlab = "Genetic value", main = "Queen effect")
hist(drones@gv[, "workersTrait"], xlab = "Genetic value", main = "Workers effect")
par(oldpar)
```

## Individual values in a colony

We continue by creating a colony from one base population virgin queen, crossing
it, and adding some workers.

```{r create_colony}
colony <- createColony(x = basePop[6])
colony <- cross(x = colony, drones = drones, checkCross = "warning")
colony <- addWorkers(x = colony, nInd = 50)
colony
```

We can access the genetic and phenotypic values of colony members with functions
`getGv()` and `getPheno()`, both of which have the `caste` argument (see more
via `help(getGv)`).

```{r getGv_and_getPheno}
getGv(colony, caste = "queen")
getGv(colony, caste = "workers") |> head(n = 4)

getPheno(colony, caste = "queen")
getPheno(colony, caste = "workers") |> head(n = 4)
```

For convenience, there are also alias functions for accessing the genetic and
phenotypic values of each caste directly.

```{r getGv_and_getPheno_caste}
getQueenGv(colony)
getWorkersGv(colony) |> head(n = 4)

getQueenPheno(colony)
getWorkersPheno(colony) |> head(n = 4)
```

Some phenotypes, such as honey yield, are only expressed if colony is at full
size. This is achieved by the `buildUp()` colony event function that adds worker
and drones and hence turns on the `production` status of the colony (to `TRUE`).
SIMplyBee includes a function `ìsProductive()` to check the production status of
a colony.

```{r build_up_colony}
# Check if colony is productive
isProductive(colony)

# Build-up the colony and check the production status again
colony <- buildUp(colony)
colony
isProductive(colony)
```

For the ease of further demonstration, we now combine workers' values into a
single data.frame.

```{r data.frame}
# Collate genetic and phenotypic values of workers
df <- data.frame(id = colony@workers@id,
                 mother = colony@workers@mother,
                 father = colony@workers@father,
                 gvQueenTrait = colony@workers@gv[, "queenTrait"],
                 gvWorkersTrait = colony@workers@gv[, "workersTrait"],
                 pvQueenTrait =  colony@workers@pheno[, "queenTrait"],
                 pvWorkersTrait = colony@workers@pheno[, "workersTrait"])
head(df)
```

To visualise correlation between queen and workers effects in workers, we plot
these effect values against each other.

```{r plot_queen_vs_worker_values}
# Covariation between queen and workers effect genetic values in workers
p <- ggplot(data = df, aes(x = gvQueenTrait, y = gvWorkersTrait)) +
  xlab("Genetic value for the queen effect") +
  ylab("Genetic value for the workers effect") +
  geom_point() +
  theme_classic()
print(p)
```

In SIMplyBee, we know genetic values of all individuals, including drones that
the queen mated with (=fathers in a colony)!

```{r fathers_values}
# Variation in patriline genetic values
getFathersGv(colony)
```

Knowing the father of each worker, we inspect variation in the distribution of
genetic values of worker by the patriline (workers from a single father drone)
for the workers effect.

```{r distribution_by_patriline, echo = FALSE, fig.height = 4.5, fig.width = 6}
# Variation in workers effect genetic values by patriline in workers
p <- ggplot(data = df, aes(x = gvWorkersTrait, colour = father)) +
  xlab("Genetic value for the workers effect") +
  geom_density() +
  theme_classic()
print(p)
```

## Colony value

However, in honeybees we usually don't observe values on individuals, but on a
colony. SIMplyBee provides functions for mapping individual values to a colony
value. The general function for this is `calcColonyValue()`, which can combine
any value and trait from any caste. There are also aliases `calcColonyGv()` and
`calcColonyPheno()`. These functions require users to specify the so-called
mapping function (via the `FUN` argument). The mapping function specifies queen
and workers traits (potentially also drone traits) and what function we want to
apply to each of them before mapping them to the colony value(s). We can also
specify whether the colony value(s) depend on the production status. For
example, if a colony is not productive, its honey yield would be 0 or
unobserved. SIMplyBee provides a general mapping function
`mapCasteToColonyValue()` and aliases `mapCasteToColonyGv()` and
`mapCasteToColonyPheno()`. These functions have arguments to cater for various
situations. By default, they first calculate caste values: leave the queen's
value as it is, sum workers' values, potentially sum drones' values, and lastly
sum all these caste values together into a colony value. Users can provide their
own mapping function(s) too!

We now calculate honey yield for our colony - a single value for the colony.

```{r colony_pheno}
# Colony phenotype value
calcColonyPheno(colony, queenTrait = "queenTrait", workersTrait = "workersTrait")
help(calcColonyPheno)
help(mapCasteToColonyPheno)
```

These colony values are not stored in a colony, because they change as colony
changes due to various events. For example, reducing the number of workers will
reduce the colony honey yield.

```{r colony_pheno_change}
# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.5) |>
  calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")
```

Please note that we assumed that the queen contributes half to colony honey
yield and workers contribute the other half. This means that removing workers
will still give a non-zero honey yield! This shows that we have to design the
mapping between individual, caste, and colony values with care!

```{r colony_pheno_change1}
# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.99) |>
  calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")
```

Finally, note that SIMplyBee currently does not provide functionality for
breeding values, dominance deviations, and epistatic deviations at caste and
colony levels, despite the availabiliy of AlphaSimR `bv()`, `dd()`, and `aa()`
functions. This is because we have to check or develop theory on how to
calculate these values across active colonies and hence we currently advise
against the use of AlphaSimR `bv()`, `dd()`, and `aa()` functions with SIMplyBee
as the output of these functions could be easily misinterpreted.

## MultiColony values

The same functions can be used on a `MultiColony` class object. Let's create an
apiary.

```{r multicolony}
apiary <- createMultiColony(basePop[7:20])
drones <- createDrones(basePop[1:5], nInd = 100)
droneGroups <- pullDroneGroupsFromDCA(drones, n = nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
apiary <- buildUp(apiary)
```

We can extract the genetic and phenotypic values from multiple colonies in the
same manner as from a single colony, by using `get*Gv()` and `get*Pheno()`
functions. The output of these function is a named list with values for each
colony or a single matrix if we set the `collapse` argument to `TRUE`.

```{r multicolony_gv}
getQueenGv(apiary) |> head(n = 4)
getQueenGv(apiary, collapse = TRUE) |> head(n = 4)
```

In a similar manner, we can calculate colony value for all the colonies in our
apiary, where the row names of the output represent colony IDs.

```{r multicolony_pheno}
colonyGv <- calcColonyGv(apiary)
colonyPheno <- calcColonyPheno(apiary)
data.frame(colonyGv, colonyPheno)
```

## Selection on colony values

Since the aim of selection is to select the best individuals or colonies for the
reproduction, we could select the best colony in our apiary based on either
genetic or phenotypic value for grafting the new generation of virgin queens. We
can use the function `selectColonies()` that takes a matrix of colony values
(the output of `calcColonyValue()` function). The default behavior is to select
the colonies with the highest value (argument `selectTop` set to `TRUE`), but
you can also select the colonies with the lowest values (argument `selectTop`
set to `FALSE`).

```{r multicolony_selection}
# Select the best colony based on gv
selectColonies(apiary, n = 1, by = colonyGv)
# Select the best colony based on phenotype
selectColonies(apiary, n = 1, by = colonyPheno)
```

The same functionality is implemented in `pullColonies()` and
`removeColonies()`.

# Honey yield and Calmness

In this section we expand simulation to two uncorrelated colony traits with
queen and workers effects, honey yield and calmness. We follow the same recipe
as in the previous section where we simulated only one colony trait.

We first reinitialize the global simulation parameters because we will define
new traits. For honey yield we will use the same parameters as before, while for
calmness trait we will assume that the trait is scored continuously in such a
way that negative values are undesirable and positive values are desirable with
zero being population mean. We will further assume the same variances for
calmness as for honey yield, and a genetic (and environmental) correlation
between the queen and workers effects of -0.4 (and 0.2) for calmness. We assume
no genetic or environmental correlation between honey yield and calmness.
Beware, this is just an example to show you how to simulate multiple colony
traits - we have made up these parameters - please use literature estimates in
your simulations!

```{r SimParamBee_2}
# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

nQtlPerChr <- 100

# Quantitative genetic parameters - for two traits, each with the queen and workers effects
meanP <- c(10, 10 / SP$nWorkers, 0, 0)
varA <- c(1, 1 / SP$nWorkers, 1, 1 / SP$nWorkers)
corA <- matrix(data = c( 1.0, -0.5,  0.0,  0.0, 
                        -0.5,  1.0,  0.0,  0.0,
                         0.0,  0.0,  1.0, -0.4, 
                         0.0,  0.0, -0.4,  1.0), nrow = 4, byrow = TRUE)
SP$addTraitA(nQtlPerChr = 100, mean = meanP, var = varA, corA = corA,
             name = c("yieldQueenTrait", "yieldWorkersTrait",
                      "calmQueenTrait", "calmWorkersTrait"))

varE <- c(3, 3 / SP$nWorkers, 3, 3 / SP$nWorkers)
corE <- matrix(data = c(1.0, 0.3, 0.0, 0.0,
                        0.3, 1.0, 0.0, 0.0,
                        0.0, 0.0, 1.0, 0.2,
                        0.0, 0.0, 0.2, 1.0), nrow = 4, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)
```

We continue by creating a base population of virgin queens and from them an
apiary with 10 full-sized colonies.

```{r base_pop_and_colony}
basePop <- createVirginQueens(founderGenomes)
drones <- createDrones(x = basePop[1:5], nInd = 100)
apiary <- createMultiColony(basePop[6:20])
droneGroups <- pullDroneGroupsFromDCA(drones, nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
apiary <- buildUp(apiary)
apiary
```

We can again inspect the genetic (and phenotypic) values of all individuals in
each colony and whole apiary with `get*Gv()` and `get*Pheno()` functions. Now,
the output contains four traits representing the queen and workers effect for
honey yield and calmness. These functions also take an `nInd` argument to sample
a number of individuals along with their values.

```{r}
getQueenGv(apiary) |> head(n = 4)
getWorkersPheno(apiary, nInd = 3) |> head(n = 4)
```

Now, we calculate colony genetic and phenotypic values for all colonies in the
apiary. Since we are simulating two traits, honey yield and calmness, we have
two ways to calculate corresponding colony values. The first way is to use the
default `mapCasteToColony*()` function in `calcColony*()` and only define
additional arguments as shown here:

```{r colony_pheno_2a}
colonyValues <- calcColonyPheno(apiary,
                                queenTrait = c("yieldQueenTrait", "calmQueenTrait"),
                                workersTrait = c("yieldWorkersTrait", "calmWorkersTrait"),
                                traitName = c("yield", "calmness"),
                                checkProduction = c(TRUE, FALSE)) |> as.data.frame()
colonyValues
```

The second way is to create our own mapping function. An equivalent outcome to
the above is shown below just to demonstrate use of your own function, but we
are simply just reusing `mapCasteToColonyPheno()` twice;)

```{r colony_pheno_2b}
myMapCasteToColonyPheno <- function(colony) {
  yield <- mapCasteToColonyPheno(colony,
                                 queenTrait = "yieldQueenTrait",
                                 workersTrait = "yieldWorkersTrait",
                                 traitName = "yield",
                                 checkProduction = TRUE)
  calmness <- mapCasteToColonyPheno(colony,
                                    queenTrait = "calmQueenTrait",
                                    workersTrait = "calmWorkersTrait",
                                    traitName = "calmness",
                                    checkProduction = FALSE)
  return(cbind(yield, calmness))
}
colonyValues <- calcColonyPheno(apiary, FUN = myMapCasteToColonyPheno) |> as.data.frame()
colonyValues
```

Again, we can now select the best colony based on the best phenotypic value for
either yield, calmness, or an index of both. Let's say that both traits are
equally important so we select on a weighted sum of both of them - we will use
the AlphaSimR `selIndex()` function that enables this calculation along with
scaling. We will represent the index such that it has a mean of 100 and standard
deviation of 10 units.

```{r}
colonyValues$Index <- selIndex(Y = colonyValues, b = c(0.5, 0.5), scale = TRUE) * 10 + 100
bestColony <- selectColonies(apiary, n = 1, by = colonyValues$Index)
getId(bestColony)
```

We see that we selected colony with ID "4", but we would be selecting a
different colony based on different selection criteria (yield, calmness, or
index).

# Strength and honey yield

In this section we change simulation to two traits where the phenotype
realisation of the first trait affects the phenotype realisation of the second
trait. Specifically, we will assume that queen's fecundity, and hence the number
of workers, is under the genetic affect of the queen and her environment.
Furthermore, we will assume as before that colony honey yield is due to the
queen effect and workers effect. Since the value of the workers effect depends
on then number of workers, we obtain correlation between fecundity and honey
yield, even if these traits would be uncorrelated on the queen level. We
emphasise that this is just an example and the biology of these traits might be
different.

We follow the same logic as before and simulate three traits that will
contribute to two colony traits, queen's fecundity, that is colony strength, and
honey yield. We assume that fecundity is only due to the queen (and not the
workers), hence we simulate only the queen effect for this trait. For honey
yield we again assume that both the queen and workers contribute to the colony
value. For speed of simulation we only simulate 100 workers per colony on
average and split honey yield mean between the queen and workers. We measure
fecundity with the number of workers, which is a count variable and for such
variables Poisson distribution is a good model. This distribution has just one
parameter (lambda) that represents both the mean and variance of the variable.
To this end we set phenotypic variance to 100 and split it into 25 for genetic
and 65 for environmental variance. As before we warn that these are just
exemplary values to demonstrate the code functionality and do not necessarily
reflect published values!

```{r SimParamBee_3}
# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

# Quantitative genetic parameters
# - the first trait has only the queen effect
# - the second trait has both the queen and workers effects
nWorkers <- 100
mean <- c(nWorkers, 10, 10 / nWorkers)
varA <- c(25, 1, 1 / nWorkers)
corA <- matrix(data = c(1.0,  0.0,  0.0,
                        0.0,  1.0, -0.5, 
                        0.0, -0.5,  1.0), nrow = 3, byrow = TRUE)
SP$addTraitA(nQtlPerChr = 100, mean = mean, var = varA, corA = corA,
             name = c("fecundityQueenTrait", "yieldQueenTrait", "yieldWorkersTrait"))

varE <- c(75, 3, 3 / nWorkers)
corE <- matrix(data = c(1.0, 0.0, 0.0,
                        0.0, 1.0, 0.3,
                        0.0, 0.3, 1.0), nrow = 3, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)
```

We continue by creating an apiary with 10 colonies.

```{r base_pop_and_colony_2}
basePop <- createVirginQueens(founderGenomes)
drones <- createDrones(x = basePop[1:5], nInd = 100)
apiary <- createMultiColony(basePop[6:20])
droneGroups <- pullDroneGroupsFromDCA(drones, nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
```

Let's explore queen's genetic and phenotypic values for fecundity and honey
yield. The below printouts show quite some variation in fecundity between queens
at the genetic, but particularly phenotypic level. This is a small example, so
we should not put too much into correlations between these three variables.
However, if you restart this simulation many times, you will notice zero
correlation on average between `fecundityQueenTrait` and the other two traits
and negative correlation on average between `yieldQueenTrait` and
`yieldWorkersTrait.` Just like we defined in the global simulation parameters.

```{r queen_values, echo = FALSE, fig.height = 5, fig.width = 6}
getQueenGv(apiary, collapse = TRUE)
queenPheno <- getQueenPheno(apiary, collapse = TRUE) |> as.data.frame()
cor(queenPheno)
plot(queenPheno)
```

We next build-up colonies in the apiary. But instead of building them all up to
the same fixed number of workers, we build them up according to queen's
fecundity. For that we use the sampling function `nWorkersColonyPhenotype()`,
that samples the number of workers based on phenotypes of colony members, in our
case `fecundityQueenTrait` in queens. Correspondingly, each colony will have a
different number of workers. Read more about this function in it's help page.

```{r colony_strength}
apiary <- buildUp(apiary, nWorkers = nWorkersColonyPhenotype,
                  queenTrait = "fecundityQueenTrait")
cbind(nWorkers = nWorkers(apiary), queenPheno)
help(nWorkersColonyPhenotype)
```

To compute the colony value for honey yield, we again employ the
`calcColonyPheno()` function. Correlating the queen and colony values we will
now see a positive correlation because our individual to colony mapping function
sums workers effect across all workers and the more workers there are the larger
the sum.

```{r colony_pheno_3, echo = FALSE, fig.height = 5.5, fig.width = 6.5}
colonyValuesPheno <- calcColonyPheno(apiary,
                                     queenTrait = "yieldQueenTrait",
                                     workersTrait = "yieldWorkersTrait")
pheno <- cbind(nWorkers = nWorkers(apiary), queenPheno, yield = colonyValuesPheno)
cor(pheno)
plot(pheno)
```