Exercise 7.8

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from patsy import dmatrix

%matplotlib inline

df = pd.read_csv('../data/Auto.csv')

# Visualize dataset
df.head()

Pre-processing data

# 'horsepower' is a string
# We would like it to be a integer, in order to appear in scatter plots.
type(df['horsepower'][0])

str

# 'horsepower' is a string because one of its values is '?'
df['horsepower'].unique()

array(['130', '165', '150', '140', '198', '220', '215', '225', '190',
       '170', '160', '95', '97', '85', '88', '46', '87', '90', '113',
       '200', '210', '193', '?', '100', '105', '175', '153', '180', '110',
       '72', '86', '70', '76', '65', '69', '60', '80', '54', '208', '155',
       '112', '92', '145', '137', '158', '167', '94', '107', '230', '49',
       '75', '91', '122', '67', '83', '78', '52', '61', '93', '148', '129',
       '96', '71', '98', '115', '53', '81', '79', '120', '152', '102',
       '108', '68', '58', '149', '89', '63', '48', '66', '139', '103',
       '125', '133', '138', '135', '142', '77', '62', '132', '84', '64',
       '74', '116', '82'], dtype=object)

# Replace '?' value by zero
df['horsepower'] = df['horsepower'].replace(to_replace = '?', value = '0')
df['horsepower'].unique()  # Check if it's ok

array(['130', '165', '150', '140', '198', '220', '215', '225', '190',
       '170', '160', '95', '97', '85', '88', '46', '87', '90', '113',
       '200', '210', '193', '0', '100', '105', '175', '153', '180', '110',
       '72', '86', '70', '76', '65', '69', '60', '80', '54', '208', '155',
       '112', '92', '145', '137', '158', '167', '94', '107', '230', '49',
       '75', '91', '122', '67', '83', '78', '52', '61', '93', '148', '129',
       '96', '71', '98', '115', '53', '81', '79', '120', '152', '102',
       '108', '68', '58', '149', '89', '63', '48', '66', '139', '103',
       '125', '133', '138', '135', '142', '77', '62', '132', '84', '64',
       '74', '116', '82'], dtype=object)

# Convert 'horsepower' to int
df['horsepower'] = df['horsepower'].astype(int)
type(df['horsepower'][0])  # Check if it's ok

numpy.int32

Relationships analysis

# Scatterplot matrix using Seaborn
sns.pairplot(df);

mpg seems to have a non-linear relationship with: displacement, horsepower and weight. Since in the book they compare mpg to horsepower, we will only analyse that relationship in this exercise.

Non-linear models

Non-linear models investigated in Chapter 7: Polynomial regression Step function Regression splines Smoothing splines Local regression Generalized additive models (GAM)

We will develop models for the first three because we didn't find functions for smoothing splines, local regression and GAM.

# Dataset
X = df['horsepower'][:,np.newaxis]
y = df['mpg']

Polynomial regression

for i in range (1,11):
    model = Pipeline([('poly', PolynomialFeatures(degree=i)), 
                      ('linear', LinearRegression())])
    model.fit(X,y)

    score = cross_val_score(model, X, y, cv=5, scoring='neg_mean_squared_error')

    print("Degree: %i  CV mean squared error: %.3f" % (i, np.mean(np.abs(score))))

Degree: 1  CV mean squared error: 32.573
Degree: 2  CV mean squared error: 30.515
Degree: 3  CV mean squared error: 27.440
Degree: 4  CV mean squared error: 25.396
Degree: 5  CV mean squared error: 24.986
Degree: 6  CV mean squared error: 24.937
Degree: 7  CV mean squared error: 24.831
Degree: 8  CV mean squared error: 24.907
Degree: 9  CV mean squared error: 25.358
Degree: 10  CV mean squared error: 25.270

The polynomial degree that best fits the relationship between mpg and horsepower is 7.

Step function

for i in range(1,11):
    groups = pd.cut(df['horsepower'], i)
    df_dummies = pd.get_dummies(groups)

    X_step = df_dummies
    y_step = df['mpg']

    model.fit(X_step, y_step)
    score = cross_val_score(model, X_step, y_step, cv=5, scoring='neg_mean_squared_error')

    print('Number of cuts: %i  CV mean squared error: %.3f' %(i, np.mean(np.abs(score))))

Number of cuts: 1  CV mean squared error: 74.402
Number of cuts: 2  CV mean squared error: 43.120
Number of cuts: 3  CV mean squared error: 36.934
Number of cuts: 4  CV mean squared error: 41.236
Number of cuts: 5  CV mean squared error: 30.880
Number of cuts: 6  CV mean squared error: 29.083
Number of cuts: 7  CV mean squared error: 29.269
Number of cuts: 8  CV mean squared error: 28.562
Number of cuts: 9  CV mean squared error: 27.776
Number of cuts: 10  CV mean squared error: 24.407

The step function that best fits the relationship between mpg and horsepower is the one with 10 steps.

Regression splines

There are two types of regression splines: splines and natural splines. The difference between these two groups is that a natural spline is a regression spline with additional boundary constraints: the natural function is required to be linear at the boundary. This additional constraint means that natural splines generally produce more stable estimates at the boundaries.

Usually, cubic splines are used. These splines are popular because most human eyes cannot detect the discontinuity at the knots. We use both cubic and natural cubic splines in this exercise.

To generate split basis representation, the package patsy is used. patsy is a Python package for describing statistical models (especially linear models, or models that have a linear component). Check the References for more details.

Cubic splines

for i in range(3,11):  # The degrees of freedom can't be less than 3 in a cubic spline
    transformed = dmatrix("bs(df.horsepower, df=%i, degree=3)" % i,
                         {"df.horsepower":df.horsepower},
                          return_type='dataframe')  # Cubic spline basis representation
    lin = LinearRegression()
    lin.fit(transformed, y)

    score = cross_val_score(lin, transformed, y, cv=10, scoring='neg_mean_squared_error')

    print('Number of degrees of freedom: %i  CV mean squared error: %.3f' %(i, np.mean(np.abs(score))))

Number of degrees of freedom: 3  CV mean squared error: 24.268
Number of degrees of freedom: 4  CV mean squared error: 22.195
Number of degrees of freedom: 5  CV mean squared error: 22.319
Number of degrees of freedom: 6  CV mean squared error: 22.028
Number of degrees of freedom: 7  CV mean squared error: 22.099
Number of degrees of freedom: 8  CV mean squared error: 22.145
Number of degrees of freedom: 9  CV mean squared error: 21.783
Number of degrees of freedom: 10  CV mean squared error: 21.965

for i in range(3,11):  # The degrees of freedom can't be less than 3 in a cubic spline
    transformed = dmatrix("cr(df.horsepower, df=%i)" % i,
                         {"df.horsepower":df.horsepower},
                          return_type='dataframe')  # Cubic spline basis representation
    lin = LinearRegression()
    lin.fit(transformed, y)

    score = cross_val_score(lin, transformed, y, cv=10, scoring='neg_mean_squared_error')

    print('Number of degrees of freedom: %i  CV mean squared error: %.3f' %(i, np.mean(np.abs(score))))

Number of degrees of freedom: 3  CV mean squared error: 27.132
Number of degrees of freedom: 4  CV mean squared error: 24.297
Number of degrees of freedom: 5  CV mean squared error: 22.317
Number of degrees of freedom: 6  CV mean squared error: 21.721
Number of degrees of freedom: 7  CV mean squared error: 21.904
Number of degrees of freedom: 8  CV mean squared error: 22.019
Number of degrees of freedom: 9  CV mean squared error: 21.926
Number of degrees of freedom: 10  CV mean squared error: 22.201