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?
 A: 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.
A: 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.
A: 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)

