Question: baySeq missing data
0
k.w.lau0 wrote:

Hi everyone,

I am trying to use baySeq with missing data. I followed the example from vignette as follows:

library(baySeq)
if(require("parallel")) cl <- makeCluster(4) else cl <- NULL
data(simData)

replicates <- c("simA", "simA", "simA", "simA", "simA", "simB", "simB", "simB", "simB", "simB")
groups <- list(NDE = c(1,1,1,1,1,1,1,1,1,1), DE = c(1,1,1,1,1,2,2,2,2,2))
CD <- new("countData", data = simData, replicates = replicates, groups = groups)
libsizes(CD) <- getLibsizes(CD)

CD <- getPriors.NB(CD, samplesize = 1000, estimation = "QL", cl = cl)
CD <- getLikelihoods(CD, cl = cl, bootStraps = 3, verbose = FALSE)

The result is the same as in the vignette. Then I introduce missing value into the simData matrix as follows:

n<-dim(simData)*dim(simData)*0.001
pos<-cbind(round(runif(n, min = 0, max = 1000)),round(runif(n, min = 0, max = 9)))
for (jj in 1:n) {
simData[pos[jj,1],pos[jj,2]]<-NA
}

Then I rerun the following steps:

CD <- new("countData", data = simData, replicates = replicates, groups = groups)
libsizes(CD) <- getLibsizes(CD)

CD <- getPriors.NB(CD, samplesize = 1000, estimation = "QL", cl = cl)
CD <- getLikelihoods(CD, cl = cl, bootStraps = 3, verbose = FALSE)

Here, when n == 10, baySeq works. However when n==100, it returns the following error

Finding priors...Error in optimise(mualt, interval = c(0, max(cts[replicates == rep]/libsizes[replicates ==  :
invalid 'xmin' value
1: In getPriors.NB(CD, samplesize = 1000, estimation = "QL", cl = cl) :
no library sizes available; inferring libsizes using default settings
2: In max(cts[replicates == rep]/libsizes[replicates == rep], na.rm = TRUE) :
no non-missing arguments to max; returning -Inf