# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="2oa_zJH8wBkj"
# http://www.science.smith.edu/~jcrouser/SDS293/labs/lab11-py.html#6.7.2-Partial-Least-Squares 

# + id="IGTdAKOBqHF4" executionInfo={"status": "ok", "timestamp": 1627152699911, "user_tz": -180, "elapsed": 273, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import scale 
from sklearn import model_selection
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.cross_decomposition import PLSRegression, PLSSVD
from sklearn.metrics import mean_squared_error

# + id="obTOufHRqpQQ" executionInfo={"status": "ok", "timestamp": 1627152700200, "user_tz": -180, "elapsed": 15, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}
df = pd.read_csv('https://raw.githubusercontent.com/borisgarbuzov/schulich_data_science_1/master/Data/Hitters.csv').dropna()
# df = df.drop(df.columns[0], axis=1)
df = df.rename(columns = {"Unnamed: 0" : "Player"})

dummies = pd.get_dummies(df[['League', 'Division', 'NewLeague']])

# + colab={"base_uri": "https://localhost:8080/", "height": 197} id="CGHtohq-qpTg" executionInfo={"status": "ok", "timestamp": 1627152700516, "user_tz": -180, "elapsed": 52, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="593888da-0ebc-4bea-ec32-8dddb2225131"
df.head()

# + colab={"base_uri": "https://localhost:8080/", "height": 197} id="l9r7e_0AqpWb" executionInfo={"status": "ok", "timestamp": 1627152700518, "user_tz": -180, "elapsed": 42, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="7d7a36fc-3068-490e-ddb9-d24d53386c5f"
dummies.head()

# + id="N4gISZq6qpZs" executionInfo={"status": "ok", "timestamp": 1627152700521, "user_tz": -180, "elapsed": 41, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}
y = df.Salary

# Drop the column with the independent variable (Salary), and columns for which we created dummy variables
X_ = df.drop(['Player','Salary', 'League', 'Division', 'NewLeague'], axis=1).astype('float64')

# Define the feature set X.
X = pd.concat([X_, dummies[['League_N', 'Division_W', 'NewLeague_N']]], axis=1)

# + id="i3PsoXDlqpcS" executionInfo={"status": "ok", "timestamp": 1627152700524, "user_tz": -180, "elapsed": 42, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}
pca = PCA()
X_reduced = pca.fit_transform(scale(X))

# + colab={"base_uri": "https://localhost:8080/", "height": 197} id="i7AYNy3hqpe6" executionInfo={"status": "ok", "timestamp": 1627152700525, "user_tz": -180, "elapsed": 42, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="4fe9122e-0340-4a72-d21c-dcb160f00f73"
pd.DataFrame(pca.components_.T).loc[:4,:5]

# + colab={"base_uri": "https://localhost:8080/", "height": 294} id="YsBSedosqpjq" executionInfo={"status": "ok", "timestamp": 1627152701382, "user_tz": -180, "elapsed": 893, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="279a9db4-c9a8-4bbf-f766-043e3827892d"
# 10-fold CV, with shuffle
n = len(X_reduced)
kf_10 = model_selection.KFold( n_splits=10, shuffle=True, random_state=1)

regr = LinearRegression()
mse = []

# Calculate MSE with only the intercept (no principal components in regression)
score = -1*model_selection.cross_val_score(regr, np.ones((n,1)), y.ravel(), cv=kf_10, scoring='neg_mean_squared_error').mean()    
mse.append(score)

# Calculate MSE using CV for the 19 principle components, adding one component at the time.
for i in np.arange(1, 20):
    score = -1*model_selection.cross_val_score(regr, X_reduced[:,:i], y.ravel(), cv=kf_10, scoring='neg_mean_squared_error').mean()
    mse.append(score)
    
# Plot results    
plt.plot(mse, '-v')
plt.xlabel('Number of principal components in regression')
plt.ylabel('MSE')
plt.title('Salary')
plt.xlim(xmin=-1);

# + [markdown] id="EK8Ydvmxvgv6"
# We see that the smallest cross-validation error occurs when M=18 components are used. This is barely fewer than M=19
#
# , which amounts to simply performing least squares, because when all of the components are used in PCR no dimension reduction occurs. However, from the plot we also see that the cross-validation error is roughly the same when only one component is included in the model. This suggests that a model that uses just a small number of components might suffice.
#
# We'll do a little math to get the amount of variance explained by adding each consecutive principal component:
#

