[R] Zero inflated negative binomial

[Achim Zeileis]

On Thu, 4 Mar 2010, Rebecca Lawrence wrote:

> Hi all,
>
> I am running the following model:
>
>> glm89.nb <- glm.nb(AvGUD ~ Year*Trt*Micro)
> where Year has 3 levels, Trt has 2 levels and Micro has 3 levels.
>
> However when I run it has a zero inflated negative binomial (as I have lots
> of zeros) I get the below error message:
>
>> Zinb <- zeroinfl(AvGUD ~ Year*Trt*Micro |1, data = AvGUD89, dist =
> "negbin")
> Error in optim(fn = loglikfun, gr = gradfun, par = c(start$count,
> start$zero,  :
>  non-finite value supplied by optim
>
> For what I have read I think the problem is that for Year level 3 there is
> no Trt 1 and for Year level 1 there is no Micro level 3.

Yes, in that case a number of coefficients are not identified. glm.nb() 
handles this more gracefully and sets the non-identified parameters to NA. 
I'll think about whether I can do something similar for zeroinfl/hurdle.

For the moment, here's a workaround: You can set up a data frame that 
contains only the identified regressors and call zeroinfl() with that. For 
illustration I use the quine data from the "MASS" package:

## packages and data
library("MASS")
library("pscl")
data("quine", package = "MASS")

## subset of data with some combination missing
quine2 <- subset(quine, !(Age == "F3" & Sex == "M"))

## glm.nb() works and yields some NA coefficients
fm1 <- glm.nb(Days ~ Eth * Sex * Age, data = quine2)

## zeroinfl() fails
fm2 <- zeroinfl(Days ~ Eth * Sex * Age | 1, data = quine2,
   dist = "negbin")

## set up version of data with non-identified regressors omitted
quine3 <- as.data.frame(model.matrix(fm1)) ## all regressors
quine3 <- quine3[, !is.na(coef(fm1))]      ## only identified
quine3 <- quine3[, -1]                     ## omit intercept
quine3$Days <- quine2$Days                 ## add response

## re-fit glm.nb()
fm1a <- glm.nb(Days ~ ., data = quine3)
## equivalent to previous fit
logLik(fm1a) - logLik(fm1)
coef(fm1a) - na.omit(coef(fm1))

## fit zeroinfl(), now works
fm2a <- zeroinfl(Days ~ . | 1, data = quine3, dist = "negbin")

hth,
Z

> I cannot find a solution to this problem, is there any way I can solve this
> so I can run the zero inflated model?
>
> Thank you for your time
> Rebecca
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.
>