How to deal with multicollinearity when performing variable selection?

Question

I have a dataset with 9 continuous independent variables. I'm trying to select amongst these variables to fit a model to a single percentage (dependent) variable, Score. Unfortunately, I know there will be serious collinearity between several of the variables.

I've tried using the stepAIC() function in R for variable selection, but that method, oddly, seems sensitive to the order in which the variables are listed in the equation...

Here's my R code (because it's percentage data, I use a logit transformation for Score):

library(MASS)
library(car)

data.tst = read.table("data.txt",header=T)
data.lm = lm(logit(Score) ~ Var1 + Var2 + Var3 + Var4 + Var5 + Var6 + Var7 +
             Var8 + Var9, data = data.tst)

step = stepAIC(data.lm, direction="both")
summary(step)

For some reason, I found that the variables listed at the beginning of the equation end up being selected by the stepAIC() function, and the outcome can be manipulated by listing, e.g., Var9 first (following the tilde).

What is a more effective (and less controversial) way of fitting a model here? I'm not actually dead-set on using linear regression: the only thing I want is to be able to understand which of the 9 variables is truly driving the variation in the Score variable. Preferably, this would be some method that takes the strong potential for collinearity in these 9 variables into account.

Answer 1

First off, a very good resource for this problem is T. Keith, Multiple Regression and Beyond. There is a lot of material in the book about path modeling and variables selection and I think you will find exhaustive answers to your questions there.

One way to address multicollinearity is to center the predictors, that is substract the mean of one series from each value. Ridge regression can also be used when data is highly collinear. Finally sequential regression can help in understanding cause-effect relationships between the predictors, in conjunction with analyzing the time sequence of the predictor events.

Do all 9 variables show collinearity? For diagnosis you can use Cohen 2003 variance inflation factor. A VIF value >= 10 indicates high collinearity and inflated standard errors. I understand you are more interested in the cause-effect relationship between predictors and outcomes. If not, multicollinearity is not considered a serious problem for prediction, as you can confirm by checking the MAE of out of sample data against models built adding your predictors one at the time. If your predictors have marginal prediction power, you will find that the MAE decreases even in the presence of model multicollinearity.

Answer 2

Because it is so hard to determine which variables to drop, it is often better not to drop variables. Two ways to proceed along this line are (1) use a data reduction method (e.g., variable clustering or principal components) and put summary scores into the model instead of individual variables and (2) put all variables in the model but do not test for the effect of one variable adjusted for the effects of competing variables. For (2), chunk tests of competing variables are powerful because collinear variables join forces in the overall multiple degree of freedom association test, instead of competing against each other as when you test variables individually.

Answer 3

If you would like to carry out variable selection in the presence of high collinearity I can recommend the l0ara package, which fits L0 penalized GLMs using an iterative adaptive ridge procedure. As this method is ultimately based on ridge regularized regression, it can deal very well with collinearity, and in my simulations it produced much less false positives whilst still giving great prediction performance as well compared to e.g. LASSO, elastic net or adaptive LASSO. Alternatively, you could also try the L0Learn package with a combination of an L0 and L2 penalty. The L0 penalty then favours sparsity (ie small models) whilst the L2 penalty regularizes collinearity. Elastic net (which uses a combination of an L1 and L2 penalty) is also often suggested, but in my tests this produced way more false positives, plus the coefficients will be heavily biased. This bias you can get rid off if you use L0 penalized methods instead (aka best subset) - it's a so-called oracle estimator, that simultaneously obtains consistent and unbiased parameter coefficients. The regularization parameters in all of these methods need to be optimized via cross validation to give optimal out of sample prediction performance. If you would also like to obtain significance levels and confidence intervals on your parameters then you can also do this via nonparametric bootstrapping.

The iterative adaptive ridge algorithm of l0ara (sometimes referred to as broken adaptive ridge), like elastic net, possesses a grouping effect, which will cause it to select highly correlated variables in groups as soon as they would enter your model. This makes sense - e.g. if you had two near-collinear variables in your model it would divide the effect equally over both.

If you're analysing proportion data you are better off using a logistic regression model btw - the l0ara package allows you to do that in combination with an L0 penalty; for the L0Learn package this will be supported shortly.