# Sane stepwise regression?

Suppose I want to build a binary classifier. I have several thousand features and only a few 10s of samples. From domain knowledge, I have a good reason to believe that the class label can be accurately predicted using only a few features, but I have no idea which ones. I also want the final decision rule to be easy to interpret/explain, further necessitating a small number of features. Certain subsets of my features are highly correlated, so selecting the most predictive few independently wouldn't work. I also want to be able to meaningfully do hypothesis testing on my features.

Is the following stepwise regression procedure reasonable under these conditions:

1. Given the features already in the model (or just the intercept on the first iteration), select the feature that produces the largest log likelihood ratio when added to the model. Use the likelihood ratio chi-square test to calculate a nominal P-value for each hypothesis test performed in this selection. The null here is that adding the extra variable to the model provides no additional predictive ability. The alternative is that it does increase predictive abilityl

2. Treat the hypotheses tested in Step 1 of each iteration as a family and calculate the false discovery rate for the smallest P-value (for the feature selected) using something like Benjamini-Hochberg.

3. Goto 1 unless some stopping criteria are met.

4. Report the false discovery rates for the individual features, but not the P-value for the model as a whole (since this will be massively inflated). Each of these multiple testing corrected P-values represents the statistical significance of that feature given all of the features previously added to the model.

Does doing something like this under these circumstances successfully avoid all of the typical criticisms of stepwise regression? Are the false discovery rates calculated in this way reasonable?

-
Is there a reason not to go with a penalized regression approach (lasso, elasticnet, etc.)? –  Ben Bolker Jan 29 '11 at 18:36

For the purposes of my answer, I will denote the binary variable of interest as $Y_i \text{ ;}(i=1,\dots,n)$ and the predictors $X_{ij} \text{ ;} (j=1,\dots,p)$ and assume that $Y$ has values of $Y=0$ and $Y=1$. It will also be convenient to define $\gamma_m$ to indicate the model $m \text{ ;}(m=1,..,M)$, such that $\gamma_m^TX_{ij}$ is equal to $X_{ij}$ if the jth variable is in the mth model, and $0$ otherwise.

I would make a modification to your method, and give a rationale. You are using a classifier model, which means you want to predict the value of a categorical variable into the future - so you should really be defining a prediction rule (given a new set of predictors $X_{j}$, how will you predict whether $Y=1$ or $Y=0$).

So I would suggest evaluating the prediction directly, rather than the likelihood ratio. However, the observation predicted should not be included in the estimation of the model (because this is exactly the situation you will face when actually using your model). So have a new step 1) (bold is my suggested change). 1) Given the features already in the model (or just the intercept on the first iteration), select the feature that produces the best predictions when added to the model.

Now you need to decide

1. what you want "best" to mean mathematically
2. how to split your data into "fitting" and "predicting" parts

I will make a suggestion for each:

1. An intuitive definition for a "good" classifier (and also computationally simple) is the proportion of correct classifications it makes. However, you may have some additional knowledge of the specific consequences of making a correct or incorrect classification (e.g. predicting correctly when $Y=1$ is twice as important than when $Y=0$). In this case you should incorporate this knowledge into the definition of "good". But for the equations in my answer I will use $F=\frac{C}{C+I}$ as the criterion ($F$="fraction" or "frequency" $C$="correct" $I$="incorrect")
2. Because you don't have a lot of data, you need as much as possible to fit the model, so a simple drop one jacknife procedure can be used. You leave observation $1$ out, fit the model with observations $2,\dots,n$, and use this to predict observation $1$. Then you leave observation $2$ out, fit the model with observations $1,3,\dots,n$, and use this to predict observation $2$; and so on until each observation has been "left out" and predicted. You will then have $n$ predictions, and you can now calculate $F=\frac{C}{n}$, the fraction of correctly predicted values for the particular model. Subscript this for the particular model $F_m$.

