Moving Beyond Linearity: Applied Exercises

import pandas as pd
import numpy as np

wage = pd.read_csv("data/Wage.csv")
wage.head()

from sklearn import preprocessing
from sklearn.model_selection import LeaveOneOut
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

np.random.seed(5)

def polynomial_regression(X_train, Y_train, X_test, Y_test, M):
    validation_MSE = {}
    test_MSE = {}

    for m in M:
        poly = preprocessing.PolynomialFeatures(degree=m)
        X_ = poly.fit_transform(X_train)

        mse = 0
        loo = LeaveOneOut() # Leave one out cross-validation
        for train_index, test_index in loo.split(X_):
            X, X_CV = X_[train_index], X_[test_index]
            Y, Y_CV = Y_train[train_index], Y_train[test_index]
            # Linear Regression (including higher order predictors)
            model = LinearRegression(fit_intercept=True)
            model.fit(X, Y)
            p = model.predict(X_CV)
            mse += mean_squared_error(p, Y_CV)
        validation_MSE[m] = mse/len(X_)

        # Compute test MSE for the model
        model = LinearRegression(fit_intercept=True)
        model.fit(X_, Y_train)
        p = model.predict(poly.fit_transform(X_test))
        test_MSE[m] = mean_squared_error(p, Y_test)

    # Plot validation MSE
    lists = sorted(validation_MSE.items()) # sorted by key, return a list of tuples
    x, y = zip(*lists) # unpack a list of pairs into two tuples
    fig = plt.figure(figsize=(15, 8))
    ax = fig.add_subplot(121)
    plt.plot(x, y, color='r')
    plt.grid()
    ax.set_xlabel('Degree of Polynomial')
    ax.set_ylabel('Validation MSE')
    ax.set_title('Validation MSE vs Degree of Polynomial')

    lists = sorted(test_MSE.items()) # sorted by key, return a list of tuples
    x, y = zip(*lists) # unpack a list of pairs into two tuples
    ax = fig.add_subplot(122)
    plt.plot(x, y, color='g')
    plt.grid()
    ax.set_xlabel('Degree of Polynomial')
    ax.set_ylabel('Test MSE')
    ax.set_title('Test MSE vs Degree of Polynomial')
    plt.show()


X_train, X_test, y_train, y_test = train_test_split(wage[['age']], wage[['wage']], test_size=0.1)
polynomial_regression(X_train, y_train.values, X_test, y_test, [1,2,3,4,5,6,8,9,10,11,12,13,14,15])

import statsmodels.api as sm

poly = preprocessing.PolynomialFeatures(degree=1)
X_ = poly.fit_transform(wage[['age']])
model1 = sm.OLS(wage[['wage']], X_)
model1 = model1.fit()

poly = preprocessing.PolynomialFeatures(degree=2)
X_ = poly.fit_transform(wage[['age']])
model2 = sm.OLS(wage[['wage']], X_)
model2 = model2.fit()

poly = preprocessing.PolynomialFeatures(degree=3)
X_ = poly.fit_transform(wage[['age']])
model3 = sm.OLS(wage[['wage']], X_)
model3 = model3.fit()

poly = preprocessing.PolynomialFeatures(degree=4)
X_ = poly.fit_transform(wage[['age']])
model4 = sm.OLS(wage[['wage']], X_)
model4 = model4.fit()

poly = preprocessing.PolynomialFeatures(degree=5)
X_ = poly.fit_transform(wage[['age']])
model5 = sm.OLS(wage[['wage']], X_)
model5 = model5.fit()

sm.stats.anova_lm(model1, model2, model3, model4, model5)

boston = pd.read_csv("data/Boston.csv")
boston.dropna(inplace=True)

poly = preprocessing.PolynomialFeatures(degree=3)
X_ = poly.fit_transform(boston[['dis']])
model = sm.OLS(boston[['nox']], X_)
result = model.fit()
print(result.summary())

x = np.linspace(boston['dis'].min()-1, boston['dis'].max()+1, 256, endpoint = True)
y = result.params.const + result.params.x1*x + result.params.x2*(x**2) + result.params.x3*(x**3)

fig = plt.figure(figsize=(15, 8))
plt.plot(x, y, 'red', label="Cubic Fit")
plt.scatter(boston['dis'], boston['nox'], alpha=1, color='green')

