######################## # variance componentes ####################### #loading phenotypes pheno <- readRDS("pheno2") head(pheno) dim(pheno) pheno$GN <- as.factor(paste(pheno$gid, pheno$N, sep = "")) head(pheno) #source("http://famuvie.github.io/breedR/src/setup_repo.R") #install.packages('breedR') library(breedR) args(remlf90) model0 <- remlf90(fixed = MPSadj2 ~ 1, random = ~ rep + gid, method = "em", data = pheno) model1 <- remlf90(fixed = MPSadj2 ~ N, random = ~ rep + gid + GN, method = "em", data = pheno) # components of variance (var.comp0 <- model0$var) (var.comp <- model1$var) #heritability without pedigree or kinship information (I = identity relationship) (r <- length(unique(pheno$rep))) (a <- length(unique(pheno$N))) hg.rep <- var.comp0[2] / (var.comp0[2] + var.comp0[3]/r) hg.rep # based on replicates hg.env <- var.comp[2] / (var.comp[2] + var.comp[3]/(a) + var.comp[4]/(a*r)) hg.env # based on enviroments ############################## # mixed model equations - BLUP ############################## # fixed effects model1$fixed #random effects model1$ranef # BLUPS blup <- model1$ranef$gid[[1]] blup # reliability (r2) and accuracy r2 <- 1 - (blup[,2])^2/(diag(diag(length(blup[,1])))*var.comp[2]) (mean(r2)) (Accuracy.blup <- mean(sqrt(r2))) (Accuracy.pheno <- sqrt(hg.env)) # confidence interval for BLUPS DMS <- blup[,2]*1.96 blups2 <- data.frame("BLUP" = blup, "DMS" = DMS) blups2$gid <- rownames(blups2) library(ggplot2) limits <- aes(ymax = blups2$BLUP.value + blups2$DMS, ymin = blups2$BLUP.value - blups2$DMS) p <- ggplot(data = blups2, aes(x = factor(gid), y = BLUP.value)) p + geom_jitter(stat = "identity", colour = "red") + geom_errorbar(limits, position = position_dodge(0.7), width = 0.15) + labs(x = "gid", y = "BLUP") # correlation between phenotypes and BLUPS means <- tapply(pheno$MPS, pheno$gid, mean) means <- means[rownames(blups2)] all(blups2$gid == names(means)) blups2$mean <- means plot(blups2$mean, blups2$BLUP.value, main = "mean x BLUP", xlab = "mean", ylab = "BLUP") abline(lm(blups2$BLUP.value ~ blups2$mean), col = "red") text(0.2,0, round(cor(blups2$mean, blups2$BLUP.value),2), col = "gray60") # selecting the best genotypes # first, define a qualile (trshld <- quantile(blup[,1], probs = 0.90)) # Then identify the best ones (selected <- blups2[blups2$BLUP.value >= trshld,]) # and estimate the response to selection (RS <- mean(selected$BLUP.value))