r 计算遗传力

Posted 2021-05-24

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
}