Ridge, Lasso, and PCR

import itertools
import pandas as pd
import numpy as np
import copy

import statsmodels.api as sm
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.decomposition import PCA
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import train_test_split, GridSearchCV, KFold
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error

import matplotlib.pyplot as plt
import seaborn as sns
from IPython.display import Markdown, display

# Reset all styles to the default:
plt.rcParams.update(plt.rcParamsDefault)
# Then make graphs inline:
%matplotlib inline

# Useful function for Jupyter to display text in bold:
def displaybd(text):
    display(Markdown("**" + text + "**"))

plt.rcParams['figure.figsize'] = (7, 6)
plt.rcParams['font.size'] = 24
plt.rcParams['legend.fontsize'] = 'large'
plt.rcParams['figure.titlesize'] = 'large'
plt.rcParams['lines.markersize'] = 10

# LOAD DATA:
hittersDF = pd.read_csv('Hitters.csv', na_values=[""])
# The first column has no name in the csv file:
hittersDF.rename(columns={hittersDF.columns[0] : "Name"}, inplace=True, copy=False)
hittersDF.set_index('Name', inplace=True)
hittersDF.dropna(inplace=True)

# CREATE X and y:

# Convert categorical variables into dummies:
dummies = pd.get_dummies(hittersDF[['League', 'Division', 'NewLeague']])
X = hittersDF.drop(['Salary', 'League', 'Division', 'NewLeague'], axis=1)
cat_features = ['League_N', 'Division_W', 'NewLeague_N']
num_features = list(X.columns)

X = pd.concat([X, dummies[cat_features]], axis=1).astype('float64')
y = hittersDF.Salary

# From here: https://kiwidamien.github.io/introducing-the-column-transformer.html
# and from "Hands-On Machine Learning with Scikit-Learn & TensorFlow" by Aurélien Géron.

# Feature scaling is super important for Ridge and Lasso.

# Pipeline for numerical columns:
num_pipeline = Pipeline([
    # Replace missing numerical values with their median values.
    # Actually, our dataset does not have missing values (we dropped them all),
    # but in case we decide not to drop the missing values, we can do the imputation.
    ('impute', SimpleImputer(strategy='median')),
    ('std_scaler', StandardScaler()),
])

# Pipeline for categorical columns:
cat_pipeline = Pipeline([
    # Do not scale data. Replace missing values with the most frequent values:
    ('impute', SimpleImputer(strategy='most_frequent'))
])

full_preprocess_pipeline = ColumnTransformer([
    ('numerical', num_pipeline, num_features),
    ('categorical', cat_pipeline, cat_features)
])

train_preprocess_pipeline = copy.copy(full_preprocess_pipeline)

X_processed = full_preprocess_pipeline.fit_transform(X)

#-- SPLIT INTO TWO TRAIN AND TEST DATASETS:

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=1)
X_train_processed = train_preprocess_pipeline.fit_transform(X_train)
# Note that we scale test data using the train data's mean and variance. This is important.
X_test_processed = train_preprocess_pipeline.transform(X_test)
print(X_train_processed[:4])

def printCoefs(regr_res, coef_names, num_features_idx, mean, var):
    scaled_coefs = np.hstack((regr_res.intercept_, regr_res.coef_))
    
    # Only numerical features are scaled. So, we need to extract them, unscale
    # the coefficients, then merge coefficients for numerical and categorical
    # features.

    cat_features_idx = [x for x in range(len(regr_res.coef_))\
                        if x not in num_features_idx]

    unscaled_num_coefs = regr_res.coef_[num_features_idx] / np.sqrt(var)
    unscaled_intercept = regr_res.intercept_ - np.sum(unscaled_num_coefs * mean)

    unscaled_coefs = np.array([None] * (len(regr_res.coef_) + 1))
    # Intercept goes first, then the rest of coefficients:
    unscaled_coefs[0] = unscaled_intercept
    unscaled_coefs[[i+1 for i in num_features_idx]] = unscaled_num_coefs
    unscaled_coefs[[i+1 for i in cat_features_idx]] = ridge.coef_[cat_features_idx]

    results = pd.DataFrame(data={'Scaled' : scaled_coefs,
                                  'Unscaled' : unscaled_coefs},
                            index=['Intercept']+coef_names)
    print(results)

# NB: glmnet uses (1/N)*||y - Xw||^2_2 + lambda/2 * ||w||^2_2, while
# Ridge uses ||y - Xw||^2_2 + alpha * ||w||^2_2. So, there is a difference.

alphas = 10**np.linspace(5, -2, 300)

# NB: Do not use the normalize option of Ridge, use scale or StandardScaler from sklearn.preprocessing
# See here: https://stackoverflow.com/questions/24098212/what-does-the-option-normalize-true-in-lasso-sklearn-do
ridge = Ridge(fit_intercept=True, normalize=False, tol=1e-12)

