<jupyter_start><jupyter_text>Exercise 7.8    Fit some of the non-linear models investigated in this chapter to the Auto data set. Is there evidence for non-linear relationships in this data set? Create some informative plots to justify your answer.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

%run ../../customModules/usefulFunctions.ipynb
# https://stackoverflow.com/questions/34398054/ipython-notebook-cell-multiple-outputs
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

from itertools import permutations
from sklearn.model_selection import LeaveOneOut
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from pygam import LinearGAM, s
df = pd.read_csv("../../DataSets/Auto/Auto.csv")
df['horsepower'] = pd.to_numeric(df['horsepower'], errors='coerce')
df = df.set_index('name')
df.dropna(inplace=True)
df.head()
_ = pd.plotting.scatter_matrix(df, figsize=(12, 12))<jupyter_output><empty_output><jupyter_text>Based on these plots we will further investigate the $\mathrm{displacement}$, $\mathrm{horsepower}$ and $\mathrm{weight}$ variables.<jupyter_code>df_y = df[['mpg']].copy()
df_x = df[['displacement', 'horsepower', 'weight']].copy()<jupyter_output><empty_output><jupyter_text>Let us first test the polynomial models.<jupyter_code>polynomialRegressionChooseOptimalDegree(df_x[['displacement']], df_y, total_degrees=5)
polynomialRegressionChooseOptimalDegree(df_x[['horsepower']], df_y, total_degrees=5)
polynomialRegressionChooseOptimalDegree(df_x[['weight']], df_y, total_degrees=5)
total_degrees = 2
columns = df_x.columns
for independent in columns:
    for i in range(2, total_degrees + 1):
        variable_name = independent + str(i)
        df_x[variable_name] = df_x[independent]**i

df_x = df_x.reindex(sorted(df_x.columns), axis=1)
df_x.head()
loocv = LeaveOneOut() # LeaveOneOut
MSE = 0
for train_index, test_index in loocv.split(df_x):
    df_x_train, df_x_test = df_x.iloc[train_index], df_x.iloc[test_index]
    df_y_train, df_y_test = df_y.iloc[train_index], df_y.iloc[test_index]

    model = LinearRegression()
    _ = model.fit(df_x_train, df_y_train)

    Y_pred = model.predict(df_x_test)
    MSE += (df_y_test.iloc[0, 0] - Y_pred[0, 0])**2        

print(f'MSE = {MSE / df_x.shape[0]}')<jupyter_output>MSE = 15.547346629694578<jupyter_text>Now we will try the cubic splines with various amounts of knots.<jupyter_code>df_x = df_x[['displacement', 'horsepower', 'weight']]
df_x_train, df_x_test, df_y_train, df_y_test = train_test_split(df_x, df_y, test_size=0.1, random_state = 1)
plotWrapperCubicSplines(df_x[['displacement']], df_x_train[['displacement']], df_x_test[['displacement']], df_y, df_y_train, df_y_test, natural=False)
plotWrapperCubicSplines(df_x[['displacement']], df_x_train[['displacement']], df_x_test[['displacement']], df_y, df_y_train, df_y_test, natural=True)
plotWrapperCubicSplines(df_x[['horsepower']], df_x_train[['horsepower']], df_x_test[['horsepower']], df_y, df_y_train, df_y_test, natural=False)
plotWrapperCubicSplines(df_x[['horsepower']], df_x_train[['horsepower']], df_x_test[['horsepower']], df_y, df_y_train, df_y_test, natural=True)
plotWrapperCubicSplines(df_x[['weight']], df_x_train[['weight']], df_x_test[['weight']], df_y, df_y_train, df_y_test, natural=False)
plotWrapperCubicSplines(df_x[['weight']], df_x_train[['weight']], df_x_test[['weight']], df_y, df_y_train, df_y_test, natural=True)<jupyter_output><empty_output><jupyter_text>Let us now extend the cubic splines as a GAM.<jupyter_code>results = {}
perm_splines = permutations(list(range(4, 10)), 3) 
for n_splines in perm_splines:
    # https://codeburst.io/pygam-getting-started-with-generalized-additive-models-in-python-457df5b4705f
    # fit a spline term on feature 0, 1 and 2
    gam = LinearGAM(s(0) + s(1) + s(2), n_splines=n_splines, fit_intercept=True,
                    lam=0, constraints=None).fit(df_x_train.to_numpy(), df_y_train.to_numpy())
    Y_pred = gam.predict(df_x_test)
    MSE = mean_squared_error(df_y_test, Y_pred)
    n_knots = tuple(splines-4 for splines in n_splines)
    results[str(n_knots)] = MSE


n_knots_arr = np.array(list(results.keys()))
MSE_arr = np.array(list(results.values()))

# plot best 'k' models
k = 10
idx = np.argpartition(MSE_arr, k)
n_knots_arr = n_knots_arr[idx[:k]]
MSE_arr = MSE_arr[idx[:k]]
std_dev = np.std(MSE_arr)

id_min = np.argmin(MSE_arr)
min_mse = n_knots_arr[id_min]  # we plot these 'best results' with a different marker
min_val = MSE_arr[id_min]  # we plot these 'best results' with a different marker
n_knots_arr = np.delete(n_knots_arr, id_min)
MSE_arr = np.delete(MSE_arr, id_min)

