# Additive genetic variance components from LMER in R

Biology Asked by rg255 on July 9, 2021

I’ve set up some dummy data in R which makes 40 genetically related lines, they are all siblings within a line so are genetically related by a factor of ½ thus additive genetic variance should be twice the variance explained by line. For the lines there are 200 individuals being measured, for three characters/traits. The first trait has low phenotypic variance, the second has high environmental variance, and the third has high genetic variance.

rm(list=ls())
re = 200 # replicate individuals per line
li = 40  # lines

# setup
set.seed(5)
data = data.frame(rep(1:li, each = re))
colnames(data)="Line"
library(nlme)
library(lme4)
par(mfrow=c(1,3))

# trait 1: little variance (Va or Ve)
data$Trait = rnorm(li*re,10,1) boxplot(data$Trait~data$Line, ylim=c(0,20), main = expression("Low V"[A]*"& Low V"[E])) var1 = var(data$Trait); mean1 = mean(data$Trait); m1 = lmer(data$Trait ~ (1|data$Line)) # trait 2: high enivronmental variance, little Va data$Trait = rnorm(li*re,10,1); data$Trait = data$Trait + rnorm(re*li,0,3)
boxplot(data$Trait~data$Line, ylim=c(0,20),main = expression("Low V"[A]*"& High V"[E]))
var2 = var(data$Trait); mean2 = mean(data$Trait); m2 = lmer(data$Trait ~ (1|data$Line))

# trait 3: high additive genetic variance, little Ve
data$Trait = rnorm(li*re,10,1); data$Trait = data$Trait+rep(rnorm(li,0,3),each=re) boxplot(data$Trait~data$Line, ylim=c(0,20), main = expression("High V"[A]*"& Low V"[E])) var3 = var(data$Trait); mean3 = mean(data$Trait); m3 = lmer(data$Trait ~ (1|data$Line))  Plots of the three traits, Then to estimate additive variance ($V_A$) I extract the line variance ($V_L$) from the lmer model and double it, # line variances (variance in additive effect of each haploid genome) m1_line = unlist(VarCorr(m1))[[1]]; m2_line = unlist(VarCorr(m2))[[1]]; m3_line = unlist(VarCorr(m3))[[1]] # additive variance (double the line variance because it is a hemiclone) m1_add = 2*m1_line; m2_add = 2*m2_line; m3_add = 2*m3_line  Residual variances should be (assuming perfect experimental design, no measurement error etc.) the estimate of environmental variance ($V_E$) and phenotypic variance ($V_P$) should be the sum of$V_L$and$V_E$, # residual variance m1_res = attr(VarCorr(m1), "sc")^2 m2_res = attr(VarCorr(m2), "sc")^2 m3_res = attr(VarCorr(m3), "sc")^2 # phenotypic variance m1_phe = m1_line + m1_res m2_phe = m2_line + m2_res m3_phe = m3_line + m3_res  Heritability is additive variance divided by the phenotypic variance$h^2 = V_A / V_P$But I think it is correct in this case to use line variance rather than additive variance (if someone could explain in an answer that would be useful), so I’ve done, # heritability (line/ (line+ residual)) m1_h2 = m1_line/ m1_phe m2_h2 = m2_line/ m2_phe m3_h2 = m3_line/ m3_phe m1_h2; m2_h2; m3_h2  My question(s): Is it appropriate to use the lmer function in R to extract variance components in this manner? Have I calculated$V_A$,$V_E$,$V_P$,$h^2$correctly? I think$V_A$and$V_E$are correct,$V_P$could perhaps be the sum of$V_A$and$V_E$rather than$V_L$, and subsequently$h^2$may be$V_A/V_P$. ## One Answer From what I understood of your code and what you are asking I am guessing that you do the following: Generating a virtual set of 40 individuals (lines) of which you have 200 measurements (repetitions). You say that they are full siblings, so they share both parents. Then you use the lmer function (which I am not familiar with) to give you the total variance, within-group variance and between-group variance, ($V_{T}, V_W, V_B$respectively). What you call$V_L$would be the$V_W$. We know that$V_P=V_A+V_D+V_E$In full siblings we also know that$V(A,A')=frac{1}2V_A+frac{1}4V_D$where$A$is the additive value and$D$is the dominance value because they share a quarter of their gene combinations. If you are treating with full siblings the$V_W$should be then the$V_E$and the$V_B$would be the$frac{1}2V_A+frac{1}4V_D$as you are seeing the variance between siblings. So no, you are not calculating correctly the$V_A$because you are missing your$V_D\$ on your counting. The code itself looks fine for what I understood, but maybe you should take that to corssvalidated with the R tag

Answered by Athe on July 9, 2021