<jupyter_start><jupyter_text>MLR Predictive Model SelectionSo now we know how to build a wide array of linear regression models. For instance if our data set contained $n$ observations of $m$ features we could build $2^m$ models without even considering interactions, polynomials, or other nonlinear transformations of the data. That's a lot of models to choose from, so in this notebook we'll introduce how you might go about selecting the best multiple linear regression model. What We'll Accomplish in This NotebookIn this notebook we'll:- Have a discussion of generalization error in predictive models- Introduce the concept of cross-validation- Use best subset selection for the Advertising data set- Practice building the best predictive model, using the tools we've learned so farLet's go!<jupyter_code># import the packages we'll use
## For data handling
import pandas as pd
import numpy as np

# We'll use this later
from numpy import meshgrid

## For plotting
import matplotlib.pyplot as plt
import seaborn as sns


## This sets the plot style
## to have a grid on a white background
sns.set_style("whitegrid")<jupyter_output><empty_output><jupyter_text>Generalization ErrorAs we've discussed before we can't just choose the model that performs the best on our training data because we could just arbitrarily make a model that fits each training point to its target value.We may also be tempted to see which model performs best on the test set we made, but this has a similar problem. We'll just be producing models that perform really well on the test set.We want to get a sense of how well our model will perform on any new random data pull, this is known as the generalization error.What can we do to get a sense for the generalization error of our model? Cross-ValidationEnter $k$-fold cross-validation.The idea behind this technique is pretty clever. We'll introduce the technique in a way that's generalizes beyond linear regression, which will come in handy down the line.When building a predictive model, the model estimate (what we've been calling $\hat{f}$) is found from minimizing a loss function. For linear regression this loss function was the MSE of the training data. Let's consider the case where we will randomly draw a new set of data (think test set) and see how well our estimate performs. Because the data was randomly drawn, the algorithm we use is a random process. So the value of the loss function (the generalization error) on this new data is an example of a random variable, let's call this random variable $G$.It would be nice to know something about the distribution of $G$, but this is tricky with a finite amount of data. However, we can leverage a popular theorem from probability theory called the law of large numbers (see the probability theory and statistics cheat sheet).If we were able to generate a bunch of random draws of $G$, say $k$ random draws, then the law of large numbers says that:$$\frac{1}{k}\sum_{i=1}^k G_i \approx E(G),$$assuming our random draws were independent.So in $k$-fold cross-validation we take our training set, and randomly split it into $k$ equally sized chunks. For each chunk we train the algortihm on the $k-1$ other chunks and then calculate the testing loss using the chunk we left out. Then we take the arithmetic mean of all $k$ testing errors. This is an approximation of the expected value of the true generalization error of the algorithm.Here's a picture to help illustrate the idea for $5$-fold cross-validation:Let's see how this works in `sklearn` to help us choose the best model for the `Advertising` data.<jupyter_code>ads = pd.read_csv("Advertising.csv")
ads_copy = ads.copy()

ads_train = ads_copy.sample(frac=.75, random_state=614)
ads_test = ads_copy.drop(ads_train.index)

## remember in notebook 3 we found these
## to be useful
ads_train['sqrt_TV'] = np.sqrt(ads_train['TV'])
ads_train['sqrtTV_radio'] = ads_train['sqrt_TV']*ads_train['radio']
## import the KFold object from sklearn
from sklearn.model_selection import KFold

## We'll need this when we fit models
from sklearn.base import clone
## Now we make a kfold object
## we'll use 5 splits
## and shuffle the data before making the splits
kfold = KFold(n_splits = 5, shuffle = True, random_state = 440)
## To make this simpler I make my data
## into an array
X = np.array(ads_train[['TV','radio']])
y = np.array(ads_train['sales'])

# You can loop through all the splits like so
for train_index, test_index in kfold.split(X,y):
    print("TRAIN:", train_index, "TEST:", test_index)
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]
## Now let's put it all together.
## It will be easier for us to make a couple functions
## to loop through all 8 possible models

