How to calculate the BLUP(Best Linear Unbiased Prediction)

In statistics, best linear unbiased prediction (BLUP) is used in linear mixed models for the estimation of random effects.
wiki link for BLUP.

Below is a R function to do it by using an add on package nlme.
*data: is a data.frame formatted input data;
*model: used to account for the fixed effect;
*random: used to account for random effect;
*trait: phenotypic name for the output data

library(nlme)
BLUP <- function(data = krn, model = KRN ~ Genotype, random = ~1 | Farm, trait = "KRN") {
    lmeout1 <- lme(model, data = data, random = random)
    ped.hat1 <- lmeout1$coef$fixed
    ped.hat1[-1] <- ped.hat1[-1] + ped.hat1[1]
    names(ped.hat1)[1] = "B73xZ001"
    names(ped.hat1) <- gsub("Genotype2", "", names(ped.hat1))
    tped <- data.frame(Genotype = names(ped.hat1), trait = ped.hat1)
    names(tped)[2] <- trait
    return(tped)
}

Example on how to use the function:
Phenotype measured in this example is KRN. The only fixed effect I believe to account for the trait variation is Genotype. That means, idealy, same genotype should have the same KRN phentoype. However, some random effect we believe could also affect the phenotype, like which field (Farm in the table) or which row (row in the table) they been planted.

ex <- read.csv("~/Documents/Test/Off-PVP_example.csv")
head(ex)

##    Row       Barcode KRN Note.x Pedigree Genotype  Note.y    Farm
## 1 1549 11-1549-22 OP  16          Ac4192      792 Off-PVP Curtiss
## 2 1549 11-1549-21 OP  16          Ac4192      792 Off-PVP Curtiss
## 3 1549 11-1549-23 OP  14          Ac4192      792 Off-PVP Curtiss
## 4 1550 11-1550-23 OP  16          Ac4226    LH195 Off-PVP Curtiss
## 5 1550 11-1550-22 OP  16          Ac4226    LH195 Off-PVP Curtiss
## 6 1550 11-1550-21 OP  16          Ac4226    LH195 Off-PVP Curtiss

krn <- BLUP(data = ex, model = KRN ~ Genotype, random = ~1 | Farm/Row, trait = "kernelrow")

## Warning: variable 'Genotype' converted to a factor

hist(krn$kernelrow, breaks = 30, col = "red")