Bug with new residual covariances

[gregorgorjanc]

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

library("AlphaSimR")
quickHaplo(nInd=10000, nChr=10, segSites=1000)

# Genetic variances
genvar <- c(1,1,1)

# Genetic correlations
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)

pop_gen0 <- newPop(founderPop)
var(pop_gen0@gv)
var(pop_gen0@gv) - CovG
var(pop_gen0@pheno)
var(pop_gen0@pheno) - CovP

Expected behavior

> CovP
     [,1] [,2] [,3]
[1,]  2.0  1.0  0.1
[2,]  1.0  2.0  0.6
[3,]  0.1  0.6  2.0

> round(var(pop_gen0@pheno), 2)
       Trait1 Trait2 Trait3
Trait1   2.01   0.48   0.10
Trait2   0.48   1.50   0.28
Trait3   0.10   0.28   1.00

> round(var(pop_gen0@pheno) - CovP, 2)
       Trait1 Trait2 Trait3
Trait1   0.01  -0.52   0.00
Trait2  -0.52  -0.50  -0.32
Trait3   0.00  -0.32  -1.00

Need to check the code in SimParam()

[gregorgorjanc]

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.

[gregorgorjanc]

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

[gaynorr]

Thanks for looking into this.