## This gets our data for us
def get_X_y(df,features,target):
    # Returns X then y
    return np.array(df[features]), np.array(df[target])

## this calculates the mse
def get_mse(model, X, y):
    # get the prediction
    pred = model.predict(X)
    
    # Returns the mse
    return np.sum(np.power(pred-y,2))/len(y)
# This function was modified from stackexchange user hughdbrown 
# at this link, 
# https://stackoverflow.com/questions/1482308/how-to-get-all-subsets-of-a-set-powerset

# This returns the power set of a set minus the empty set
def powerset_no_empty(s):
    power_set = []
    x = len(s)
    for i in range(1 << x):
        power_set.append([s[j] for j in range(x) if (i & (1 << j))])
            
    return power_set[1:]

powerset_no_empty(['TV','radio','newspaper','sqrt_TV','sqrtTV_radio'])
possible_features = powerset_no_empty(['TV','radio','newspaper','sqrt_TV','sqrtTV_radio'])
## Now make an array that will hold the mses
## for all the models
## the columns represent each possible model
MSEs = np.empty((5,len(possible_features)))
from sklearn.linear_model import LinearRegression
## Make a regression model
reg = LinearRegression(copy_X = True)
## keep track of what split we're on
i = 0

## This is for the initial input into the kfold object
X, y = get_X_y(ads_train, possible_features[0], 'sales')

## Perform CV
for train_index, test_index in kfold.split(X):
    ## For each possible model
    for j in range(len(possible_features)):
        ## get X and y
        X, y = get_X_y(ads_train, possible_features[j], 'sales')


        # Get the cv train test split
        X_train, X_test = X[train_index], X[test_index]
        y_train, y_test = y[train_index], y[test_index]
        
        # Cloning the regression makes a fresh regression 
        # model for each run
        clone_reg = clone(reg)
        
        # fit the model
        clone_reg.fit(X_train,y_train)
        
        MSEs[i,j] = get_mse(clone_reg, X_test, y_test)
    
    ## We'll now move to the next split
    i = i + 1
# Here are the MSEs
MSEs
## We can get the mean MSE for each model 
## across the CV splits like so
np.mean(MSEs, axis=0)
## We can get the mean MSEs using np.mean
## axis = 0 tells np.mean to take the column mean
## we can get where the min occurs with argmin
np.argmin(np.mean(MSEs, axis = 0))
print("The model with the lowest mean CV MSE",
     "was the one with", possible_features[np.argmin(np.mean(MSEs, axis = 0))],
     "as the features. This model had a mean CV MSE of",
     np.round(np.min(np.mean(MSEs, axis=0)),5))
## Another popular measure is root mse
## this is because it has the same dimension as 
## the target variable

## It can be interpreted as how far off we were from the
## true value on average
np.round(np.min(np.sqrt(np.mean(MSEs, axis = 0))),5)<jupyter_output><empty_output><jupyter_text>Here we brute forced our way through all the models because we had a small number of predictors and a decent sample size. This really doesn't work when your data has a large number of predictors (there are too many models to check, and you need a large quantity of samples to fit models with a lot of predictors this is known as the curse of dimensionality), or if you have a small sample size (makes it difficult to split your data even further). We'll learn other techniques for model selection in those cases later on. Also in many cases it doesn't make sense to even include a predictor in the model, for example because it has no association with your target.There are two other algorithms we touch on in the HW called backwards and forwards selection that are greedy algorithms. The benefit of these approaches is that they run more quickly than brute force, the problem is that they might not give you the "best model".For now we'll stick to examining correlations, scatter plots to provide plausible features, and then use cross validation to pick the best subset of those plausible features.  You CodeGo ahead and try to build the best model you can to predict `carseats` `Sales`. Everyone's model may be different so do your best :)<jupyter_code>carseats = pd.read_csv("carseats.csv")
## Code here









## Code here









## Code here<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup