baySeq missing data
0
0
Entering edit mode
k.w.lau • 0
@kwlau-7687
Last seen 9.5 years ago
United Kingdom

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.

 

bayseq missing data • 1.5k views
ADD COMMENT

Login before adding your answer.

Traffic: 772 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6