You would then calculate $F_m$ for each model $(m=1,\dots,M)$, and pick the model which predicts the best $m=\text{argmax}_{m\in M} F_m$. Note that the good thing about the above method is that you do not need to worry about how many variables are in your model or how correlated these variables are (unless it makes it impossible to actually fit the model). This is because the model is fit separately to the prediction, so bias due to over-fitting, or degradation due to numerical instability will show up in the poorer predictions.

In a step-wise situation it is done sequentially, so at the $sth$ you have $M_s=p+1$ models to choose between: one each for "removing" each $X_{j}$ which is in the model, one for "adding" each $X_{j}$ which is not in the model, and one for keeping the model unchanged (you stop the procedure when you choose this model, and this is your final model). If there is a tie, you need an additional criteria to split the winners (or you could let you algorithm "branch" off, and see where each "branch" ends up, then take the "branch" which had the best predictions at its final step)

Step-wise can be risky because you may find "local maximums" instead of "global maximums", especially because you have such a large number of predictors (this is a big "space" to optimise over, and is probably multi-modal - meaning there are many "best" models)

The good thing about this is that the model you choose has a clear, directly relevant interpretation: The model which predicted the highest proportion of results correctly, out of the alternatives considered. And you have a clear measure of exactly how good your binary classifier is (it classified $100F$ percent correctly).

I think you will find this a lot easier to justify your choice of final model to a non-statistician, rather than trying to explain why the p-value indicates the model is good.

And for hypothesis testing, you can declare any effect remaining in your final model as "significant" in that the relationships contained in this model we able to re-produce the data ($Y$) the most effectively.

Two final remarks:

1. You could also use this machinery to decide if step-wise is better than forward selection (only add variables) or backward selection (start from full model, and only remove variables).
2. You could fit the full model (or any model with $p\geq n$) by "ridging" the model, which amounts to adding a small number to the diagonal elements of the $X^TX$ matrix, or $X^TWX$ for GLMs before inverting when calculating your betas, to give $(X^TX+\lambda I)^{-1}X^TY$ or $(X^TWX+\lambda I)^{-1}X^TWY$. Basically, $\lambda$ constrains the sum of squares of the betas to be less than a particular value, increasing the value of $\lambda$ decreases this constraint (which is a "smooth" model selection procedure in its own right, if you think about it).
-
Thanks for the advice. There are three problems with this, though: 1. I care about quantifying the uncertainty in my predictions and the contribution of each variable, not just binary prediction accuracy. 2. Given the nature of my dataset, it's way too computationally intensive. 3. From domain knowledge, I believe that the local optima issue is not important. –  dsimcha Jan 29 '11 at 17:18
How can it be computationally intensive, when you have less than $100$ observations? The drop-one-jacknife is very quick for small $n$. And the fraction $F$ is a measure of the uncertainty of your predictions, namely the probability that you will classify correctly given your training data, and your selected model. –  probabilityislogic Jan 29 '11 at 18:16
Actually you're right. This procedure would be part of a larger codebase and I forgot that some of the rest of the code wouldn't need to be re-run for every jackknife iteration. The other two points still apply, though. –  dsimcha Jan 29 '11 at 18:31
@dsimcha - if all you're doing is building a binary classifier, why would you care about anything other than the predictions it makes? Surely this has to be the number 1 priority, with the contribution from each variable only needed to help understand why it is the best classifier. A way you can get the contribution for each variable, is to re-calculate $F$ removing each variable, one at a time, and check how this compares to the $F$ for the chosen model ($F_{chosen}-F_{(-j)}$ is like a marginal effect for the jth variable on the predictive accuracy due to each variable) –  probabilityislogic Jan 30 '11 at 1:10
The other issue is that my answer is generic in the criterion, because I don't know what the consequences of using this classifier are. If you know more background information, then $F$ should be changed to reflect what is important when you will actually be using the classifier and not what's important when you are fitting the model. I believe I made this clear in my answer. –  probabilityislogic Jan 30 '11 at 7:10