plt.xlabel('dis')
plt.ylabel('nox')
plt.title('Cubic Fit')
plt.legend(loc='best')

plt.grid()
plt.show()

def add_polt(ax, x, y, RSS, degree):
    ax.plot(x, y, 'red', label="Degree: " +str(degree))
    ax.scatter(boston['dis'], boston['nox'], alpha=1, color='green')
    ax.set_ylim(boston['nox'].min()-0.1, boston['nox'].max()+0.1)
    ax.set_xlabel('dis')
    ax.set_ylabel('nox')
    ax.set_title('Degree: ' +str(degree) + ' Residuals: ' + str(RSS))
    ax.legend(loc='best')
    ax.grid()

fig = plt.figure(figsize=(15, 40))
x = np.linspace(boston['dis'].min()-1, boston['dis'].max()+1, 256, endpoint = True)

degree = 1
poly = preprocessing.PolynomialFeatures(degree=degree)
X_ = poly.fit_transform(boston[['dis']])
model = sm.OLS(boston[['nox']], X_)
result = model.fit()
y = result.params.const + result.params.x1*x
ax = fig.add_subplot(5, 2, degree)
add_polt(ax, x, y, result.ssr, degree)

degree = 2
poly = preprocessing.PolynomialFeatures(degree=degree)
X_ = poly.fit_transform(boston[['dis']])
model = sm.OLS(boston[['nox']], X_)
result = model.fit()
y = result.params.const + result.params.x1*x + result.params.x2*(x**2)
ax = fig.add_subplot(5, 2, degree)
add_polt(ax, x, y, result.ssr, degree)

degree = 3
poly = preprocessing.PolynomialFeatures(degree=degree)
X_ = poly.fit_transform(boston[['dis']])
model = sm.OLS(boston[['nox']], X_)
result = model.fit()
y = result.params.const + result.params.x1*x + result.params.x2*(x**2) + result.params.x3*(x**3)
ax = fig.add_subplot(5, 2, degree)
add_polt(ax, x, y, result.ssr, degree)

degree = 4
poly = preprocessing.PolynomialFeatures(degree=degree)
X_ = poly.fit_transform(boston[['dis']])
model = sm.OLS(boston[['nox']], X_)
result = model.fit()
y = result.params.const + result.params.x1*x + result.params.x2*(x**2) + result.params.x3*(x**3)
+ result.params.x4*(x**4)
ax = fig.add_subplot(5, 2, degree)
add_polt(ax, x, y, result.ssr, degree)

plt.show()

np.random.seed(5)
X_train, X_test, y_train, y_test = train_test_split(boston[['dis']], boston[['nox']], test_size=0.2)
polynomial_regression(X_train, y_train.values, X_test, y_test, [1,2,3,4,5,6,8,9,10])

college = pd.read_csv("data/College.csv")
college = college.rename(columns={'Unnamed: 0': 'Name'})
college['Private'] = college['Private'].map({'Yes': 1, 'No': 0})
college.head()

from sklearn.linear_model import LinearRegression
from sklearn.feature_selection import RFE
from sklearn.metrics import mean_squared_error

np.random.seed(1)
msk = np.random.rand(len(college)) < 0.8
train = college[msk]
test = college[~msk]

total_features = 17
validation_MSE = {}
for m in range(1, total_features+1):
    estimator = LinearRegression(fit_intercept=True)
    selector = RFE(estimator, m, step=1)
    selector.fit(train.drop(['Name', 'Outstate'], axis=1), train['Outstate'])
    predictions = selector.predict(test.drop(['Name', 'Outstate'], axis=1))
    mse = mean_squared_error(predictions, test['Outstate'])
    validation_MSE[m] = mse/len(test)

# Plot validation MSE
lists = sorted(validation_MSE.items()) # sorted by key, return a list of tuples
x, y = zip(*lists) # unpack a list of pairs into two tuples
fig = plt.figure(figsize=(15, 8))
ax = fig.add_subplot(111)
plt.plot(x, y, color='r')
plt.grid()
ax.set_xlabel('Number of features')
ax.set_ylabel('Validation MSE')
ax.set_title('Validation MSE vs Number of features')