Pick the best combination of n choose k models
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]) --&gt; () (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