glm.nb error

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

dat =
read.table("test",col.names=c("gene","b1","b2","b3","b4","c1","c2","c3","c4"))

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


On Fri, Jun 7, 2013 at 10:15 AM, Marc Schwartz <marc_schwartz at me.com> wrote:

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

Dear R Community,

I have encountered a problem while using the R function glm.nb.
The code that produce the error was following two lines:

group=c(1,1,1,1,0,0,0,0)
fit=glm.nb(y~group)

While the y contains 8 sets of number like:
gene275        0       1       0       0       1       5       1       0

Error message:

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: There were 50 or more warnings (use warnings() to see the
first 50)
Execution halted


Information of my system:

sessionInfo()

R version 3.0.1 (2013-05-16)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C                 LC_NAME=C
[9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

Does anyone happen to have some hit on how to solve this?
Appreciate for any response.

Thanks in advance,

Daofeng

There is something wrong with your actual 'y' or 'group' that is not
evident from the above info:


group <- c(1, 1, 1, 1, 0, 0, 0, 0)
y <- c(0, 1, 0, 0, 1, 5, 1, 0)

require(MASS)

Loading required package: MASS

glm.nb(y ~ group)

Call:  glm.nb(formula = y ~ group, init.theta = 1.711564307, link = log)

Coefficients:
(Intercept)        group
     0.5596      -1.9459

Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
Null Deviance:      10.23
Residual Deviance: 6.848        AIC: 25.25


Check str(y) and str(group)

You should also be sure to note in your posts when you are using a
function from a non-base package, in this case MASS, which is not indicated
in your sessionInfo() above, so something is amiss there as well.

Regards,

Marc Schwartz