gaynorr / AlphaSimR

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

Question about pop@gxe #41

Closed marcio-resende closed 2 years ago

marcio-resende commented 2 years ago

I don't think this is necessarily and issue, but I am trying to understand the gxe values saved within pop@gxe when a trait is simulated with addTraitADG. See an example below. The values stored in F1@gxe are in a different scale compared to the values of F1@pheno and F1@gv. I tried different values of the environmental p-value (p) when setting the phenotype, but that leads to small changes.

Any thoughts?

#########----------------
require(AlphaSimR)
#Create founder haplotypes
founderPop <- runMacs(nInd = 10, nChr = 10, species = "MAIZE", ploidy = 2L)

#Set simulation parameters
MeanG = 150
VarG =  30 

ddVar = 0.93 
ddMean = 0.3 

H2 = 0.30
VarGE = 30

SP <- SimParam$new(founderPop)

SP$addTraitADG(300,mean=MeanG,var=VarG,
               varGxE=VarGE,meanDD=ddMean,varDD=ddVar)
SP$setVarE(H2=H2)

pop = newPop(founderPop, simParam=SP)
F1 = self(pop)
F1_pop = setPheno(F1, reps = 2)

#Assessing the gv values of the 10 individuals
F1_pop@gv

#Assessing the phenotype values
F1_pop@pheno

#Assessing the GxE interaction values
F1_pop@gxe
gaynorr commented 2 years ago

It looks like the strange magnitude of the GxE effects is due to environmental variance being set to a very small number by default (1e-06) instead of it being set to zero as intended. I missed changing the defaults for two of the trait types. This shouldn't have any real effect on the simulation, but it does look strange.

The underlying model is described in the traits vignette. Below is some code following on from your above example to show what's happening behind the scenes and how the GxE effects get used.

# Create phenotypes with no residual error for two locations
# Default functionality is equivalent p=runif(1)
Loc1 = setPheno(F1, varE=0, p=0.25, onlyPheno=TRUE)
Loc2  = setPheno(F1, varE=0, p=0.75, onlyPheno=TRUE)
cor(Loc1, Loc2)

# Reconstruct Loc1 yield manually
ManLoc1 = F1@gv + F1@gxe[[1]]*qnorm(0.25, sd=sqrt(SP$traits[[1]]@envVar))
cor(Loc1, ManLoc1)

# Examine underlying function calls
# The above calculation takes place within calcPheno, which is nested in setPheno
setPheno  
AlphaSimR:::calcPheno