coefs = []
for a in alphas:
    ridge.set_params(alpha=a)
    ridge.fit(X_processed, y)
    coefs.append(ridge.coef_)

ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1])  # reverse axis
plt.axis('tight')
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Ridge coefficients as a function of the regularization');

ridge2 = Ridge(alpha=4, normalize=False, tol=1e-12)
ridge2.fit(X_train_processed, y_train)
yhat = ridge2.predict(X_test_processed)
print("Ridge regression with alpha =", 4)
print("MSE:", mean_squared_error(y_test, yhat))

# NB: Don't use 'RidgeCV'!
# See here: https://stackoverflow.com/questions/39466671/use-of-scaler-with-lassocv-ridgecv

ridge_pipeline = Pipeline([
    ('scaler', train_preprocess_pipeline),
    ('regressor', Ridge(fit_intercept=True, normalize=False, tol=1e-12))
])

# We need to supply parameter alpha to Ridge. Note the convention adopted to
# name the parameters: name of the pipeline step, followed by a double underscore (__),
# then finally the name of the parameter within the step
pipeline_params = {'regressor__alpha': alphas}

# Using our pipline with GridSearchCV allows us to scale the train part of X_train in each fold
# and then apply that scaling to the test part of X_train in each fold. This is the correct way to
# K-fold cross-validation. This way there is no information leakage from the train part of a fold to
# the test part of a fold.
kf = KFold(n_splits=10, shuffle=True, random_state=1)
gridsearch = GridSearchCV(ridge_pipeline, param_grid=pipeline_params, cv=kf,
                          scoring='neg_mean_squared_error',
                          return_train_score=False, # We don't need them
                          n_jobs=6, # Use 6 CPU cores
                          verbose=1, iid=False).fit(X_train, y_train)

best_ridge_alpha = gridsearch.best_params_['regressor__alpha']
print("Best alpha:", best_ridge_alpha)
print('MSE on the test dataset is: ', -gridsearch.score(X_test, y_test))

ridge.set_params(alpha=best_ridge_alpha)
fitted_scaler = full_preprocess_pipeline.transformers_[0][1].steps[1][1]
displaybd("Regression coefficients of the optimal ridge regression on the full dataset:")
printCoefs(ridge, num_features + cat_features,
                list(range(len(num_features))),
                fitted_scaler.mean_, fitted_scaler.var_)

lasso = Lasso(fit_intercept=True, normalize=False, tol=1e-12, max_iter=100000)

coefs = []
for a in alphas:
    lasso.set_params(alpha=a)
    lasso.fit(X_processed, y)
    coefs.append(lasso.coef_)

ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1])  # reverse axis
plt.axis('tight')
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Lasso coefficients as a function of the regularization');

# NB: Again, don't use 'LassoCV'!
# See here: https://stackoverflow.com/questions/39466671/use-of-scaler-with-lassocv-ridgecv

lasso_pipeline = Pipeline([
    ('scaler', train_preprocess_pipeline),
    ('regressor', Lasso(fit_intercept=True, normalize=False, tol=1e-12, max_iter=10000))
])

# We need a much lower alpha for Lasso (it means a different thing in Lasso versus Ridge).
# I just divide it by 10:
pipeline_params = {'regressor__alpha': alphas / 10}

# Using our pipline with GridSearchCV allows us to scale the train part of X_train in each fold
# and then apply that scaling to the test part of X_train in each fold. This is the correct way to
# K-fold cross-validation. This way there is no information leakage from the train part of a fold to
# the test part of a fold.
kf = KFold(n_splits=10, shuffle=True, random_state=1)
gridsearch = GridSearchCV(lasso_pipeline, param_grid=pipeline_params, cv=kf,
                          scoring='neg_mean_squared_error',
                          return_train_score=False, # We don't need them
                          n_jobs=6, # Use 6 CPU cores
                          verbose=1, iid=False).fit(X_train, y_train)

best_lasso_alpha = gridsearch.best_params_['regressor__alpha']
print("Best alpha:", best_lasso_alpha)
print('MSE on the test dataset is: ', -gridsearch.score(X_test, y_test))

lasso.set_params(alpha=best_lasso_alpha)
fitted_scaler = full_preprocess_pipeline.transformers_[0][1].steps[1][1]
displaybd("Regression coefficients of the optimal lasso regression on the full dataset:")
printCoefs(lasso, num_features + cat_features,
                list(range(len(num_features))),
                fitted_scaler.mean_, fitted_scaler.var_)

pca_nonscaled = PCA()
pca_scaled = PCA()

