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))[]; m2_line = unlist(VarCorr(m2))[]; m3_line = unlist(VarCorr(m3))[] # 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
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$.