[R] glm.nb, anova.negbin

Prof Brian Ripley ripley at stats.ox.ac.uk
Thu Sep 27 13:39:56 CEST 2001


On Thu, 27 Sep 2001, juli g. pausas wrote:

> That's solve the problem. Thank you very much.
> I think that there is still a minor problem, in the display of the results only.
> When ommiting the intercept, the dispaly in the anova do not show the P(>|Chi|) of
> the first var. This is very minor problem but you may like to know it.

It is not a problem at all.  Look at the df. Those two models are not
comparable: one has an intercept, the other does not. Ideally this would
fit a null model, but again that's an R/S difference as to how glm
internals work.

> t.tax <- rnbinom(200,mu=1,size=1)
> t.var1 <- rnorm(200, t.tax, 5)
> t.var2 <- runif(200)
> t <- glm.nb(t.tax ~ t.var1)
> t2 <- update(t, .~ .+t.var2)
> t3 <- update(t3, .~ . -1 )
> anova(t3)
>
> Analysis of Deviance Table
>
> Model: Negative Binomial(1.3652), link: log
>
> Response: t.tax
>
> Terms added sequentially (first to last)
>
>
>         Df Deviance Resid. Df Resid. Dev P(>|Chi|)
> NULL                      199    205.098
> t.var1   0    4.597       199    200.501
> t.var2   1    1.534       198    198.967     0.215
> t.var3   1    0.055       197    198.912     0.815
> Warning messages:
> 1: tests made without re-estimating theta in: anova.negbin(t3)
> 2: NaNs produced in: pchisq(q, df, lower.tail, log.p)
>
>
>
> Prof Brian Ripley escribió:
>
> > That's an error in anova.glm: try it with glm rather than glm.nb.
> >
> > The line in anova.glm like
> >
> >             fit <- method(x = x[, varseq <= i], y = object$y,
> >
> > should be
> >
> >             fit <- method(x = x[, varseq <= i, drop = FALSE], y = object$y,
> >
> > You are really good at finding long-standing bugs!
> >
> > On Wed, 26 Sep 2001, juli g. pausas wrote:
> >
> > > > > I looked quickly, and it seems that this error comes actually from
> > > > > `anova.glm'  and from the following `get':
> > > > >
> > > > >         if (!is.function(method))
> > > > >             method <- get(method, mode = "function")
> > > > >
> > > > > `anova.negbin' sends the  `negbin' `object' to `anova.glm', but the
> > > > > negbin object does not have a method (`glm' objects have method
> > > > > "glm.fit"). There seem to be some other differences in `negbin' and
> > > > > `glm' objects, so that `anova.glm' can't handle them (e.g., glm.control
> > > > > object is missing in `negbin' object).
> > > > >
> > > > > I would like to help, but I don't have time now (got to mark exams).
> > > >
> > > > Yes, that's it: R is incompatible with S.  Although the information
> > > > is in the saved call, R uses a duplicate copy.  In glm.nb() adding
> > > >     fit$method <- method
> > > >     fit$control <- control
> > > > at the end seems to fix it, but without an example I can reproduce,
> > > > I cannot of course be sure.
> > > >
> > >
> > > Thank you very much. This fix the problem on my own dataset.
> > > I have tried to test in other reproducible exemples, following the example by
> > > Ben Bolker (many thanks!), and it works well. I've only got problems when
> > > trying to fit -1. I'm not sure if this is because it doen't make sense with
> > > this dataset or because there is still a bug there (see below). Anyway, my
> > > main problem is solved now. Many thanks.
> > >
> > > t.tax <- rnbinom(200,mu=1,size=1)
> > > t.var1 <- runif(200)
> > > t.var2 <- runif(200)
> > > t.var3 <- runif(200)
> > > t.var4 <- runif(200)
> > > t.var5 <- runif(200)
> > >
> > > ## ignore warning messages about convergence
> > > t <- glm.nb(t.tax ~ t.var1)
> > > t2 <- update(t, .~ .+t.var2)
> > > t3 <- update(t2, .~ .+t.var3)
> > > t4 <- update(t3, .~ .+t.var4)
> > > t5 <- update(t4, .~ .+t.var5)
> > > t6 <- update(t5, .~ . -1 )
> > >
> > >
> > > anova(t, t2, t3, t4, t5, t6) ## OK
> > > anova(t5)   ## OK
> > > summary(t6) ## OK
> > > anova(t6)   ## error
> > >
> > > Error in "names<-.default"(*tmp*, value = c("", "", "", "", "", "", "",  :
> > >         names attribute must be the same length as the vector
> > > In addition: Warning message:
> > > tests made without re-estimating theta in: anova.negbin(t6)
> > >
> > >
> > >
> > >
> > > > And that was my main point: if reporting problems *please* report the
> > > > results of traceback(), and if possible a small reproducible example.
> > > >
> > > > --
> > > > Brian D. Ripley,                  ripley at stats.ox.ac.uk
> > > > Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> > > > University of Oxford,             Tel:  +44 1865 272861 (self)
> > > > 1 South Parks Road,                     +44 1865 272860 (secr)
> > > > Oxford OX1 3TG, UK                Fax:  +44 1865 272595
> > >
> > > --
> > > Juli G. Pausas   (Mr)
> > > Centro de Estudios Ambientales del Mediterraneo (CEAM)
> > > C/ C.R. Darwin 14, Parc Tecnologic,
> > > 46980 Paterna, Valencia, SPAIN
> > > Tel: (+ 34) 96 131 8227; Fax: (+ 34) 96 131 8190
> > > mailto:juli at ceam.es
> > > http://www.gva.es/ceam
> > >
> > > GCTE Fire Network - http://www.gva.es/ceam/FireNetwork
> > >
> > >
> > >
> >
> > --
> > Brian D. Ripley,                  ripley at stats.ox.ac.uk
> > Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> > University of Oxford,             Tel:  +44 1865 272861 (self)
> > 1 South Parks Road,                     +44 1865 272860 (secr)
> > Oxford OX1 3TG, UK                Fax:  +44 1865 272595
>
> --
> Juli G. Pausas   (Mr)
> Centro de Estudios Ambientales del Mediterraneo (CEAM)
> C/ C.R. Darwin 14, Parc Tecnologic,
> 46980 Paterna, Valencia, SPAIN
> Tel: (+ 34) 96 131 8227; Fax: (+ 34) 96 131 8190
> mailto:juli at ceam.es
> http://www.gva.es/ceam
>
> GCTE Fire Network - http://www.gva.es/ceam/FireNetwork
>
>
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list