# How to do logistic regression subset selection?

I am fitting a binomial family glm in R, and I have a whole troupe of explanatory variables, and I need to find the best (R-squared as a measure is fine). Short of writing a script to loop through random different combinations of the explanatory variables and then recording which performs the best, I really dont know what to do. And the leaps function from package leaps does not seem to do logistic regression.

Any help or suggestions would be greatly appreciated Leendert

-
There exist functions that perform automatic search. You should have a look at the step function. Section 5.4 illustrates that point: data.princeton.edu/R/glms.html – ocram Mar 15 '11 at 8:43
I'm sorry but my post has been edited so that it no longer asks my question. I have 35 (26 significant) explanatory variables in my logistic regression model. I need the best possible combination of 8, not the best subset, and at no point was I interested in a stepwise or all subsets style approach. There is no wiggle room in this 8. I just thought someone may have know how I could fit all combinations of 8 explantory variables and it could tell me which maximises the likelihood (sorry about the R-squared brain fart but AIC isnt relevant either since I have a fixed number of parameters, 8). – Leendert Mar 16 '11 at 10:25
You can revert to the previous version of your post, or combine both edits. I'm sure @mpiktas was of good intention when trying to improve its appearance and just didn't notice the No. parameters. – chl Mar 16 '11 at 18:09
@ Everyone: Thank you very much. In the end I used many different thing in the hope they would all give similar answers. And they did. I used the BMA, bestglm and glmnet packages as well as the step function. Fitted models with all of them, and there was not discrepancy in what BMA with maxcol = 9 and step deemed the best model. All the experts in the field around me seemed very content with the variables, and felt that it was quite progressive. So thanks for all the input. I really used it all. – Leendert Mar 23 '11 at 8:22

Stepwise and "all subsets" methods are generally bad. See Stopping Stepwise: Why Stepwise Methods are Bad and what you Should Use by David Cassell and myself (we used SAS, but the lesson applies) or Frank Harrell Regression Modeling Strategies. If you need an automatic method, I recommend LASSO or LAR. A LASSO package for logistic regression is available here, another interesting article is on the iterated LASSO for logistic

-
(+1) About R packages, there're also glmnet (implementation with coordinate descent algo., Friedman and coll.) and penalized (allows to keep some var. unpenalized). Of note, F. Harrell provides penalized ML estimation for GLMs (see lrm, or his RMS textbook for further info). – chl Mar 15 '11 at 11:31
(+1) Nice article, it seems I have to start going far beyond the author states in the question (not the first time I didn't). @chl (+1) perfect alternative suggestions too. – Dmitrij Celov Mar 15 '11 at 15:04
@chl: +1 for glmnet, that's a great package. – Zach Mar 16 '11 at 16:07
@chl Thanks! One of the problems with R is keeping track of packages (there are so many!) and which are best. The task views do help – Peter Flom Mar 18 '11 at 10:42

First of all $R^2$ is not an appropriate goodness-of-fit measure for logistic regression, take an information criterion $AIC$ or $BIC$, for example, as a good alternative.

Logistic regression is estimated by maximum likelihood method, so leaps is not used directly here. An extension of leaps to glm() functions is the bestglm package (as usually recommendation follows, consult vignettes there).

You may be also interested in the article by David W. Hosmer, Borko Jovanovic and Stanley Lemeshow Best Subsets Logistic Regression // Biometrics Vol. 45, No. 4 (Dec., 1989), pp. 1265-1270 (usually accessible through the university networks).

-
+1 For commenting on $R^2$. – whuber Mar 15 '11 at 15:22
While your comment about $R^2$ being worse than $BIC,AIC$ is useful in general, it actually makes no difference unless you are comparing models of different sizes. The OP clearly states that they are only interest in $8$ variable models, so $BIC$ and $AIC$ will revert back to choosing the model with the highest likelihood. This is equivalent to $R^2$ fitting. – probabilityislogic Dec 21 '11 at 13:32
Thanks for the remark, but comments by chl below explains why fixed number of explanatory variables is dangerous. By the way the answer appeared earlier than the comment regarding $8$ (up to?) variables rstriction ;) – Dmitrij Celov Dec 22 '11 at 10:18

One idea would be to use a random forest and then use the variable importance measures it outputs to choose your best 8 variables. Another idea would be to use the "boruta" package to repeat this process a few hundred times to find the 8 variables that are consistently most important to the model.

-
@Zach Do you suggest to rely on RFs to perform feature selection, and then apply a GLM -- in this case, there's risk of overfitting or over-optimism --, or to use RFs (with standard measures of var. importance, or all-relevant selection) as a standalone tool? – chl Mar 16 '11 at 18:05
@chl: I was suggesting using the RFs to perform feature selection, and then apply the GLM. I agree there's risk of over fitting, but the OP said he needed exactly 8 variables. – Zach Mar 16 '11 at 20:47
@Zach "exactly 8 variables"... Then, you are somewhat soft-thresholding the potential variables of interest based on a measure of variable importance (which thanks to being based on permutation and double resampling is supposed to be bias-free) and then reinject them into an GLM. IMHO, you break the control on overfitting exerted through bagging. This is also described in ESLII from Hastie et al.: feature selection, if any, must be included in the cross-validation procedure (where cross-validation include the assessment of model performance). – chl Mar 16 '11 at 20:57
@chl: Doesn't the r package "boruta" cross validate by running the random forest several time? Are you basically saying we need "meta-cross-validation," where you do variable selection and and the fit your model on random subsets of the data? – Zach Mar 17 '11 at 2:49
@Zach My point was that bagging, or in particular RFs, prevent from overfitting (to a certain extent) provided you remain in the same toolchain. If you use the results from RFs and then go see how another model would perform with the same data, then you break the CV loop. But why not use directly RFs for classification? Another solution would be to use RFs on a training sample and then apply the GLM on an held-out sample (that could potentially be cross-validated too). – chl Mar 17 '11 at 7:56