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
On Fri, Jun 7, 2013 at 10:15 AM, Marc Schwartz <marc_schwa...@me.com> wrote:
>
> On Jun 7, 2013, at 9:44 AM, Daofeng Li <lid...@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
>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
