Additive genetic variance components from LMER in R

Question

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$.

Answer 1

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