Is this a legit way to make a variable/predictor/dummy selection?

(My goal is forecasting with the selected variables)

fit <- train(train.values ~ .,data=train.data, method='glmnet') # train.data includes all variables

#getting the coefficients of the final model
coefficients <- coef(fit$finalModel, fit$bestTune$lambda)

#create a list of the selected coefficients
variables <- names(coefficients[which(Coefficients != 0),])

Due to reading lots of stuff on this platform i am aware of the fact that stepAIC() is not that great of a choice for a variable selection. After the variable selection is done i would use those variables for predictions with glmnet as well as for a linear model.

  • 1
    $\begingroup$ Doesn't caret use CV as it's default for model selection ? It certainly is a choice. If one is using CV, where does AIC enter the picture ? $\endgroup$
    – meh
    Oct 10, 2015 at 15:25
  • $\begingroup$ As far as i understood here link using the AIC is not that optimal. I am looking for an automated decision to choose my variables. How would you do it $\endgroup$
    – RandomDude
    Oct 11, 2015 at 8:18

2 Answers 2


What you're trying to do here is to identify most "important" variables within glmnet and then trying to pass your features to another model. This is not optimal as Max Kuhn writes here:

In many cases, using these models with built-in feature selection will be more efficient than algorithms where the search routine for the right predictors is external to the model. Built-in feature selection typically couples the predictor search algorithm with the parameter estimation and are usually optimized with a single objective function (e.g. error rates or likelihood).

From theoretical perspective: different models allow for different degree of flexibility and protection against overfit. Why would you use set of optimal parameters from one model in another one?

The only scenario I can think of using glmnet as a feature filter if you are forced to use linear regression and you have a wide matrix of features (more features than cases). There is nothing wrong with this, simply the results (RMSE, R^2) may be suboptimal.

  • 1
    $\begingroup$ From a practical perspective you would agree that i will identify the "important" variables with this method? I want to compare the performance of different models (lm, glmnet, rlm, auto.arima()) for a given dataset. Since i have lots of variables i somehow need to separate "important" from "unimportant". $\endgroup$
    – RandomDude
    Oct 11, 2015 at 10:47
  • 1
    $\begingroup$ Yes, you would. I tried several times prefiltering list of features for most "important" -- with glmnet (as you did !=0), svm with regularization (Python), and random forest (most important) -- and then passing this variables to another model: all the time the results were inferior to having selected variables with built-in feature selection. My point -- in places where you compete for a result, like Kaggle eg, that won't do. If you decide to do it, eg for research purpose, it's your choice. But why you cant do lm on all the variables??? Multicollinearity, wide matrix? $\endgroup$ Oct 11, 2015 at 11:08
  • $\begingroup$ I am dealing with time series data where i am not sure which events (mostly moving holidays including lead and lag) have an effect on the data or not. So i am putting all possible holidays in and see what comes out. Identifying the effects before is not that easy due to multi-seasonality (daily data, several years). I am also using lm on all variables - but its always 'worse' than the one with variable selection. + for sure calculation of models is way faster with less variables (especially auto.arima). Any thoughts about this method? $\endgroup$
    – RandomDude
    Oct 11, 2015 at 11:20
  • $\begingroup$ It is for sure not only about the moving holidays, but also which day of the month, month of the year have an importance - and which not... I know that the whole method is from a theoretical perspective far away from perfect, but the results on the different datasets i have to forecast is "decent". $\endgroup$
    – RandomDude
    Oct 11, 2015 at 11:40
  • 1
    $\begingroup$ I do not have much experience with time series, but from what I know, in your case I would invest in ARIMA: with holidays and other special events as dummies, with perhaps couple of days around, plus multiseasonality, perhaps via fourier transform (as discussed in another ts post where you commented). In any case, I would pay close attention to residual diagnostics to see if something is left there in terms of unexplained patterns (spikes?) or autoregression. $\endgroup$ Oct 11, 2015 at 11:45

I have read academic papers citing the effectiveness of using Lasso for variable selection as well as actually putting it into practice myself.

The following code block identifies features from your data set.

##returns variables from lasso variable selection, use alpha=0 for ridge
  x=x[,-1] ##remove intercept


  co<-coef(glmnet1,s = "lambda.1se")
  variables<-variables[!(variables %in% '(Intercept)')];
  return( c(yvar,variables));

(I cannot take 100% credit for this code as I am sure it is adapted from some place - most likely here: Using LASSO from lars (or glmnet) package in R for variable selection )

While on the topic of variable selection, I have also found that VIF (variable inflation factor) to be effective especially when cross-validated.

#returns selected variables using VIF and kfolds cross validation 
  for(v in names(df)){
  for(i in 1:folds){   
    if(trace) message("fold ",i);
    rows=f$subsets[f$which!=i] ##leave one out 
    xdf=df[rows,names(df) != yvar]; #remove output var    
    if(trace) say("trying ",i,yvar,nrow(df),length(y)," subsize=",min(200,floor(nrow(xdf))));
    if(trace) print(names(xdf)[vifResult$select]);
    for(v in names(xdf)[vifResult$select]){
      findings[[v]]=findings[[v]]+1; #vote
  findings=(sort(unlist(findings),decreasing = T))    
  if(trace) print(findings[findings>0]); 
  return( c(yvar,names(findings[findings==findings[1]])) )  

Both of the above ez functions return an vector of variable names. The following code block converts the return values to a formula.

#converts ezvif or ezlasso results into formula
ezformula=function(v,operator=' + '){
  return(as.formula(paste(v[1],'~',paste(v[-1],collapse = operator))))

I hope this is helpful.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.