How can I select a model using subset selection and cross-validation?

Question

I try to find a model using logistic regression. More precisely, what I did so far, is using stepwise regression and subset selection (although I know, it is often a bad idea) to find the "best" model. Clearly, depending on the information criteria I used, I got different results.

Now, I found an interesting example on page 250 in the book "An Introduction to Statistical Learning". They chose among the models of different sizes using cross-validation, that is they make predictions for each model and compute the test errors. Eventually, the compute the cross validation error and choose the model corresponding to the minimal average cross-validation error.

However, the function regsubsets of the R package "leaps" is only working for linear models. How can I implement this for logistic regression or glm models in general?

My idea was, to just estimate the models within a cross-validation using the step function of the "stats" package and then kind of take the average number of features (which is determined by minimum AIC, for example). Is this a legitimate approach?

Answer 1

I would recommend moving on one more page, to page 251 of ISLR, and using ridge regression as illustrated there instead to get your "best" model.

If you repeat a search for the "best" variable-selection model on multiple bootstrap samples of the same data set, even with identical criteria for the "information criterion," you will almost certainly get a surprisingly large collection of "best" models. Try it on your data and see. So which one is really "best"? Will any generalize well? How will you correct for the bias introduced by looking for the "best" predictors in a particular data sample?

Unless you have a compelling need to minimize the number of predictor variables, ridge regression will give you an appropriately penalized model with a penalty chosen by cross-validation. Unlike variable-selection approaches, ridge regression handles collinearity among predictor variables well; it is likely to generalize well and the ridge coefficients should vary relatively little among bootstrap samples. If your interest is in prediction then the retention of all variables by ridge regression is likely to provide better performance than any variable-selection model. The glmnet package used in the exercise starting on page 251 includes facilities for ridge regression on generalized linear models, including logistic regression if you use the setting family = "binomial".

If you have a compelling need to cut down on predictor variables you can tweak the glmnet parameters to perform LASSO or the hybrid elastic net. The choice of penalty (and thus the number of variables) can be determined by cross-validation. As with any variable-elimination technique, however, there will still be some arbitrary selection among individual collinear predictors in models developed those ways.

Answer 2

Many packages provide cross validation directly. However, since you seem to be learning the principles, why not roll your own. For example like this

Divide the data into K, lets say 10, chunks. 
Round 1: Designate chunk 1 as test data, the other 9 as training data.
Build your model on the training data, test it on the test data.
Store the testing error. 

Round 2: Designate chunk 2 as test data, the other 9 as training data. 
...

And so on.

After 10 iterations, you have 10 error scores, one for each round. Average those and you have cross-validated your model. Once you have done that for many models, choose one and retrain it on the entire data set. That is then your final model.

Which error score to use? You say it's a logistic regression, is it a binary classifier? Are the two potential mis-classifications equally critical? If so, why not simply use hit-ratio? It does the job and importantly communicates well. If you have other problems, then other error functions might suit you better. AIC / BIC are well understood in the statistical community and are fine to use as well, but communicates less well to people outside the profession. You know your problem and your audience, keep both in mind.

With regards to penalized regression, which I didn't bring to mind as your question indicates you're learning CV, but @Edm is right in pointing them out. Here are some decent slides. http://statweb.stanford.edu/~tibs/sta305files/Rudyregularization.pdf

Understanding the difference of the L1 and L2 norm will help you grasp why Lasso (L1) might be better suited for model selection that Ridge (L2).

Good luck!