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)[1]*dim(simData)[2]*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
In addition: Warning messages:
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
I tried to increase the samplesize but no luck. Please helps.