[R] glm.nb error
Sarah Goslee
sarah.goslee at gmail.com
Fri Jun 7 18:25:51 CEST 2013
As Marc already pointed out, take a close look at this part of your loop:
R> i <- 6
R>
R> y <- as.numeric(data[i,-1])
R> y
[1] 3 3 3 3 4 4 4 4
R> group
[1] 1 1 1 1 0 0 0 0
R> fit <- glm.nb(y~group)
Error in while ((it <- it + 1) < limit && abs(del) > eps) { :
missing value where TRUE/FALSE needed
What do you expect to happen there?
In general, it's a better practice to pre-specify the size of result
(eg matrix(NA, nrow=n, ncol=4) ) and fill it as you go, rather than
using rbind() within a loop. Much more memory-efficient.
Sarah
On Fri, Jun 7, 2013 at 11:58 AM, Daofeng Li <lidaof at gmail.com> wrote:
> Sorry Sarah.
>
>> dput(dat)
> structure(list(gene = structure(c(1L, 3L, 4L, 5L, 6L, 7L, 8L,
> 9L, 10L, 2L), .Label = c("gene1", "gene10", "gene2", "gene3",
> "gene4", "gene5", "gene6", "gene7", "gene8", "gene9"), class = "factor"),
> b1 = c(18L, 15L, 10L, 4L, 0L, 3L, 0L, 4L, 11L, 6L), b2 = c(15L,
> 8L, 9L, 0L, 1L, 3L, 4L, 4L, 6L, 3L), b3 = c(13L, 8L, 8L,
> 4L, 0L, 3L, 0L, 7L, 13L, 6L), b4 = c(13L, 7L, 12L, 3L, 0L,
> 3L, 2L, 3L, 10L, 3L), c1 = c(16L, 0L, 9L, 0L, 1L, 4L, 2L,
> 0L, 13L, 4L), c2 = c(9L, 12L, 11L, 5L, 5L, 4L, 2L, 6L, 7L,
> 7L), c3 = c(20L, 18L, 12L, 0L, 1L, 4L, 0L, 6L, 12L, 6L),
> c4 = c(24L, 4L, 12L, 0L, 0L, 4L, 0L, 12L, 9L, 3L)), .Names = c("gene",
> "b1", "b2", "b3", "b4", "c1", "c2", "c3", "c4"), class = "data.frame",
> row.names = c(NA,
> -10L))
>
> above was the dput(dat). Thanks.
>
> Daofeng
>
>
> On Fri, Jun 7, 2013 at 10:47 AM, Sarah Goslee <sarah.goslee at gmail.com>
> wrote:
>>
>> 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
More information about the R-help
mailing list