[R] glm.nb error

[Sarah Goslee]

Hi,

On Fri, Jun 7, 2013 at 11: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

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

Is this what's in the "test" file that your code reads in? We don't
have that file, so can't run your code.

If it is, then the output of

dput(dat)

would be enormously more useful than copy and pasting your data file
into your email.

And yes, the full code is very little like the pair of lines you
originally pasted in.

Sarah

> 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
>>
>



-- 
Sarah Goslee
http://www.functionaldiversity.org