fig, ax1 = plt_func.subplots(1, 1, constrained_layout=True, figsize=(8, 4))
_ = ax1.scatter(n_knots_arr, MSE_arr)
_ = ax1.scatter(min_mse, min_val, marker='x', s=40, c=plt_func.rcParams['axes.prop_cycle'].by_key()['color'][0], label='smallest MSE')
_ = ax1.set_xlabel(r'knots')
_ = ax1.set_ylabel('Cross-validation MSE')
_ = ax1.legend()
# partial dependency plots of best model with full dataset
n_knots = [0, 1, 3]
n_splines = [el + 4 for el in n_knots]
gam = LinearGAM(s(0) + s(1) + s(2), n_splines=n_splines, fit_intercept=True,
                    lam=0, constraints=None).fit(df_x.to_numpy(), df_y.to_numpy())
columns = df_x.columns
numberOfSubplots = df_x.shape[1]
confidence_width = 0.95
fig, axes = plt_func.subplots(numberOfSubplots, 1, constrained_layout=True, figsize=(8, 4*numberOfSubplots))
for i, ax in enumerate(axes):
    X_pred = gam.generate_X_grid(term=i)
    Y_pred = gam.predict(X_pred)
    
    partial_dependence, confidence_intervals = gam.partial_dependence(term=i, X=X_pred, width=confidence_width)
    lower_conf_int = confidence_intervals[:, 0]
    upper_conf_int = confidence_intervals[:, 1]

    diff = Y_pred[0] - partial_dependence[0]
    partial_dependence += diff
    lower_conf_int += diff
    upper_conf_int += diff
    
#     _ = ax.scatter(df_x[columns[i]], df_y, label='training data')
    _ = ax.plot(X_pred[:, i], partial_dependence, '-', color='darkorchid', linewidth=2, label='prediction')
    _ = ax.fill_between(X_pred[:, i], lower_conf_int, upper_conf_int, color='#888888', alpha=0.4, label=f"confidence interval ({int(confidence_width*100)}%)")
    
    prediction_intervals = gam.prediction_intervals(X_pred, width=confidence_width)
    lower_pred_int = prediction_intervals[:, 0]
    upper_pred_int = prediction_intervals[:, 1]
    _ = ax.fill_between(X_pred[:, i], lower_pred_int, upper_pred_int, color='#888888', alpha=0.1, label=f"prediction interval ({int(confidence_width*100)}%)")
    
    _ = ax.legend()
    _ = ax.set_xlabel(columns[i])
    _ = ax.set_ylabel(r'$f_{%s}(\mathrm{%s})$' % (i+1, columns[i]))
    _ = ax.set_title(f'partial dependency plot of {columns[i]}')<jupyter_output><empty_output><jupyter_text>We can also fit a GAM using smoothing splines.<jupyter_code>results = {}
perm_lambda = permutations(np.linspace(0, 5, 21, endpoint=True), 3) 
for lambda_ in perm_lambda:
    # https://codeburst.io/pygam-getting-started-with-generalized-additive-models-in-python-457df5b4705f
    # fit a spline term on feature 0, 1 and 2
    gam = LinearGAM(s(0) + s(1) + s(2), n_splines=4, fit_intercept=True,
                    lam=lambda_, constraints=None).fit(df_x_train.to_numpy(), df_y_train.to_numpy())
    Y_pred = gam.predict(df_x_test)
    MSE = mean_squared_error(df_y_test, Y_pred)
    results[str(lambda_)] = MSE

lambda_arr = np.array(list(results.keys()))
MSE_arr = np.array(list(results.values()))

# plot best 'k' models
k = 5
idx = np.argpartition(MSE_arr, k)
lambda_arr = lambda_arr[idx[:k]]
MSE_arr = MSE_arr[idx[:k]]
std_dev = np.std(MSE_arr)

id_min = np.argmin(MSE_arr)
min_mse = lambda_arr[id_min]  # we plot these 'best results' with a different marker
min_val = MSE_arr[id_min]  # we plot these 'best results' with a different marker
lambda_arr = np.delete(lambda_arr, id_min)
MSE_arr = np.delete(MSE_arr, id_min)

fig, ax1 = plt_func.subplots(1, 1, constrained_layout=True, figsize=(8, 4))
_ = ax1.scatter(lambda_arr, MSE_arr)
_ = ax1.scatter(min_mse, min_val, marker='x', s=40, c=plt_func.rcParams['axes.prop_cycle'].by_key()['color'][0], label='smallest MSE')
_ = ax1.set_xlabel(r'$\lambda$')
_ = ax1.set_ylabel('Cross-validation MSE')
_ = ax1.legend()<jupyter_output><empty_output><jupyter_text>Finally, we use pyGAM's gridsearch() method to search for the optimal smoothing parameters.<jupyter_code>gam = LinearGAM().gridsearch(df_x_train.to_numpy(), df_y_train.to_numpy())
Y_pred = gam.predict(df_x_test)
MSE = mean_squared_error(df_y_test, Y_pred)
print(f'MSE = {MSE}\n')
gam.summary()<jupyter_output>100% (11 of 11) |########################| Elapsed Time: 0:00:00 Time:  0:00:00
{'fraction_non_alphanumeric': 0.14338773213280812, 'fraction_numerical': 0.020146314012380418, 'mean_word_length': 6.069212410501193, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '384936', 'license': 'MIT', 'n_tokens_mistral': 3607, 'n_tokens_neox': 3426, 'n_words': 630, 'path': 'Chapter07/Exercise08/7_8.ipynb', 'repo_name': 'HerrBroekman/Applied-solution-of-An-Introduction-to-Statistical-Learning-written-in-Python', 'size': '8885'}