gaynorr / AlphaSimR

R package for breeding program simulations
https://gaynorr.github.io/AlphaSimR/
Other
40 stars 17 forks source link

switchTrait() introduces NA values in VarE #188

Closed irenedecarlos closed 1 month ago

irenedecarlos commented 3 months ago

Describe the bug I am adding 4 quantitative traits to the SimParam object (SP), adding a genetic variance and genetic correlation matrix as well as the environmental variance and the environmental correlation matrix. I want to switch traits after creating the first population object and then recreate the first population with the switched traits. However, when I try the function switchTrait(), I get Na values introduced into the envorionmental variance varE that I see when I inspect the SP object.

Steps To Reproduce

library(AlphaSimR) founderGenomes<- quickHaplo(800,4,segSites = 1000) SP <- SimParam$new(founderGenomes) mean <- c(0,0,0,0) varA <- c(0.25,0.25,0.1,0.1) corA <- matrix(data = c( 1.0, 0.9, 0.0, 0.0, 0.9, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.0, 0.5, 1.0), nrow = 4, byrow = TRUE)

SP$addTraitA(nQtlPerChr = 100, mean = mean, var = varA, corA = corA, name = c("trait1", "trait2", "trait3", "trait4"))

varE <- c(0.75,0.75,0.9,0.9)

corE <- matrix(data = c( 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0 ), nrow = 4, byrow = TRUE)

SP$setVarE(varE = varE, corE = corE)

basePop = newPop(founderGenomes) SP$resetPed() traits = SP$traits SP$switchTrait(traitPos= 4,lociMap= traits[[3]]) SP$switchTrait(traitPos= 3,lociMap= traits[[4]])

SP

basePop = newPop(founderGenomes)

Expected behavior I expected to be able to switch the traits and not get NA values introduced into the environmental variance.

Additional context OUTPUT of SP after switching traits .traits: list .varA: 0.25 0.25 0.1 0.1 .varE: 0.75 0 0 0 0 0.75 0 0 0 0 NA 0 0 0 0 NA .varG: 0.25 0.25 0.1 0.1

gregorgorjanc commented 3 months ago

@irenedecarlos 16 comes from 4x4 matrix, but the real issue is that switchTrait() induced NAs into the varE. Please adapt the above example with using switchTrait() and show the ouput.

gaynorr commented 3 months ago

There is definitely a problem with the way it is currently implemented, so I've made a change. However, I might make more changes later.

There's a lot going on here with how AlphaSimR is trying to handle the various things the user may do. First, I should mention that switchTrait comes with a varE argument that is supposed to be used to pass along the residual variance for the new trait. The default value for this variance is NA, this is consistent with what is create when using a new trait function.

The problem here is that .varE has two possible states. If no corE has been specified, .varE is a vector of length nTraits. In this case, switchTrait is working as intended. If corE has been specified, .varE is a matrix and the behavior of switchTrait can be strange. I've changed the behavior to drop the correlations and revert .varE to a vector. This will also trigger a warning to the user.

These changes should prevent the code crashing with errors, but I'm not doing a very good job guarding against the user doing something unintended.