I am doing a research project as a 2nd author on a paper exploring the properties of a novel algorithm for Optimal Variable Selection where I am running the benchmark Variable Selection Methods. Each of these 3 Benchmarks has been run on the same set of 260,000 synthetic datasets with known properties. But, because this is intended for publication, the lead author has asked that I re-perform my analysis using a different statistical programming language or software package, however, I really only know R, so I told him I can just do it using a different function from a different package in R but with the same random seed beforehand.

My code for loading in the data and running my Backward Elimination & Forward Selection Stepwise Regressions (miscellaneous lines of code used for sorting, reformatting, and preprocessing the data has been omitted) is included below with the context being explained via my comments in the code:

# Extract all of the individual spreadsheet containing workbooks
# in the file folder called 'Data' which is filled
# random synthetic observations to run FS, Stepwise, and eventually EER
# on to compare the results. There are 260k spreadsheets in this folder.
folderpath <- "C:/Users/Spencer/.../datasets folder"
paths_list <- list.files(path = folderpath, full.names = T, recursive = T)

## This command reads all of the data in each of the N csv files and stores that 
## data for each of them in their own data.table (all data.tables are data.frames) 
## and stores all of N data.tables in a list object.
CL <- makeCluster(detectCores() - 2L)
clusterExport(CL, c('paths_list'))
system.time(datasets <- parLapply(cl = CL, X = paths_list, 
                                  fun = data.table::fread))

system.time(Structural_IVs <- lapply(datasets, function(j) {j[1, -1]}))
  True_Regressors <- lapply(Structural_IVs, function(i) 
    { names(i)[i == 1] }))

### Run a Backward Elimination Stepwise Regression
### on each of the 260,000 datasets using the step function.
set.seed(11)      # for reproducibility
system.time( BE.fits <- parLapply(cl = CL, X = datasets, \(X) {
  full_models <- lm(X$Y ~ ., X)
  back <- stats::step(full_models, scope = formula(full_models), 
                      direction = 'back', trace = FALSE) }) )

# extract the coefficients and their corresponding variable names
BE_Coeffs <- lapply(seq_along(BE.fits), function(i) coef(BE.fits[[i]]))

# extract the names of all IVs selected by them without their intercepts
IVs_Selected_by_BE <- lapply(seq_along(BE.fits), 
                             \(i) names(coef(BE.fits[[i]])[-1]))

And from there, to make a long story short, I just compared what is returned by True_Regressors for each dataset to whatever is returned by IVs_Selected_by_BE for that same dataset, how I do this will presumably be the same or very similar. What I need is a good suggestion for another package and function to use to do all of this over again as a great big sanity check!

For completeness, my complete code used to run all the Forward Selection Stepwise Regressions is included below as well:

### Run a Forward Selection Stepwise Regression
### function on each of the 260,000 datasets.
set.seed(11)      # for reproducibility
system.time( FS.fits <- parLapply(cl = CL, X = datasets, \(X) {
  nulls <- lm(X$Y ~ 1, X)
  full_models <- lm(X$Y ~ ., X)
  forward <- stats::step(object = nulls, direction = 'forward',
                         scope = formula(full_models), trace = FALSE) }) )

# extract the coefficients and their corresponding variable names
FS_Coeffs <- lapply(seq_along(FS.fits), function(i) coef(FS.fits[[i]]))

# assign all regressors selected by Forward Stepwise Regression,
# not including the Intercepts, to IVs_Selected_by_FS
IVs_Selected_by_FS <- lapply(seq_along(FS.fits), 
                             \(i) names(coef(FS.fits[[i]])[-1]))

I have already attempted to do this using the ols_step_backward_aic function from the olsrr package, but I cannot get it to run for the life and me and it is continuously throwing up errors that don't make much sense. So, now I am looking for an alternative alternative way of running Stepwise Regressions in R.

  • 1
    $\begingroup$ I am not too experienced with regression in R, but based on this article it seems the MASS, leaps, and caret packages have functions that can be used to perform stepwise regression. $\endgroup$
    – Oxbowerce
    Commented Jan 19, 2023 at 10:20
  • $\begingroup$ @Oxbowerce the caret package has functions or arguments to run pretty much any type of regression or classification model imaginable, but it has cross-validation and/or resampling built into it, so I cannot use it for this project unfortunately. But thank you for the other suggestions, I'll look into them for sure!! $\endgroup$
    – Marlen
    Commented Jan 19, 2023 at 13:13

1 Answer 1


Based on the comment from Oxbowerce, I looked into using the MASS library instead of olsrr, and fortunately, the inherent similarity in the syntax of the MASS::stepAIC() function and the stats::step() functions are quite high actually, at least in this application!

Here is the version of my backward elimination stepwise regression function estimated using this new function instead:

### Step 2: Run a Backward Elimination Stepwise Regression
### function on each of the 260,000 datasets.
set.seed(11)      # for reproducibility
system.time( BE.fits <- lapply(X = datasets, \(ds_i) {
  null_models <- lm(ds_i$Y ~ 1, ds_i)
  full_models <- lm(ds_i$Y ~ ., data = ds_i)
  back <- MASS::stepAIC(full_models, direction = 'backward', 
                        scope = list(upper = full_models, lower = null_models), 
                        trace = FALSE) }) ) 

As you can see, the only substantive syntactical difference is in the scope argument (because using 'backward' also works for the step function, you don't have to use 'back' as I did here) where it changes from using formula to list and only requiring the upper argument specified to requiring both as my original forward selection stepwise regression estimation code using the stat function does.

The neatest part about using this scheme to reproduce my results is that I can reuse the following lines of code which isolate and save the coefficient estimates on and then the names of the Variables selected as is, no adjustments whatsoever are necessary!

p.s. One important note for anyone trying to perform a similar reproduction of prior results when you originally used the stats::step() function is that you must use an alternative which also uses the AIC as its selection criterion because by default, the stats::step() function uses this criterion when selecting optimal candidate factors.


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.