Skip to content

question on nbinom1

6 messages · Mollie Brooks, Don Cohen, Ben Bolker

#
In response to a question about use of nbinom1 in glmTMB, specifically
about this paper by VerHoef J.M. & Boveng:
 https://digitalcommons.unl.edu/cgi/viewcontent.cgi?article=1141&context=usdeptcommercepub
which says that AIC cannot be used to choose between a quasi-Poisson model
and a negative binomial,

Ben Bolker writes:

 >     nbinom1 is *not* quasi-: see Hardin & Hilbe 2007 as suggested in 
 > ?nbinom1.  (It can't be fitted in a standard GLM framework, since this 
 > parameterization of the nbinom doesn't fit into the exponential family, 
 > but glmmTMB is more flexible than that ...)

I've been trying to figure out what Hardin & Hilbe have to say about this,
and so far failing to find the connection.
(BTW, is nbinom1 the same as what they call NB-C ?)

Does "nbinom1 is *not* quasi-" mean that the AIC reported when I use 
nbinom1 really IS comparable to the one reported when I use nbinom2?
And that loglik is actually being computed for a "real" distribution?
Does it matter whether I also use random effects, zero inflation, 
offset, etc. ?

Alternatively, if AIC for nbinom1 is really NOT comparable to AIC for the
other (real?) distributions, how is one supposed to determine which model
better fits the data?  VerHoef suggests trying to determine whether the
variance for different subpopulations with different means more closely
resembles linear or quadratic functions of the mean, but this is not at
all clear in the data I'm trying to analyze.
#
I don?t have a copy of Hardin & Hilbe 2007 on hand, but I answered a few of your questions below.
Yes
Yes
No it doesn?t matter; the AIC is comparable when fitting those types of models. 

Some examples are 
https://cran.r-project.org/web/packages/glmmTMB/vignettes/glmmTMB.pdf <https://cran.r-project.org/web/packages/glmmTMB/vignettes/glmmTMB.pdf>
and 
appendix A here https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf <https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf>

cheers,
Mollie

  
  
#
Mollie Brooks writes:

 > I don't have a copy of Hardin & Hilbe 2007 on hand, but I answered
 > a few of your questions below.

Thank you.

One more question:

How can I compute the nbinom1 distribution?
Is there a formula for the pdf or cdf ?  An R function ?
1 day later
#
I think the easiest way to get a numerical representation of the distribution from a fitted model would be using the simulate function.

There?s an example of how to do that on pages 392-393 of this pdf (including Figs 6 and 7)
https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf <https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf>

cheers,
Mollie

  
  
1 day later
#
I believe the R d/p/q/r functions corresponding to glmmTMB's 
implementation of *nbinom1 would look like this:


rnbinom1 <- function(n, mu, phi) {
     ## var = mu*(1+phi) = mu*(1+mu/k) -> k = mu/phi
     rnbinom(n, mu=mu, size=mu/phi)
}

dnbinom1 <- function(x, mu, phi, ...) {
     dnbinom(n, mu=mu, size=mu/phi, ...)
}

pnbinom1 <- function(q, mu, phi, ...) {
     pnbinom(q, mu=mu, size=mu/phi, ...)
}

qnbinom1 <- function(p, mu, phi, log=FALSE) {
     pnbinom(p, mu=mu, size=mu/phi, ...)
}


   (there would be an even more clever/inscrutable way to do this by 
transforming the body of the code, without repeating oneself so much, 
but it would probably be a bad idea)
On 10/12/20 6:34 AM, Mollie Brooks wrote:
#
Thank you.
This is perfectly clear and enables further analysis.
I hope you can arrange to put it in some (even better, much) documentation.
I wish it were in wikipedia.

Ben Bolker writes:
 > 
 >    I believe the R d/p/q/r functions corresponding to glmmTMB's 
 > implementation of *nbinom1 would look like this:
 > 
 > 
 > rnbinom1 <- function(n, mu, phi) {
 >      ## var = mu*(1+phi) = mu*(1+mu/k) -> k = mu/phi
 >      rnbinom(n, mu=mu, size=mu/phi)
 > }
 > 
 > dnbinom1 <- function(x, mu, phi, ...) {
 >      dnbinom(n, mu=mu, size=mu/phi, ...)
 > }
 > 
 > pnbinom1 <- function(q, mu, phi, ...) {
 >      pnbinom(q, mu=mu, size=mu/phi, ...)
 > }
 > 
 > qnbinom1 <- function(p, mu, phi, log=FALSE) {
 >      pnbinom(p, mu=mu, size=mu/phi, ...)
 > }
 > 
 > 
 >    (there would be an even more clever/inscrutable way to do this by 
 > transforming the body of the code, without repeating oneself so much, 
 > but it would probably be a bad idea)
 >
> On 10/12/20 6:34 AM, Mollie Brooks wrote:
> > I think the easiest way to get a numerical representation of the distribution from a fitted model would be using the simulate function.
 > > 
 > > There's an example of how to do that on pages 392-393 of this pdf (including Figs 6 and 7)
 > > https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf <https://journal.r-project.org/archive/2017/RJ-2017-066/RJ-2017-066.pdf>
 > > 
 > > cheers,
 > > Mollie
 > >
> >> On 10Oct 2020, at 14:23, Don Cohen <don-lme4 at isis.cs3-inc.com> wrote:
> >>
 > >> Mollie Brooks writes:
 > >>
 > >>> I don't have a copy of Hardin & Hilbe 2007 on hand, but I answered
 > >>> a few of your questions below.
 > >>
 > >> Thank you.
 > >>
 > >> One more question:
 > >>
 > >> How can I compute the nbinom1 distribution?
 > >> Is there a formula for the pdf or cdf ?  An R function ?
 > >>
 > > 
 > > 
 > > 	[[alternative HTML version deleted]]
 > > 
 > > _______________________________________________
 > > R-sig-mixed-models at r-project.org mailing list
 > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
 > >
 > 
 > _______________________________________________
 > R-sig-mixed-models at r-project.org mailing list
 > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models