Pick the best combination of n choose k models

Question

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)

           # Append RSS to RSS Lists
           RSS_list.append(temp_results[0])

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
    RSS_list = []
    # 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 RSS to RSS Lists
                RSS_list.append(temp_results[0])

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

                # 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,
        'RSS' : RSS_list,
        'R-Squared' : R_squared_list,
        'Features' : feature_list
    })

    # Finding the Best Subsets for each number of features

    # The smallest RSS
    df_min = df[df.groupby('No_of_Features')['RSS'].transform(min) == df['RSS']]
    # 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['min_RSS'] = df.groupby('No_of_Features')['RSS'].transform(min)
    df['max_R_Squared'] = df.groupby('No_of_Features')['R-Squared'].transform(max)

This code is taken from my IPython notebook.

Answer 1

comments

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
        RSS, R_squared = fit_linear_reg(X[list(combo)], Y)
        yield {
            "No_of_Features": len(combo),
            "RSS": RSS,
            "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)),
        index=["No_of_Features", "RSS", "R-Squared", "Features"],
    )

    summary = df.groupby("No_of_Features")["R-Squared"].agg(
        {"RSS": "min", "R-Squared": "max"}
    )
    df_min = df.loc[df.groupby("No_of_Features")["RSS"].idxmin()]
    df_max = df.loc[df.groupby("No_of_Features")["R-Squared"].idxmax()]
    return df, df_min, df_max, summary

And then you can pass these results on to the plotting function