danielecook
2/16/2017 - 5:45 PM

Calculate Heritability

Calculate Heritability

H2.fun <- function(x){
  pdata <- x
  pdata = split(pdata$pheno, pdata$isotype)
  pdata.notNAcnt = sapply(pdata, function(x){sum(!is.na(x))})
  pdata[pdata.notNAcnt<2]=NULL
  pdata.melted = melt(pdata)
  names(pdata.melted)=c('pheno', 'isotype')
  pdata.melted$strain=as.factor(pdata.melted$isotype)
  reffMod = lme4::lmer(pheno ~ 1 + (1|strain), data=pdata.melted)
  Var_Random_effect <- as.numeric(VarCorr(reffMod))
  Var_Residual <- attr(VarCorr(reffMod), "sc")^2
  H2 <- Var_Random_effect/(Var_Random_effect+Var_Residual)
  H2
}