# Pick the best combination of n choose k models

## Algorithm

Let $$\M_0\$$ denote the null model which contains no predictors. This model simply predicts the sample mean of each observation.

For $$\k=1,2,\ldots,n\$$:

• Fit all $$\n \choose k\$$ models that contain exactly $$\k\$$ predictors.
• Pick the best among these $$\n \choose 𝑘\$$ models, and call it $$\M_k\$$. Here the best is defined as having the smallest RSS or equivalent measure.

Select the single best model among $$\M_0,M_1,\ldots,M_n\$$ using cross validated prediction error, $$\C_p\$$, BIC, $$\R^2_{\mathit{adj}}\$$ or any other method.

I looked up a tutorial and created a function based upon its script. It's essentially used so I can select dependent variables that's a subset of a data frame. It runs but it is very very slow.

How would I flatten a nested for-loop such as this?

# Loop over all possible combinations of k features
for k in range(1, len(X.columns) + 1):
# Looping over all possible combinations: from 11 choose k
for combo in itertools.combinations(X.columns,k):
# Store temporary results
temp_results = fit_linear_reg(X[list(combo)],Y)



I tried implementing an enumerate version but it did not work. I'm not sure, how I can flatten the nested for loop, such that I can append the results of a function to a list.

# This function takes in a subset of a dataframe representing independent
# variables (X) and a column for dependent variable (Y). This function fits
# separate models for each possible combination of the k predictors (which is
# based on the column length of X) and then select the best subset. The
# resulting output is a dataframe.

def BestSubsetSelection(X,Y):
# number of predictors
k = len(X.columns)
# Store the RSS from a linear regression model
# Store the R-square from a linear regression model
R_squared_list = []
# Store the features for a given iteration.
feature_list = []
# Store the number of features used for a given iteration. This corresponds with the feature_list.
numb_features = []

# Loop over all possible combinations of k features
for k in range(1, len(X.columns) + 1):
# Looping over all possible combinations: from 11 choose k
for combo in itertools.combinations(X.columns,k):
# Store temporary results
temp_results = fit_linear_reg(X[list(combo)],Y)

# Append R-Squared TO R-Squared list
R_squared_list.append(temp_results)

# Append Feature/s to Feature list
feature_list.append(combo)

# Append the number of features to the number of features list
numb_features.append(len(combo))

df = pd.DataFrame({
'No_of_Features': numb_features,
'R-Squared' : R_squared_list,
'Features' : feature_list
})

# Finding the Best Subsets for each number of features

# The Largest R-Squared Value
df_max = df[df.groupby('No_of_Features')['R-Squared'].transform(min) == df['R-Squared']]
display(df_min)
display(df_max)

# Adding columns to the dataframe with RSS and R-Squared values of the best subset
df['max_R_Squared'] = df.groupby('No_of_Features')['R-Squared'].transform(max)



This code is taken from my IPython notebook.

• @Peilonrayz I hope the modifications clarifies the situation. I was also told that if all combinations are to be checked then there is no way it could be linear, but does that mean it cannot be made more efficient? Jun 12, 2019 at 1:23
• Is there a good reason for using different measures to select the best between two models with the same number of predictors vs two models with a different number of predictors? It seems unnecessarily complicated to me. Jun 12, 2019 at 7:09

comments should explain why you do something, not what you do. # Append R-Squared TO R-Squared list adds nothing of value. On the contrary, it uses vertical space, and if ever you change something of the code you will need to change the coàmments as well

# This function takes in a subset of a dataframe representing independent
# variables (X) and a column for dependent variable (Y). This function fits
# separate models for each possible combination of the k predictors (which is
# based on the column length of X) and then select the best subset. The
# resulting output is a dataframe.


could be the docstring of the method if it were correct. This function does not return a DataFrame, but only prints the results

# functions

Now you have 1 monster function that: - calls the test - aggregates the values - displays the results

Better would be to split this

# getting the results

instead of having to append the results to 4 lists, I would extract this to a generator that yields a dict and then use something like DataFrame.from_records to combine this.

## powerset

what you do here:

for k in range(1, len(X.columns) + 1):
# Looping over all possible combinations: from 11 choose k
for combo in itertools.combinations(X.columns,k):


looks a lot like the powerset itertools-recipe, so let's use that one:

from itertools import chain, combinations

def powerset(iterable):
"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
s = list(iterable)
return chain.from_iterable(combinations(s, r) for r in range(len(s) + 1))

def generate_combo_comparisons(X, Y):
for combo in powerset(X.columns):
if not combo:  # the first combo is empty
continue
yield {
"No_of_Features": len(combo),
"R-Squared": R_squared,
"Features": combo,
}


To get the maximum indices of each group, you can use groupby.idxmax

You add columns 'min_RSS' to the original DataFrame. Better here would be to generate a new summary DataFrame

def subset_results(X, Y):
df = pd.DataFrame.from_records(
data=list(generate_combo_comparisons(X, Y)),