# Automatisation of GLM analysis with negative binomial errors

I am new to R and some help would be of great use to me.

Basically, I need to perform a GLM analysis with negative binomial errors and with fixed factors, no covariates and no random effects. My factors are of type: year (1-4), site (1-3), sex (1-2), age (1-3), with a sample size of around 5000.

Currently I am fitting GLMs with negative binomial errors, using the glm.nb routine in the MASS library of R. I then need to remove all nonsignificant interactions by hand (using the R console). Can someone point me into the right direction of how can I rather automate the process, in such a way that R will choose desired factors by trying different combinations and ignoring nonsignificant ones?

(I am sorry if the question lack details, I can add more info if you need it.)

The stepAIC function in MASS can perform the kinds of variable selection you are looking for.
In addition, the leaps package appears to have similar capacities. That being said, I have not used it, so cannot speak directly on its efficacy.
glmulti provides a wrapper for glm and similar functions (glm.nb, etc.), automatically generating all possible models (under constraints set by the user) with the specified response and explanatory variables, and finding the best models in terms of some Information Criterion (AIC, AICc or BIC). It can handle very large numbers of candidate models and features a Genetic Algorithm to find the best models when an exhaustive screening of the candidates is not feasible.