Choosing one variable from each of 3 buckets of variables I have a regression model that looks like the following
glm.nb(formula = y ~ Gender + Age + x1 + x2 + x3, data = df)

In my problem, there are 20 possible choices of variables for x1, 20 possible choices for x2, and 20 possible choices for x3.  Gender and Age must be in the model.  This leaves me with 20*20*20 = 8,000 possible regressions.  I was able to create a program that ran all of these regressions and deliver me the lowest AIC, but I was wondering if there was a library that already does this.
I do not consider what I will find to be the "best" model in any statistical manner, but I do find this exercise useful for exploring my data.
I have already attempted using bestglm and leaps.  I do not believe these programs allow for specifying the choice of variable from multiple bucket of variables.
 A: Since you have three categorical variables with $20$ categories each, plus gender and age, that gives you a total of $3 \times 19+1 = 58$ binary variables and one continuous variable.  If you are willing to proceed without interaction effects, that gives you a model with $60$ coefficients (including an intercept term).  That is a relatively manageable number of terms, and with a reasonable-sized data set, you should have an adequate number of residual degrees-of-freedom, and get reasonable estimates of your parameters.
Your problem emerges if you decide you want to include interaction effects between the categorical variables.  Adding two- and three-way interactions between your three categorical variables, gives you $20^3 = 8000$ parameters in your model instead of $57$.  (If you also interact with gender and age this quadruples again.)  That is a large number of parameters, and you would need a large amount of data to ensure that you have an adequate number of residual degrees-of-freedom.  Even if you have sufficient data for this, it is dubious to cherry-pick interaction terms from categorical variables in the manner you are describing - that is a classic example of data-dredging.  Instead of doing this, you should either include the whole set of interaction terms for a categorical interaction, or remove the whole set of interaction terms.  It is not legitimate to cherry-pick interaction terms from within indicator values for a broad categorical variable. 
A: If you have 60 possible covariates, and just want to be able to use the model to build predictions and are not that concerned with interpretability, you might build a random forest on a training set of your data and see what kind of predictive power you could get from the model it builds.  The package randomForest in r can help you with this.
