######################## # variance components and BLUPS ####################### #loading phenotypes pheno <- readRDS("pheno2") head(pheno) dim(pheno) # lets create a new factor, genotype x nitrogen (GxE) 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) ?remlf90 args(remlf90) model0 <- remlf90(fixed = MPSadj ~ 1, random = ~ rep + gid, method = "em", data = pheno) model1 <- remlf90(fixed = MPSadj ~ 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)) # based on replicates (hg.env <- var.comp[2] / (var.comp[2] + var.comp[3]/(a) + var.comp[4]/(a*r))) # based on enviroments # test for random effects - LRT model1 <- remlf90(fixed = MPSadj ~ N, random = ~ rep + gid + GN, method = "em", data = pheno) model1.red1 <- remlf90(fixed = MPSadj ~ N, random = ~ rep + gid, method = "em", data = pheno) model1.red2 <- remlf90(fixed = MPSadj ~ N, random = ~ rep, method = "em", data = pheno) (LRT.GN <- -(model1$fit$`-2logL` - model1.red1$fit$`-2logL`)) dchisq(LRT.GN, 1) #p.value (LRT.gid <- -(model1.red1$fit$`-2logL` - model1.red2$fit$`-2logL`)) dchisq(LRT.gid, 1) #p.value # test for fixed effects - Wald output <- data.frame() # function to estimate by Wald t-test wald.breedr <- function(model, alpha = 0.05) { fctrs <- names(model$fixed) for (i in 1:length(fctrs)) { dataset <- na.omit(model$fixed[[fctrs[i]]][[1]]) output <- rbind(output, data.frame( Bj = diff(range(dataset[,1])), epBj = dataset[which(abs(dataset[,1]) == max(abs(dataset[,1]))), 2], Zj = diff(range(dataset[,1])) / dataset[which(abs(dataset[,1]) == max(abs(dataset[,1]))), 2], df = length(dataset[,1]) - 1, t = qt((1 - alpha / 2), length(dataset[,1]) - 1), p.value = pt(diff(range(dataset[,1])) / dataset[which(abs(dataset[,1]) == max(abs(dataset[,1]))), 2], length(dataset[,1]) - 1, lower.tail = F), alpha = alpha, Wald = ifelse((diff(range(dataset[,1])) / dataset[which(abs(dataset[,1]) == max(abs(dataset[,1]))), 2]) > qt((1 - alpha / 2), length(dataset[,1]) - 1), "*", "ns") )) } rownames(output) <- fctrs return(output) } (wald <- wald.breedr(model1, alpha = 0.01)) ############################## # mixed model equations - BLUP ############################## # fixed effects model1$fixed fixef(model1) #random effects model1$ranef # BLUPS blup <- model1$ranef$gid[[1]] head(blup) summary(model1) # reliability (r2) = 1 - PEV / (Vg + Vg*Fii) # near to heritability # If F == 0, then r = 1 - PEV/Vg # If F == 1, then r = 1 - PEV/2Vg r2 <- 1 - (blup[,2]) ^ 2 / (diag(diag(length(blup[,1]))) * var.comp[2]) (mean(r2)) # Accuracy = sqrt(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.5), width = 0.10) + labs(x = "gid", y = "BLUP") # correlation between phenotypes and BLUPS means <- tapply(pheno$MPSadj, 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(-1,0, round(cor(blups2$mean, blups2$BLUP.value),2), col = "red") ############################## # selecting the best genotypes # first, define a qualile, in this case 90% (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))