ODP instability

[Arthfael]

I have recently added both the LRT and ODP tests to the statistics which I routinely run for my proteomics workflow, however whilst the LRT is not problematic, the OPD often fails with:

Error in kl(temp.center.fitFull, temp.fitFull, center.var, sigma2, n = n) : 
  NA/NaN/Inf in foreign function call (arg 1)

which seems to be related to how it calls the kl function from package dependency splines.

Adding lambda = 0 to the arguments of the edge::odp() call (based on this link ) does help sometimes, but not always.

Disclaimer: I do filter my data before feeding it to the functions (keeping currently only rows with 2+ valid values per group. But I tend to deal with few replicates (3-5 usually), this might be the problem.
I am happy to provide an example dataset if that helps. Normally I would try to investigate deeper into the function's code, but I am unfamiliar with how to access code for standardGeneric functions.

[ajbass]

I am happy to take a look!

Can you send sample code / data so I can investigate what is going on?

[Arthfael]

Sorry that it took so long. I have attached a zip containing the de_obj and odp_fit files as .RDS files. Here's how I run ODP currently:

require(edge)
#
Coeff <- c("Genotype", "Treatment")
expMap <- data.frame(Reference = rep(c(rep(FALSE, 3), rep(TRUE, 3)), 2),
                     Genotype = c(rep("Columbia", 6), rep("cre1.12", 6)),
                     Treatment = rep(c(rep("2uM_tZ", 3), rep("DMSO", 3)), 2),
                     Group = c(rep("Exp1___Columbia___2uM_tZ", 3),
                               rep("Exp1___Columbia___DMSO", 3),
                               rep("Exp1___cre1.12___2uM_tZ", 3),
                               rep("Exp1___cre1.12___DMSO", 3)),
                     Sample_name = c(paste0("Exp1___Columbia___2uM_tZ___", 1:3),
                                     paste0("Exp1___Columbia___DMSO___", 1:3),
                                     paste0("Exp1___cre1.12___2uM_tZ___", 1:3),
                                     paste0("Exp1___cre1.12___DMSO___", 1:3)))
#
rg <- 1:6
expMap <- expMap[rg,]
# For now I am doing two-group comparisons,
# although I do understand that ODP is meant to be for multiple hypotheses testing (like an F-test) so this will go eventually.
#
kol <- paste0("log10(Expr.) - ", expMap$Sample_name[rg])
Grps <- unique(expMap$Group)
tst <- sapply(Grps, function(grp) {
  smpls <- expMap
  kl <- paste0("log10(Expr.) - ", expMap$Sample_name[which(expMap$Group == grp)])
  apply(myData[, kl], 1, function(x) {
    length(which((!is.na(x))&(is.finite(x))))
  })
})
w <- which(apply(tst, 1, function(x) {
  min(x) >= 2
}))
tmpDat <- myData[w, kol]
tmpDat2 <- as.matrix(tmpDat)
tmpDat2[which(!is.finite(tmpDat2), arr.ind = TRUE)] <- NA
tmpDat2[which(is.nan(tmpDat2), arr.ind = TRUE)] <- NA
colnames(tmpDat2) <- expMap$Sample_name
#
bioVar <- magrittr::set_rownames(as.matrix(expMap[, Coeff, drop = FALSE]), NULL)
bioVar <- bioVar[, which(apply(bioVar, 2, function(x) { length(unique(x)) }) > 1), drop = FALSE]
stopifnot(ncol(bioVar) > 0)
bsCall <- "de_obj <- edge::build_study(data = tmpDat2, bio.var = bioVar)"
#cat(bsCall)
eval(parse(text = bsCall), envir = .GlobalEnv)
#
# optimal discovery procedure
odpTmp <- try({
  
  # Fit
  odp_fit <- edge::fit_models(de_obj, stat.type = "odp")
  
  mySeed <- 1234567
  de_odp <- try(edge::odp(de_obj, odp_fit, bs.its = 50, n.mods = 50, verbose = TRUE, seed = mySeed), silent = TRUE)
  
  if ("try-error"%in% class(de_odp)) {
    de_odp <- edge::odp(de_obj, odp_fit, bs.its = 50, n.mods = 50, verbose = TRUE, seed = mySeed, lambda = 0)
  }
  
  # lambda = 0 is from troubleshooting common error:
  # "Error in smooth.spline(lambda, pi0, df = smooth.df) : missing or infinite values in inputs are not allowed",
  # cf. https://support.bioconductor.org/p/105623/
  
}, silent = TRUE)

# Code to save RDS files
fl <- file(paste0(wd, "/de_obj.RDS"), open = "wb")
saveRDS(de_obj, fl, compress = FALSE)
close.connection(fl)
fl <- file(paste0(wd, "/odp_fit.RDS"), open = "wb")
saveRDS(odp_fit, fl, compress = FALSE)
close.connection(fl)

RDS.zip

[Arthfael]

This does seem to just be a case of NAs in the data not being filtered properly. I had assumed that at least 2 non-missing values per condition was acceptable - it works for edge's LRT implementation. However, if I apply stricter filtering (no missing values at all), ODP does work. Is there somewhere a guide to how filtering should be handled for each algorithm?