# + colab={"base_uri": "https://localhost:8080/"} id="01MIPRNHqpmd" executionInfo={"status": "ok", "timestamp": 1627152701386, "user_tz": -180, "elapsed": 30, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="af9ee944-26a4-4006-863f-4c13c3820007"
np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)

# + colab={"base_uri": "https://localhost:8080/", "height": 294} id="it6NCz2_vWBa" executionInfo={"status": "ok", "timestamp": 1627152702074, "user_tz": -180, "elapsed": 708, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="0db23ddd-3bd9-4953-9602-8ea6a65b180c"
pca2 = PCA()

# Split into training and test sets
X_train, X_test , y_train, y_test = model_selection.train_test_split(X, y, test_size=0.5, random_state=1)

# Scale the data
X_reduced_train = pca2.fit_transform(scale(X_train))
n = len(X_reduced_train)

# 10-fold CV, with shuffle
kf_10 = model_selection.KFold( n_splits=10, shuffle=True, random_state=1)

mse = []

# Calculate MSE with only the intercept (no principal components in regression)
score = -1*model_selection.cross_val_score(regr, np.ones((n,1)), y_train.ravel(), cv=kf_10, scoring='neg_mean_squared_error').mean()    
mse.append(score)

# Calculate MSE using CV for the 19 principle components, adding one component at the time.
for i in np.arange(1, 20):
    score = -1*model_selection.cross_val_score(regr, X_reduced_train[:,:i], y_train.ravel(), cv=kf_10, scoring='neg_mean_squared_error').mean()
    mse.append(score)

plt.plot(np.array(mse), '-v')
plt.xlabel('Number of principal components in regression')
plt.ylabel('MSE')
plt.title('Salary')
plt.xlim(xmin=-1);

# + colab={"base_uri": "https://localhost:8080/"} id="_fMClX8fqpuo" executionInfo={"status": "ok", "timestamp": 1627152702075, "user_tz": -180, "elapsed": 18, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="d47138b2-5a81-4051-b241-a420fc81ce13"
X_reduced_test = pca2.transform(scale(X_test))[:,:7]

# Train regression model on training data 
regr = LinearRegression()
regr.fit(X_reduced_train[:,:7], y_train)

# Prediction with test data
pred = regr.predict(X_reduced_test)
mean_squared_error(y_test, pred)

# + [markdown] id="pTXXamIQv063"
# Partial Least Squares

# + colab={"base_uri": "https://localhost:8080/", "height": 312} id="AtbKXCPsqpxc" executionInfo={"status": "ok", "timestamp": 1627152702683, "user_tz": -180, "elapsed": 615, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="f7200fe8-4291-4ea8-887d-ea14fd7f1a40"
n = len(X_train)

# 10-fold CV, with shuffle
kf_10 = model_selection.KFold(n_splits=10, shuffle=True, random_state=1)

mse = []

for i in np.arange(1, 20):
    pls = PLSRegression(n_components=i)
    score = model_selection.cross_val_score(pls, scale(X_train), y_train, cv=kf_10, scoring='neg_mean_squared_error').mean()
    mse.append(-score)

# Plot results
plt.plot(np.arange(1, 20), np.array(mse), '-v')
plt.xlabel('Number of principal components in regression')
plt.ylabel('MSE')
plt.title('Salary')
plt.xlim(xmin=-1)

# + colab={"base_uri": "https://localhost:8080/"} id="sIiL_9JjvsMW" executionInfo={"status": "ok", "timestamp": 1627152702960, "user_tz": -180, "elapsed": 287, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}} outputId="07a3cb9f-bb61-41e0-cafd-3205a813ef93"
pls = PLSRegression(n_components=2)
pls.fit(scale(X_train), y_train)

mean_squared_error(y_test, pls.predict(scale(X_test)))

# + id="KfEoyOg4vsPP" executionInfo={"status": "ok", "timestamp": 1627152702962, "user_tz": -180, "elapsed": 35, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}


# + id="HO_gkGU6vsR6" executionInfo={"status": "ok", "timestamp": 1627152702964, "user_tz": -180, "elapsed": 35, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}


# + id="L5GRCaKCvsVE" executionInfo={"status": "ok", "timestamp": 1627152702965, "user_tz": -180, "elapsed": 34, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}


# + id="jbJvvZSVvsXg" executionInfo={"status": "ok", "timestamp": 1627152702980, "user_tz": -180, "elapsed": 47, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}


# + id="VFLuuuZFqp0P" executionInfo={"status": "ok", "timestamp": 1627152702981, "user_tz": -180, "elapsed": 42, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "08851466340381164023"}}


StarcoderdataJupyterScriptsDedupFiltered