Liability-Threshold Model

[gaynorr]

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