Error in rinvGauss(n = ETA[[j]]$p, nu = nu, lambda = ETA[[j]]$lambda2) : nu must be positive

[isaamael]

Dear Author
I am trying to run the following script for BGLR analysis, all the required parameters for the script have been passed in correctly, however as the number of iterations increases, the software starts to report an error that I can't handle：Error in rinvGauss(n = ETA[[j]]$p, nu = nu, lambda = ETA[[j]]$lambda2)
nu must be positive:
Script Content:

library(BGLR)
library(optparse)

option_list = list(
  make_option("--pheno", type="character", default=NULL, 
              help="Phenotype file name", metavar="phenofilename"),
  make_option("--geno", type="character", default=NULL, 
              help="Genotype file name, it is a 012 type matrix which created by vcftools", metavar="genofilename"),
  make_option("--out", type="character", default=NULL, 
              help="Output model result name", metavar="outfilename"),
  make_option("--cv", type="integer", default=10, 
              help="Number of cross-validation folds", metavar="cvnum"),
  make_option("--fold", type="integer", default=5, 
              help="Number of individuals to sample in each fold", metavar="foldnum"),
  make_option("--iter", type="integer", default=16000,
              help="Number of interation times", metavar="interation"),
  make_option("--burn", type="integer", default=5000,
              help="Number of burn times", metavar="burn"),
  make_option("--outputdir", type="character", default=NULL,
              help="providing a path to output the result and idx", metavar="outputdir")
)


opt_parser = OptionParser(option_list=option_list)
opt = parse_args(opt_parser)


set.seed(42)
options(digits = 4)

SNP=read.table(opt$geno)
META=read.table(opt$pheno,header = TRUE)

options(digits =4)

rrBLUPres = matrix(ncol = 5, nrow = nrow(META) * ncol(META) * opt$cv)
colnames(rrBLUPres) = c("trait", "heritability", "accuracy", "r2", "set")
idxlist = matrix(nrow = opt$cv, ncol = nrow(META) * 1 / opt$fold)
rownames(idxlist) = c(seq(1,opt$cv))

nIter=opt$iter
burnIn=opt$burn
set.seed(42)
n=nrow(META)
nTST=round(as.numeric(n/opt$fold),0)

X1=scale(SNP)
X1[is.na(X1)] <- 0
for (i in seq(1,opt$cv)){
  testidx = sample(as.numeric(rownames(META)),as.numeric(nrow(META) * 1 / opt$fold))
  sort_testidx = sort(testidx)
  idxlist[i,] = sort_testidx
  pheno = META
  pheno[sort_testidx,] = NA

  setwd(opt$outputdir)

  for (p in seq(2,ncol(META))){
    varfile_dir <- paste0("varfile/BLcv", i, "_META", p)
    dir.create(varfile_dir, recursive = TRUE)
    setwd(varfile_dir)
    y=pheno[,p]
    y_true = META[,c(1,p)]
    idx_pre = intersect(which(is.na(y_true[,2])==F),sort_testidx)
    idx_r2 = which(is.na(y_true[,2])==F)
    
    fmBL=BGLR(y=y,ETA=list(list(X=X1,model='BL',saveEffects=T)),nIter=nIter,burnIn=burnIn)
    prediction=fmBL$yHat
    
    varU = scan('ETA_1_lambda.dat')
    #varU = na.omit(varU)
    varE = scan('varE.dat')
    #varE = na.omit(varE)
    h2_2 = varU/(varU+varE)
    
    accuracy = cor(prediction[idx_pre], y_true[idx_pre,2],use = "complete.obs")
    r2=1 - sum((prediction[idx_r2]-y_true[idx_r2,2])^2)/
      sum((y_true[idx_r2,2]-mean(y_true[idx_r2,2]))^2)
    heritability = mean(h2_2)
    
    rrBLUPres[(i-1)*ncol(META)+p,1] = colnames(META)[p]
    rrBLUPres[(i-1)*ncol(META)+p,2] = heritability
    rrBLUPres[(i-1)*ncol(META)+p,3] = accuracy
    rrBLUPres[(i-1)*ncol(META)+p,4] = r2
    rrBLUPres[(i-1)*ncol(META)+p,5] = i
    
    rrBLUPres = rrBLUPres[order(rrBLUPres[,1]),]
    setwd(opt$outputdir)
    unlink(varfile_dir, recursive = TRUE)
  }
}
setwd(opt$outputdir)
rrBLUPres_clean <- na.omit(rrBLUPres)
write.table(rrBLUPres_clean, file = opt$out, 
            col.names = T, row.names = F, quote = F, sep="\t")
write.table(idxlist, file = "idx.BGLR.BL.out", 
            col.names = T, row.names = T, quote = F, sep="\t")

Error:

Iter=1648 Time/Iter=0.071
varE=28.741
---------------------------------------
Iter=1649 Time/Iter=0.072
varE=28.545
---------------------------------------
Iter=1650 Time/Iter=0.093
varE=28.623
Error in rinvGauss(n = ETA[[j]]$p, nu = nu, lambda = ETA[[j]]$lambda2) : 
  nu must be positive
Warning in BGLR(y = y, ETA = list(list(X = X1, model = "BL", saveEffects = T)),  :
  tau2 was not updated  in iteration 1651 due to numeric problems with beta

---------------------------------------
Iter=1651 Time/Iter=0.086
varE=NaN
Error in if (any(nu <= 0)) stop("nu must be positive") : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In rnorm(n = nNa, sd = sdE) : NAs produced
2: In rnorm(n = 1, sd = sqrt(1/C)) : NAs produced
Warning in BGLR(y = y, ETA = list(list(X = X1, model = "BL", saveEffects = T)),  :
  tau2 was not updated  in iteration 1652 due to numeric problems with beta

Can I make any attempts to troubleshoot this issue or are you aware of the possible cause or location of the problem?
Thank you for your help, I would appreciate it!