At a minimum, make script liability-threshold model available in AlphaSimR discussions thread. Better yet, role this functionality into AlphaSimR directly. Script is below.
# Based on the liability-threshold model
library(AlphaSimR)
prevelance = 0.95
varG = 1
H2 = 0.2
# Calculate varE
varE = varG/H2 - varG
# Calculate truncation point
t = qnorm(prevelance, sd=sqrt(varG+varE),lower.tail=F)
phenoTL = function(Y,whichTrait=1,trunc=t){
return(Y[,whichTrait,drop=FALSE]>trunc)
}
founderPop = quickHaplo(1000, 10, 1000)
SP = SimParam$
new(founderPop)$
addTraitA(1000,var=varG)$
setVarE(H2=H2)
pop = newPop(founderPop)
meanO = meanL = numeric(20)
for(i in 1:20){
pop = selectCross(pop,100,trait=phenoTL,nCrosses=1000)
meanO[i] = mean(phenoTL(pop@pheno))
meanL[i] = meanG(pop)
}
meanO
meanL
At a minimum, make script liability-threshold model available in AlphaSimR discussions thread. Better yet, role this functionality into AlphaSimR directly. Script is below.