Lovell — Aug 7, 2014, 1:35 PM
#H2 of variance
rm(list=ls())
setwd("~/Desktop")
varin<-read.csv("si_2009_phenos_allVC.csv")
library(lme4)
Warning: package 'lme4' was built under R version 3.1.1
Loading required package: Matrix
Loading required package: Rcpp
head(varin)
Site Genotype N.Rows CV.NFruits. Variance.NFruits.
1 Italy 3 16 55.64 50.83
2 Italy 6 17 41.45 90.88
3 Italy 7 17 61.77 95.53
4 Italy 8 16 67.10 198.53
5 Italy 11 18 46.51 34.71
6 Italy 13 17 49.63 53.69
# there is one measure of variance / genotype / site. 398 genotypes
mod<-lmer(CV.NFruits. ~ Site + (1|Genotype), data=varin)
summary(mod)
Linear mixed model fit by REML ['lmerMod']
Formula: CV.NFruits. ~ Site + (1 | Genotype)
Data: varin
REML criterion at convergence: 6502
Scaled residuals:
Min 1Q Median 3Q Max
-2.510 -0.642 -0.108 0.537 4.397
Random effects:
Groups Name Variance Std.Dev.
Genotype (Intercept) 16.7 4.08
Residual 191.7 13.84
Number of obs: 796, groups: Genotype, 398
Fixed effects:
Estimate Std. Error t value
(Intercept) 48.133 0.724 66.5
SiteSweden 0.342 0.981 0.3
Correlation of Fixed Effects:
(Intr)
SiteSweden -0.678
Var_Random_effect <- as.numeric(VarCorr(mod)) #we extract the variance of the random effect (which capture individual variation not captured by fixed effects)
Var_Residual <- attr(VarCorr(mod), "sc")^2 #we extract the residual variance (which corresponds to within individual variation in this case where no fixed effects differ within indiv)
H2<-Var_Random_effect/(Var_Random_effect+Var_Residual)*100
print(H2)
[1] 8.007