R/bGREG.R

In gcostaneto/envirotype: GIS platforms to support envirorytping in breeding programs

# Bayesian genotype-regression with environmental covariables
bGREG = function(df, W, y,
                 Kernel.G = NULL, 
                 Kernel.E = NULL, 
                 gid,env,
                 g.model  = "FIXED",
                 w.model  = "FIXED",
                 e.model  = "FIXED", 
                 ge.model = "FIXED",
                 nIter    = 15E3,
                 burnIn   = 5E3,
                 thin     = 10,
                 R2       = .6,
                 verbose  = FALSE){
  require(BGLR);
  require(BMTME);
  W             = as.matrix(W);
  .Y            = data.frame(gid=df[,gid],env=df[,env],value=df[,y]);
  .Zg           = model.matrix(~0+gid,.Y);
  colnames(.Zg) = gsub(colnames(.Zg),pattern = "gid",replacement = "");
  .Ze           = model.matrix(~0+env,.Y);
  colnames(.Ze) = gsub(colnames(.Ze),pattern = "env",replacement = "");
  cat("#'-------------------------------------- \n")
  cat(paste0("starting: ",Sys.time(),"\n"));
  cat("#'-------------------------------------- \n")
  cat("Bayesian Model with: \n");
  cat(paste0(e.model,"   Environmental effect \n"));
  cat(paste0(w.model,"   Covariables effect \n"));
  cat(paste0(g.model,"   Genotypic  effect \n"));
  cat(paste0(ge.model,"  GxE effect \n"));
  cat("#'-------------------------------------- \n")
  #'-------------------------------------------------------------------------------------#
  # Environmental covariables (W)
  #'-------------------------------------------------------------------------------------#
  if( w.model == "FIXED"){Wn = W}
  if(!w.model == "FIXED"){
    # var-cov matrix for w
    .WtW  = crossprod(W)/nrow(W);
    .WtW  = .WtW  + diag(1E-3,nrow=nrow(.WtW))
    .LC = t(chol(.WtW))
    Wn = W %*% .LC
  }
  #'-------------------------------------------------------------------------------------#
  #' Genotype effect (G)
  #'-------------------------------------------------------------------------------------#
  if( g.model == "FIXED"){Gn = .Zg}
  if(!g.model == "FIXED"){
  if(is.null(Kernel.G)) {Kernel.G =diag(1,nrow=ncol(.Zg))}
  .LG = t(chol(Kernel.G));
  # Var-cov genetic (ZGZ')
  .ZGZt = tcrossprod(.Zg%*%.LG);
  # Var-cov gennotype x environment (ZGZ')
  .ZWZt = tcrossprod(Wn)*.ZGZt;
  Gn = .Zg %*% .LG}
  
  #'-------------------------------------------------------------------------------------#
  #' Environmental Effect (E)
  #'-------------------------------------------------------------------------------------#
  if( e.model == "FIXED"){En = .Ze}
  if(!e.model == "FIXED"){
    if(is.null(Kernel.E)) {Kernel.E = diag(1,nrow=ncol(.Ze))}
    .LE = t(chol(Kernel.E));
    # Var-cov genetic (ZEZ')
    .ZEZt = tcrossprod(.Ze%*%.LE);
    En = .Ze %*% .LE}
  
  if(w.model == "FIXED" || g.model=="FIXED"){
    ge.model = 'FIXED';
    .ZWZt = (.Zg%*%t(.Zg))%*%(W%*%t(W))
  }
  #'-------------------------------------------------------------------------------------#
  .ETA =  list(
           W = list(X =  Wn   , model = w.model),
           E = list(X =  En   , model = e.model),
           G = list(X  = Gn   , model = g.model),
          GE = list(X  = .ZWZt, model = ge.model));
  #'-------------------------------------------------------------------------------------#
  cat(paste0("Bayesian GREG with nIter= ",nIter," and burnIn ",burnIn,"\n"));
  #'-------------------------------------------------------------------------------------#
  .mixed = BGLR(y = .Y$value,ETA=.ETA,nIter = nIter,
                burnIn = burnIn,
                thin = thin,
                verbose = verbose,R2=R2);
  #'-------------------------------------------------------------------------------------#
  #' Predicted yhat
  .Y = data.frame(.Y, yhat = .mixed$yHat, sd.yhat = .mixed$SD.yHat);
  #'-------------------------------------------------------------------------------------#
  cat(paste0("ending: ",Sys.time(),"\n"));
  #'-------------------------------------------------------------------------------------#
  return(list(mixed.solve = .mixed,df = .Y));
}