# In ridge and lasso we were not scaling categorical variables.
# Similarly, we can do PCA without scaling categorical variables.
# Moreover, it is debatable, whether it is appropriate at all to
# include binary variables into PCA. Probably, it is OK.
# See here: https://stats.stackexchange.com/questions/16331/doing-principal-component-analysis-or-factor-analysis-on-binary-data
X_reduced_nonscaled = pca_nonscaled.fit_transform(X_processed)

scaler = StandardScaler()
X_reduced_scaled = pca_scaled.fit_transform(scaler.fit_transform(X))

#-- PLOT THE SHARE OF EXPLAINED VARIANCE:

ax = plt.gca()
nonscaled_cat_line, = \
    ax.plot(range(1, X_reduced_nonscaled.shape[1] + 1),
            np.cumsum(pca_nonscaled.explained_variance_ratio_), 
            '-o', markersize=4, color='blue')
scaled_cat_line,  = \
    ax.plot(range(1, X_reduced_scaled.shape[1] + 1),
            np.cumsum(pca_scaled.explained_variance_ratio_),
            '-o', markersize=4, color='red')
ax.set_xticks(range(1, X_reduced_scaled.shape[1] + 1, 2));
plt.xlabel('Number of components');
plt.ylabel('Share of explained variance');
plt.legend([nonscaled_cat_line, scaled_cat_line],
           ['Categorical scaled', 'Categorical unscaled'],
           fontsize=18);

n_features_to_test = np.arange(1, X.shape[1] + 1)
pipeline_params = {'reduce_dim__n_components': n_features_to_test}
kf = KFold(n_splits=10, shuffle=True, random_state=90)

# Pipeline for unscaled categorical variables:
pca_regr_pipe_unscaled = Pipeline([
    ('scaler', full_preprocess_pipeline),
    ('reduce_dim', PCA()),
    ('regressor', LinearRegression())
])

# Pipeline for scaled categorical variables:
pca_regr_pipe_scaled = Pipeline([
    ('scaler', StandardScaler()),
    ('reduce_dim', PCA()),
    ('regressor', LinearRegression())
])

gridsearch_unscaled = GridSearchCV(pca_regr_pipe_unscaled,
                                   param_grid=pipeline_params, cv=kf,
                                   scoring='neg_mean_squared_error',
                                   return_train_score=False,
                                   n_jobs=6, # Use 6 CPU cores
                                   verbose=1, iid=False).fit(X, y)

gridsearch_scaled = GridSearchCV(pca_regr_pipe_scaled,
                                 param_grid=pipeline_params, cv=kf,
                                 scoring='neg_mean_squared_error',
                                 return_train_score=False,
                                 n_jobs=6, # Use 6 CPU cores
                                 verbose=1, iid=False).fit(X, y)

rmse_unscaled = np.sqrt(-gridsearch_unscaled.cv_results_["mean_test_score"])
rmse_scaled = np.sqrt(-gridsearch_scaled.cv_results_["mean_test_score"])

ax = plt.gca()
nonscaled_cat_line, = \
    ax.plot(range(1, len(rmse_unscaled) + 1), rmse_unscaled, 
            '-o', markersize=4, color='blue')
scaled_cat_line,  = \
    ax.plot(range(1, len(rmse_scaled) + 1), rmse_scaled, 
            '-o', markersize=4, color='red')
ax.set_xticks(range(1, len(rmse_unscaled) + 1, 2));
plt.xlabel('Number of components');
plt.ylabel('Root of CV MSE');
plt.legend([nonscaled_cat_line, scaled_cat_line],
           ['Categorical scaled', 'Categorical unscaled']);

gridsearch_scaled = GridSearchCV(pca_regr_pipe_scaled,
                                 param_grid=pipeline_params, cv=kf,
                                 scoring='neg_mean_squared_error',
                                 return_train_score=False,
                                 n_jobs=6, # Use 6 CPU cores
                                 verbose=1, iid=False).fit(X_train, y_train)

#-- PLOT CROSS-VALIDATED (ROOT OF) MEAN SQUARED ERROR:

rmse_scaled = np.sqrt(-gridsearch_scaled.cv_results_["mean_test_score"])

ax = plt.gca()
ax.plot(range(1, len(rmse_scaled) + 1), rmse_scaled, 
        '-o', markersize=4, color='red')
ax.set_xticks(range(1, len(rmse_scaled) + 1, 2));
plt.xlabel('Number of components');
plt.ylabel('Root of CV MSE');
plt.title('Scaled categorical variables')

pca_regr_pipe_scaled.set_params(reduce_dim__n_components=7)
pca_regr_pipe_scaled.fit(X_train, y_train)
yhat = pca_regr_pipe_scaled.predict(X_test)
print("PCA regression with 7 components")
print("Test MSE:", mean_squared_error(y_test, yhat))