Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.
up vote 5 down vote favorite

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?

share|improve this question

3 Answers

up vote 8 down vote accepted

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)
share|improve this answer
2  
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
up vote 8 down vote

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.

share|improve this answer
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
up vote 3 down vote

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.

share|improve this answer
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.