R/gBLUP.eps.R

In YaoZhou89/GST: Genetic Selection Tools

`gBLUP.eps`<-function(X=NULL,K=NULL,mm=NULL,y=NULL){
  # Objects: prediciton using BLUP
  
  # Input: 
  #     K: list, Kinship with inference and reference
  #     mm: list, generated by calcAIvaf
  #     y: vector, inference individuals values are NA
  #     X: n by p, p is the number of covariates
  
  # Output:
  #     BLUP: blup values for all individuals; use the same function as GAPIT
  #     phe: predicted phenotype for all individuals
  
  ## evaluate the y value
  
  index.inf= which(is.na(y))
  index.ref = which(!is.na(y))
  n.inf = length(index.inf)
  n.ref = sum(!is.na(y))
  n = length(y)
  cat("The individual number of Inference population is:",length(index.inf),sep=" ")
  n.random = length(K) # Kinship number
  if(n.random==0) stop("This no random effect! Kinship Matrix is needed.")
  if(is.null(X)) X = matrix(rep(1,n),n)
  Ind = matrix(diag(n.ref),n.ref,n.ref)
  BLUP.left = matrix(0,n,n.ref)
  Var = mm$Var
  for(i in 1:n.random){
    BLUP.left = BLUP.left + Var[i]* K[[i]][,index.ref] %*% solve(K[[i]][index.ref,index.ref]*Var[1]+Ind*Var[length(Var)])
  }
  BLUP.right =as.matrix(y[index.ref]) - X[index.ref,]%*%mm$Bhat
  BLUP = BLUP.left %*% BLUP.right
  phe = X %*% mm$Bhat + BLUP
  LL = mm$llike
  return(list(BLUP = BLUP, phe = phe, Var = Var, LL = LL))
}