# No valid coefficients for NegBin regression

I am doing multiple regression with some data (5 predictors, 1 response). Since the response is discrete and non-negative, I thought I would try Poisson regression. However, the data are significantly overdispersed (variance > mean), so I am now trying negative binomial regression.

I was able to fit the model with this code.

library(MASS)
model.nb <- glm.nb(Response ~ Pred1 + Pred2 + Pred3 + Pred4 + Pred5 - 1, data=d)


Now I would like to see if I can get a better fit by including interactions between the predictors. However, when I try to do so, I get the following error.

> model.nb.intr <- glm.nb(Response ~ Pred1 * Pred2 * Pred3 * Pred4 * Pred5 - 1, data=d)
Error: no valid set of coefficients has been found: please supply starting values


Any ideas what may be causing this?

-

Before jumping to a model that includes all interactions, you can try adding only the 2-way interactions:

model.nb.intr <- glm.nb(Response ~ (Pred1 + Pred2 + Pred3 + Pred4 + Pred5)^2 - 1, data=d)

-
slightly simpler would be: Response ~ .^2 -1, data = DF, but excellent point! – Gavin Simpson Dec 13 '10 at 18:00
Thanks for this suggestion! – Daniel Standage Dec 13 '10 at 20:45
Ok, ran the code based on this suggestion and now I have the error "Warning while fitting theta: alternation limit reached". What is the alternation limit and how do I address this? – Daniel Standage Dec 13 '10 at 20:55
As described on the ?glm.nb help page: "An alternating iteration process is used. For given theta the GLM is fitted using the same process as used by glm(). For fixed means the theta parameter is estimated using score and information iterations. The two are alternated until convergence of both." Your model does not converge (you reached the maximum of the theta iterations). You could try setting the trace parameter to see what's going on. Or you could use a different package - I believe VGAM might use a different algorithm. – Aniko Dec 13 '10 at 21:04

Your model is too complex for the computer to work out some reasonable starting values that do not lead to infinite deviance when doing the glm.fit iterations.

Have you got enough data to estimate all these interactions? Do you think it is plausible for all predictors to interact with each other? If not, think about which predictors might interact and include only those terms.

The error asks you to supply some starting values for it to work from. For this, you need to supply a vector of parameter values as argument start; from ?glm:

   start: starting values for the parameters in the linear predictor.


You need to supply 31 model parameters (I hope you have many 1000s of data points?) to start, in this order:

> colnames(model.matrix(Y ~ Pred1*Pred2*Pred3*Pred4*Pred5 -1, data = DF))
[1] "Pred1"                         "Pred2"
[3] "Pred3"                         "Pred4"
[5] "Pred5"                         "Pred1:Pred2"
[7] "Pred1:Pred3"                   "Pred2:Pred3"
[9] "Pred1:Pred4"                   "Pred2:Pred4"
[11] "Pred3:Pred4"                   "Pred1:Pred5"
[13] "Pred2:Pred5"                   "Pred3:Pred5"
[15] "Pred4:Pred5"                   "Pred1:Pred2:Pred3"
[17] "Pred1:Pred2:Pred4"             "Pred1:Pred3:Pred4"
[19] "Pred2:Pred3:Pred4"             "Pred1:Pred2:Pred5"
[21] "Pred1:Pred3:Pred5"             "Pred2:Pred3:Pred5"
[23] "Pred1:Pred4:Pred5"             "Pred2:Pred4:Pred5"
[25] "Pred3:Pred4:Pred5"             "Pred1:Pred2:Pred3:Pred4"
[27] "Pred1:Pred2:Pred3:Pred5"       "Pred1:Pred2:Pred4:Pred5"
[29] "Pred1:Pred3:Pred4:Pred5"       "Pred2:Pred3:Pred4:Pred5"
[31] "Pred1:Pred2:Pred3:Pred4:Pred5"


I would use the coefs from the first model to fill in the first 5 starting values and then what you do about the others is up to you. You could try starting them all off at 1 and see if that will get the model to fit?

You might also benefit from code in the pscl package which can fit hurdle and zero-inflated models to count data.

-
 I've got plenty of data (almost 300k data points), but I am not convinced there are significant interactions between each predictor. I wanted to explore this possibility, but I guess I bit off more than I could chew. – Daniel Standage Dec 13 '10 at 17:39 @Daniel by all means then give it a go, but you will need to supply some starting points for the algorithm – Gavin Simpson Dec 13 '10 at 17:47

If you can't get satisfaction with R you can fit this model and more complicated ones with AD Model Builder which is free software available at http://admb-project.org. ADMB permits you to model the over dispersion in a variety of ways, rather than being confined to the GLM paradigm. I can advise you if you are interested.

-
 Dave, wouldn't glmmADMB work for this problem (as a start)? – Ben Bolker Dec 14 '10 at 22:42 It might require some tweaking. glmmADMB is intended to fit a Negative Binomial mixed model. It would be fairly simple to modify it for just a Negative Binomial regression. – dave fournier Dec 15 '10 at 14:16