Closed gregorgorjanc closed 1 year ago
Just to add, that the SP$varE
is stored well
> CovE
[,1] [,2] [,3]
[1,] 1.0 0.5 0.0
[2,] 0.5 1.0 0.3
[3,] 0.0 0.3 1.0
> SP$varE
[,1] [,2] [,3]
[1,] 1.0 0.5 0.0
[2,] 0.5 1.0 0.3
[3,] 0.0 0.3 1.0
> SP$varE - CovE
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
[3,] 0 0 0
So, the error must be in pheno generation.
The above pull request should fix this issue - I get this behaviour now
founderPop <- quickHaplo(nInd=10000, nChr=10, segSites=1000)
genvar <- c(1,1,1)
gencor <- matrix(c(1.0,0.5,0.1,
0.5,1.0,0.3,
0.1,0.3,1.0), nrow = 3, ncol = 3)
CovG <- diag(sqrt(genvar)) %*% gencor %*% t(diag(sqrt(genvar)))
resvar <- c(1,1,1) # h2 = 0.5 for all three traits
rescor <- matrix(c(1.0,0.5,0.0,
0.5,1.0,0.3,
0.0,0.3,1.0), nrow = 3, ncol = 3)
CovE <- diag(sqrt(resvar)) %*% rescor %*% t(diag(sqrt(resvar)))
CovP <- CovG + CovE
SP = SimParam$new(founderPop)
SP$addTraitA(nQtlPerChr = 500, var = genvar, mean = rep(0, length(genvar)), corA = gencor)
SP$setVarE(varE = resvar, corE = rescor)
CovE
SP$varE
SP$varE - CovE
pop_gen0 <- newPop(founderPop)
round(var(pop_gen0@gv), digits = 2)
round(var(pop_gen0@gv) - CovG, digits = 2)
round(var(pop_gen0@pheno), digits = 2)
round(var(pop_gen0@pheno) - CovP, digits = 2)
> round(var(pop_gen0@pheno), digits = 2)
Trait1 Trait2 Trait3
Trait1 1.99 0.96 0.12
Trait2 0.96 1.97 0.60
Trait3 0.12 0.60 1.97
> round(var(pop_gen0@pheno) - CovP, digits = 2)
Trait1 Trait2 Trait3
Trait1 -0.01 -0.04 0.02
Trait2 -0.04 -0.03 0.00
Trait3 0.02 0.00 -0.03
Thanks for looking into this.
Describe the bug There seems something wrong with how new residual covariances are used - in AlphaSimR_1.3.2
Steps To Reproduce Use this script
Expected behavior
Need to check the code in
SimParam()