Quantifying the performance of Stepwise Regression ran on Monte Carlo generated datasets & comparing them to your method of interest

Question

The source data files and scripts referenced here and from whom lines of code are included here can be found in my GitHub Repository for this collaborative research project exploring the properties of a newly proposes Statistical Learning Algorithm for basic Optimal Variable Selection purposes. All code in this question comes from the "Benchmarks 2 & 3 - BE & FS Stepwise Regression (manually on a single dataset)" Rscript and the data loaded into that script is the one in it called "0-11-3-462", which can be found in the "last 40" folder in my aforementioned GitHub Repo for this project.

All of the synthetic datasets used in this project were generated via Monte Carlo Simulation by way of a custom Microsoft Excel Macro in such a way that each of them contains 30 candidate regressor columns (X1:X30), with 500 observations on each; and 500 such observations on a dependent variable column Y. Plus, three rows above, one header with each column name as just specified, and two more rows above that which together indicate the true underlying structural regression equation parameters characterizing that dataset.

The goal here is to (in RStudio) run my 2nd Benchmark Comparison Algorithm, namely Backward Stepwise Regression, on one of the datasets, namely, the one called '0-11-3-462' in my "top 40" file-folder, and write a function which quantifies how many of the variables selected by my Stepwise Regression on that dataset are correct.

After transforming, manipulating, and pre-processing data in R so that all the observations on each of the 30 candidate regressors stored in a 500 by 30 matrix called X_obs_matrix, and all 500 observations on the DV in a vector called Y_obs_vector, I the run my Regression using the following code:

set.seed(11)      # for reproducibility
full_model <- lm(formula = Y ~ ., data = All_sample_obs)
BE.fit <- step(full_model, direction = "backward", trace = 0)

Then store its results (the candidate VSs selected by it) using:

BE_Coeffs <- coef(BE.fit)

IVs_Selected_by_BE <- names(BE_Coeffs[-1])

Where the results for this dataset are:

> IVs_Selected_by_BE
 [1] "X2"  "X3"  "X4"  "X5"  "X8"  "X10" "X13" "X15" "X17" "X19" "X20"
[12] "X22" "X25" "X26" "X27" "X28" "X30"

And the true underlying (structural) parameters are by:

> True_Regressors
  X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 X13 X14 X15 X16 X17 X18 X19 X20
1  0  0  0  1  0  0  0  0  0   1   0   0   1   0   0   0   0   0   1   1
  X21 X22 X23 X24 X25 X26 X27 X28 X29 X30
1   0   1   0   0   1   1   1   1   0   1

Note: The 1s at the start of each line do not indicate included candidate variables/regressors, only the 1s and 0s directly beneath the last character in each candidate regressors name printed out directly above it.

So the 1 above the X4 in row 3 and the 4 in row 2 indicates that the X4 regressor is a true population parameter; along with X10, X13, X19, X20, X22, X25, X26, X27, X28, & X30 as well.

Each of these 260k datasets are named in the following format: n1-n2-n3-n4, here, what n1, n2, and n4 mean is not relevant, but, n3 tells you how many parameters are included in the true underlying population regression specification. So, one possible model performance metric generating function could be a function which returns the # of correctly selected regressor candidates by Backward Stepwise Regression (which we'll call TSRs) out of/divided by the # of regressors included in the structural model (which we'll call SRs), which means it should return TSR/SR for each dataset.

I need a create a function (or functions) which calculates the Sensitivity (aka the TPR) for each candidate regressor, and also for the overall selected regression specification.

p.s. This function will be used to assess and compare the accuracy of Backward Elimination Stepwise Regressions, that are trained on 260k different Monte Carlo generated random (csv file-formatted) datasets. Their performance will then be compared with that of at least 2 other benchmark variable selection methods (only LASSO Regression and FS Stepwise Regression for now) as well as the new Variable Selection algorithm being studied in this collaborative machine/statistical learning research project known as Estimated Exhaustive Regression.

Answer 1

The first step is to quantify the total number of 'Positives', i.e., the total number of structural factors explaining each dataset, and this is straight-forward since you already have that number stored implicitly in your True_Regressors object. All you will need to do is use:

num_of_Positives <- lapply(True_Regressors, function(i) { length(i) })

Then, you need write a function which determines and returns you the names of each of the 'True Positives', that is, the names of all the candidate regressors/independent variables your stepwise regression selected which are actually structural variables. One way to do that would be this:

True_Positives <- lapply(seq_along(All_sample_obs), \(k)
                              sum(IVs_Selected_by_BE[[k]] %in% 
                                    True_Regressors[[k]]))

And because the definition of the True Positive Rate is the number of True Positives over the total number of all Positives, you can calculate this using:

TPRs = lapply(seq_along(All_sample_obs), \(j)
                  j <- (BE_TPs[[j]]/num_of_Positives[[j]]))