# 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](https://www.python.org/dev/peps/pep-0257/) 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 `yield`s a `dict` and then use something like  [`DataFrame.from_records`][1] 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][2], 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`](https://stackoverflow.com/a/51496653/1562285)

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

  [1]: http://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.from_records.html#pandas.DataFrame.from_records
  [2]: https://docs.python.org/3/library/itertools.html#itertools-recipes