[R] glm.nb error

[Marc Schwartz]

On Jun 7, 2013, at 10:36 AM, Daofeng Li <lidaof at gmail.com> wrote:

> Thank you Sarah and Marc for your fast and nice response.
> Apology for didn't include all information.
> 
> I have a input file like following:
> 
> gene1	18	15	13	13	16	9	20	24
> gene2	15	8	8	7	0	12	18	4
> gene3	10	9	8	12	9	11	12	12
> gene4	4	0	4	3	0	5	0	0
> gene5	0	1	0	0	1	5	1	0
> gene6	3	3	3	3	4	4	4	4
> gene7	0	4	0	2	2	2	0	0
> gene8	4	4	7	3	0	6	6	12
> gene9	11	6	13	10	13	7	12	9
> gene10	6	3	6	3	4	7	6	3
> 
> I am using following R code to compare the difference between the 1st 4 numbers against 2nd 4 numbers:
> 
> library(MASS)
> library(coin)
> suppressPackageStartupMessages(suppressWarnings(library(tcltk)))
> library(qvalue)
> library(plyr)
> dat = read.table("test",col.names=c("gene","b1","b2","b3","b4","c1","c2","c3","c4"))
> index=(apply(dat[,-1],1,sum)>0)
> data = dat[index,]
> group=c(1,1,1,1,0,0,0,0)
> n=nrow(data)
> result=NULL
> for (i in 1:n){
>   print(i)
>   y=as.numeric(data[i,-1])
>   if (all((y-mean(y))==0))
>     result=rbind(result,c(0,0,0,1))
>   else {
>     fit=glm.nb(y~group)
>     result=rbind(result,summary(fit)$coef[2,])
>   }
> }
> qval = qvalue(result[,4])
> fdr=0.1
> index=(qval$qvalues<fdr)
> dat.result = data[index,]
> write.table(dat.result,file="test_result",row.names=F,quote=F)
> 
> If you use this input file and code, would reproduce the same error:
> 
> Loading required package: methods
> Loading required package: survival
> Loading required package: splines
> Loading required package: mvtnorm
> Loading required package: modeltools
> Loading required package: stats4
> 
> Attaching package: ‘plyr’
> 
> The following object is masked from ‘package:modeltools’:
> 
>     empty
> 
> [1] 1
> [1] 2
> [1] 3
> [1] 4
> [1] 5
> [1] 6
> Error in while ((it <- it + 1) < limit && abs(del) > eps) { : 
>   missing value where TRUE/FALSE needed
> Calls: glm.nb -> as.vector -> theta.ml
> In addition: Warning messages:
> 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
>   iteration limit reached
> 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace >  :
>   iteration limit reached
> Execution halted
> 
> So might be the error was in 6th line, not the line I pasted before (5th line)? Sorry about that.
> 
> Thanks.
> 
> Daofeng


Your 'y' at that point in the loop is:
  
> y
[1] 3 3 3 3 4 4 4 4

and 'group' is:

> group
[1] 1 1 1 1 0 0 0 0


> glm.nb(y ~ group)
Error in while ((it <- it + 1) < limit && abs(del) > eps) { : 
  missing value where TRUE/FALSE needed


Think about it...  :-)

Regards,

Marc Schwartz