Simulating Discrete/Categorical Phenotypes

[MattKinghorn]

I am trying to use AlphaSimR to simulate a breeding program.

I am wanting to simulate a discrete trait such as heifer pregnancy where 1 = success and 0 = failure. I tried to do this using the setPheno function and then modifying the phenotype afterwards, but I am not sure of this is the correct procedure. It feels like I might be distorting phenotypic and genotypic relationship by doing this.

basePopF = setPheno(basePopF, h2 = 0.10)

threshold = 0
basePopF@pheno[,1] = ifelse(basePopF@pheno[,1] > threshold, 1,0)

Is there an option in the setPheno function to generate traits that are not continuous, such as heifer pregnancy?

[gaynorr]

@MattKinghorn, you can pretty freely change the 'pheno' slot in populations, so your method should work just fine. I have some code from doing something similar and I just never got around to generalizing it enough and folding it into AlphaSimR. Here's the function I used and some example code. I was modifying the trait at the time of selection, but you could just as easily do it in the pheno slot like you're doing.

# 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

[MattKinghorn]

@gaynorr thank you so much for your response. What assumptions are made in the setPheno function in terms of distribution? The underlying genetic effects are simulated using SP$addtraitA and are draws from a standard normal distribution, right? Does the setPheno also assume a standard normal distribution as well, just with a different variance based on the heritability that you assign?

[gaynorr]

@MattKinghorn , yes the default sampling for additive QTL effects is standard normal distribution. However, you can also specify a gamma distribution that has a random sign (+/-) assigned to the effects (this could equally be considered a t distribution). The variance of the effects then gets scaled to the user specified variance (default is 1). This variance can either be additive or total genetic, which are equivalent if you only have additive effects.

The setPheno function really just takes the genetic values in the 'gv' slot and adds noise to them. This noise comes from a Gaussian distribution with variance equal to varE in the function argument or SP$varE if the argument is left at its default value of NULL. The h2 and H2 arguments are alternative methods for setting varE. They calculated varE needed to hit the desired heritability in the founder population. Note that the founder population doesn't change, but your current population will. Thus, you're not assigning the heritability for a specific population in setPheno.

If the trait has GxE, there is an extra step of pulling the GxE effects into the phenotype. I'd suggest checking out the Traits vignette for details, if needed.

[MattKinghorn]

@gaynorr a very informative response, that helps a lot. Thank you so much!