Caution about trait means for queen and worker traits and mapping them to colony trait/value

[gregorgorjanc]

Here is an example that shows we need to be careful how we set trait means - setting trait population mean only to, say, queen trait, might not be a correct approach in every case - it will depend on how we want to map individual values to colony values. So, think how various colony events will interact with the mapping function when calculating colony values.

First, what triggered this is working on the Quantitative genetics vignette:

# ... see the vignette for the complete code - the key for this demonstration is this trait mean setting (which is not what I recommend anymore!
mean <- c(20, 0)
varA <- c(1, 1 / SP$nWorkers)
# ... much code skipped (see the vignette)!

# Calculate colony phenotype value
> calcColonyPheno(colony, queenTrait = "queenTrait", workersTrait = "workersTrait")
         [,1]
[1,] 19.75156

# Calculate colony phenotype value on a reduced colony
> calcColonyPheno(removeWorkers(colony, p = 0.5), queenTrait = "queenTrait", workersTrait = "workersTrait")
        [,1]
[1,] 21.5535

What just happened here? I removed ~50% of workers and got more honey? This isn't what I want!

The calcColonyPheno() essentially takes queen effect pheno of the queen and adds the sum of workers effect values from workers. Let's see what is going on:

> # Colony phenotype value
> (a <- getQueenPheno(colony)[, "queenTrait"])
queenTrait 
  20.49738 
> (b <- sum(getWorkersPheno(colony)[, "workersTrait"]))
[1] -0.7458214
> a + b
queenTrait 
  19.75156 

So, this is OK. Now for the reduced colony:

> # Colony phenotype value from a reduced colony
> (c <- sum(getWorkersPheno(removeWorkers(colony, p = 0.5))[, "workersTrait"]))
[1] 1.277538
> a + c
queenTrait 
  21.77492 

OK, the sums are correct, but I think this summing is not what we want - bigger colony should give more honey!

If I now change the trait means and variances from

mean <- c(20, 0)
varA <- c(1, 1 / SP$nWorkers)

to

mean <- c(10, 10 / SP$nWorkers)
varA <- c(1, 1 / SP$nWorkers)

we get this correct behaviour

> # Colony phenotype value
> calcColonyPheno(colony, queenTrait = "queenTrait", workersTrait = "workersTrait")
         [,1]
[1,] 29.19563

> # Colony phenotype value from a reduced colony
> calcColonyPheno(removeWorkers(colony, p = 0.5), queenTrait = "queenTrait", workersTrait = "workersTrait")
         [,1]
[1,] 20.54228

Here is the breakdown as before:

> # Colony phenotype value
> (a <- getQueenPheno(colony)[, "queenTrait"])
queenTrait 
  10.31003 
> (b <- sum(getWorkersPheno(colony)[, "workersTrait"]))
[1] 18.8856
> a + b
queenTrait 
  29.19563 

> # Colony phenotype value from a reduced colony
> (c <- sum(getWorkersPheno(removeWorkers(colony, p = 0.5))[, "workersTrait"]))
[1] 8.481791
> a + c
queenTrait 
  18.79182