Chapter_6_sec_6.5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.5 Lab: Linear Models and Regularization Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5.1 Subset Selection Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This notebook contains the code for best subset selection, \n",
    "so this notebook may take longer time to run, for faster run, make\n",
    "max_feature into a smaller number\n",
    "\"\"\"\n",
    "max_feature = 3\n",
    "\"\"\"\n",
    "As we move to more ambgious topics, I would try to add in more comments on the approaches and results.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import pandas as pd \n",
    "import math\n",
    "import random\n",
    "import itertools\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.graphics.regressionplots import *\n",
    "from sklearn import datasets, linear_model\n",
    "from sklearn.decomposition import PCA \n",
    "from sklearn.cross_decomposition import PLSRegression, PLSSVD\n",
    "from sklearn.preprocessing import StandardScaler, scale \n",
    "from sklearn import linear_model \n",
    "from sklearn.model_selection import cross_val_predict \n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import Ridge, RidgeCV, Lasso, LassoCV\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# read the data and take a look at the data\n",
    "Hitters = pd.read_csv('data/Hitters.csv', header=0, na_values='NA')\n",
    "print(list(Hitters)) # get the header of this data\n",
    "print(Hitters.shape) # get the dimension of this \n",
    "Hitters.head() # pull a sample of this dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# clean the data a bit to remove NAs\n",
    "print(np.sum(pd.isnull(Hitters['Salary']))) # number of NAs in Salary column'\n",
    "print(Hitters['Salary'].isnull().sum())\n",
    "\n",
    "Hitters = Hitters.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "print(Hitters.shape)\n",
    "print(Hitters['Salary'].isnull().sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# do some feature engineering to prepare for the model training\n",
    "y = Hitters.Salary  \n",
    "\n",
    "\"\"\"\n",
    "take care of the features \n",
    "1. change category into dummy variables \n",
    "2. Choose (n-1) dummy variable into the feature set: n is the unique values of each categorical variable.\n",
    "\"\"\"\n",
    "\n",
    "dummies = pd.get_dummies(Hitters[['League', 'Division', 'NewLeague']])\n",
    "print(dummies.head())\n",
    "X_prep = Hitters.drop (['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n",
    "X = pd.concat([X_prep,  dummies[['League_A', 'Division_E', 'NewLeague_A']]], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let us get it ready for best subset selection \n",
    "\"\"\"\n",
    "Since in Python there is no well-defined function for best subset selection, \n",
    "we will need to define some functions ourselves.\n",
    "1. Define a function to run on a subset of feature and extract RSS\n",
    "2. Select the best model (models) for a fix number of features\n",
    "\"\"\"\n",
    "def getRSS(y, X, feature_list):\n",
    "    model = sm.OLS(y, X[list(feature_list)]).fit()\n",
    "    RSS = ((model.predict(X[list(feature_list)]) - y) ** 2).sum()\n",
    "    return {'Model':model, \"RSS\":RSS}\n",
    "\n",
    "def bestModel(y, X, K):\n",
    "    results = []\n",
    "    for c in itertools.combinations(X.columns, K):\n",
    "        results.append(getRSS(y, X, c))     \n",
    "    model_all =  pd.DataFrame(results)\n",
    "    \n",
    "    best_model = model_all.loc[model_all[\"RSS\"].idxmin()] ## this could be modified to have the top several models\n",
    "    return best_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# execute the best model selection \n",
    "models = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "for i in range(1,(max_feature+1)):  # for illustration purpuse, I just run for 1 - max_fearure features \n",
    "    models.loc[i] = bestModel(y, X, i)\n",
    "    \n",
    "print(models.loc[2, 'Model'].summary())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this show an example to plot the RSS of best models with different number of parameters\n",
    "plt.figure()\n",
    "plt.plot(models[\"RSS\"])\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('RSS')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# find the adjust R^2, use dir() to identify all available attributes\n",
    "rsquared_adj = models.apply(lambda row: row[1].rsquared_adj, axis=1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "The following graph shows the adj R^2 is still increasing, \n",
    "in this case, it is a good idea trying models with more features. \n",
    "\"\"\"\n",
    "plt.figure()\n",
    "plt.plot(rsquared_adj)\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('Adjust R^2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5.2 Forward and Backward Stepwise Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward Stepwise Selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we can use the previous user defined function 'def getRSS(y, X, feature_list)' \n",
    "# to add 1 feature at a time (start from 0 feature) for forward stepwise selection\n",
    "# or delete 1 feature at a time(start from all the features) for backward stepwise selection. \n",
    "def forward_select(y, X, feature_list):\n",
    "    remaining_predictors = [p for p in X.columns if p not in feature_list]\n",
    "    results = []\n",
    "    for p in remaining_predictors:\n",
    "        results.append(getRSS(y, X, feature_list+[p]))\n",
    "\n",
    "    models = pd.DataFrame(results)\n",
    "    best_model = models.loc[models['RSS'].idxmin()]\n",
    "    return best_model  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "models2 = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "feature_list = []\n",
    "for i in range(1,len(X.columns)+1):\n",
    "    models2.loc[i] = forward_select(y, X, feature_list)\n",
    "    feature_list = models2.loc[i][\"Model\"].model.exog_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we can compare the results of best subset selection and the forward selection\n",
    "print('Best max_feature variable from best subset selection on tranining')\n",
    "print(models.loc[max_feature, 'Model'].params)\n",
    "print('\\n---------------------------------------------')\n",
    "print('Best max_feature variable from forward selection on tranining')\n",
    "print(models2.loc[max_feature, 'Model'].params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backward Stepwise Selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def backward_select(y, X, feature_list):\n",
    "    results = []\n",
    "    for combo in itertools.combinations(feature_list, len(feature_list)-1):\n",
    "        results.append(getRSS(y, X, combo))\n",
    "\n",
    "    models = pd.DataFrame(results)\n",
    "    best_model = models.loc[models['RSS'].idxmin()]\n",
    "    return best_model  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the backward selection starts from all the variables of features\n",
    "models3 = pd.DataFrame(columns=[\"RSS\", \"Model\"], index = range(1,len(X.columns)))\n",
    "feature_list = X.columns\n",
    "\n",
    "while(len(feature_list) > 1):\n",
    "    models3.loc[len(feature_list)-1] = backward_select(y, X, feature_list)\n",
    "    feature_list = models3.loc[len(feature_list)-1][\"Model\"].model.exog_names\n",
    "\n",
    "print(models3.loc[max_feature, \"Model\"].params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5.3 Choosing Among Models Using the Validation Set Approach and Cross-Validation\n",
    "In previous sections, we defined the 'best' model based on some statistics (R^2, adj R^2, AIC, BIC, etc) of the training dataset. This may cause 'overfitting' problemm which means the best model on training data can not generalize well to new data. In this section, validation approach will be discussed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# randomly split the data into traning dataset and validation dateset\n",
    "np.random.seed(seed = 21)\n",
    "train_index = np.random.choice([True, False], size = len(y), replace = True, p = [0.7, 0.3]) \n",
    "# random select ~70% of data into traning sample\n",
    "# the rest of the samples will be in testing set.\n",
    "test_index = np.invert(train_index)\n",
    "X_train= X[train_index]\n",
    "y_train = y[train_index]\n",
    "X_test = X[test_index]\n",
    "y_test = y[test_index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "We can recyle the old functions. Modification is needed to compute the RSS for the testing data. \n",
    "So we need to add both train and test into the function input (Implement)\n",
    "-OR-: we can wrap the train and test split step into the function(Not Implemented)\n",
    "\"\"\"\n",
    "def getRSS_validation(y_train, X_train, y_test, X_test,  feature_list):\n",
    "    model = sm.OLS(y_train, X_train[list(feature_list)]).fit()\n",
    "    RSS = ((model.predict(X_test[list(feature_list)]) - y_test) ** 2).sum()\n",
    "    return {'Model':model, \"RSS\":RSS}\n",
    "\n",
    "def bestModel_validation(y_train, X_train, y_test, X_test, K):\n",
    "    results = []\n",
    "    for c in itertools.combinations(X_train.columns, K):\n",
    "        results.append(getRSS_validation(y_train, X_train, y_test, X_test, c))     \n",
    "    model_all =  pd.DataFrame(results)\n",
    "    \n",
    "    best_model = model_all.loc[model_all[\"RSS\"].idxmin()] ## this could be modified to have the top several models\n",
    "    return best_model\n",
    "\n",
    "\n",
    "def forward_select_validation(y_train, X_train, y_test, X_test,  feature_list):\n",
    "    remaining_predictors = [p for p in X_train.columns if p not in feature_list]\n",
    "    results = []\n",
    "    for p in remaining_predictors:\n",
    "        results.append(getRSS_validation(y_train, X_train, y_test, X_test, feature_list+[p]))\n",
    "\n",
    "    models = pd.DataFrame(results)\n",
    "    best_model = models.loc[models['RSS'].idxmin()]\n",
    "    return best_model\n",
    "\n",
    "def backward_select_validation(y_train, X_train, y_test, X_test,  feature_list):\n",
    "    results = []\n",
    "    for combo in itertools.combinations(feature_list, len(feature_list)-1):\n",
    "        results.append(getRSS_validation(y_train, X_train, y_test, X_test,  combo))\n",
    "\n",
    "    models = pd.DataFrame(results)\n",
    "    best_model = models.loc[models['RSS'].idxmin()]\n",
    "    return best_model  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "models_validation = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "for i in range(1,(max_feature+1)):  # for illustration purpuse, I just run for 1 - max_fearure features \n",
    "    models_validation.loc[i] = bestModel_validation(y_train, X_train, y_test, X_test, i) \n",
    "    \n",
    "    \n",
    "\"\"\"\n",
    "change the function to  forward_select_validation (.) or backward_select_validation(.) \n",
    "for forward selection or backward selection\n",
    "\"\"\" \n",
    "    \n",
    "models2_forward = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "feature_list = []\n",
    "for i in range(1,len(X.columns)+1):\n",
    "    models2_forward.loc[i] = forward_select_validation(y_train, X_train, y_test, X_test,  feature_list)\n",
    "    feature_list = models2_forward.loc[i][\"Model\"].model.exog_names    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print('Best max_feature variable from best subset selection on tranining')\n",
    "print(models.loc[max_feature, 'Model'].params)\n",
    "print('\\n---------------------------------------------')\n",
    "print('Best max_feature variable from forward selection on tranining')\n",
    "print(models2.loc[max_feature, 'Model'].params)\n",
    "print('\\n---------------------------------------------')\n",
    "print('Best max_feature variable from backward selection on tranining')\n",
    "print(models3.loc[max_feature, 'Model'].params)\n",
    "print('\\n---------------------------------------------')\n",
    "print('Best max_feature variable from best subset selection on traning and validation split')\n",
    "print(models_validation.loc[max_feature, 'Model'].params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this show an example to plot the RSS of best models with different number of parameters for best subset with validation\n",
    "plt.figure()\n",
    "plt.plot(models_validation[\"RSS\"])\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('RSS')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# this show an example to plot the RSS of best models with different number of parameters for forward selection with validation\n",
    "plt.figure()\n",
    "plt.plot(models2_forward[\"RSS\"])\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('RSS')\n",
    "plt.show()\n",
    "\"\"\"\n",
    "Note:from above graph, that 6 variables model gives us the best RSS under forward selection. \n",
    "To learn the final model, it is also recommendated to re-train the model on entire data (train + validation). \n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This approach is similar to the previous validation idea. \n",
    "The difference is that we break the entire dataset into K different folds. \n",
    "Each run, the model will be trained on all the data from K-1 folds and tested on the remaining fold.\n",
    "Advantages: \n",
    "1. There will be multiple metrics out of testing results => distribution of testing RSS, etc\n",
    "2. The size of traning dataset is much closer to the size of original dateset. This will remove some biases caused by\n",
    "the size difference.\n",
    "\"\"\"\n",
    "\n",
    "k = 10\n",
    "np.random.seed(seed = 21)\n",
    "train_index = np.random.choice(k, size = len(y), replace = True)  # Randomly assign each observations into folds\n",
    "cv_errors = pd.DataFrame(columns=range(1,k+1), index=range(1,len(X.columns) + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# train and test on each fold\n",
    "models_cv = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "for j in range(1,k+1):\n",
    "    feature_list = []\n",
    "    for i in range(1,len(X.columns)+1):\n",
    "        models_cv.loc[i] = forward_select_validation(y[train_index!= (j-1)], X[train_index != (j-1)], \n",
    "                                                     y[train_index == (j-1)],X[train_index == (j-1)], \n",
    "                                                     feature_list)\n",
    "        \n",
    "        cv_errors[j][i] = models_cv.loc[i][\"RSS\"]\n",
    "        feature_list = models_cv.loc[i][\"Model\"].model.exog_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# similar to above, let us plot the performance of each model \n",
    "cv_errors_mean = cv_errors.mean(axis = 1)\n",
    "plt.figure()\n",
    "plt.plot(cv_errors_mean)\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('RSS')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "From the above plot, we can see that the model with 5 variables yielded the smallest RSS.\n",
    "We can take a closer look at that model summary. \n",
    "We can also see that the model performance for variables 4 - 12 are similar.\n",
    "\"\"\"\n",
    "print(models_cv.loc[5, \"Model\"].summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5.2 Ridge Regression and the Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Ridge and Lasso are the two regularization methods to deal with overfitting. \n",
    "We do not to run multiple backward/forward selection. \n",
    "In Ridge and Lasso , we would need to run the regularization strength.\n",
    "\"\"\"\n",
    "\n",
    "# similar to before, we need to prepare the features(indepedent) variables and response(depedent) varisble.\n",
    "Hitters = pd.read_csv('data/Hitters.csv', header=0, na_values='NA')\n",
    "Hitters = Hitters.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "dummies = pd.get_dummies(Hitters[['League', 'Division', 'NewLeague']])\n",
    "\n",
    "y = Hitters.Salary  # the response variable \n",
    "X_prep = Hitters.drop (['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n",
    "X = pd.concat([X_prep,  dummies[['League_A', 'Division_E', 'NewLeague_A']]], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ridge Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# next, we will generate a few canadidates regularization strength(in sklearn, the keyword is alphas) \n",
    "# for our Ridge regression. In R, alpha is a switch for Ridge and Lasso methods.\n",
    "alphas = 10**np.linspace(10,-2,100)\n",
    "print(alphas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# associated with each value of alpha is a vector of ridge regression coefficients, \n",
    "# stored in a matrix that can be accessed by coeffs. In this case, it is a 19×100, \n",
    "# 19 is the dimension of the features + (intercept needs to call separately) and 100 \n",
    "# is the len of the alphas. The result is a numpy series with len 100 and len(coffes[0]) is 19. \n",
    "# In this specific implementation, the default is no intercept. \n",
    "# here I used list to store them, we could also use dictionary to store them as we did before.\n",
    "ridge = Ridge(fit_intercept=True, normalize=True)\n",
    "coeffs = []\n",
    "intercepts = []\n",
    "\n",
    "for a in alphas:\n",
    "    ridge.set_params(alpha=a)\n",
    "    ridge.fit(X, y)\n",
    "    coeffs.append(ridge.coef_)\n",
    "    intercepts.append(ridge.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(len(coeffs))\n",
    "print(len(coeffs[0]))\n",
    "print(len(intercepts))\n",
    "print(intercepts[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# as design, as alphas get bigger, magnitude of coefficients turn to be closer to zero. \n",
    "# one thing to remember is that the decay is quite smooth\n",
    "ax = plt.gca()\n",
    "ax.plot(alphas, coeffs)\n",
    "\"\"\"\n",
    "Typically, we scale the alpha in log scale.\n",
    "Try to plot it without. \n",
    "\"\"\"\n",
    "ax.set_xscale('log') \n",
    "plt.axis('tight')\n",
    "plt.xlabel('alpha')\n",
    "plt.ylabel('weights')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let us take one example\n",
    "\"\"\"\n",
    "We may notice the coefficients l2 norm is different from R output \n",
    "I tried a few different normalization methods but still did not get the exact same output \n",
    "\"\"\"\n",
    "sample = 49\n",
    "print(alphas[sample])\n",
    "print(math.sqrt(sum(map(lambda x:x*x, coeffs[sample]))))\n",
    "print(coeffs[sample])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we now split the samples into a training set and a test set in order to estimate \n",
    "# the test error of ridge regression and the lasso. \n",
    "# Python provides a built-in function to produce training and test data set.\n",
    "X_train, X_test , y_train, y_test = train_test_split(X, y, test_size=0.66)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# train the model and do prediction on test dataset\n",
    "ridge = Ridge(fit_intercept=True, normalize=True, alpha=4)\n",
    "ridge.fit(X_train, y_train)             # Fit a ridge regression on the training data\n",
    "pred = ridge.predict(X_test)           # Use this model to predict the test data\n",
    "print(pd.Series(ridge.coef_, index=X.columns)) # Print coefficients\n",
    "print(mean_squared_error(y_test, pred))        # Calculate the test MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# to select best alpha, we will use cross validation. \n",
    "# as standard, we will report test set performance as the final performance metric\n",
    "ridgecv =  RidgeCV(alphas, scoring='neg_mean_squared_error', normalize = True)\n",
    "ridgecv.fit(X_train, y_train)\n",
    "ridgecv.alpha_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# produce the mean squred error of the selected alpha on test dataset\n",
    "ridge_best = Ridge(alpha=ridgecv.alpha_, normalize=True)\n",
    "ridge_best.fit(X_train, y_train)\n",
    "mean_squared_error(y_test, ridge_best.predict(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we can exame the values of the coefficients \n",
    "\"\"\"\n",
    "If we exam the values of the coefficients, most of them are tiny, but none of them is zero.\n",
    "This is a diff between ridge and lasso. Ridge shrinks the coefficients proportionally.\n",
    "\"\"\"\n",
    "pd.Series(ridge_best.coef_, index=X.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  The Lasso\n",
    "#### We saw that ridge regression with a wise choice of λ can outperform least squares as well as the null model on the Hitters data set. We now ask whether the lasso can yield either a more accurate or a more interpretable model than ridge regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# same precedure as before, we will use cross validation to select the best alpha\n",
    "lasso= Lasso(normalize=True, max_iter=1e5) \n",
    "coeffs = []\n",
    "\n",
    "for a in alphas:\n",
    "    lasso.set_params(alpha=a)\n",
    "    lasso.fit(X_train, y_train)\n",
    "    coeffs.append(lasso.coef_)\n",
    "    \n",
    "ax = plt.gca()\n",
    "ax.plot(alphas, coeffs)\n",
    "ax.set_xscale('log')\n",
    "plt.axis('tight')\n",
    "plt.xlabel('alpha')\n",
    "plt.ylabel('weights')\n",
    "plt.show()\n",
    "\n",
    "lassocv = LassoCV(alphas=None, cv=10, max_iter=1e5, normalize=True)\n",
    "lassocv.fit(X_train, y_train)\n",
    "\n",
    "lasso.set_params(alpha=lassocv.alpha_)\n",
    "lasso.fit(X_train, y_train)\n",
    "mean_squared_error(y_test, lasso.predict(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# some of the coefficients should reduce to exact zero\n",
    "pd.Series(lasso.coef_, index=X.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.5.3 PCR and PLS Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Principal Components Regression\n",
    "Principal components regression (PCR) can be performed using the PCA function + regression afterwards. The PCA function is part of the scikit-learn module. In this section, we will continue using Hitters data, in order to predict Salary. Again, we will drop NA and deal with the categorical variables from the data set as what we did in Section 6.5. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Hitters = pd.read_csv('data/Hitters.csv', header=0, na_values='NA')\n",
    "Hitters = Hitters.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "y = Hitters.Salary  # the response variable \n",
    "dummies = pd.get_dummies(Hitters[['League', 'Division', 'NewLeague']])\n",
    "X_prep = Hitters.drop (['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n",
    "X = pd.concat([X_prep,  dummies[['League_A', 'Division_E', 'NewLeague_A']]], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let us do PCA on the input dataset. \n",
    "# since the units of the variables are different, it is always recommended to scale the variables.\n",
    "pca = PCA()\n",
    "X_pca = pca.fit_transform(scale(X))\n",
    "regr = linear_model.LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Here comes a problem, how to choose the number of PCs. We still use cross-validation. \n",
    "We compute MSE for validation set at different number of PCs, \n",
    "and choose the one with lowest validation MSE as the optimal number for PCA dimension reduction. \n",
    "\n",
    "[NOT fully covered here, covered in Chapter 12] We may decide on the number of PCS by the explained variance ratio.\n",
    "\"\"\"\n",
    "\n",
    "def pcr(X,y,pc):\n",
    "    ''' Principal Component Regression in Python'''\n",
    "    ''' Step 1: PCA on input data'''\n",
    "\n",
    "    # Define the PCA object\n",
    "    pca = PCA()\n",
    "\n",
    "    # Preprocessing (1): first derivative\n",
    "    X_pca = pca.fit_transform(scale(X))[:,:pc]\n",
    "\n",
    "    ''' Step 2: regression on selected principal components'''\n",
    "\n",
    "    # Create linear regression object\n",
    "    regr = linear_model.LinearRegression()\n",
    "    \n",
    "    # Fit\n",
    "    regr.fit(X_pca, y)\n",
    "\n",
    "    # Calibration\n",
    "    y_train = regr.predict(X_pca)\n",
    "\n",
    "    # Cross-validation\n",
    "    y_cv = cross_val_predict(regr, X_pca, y, cv=20)\n",
    "\n",
    "    # Calculate scores for training and cross-validation\n",
    "    score_train = r2_score(y, y_train)\n",
    "    score_cv = r2_score(y, y_cv)\n",
    "\n",
    "    # Calculate mean square error for training and cross validation\n",
    "    mse_train = mean_squared_error(y, y_train)\n",
    "    mse_cv = mean_squared_error(y, y_cv)\n",
    "\n",
    "    return(y_cv, score_train, score_cv, mse_train, mse_cv)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mse_train = []\n",
    "mse_cv = []\n",
    "\n",
    "# calculate MSE using CV for the 19 principle components, adding one component at the time.\n",
    "for i in np.arange(1, 20):\n",
    "    results =  pcr(X,y,i)\n",
    "    mse_train.append(results[3])\n",
    "    mse_cv.append(results[4])\n",
    "    \n",
    "# plot results    \n",
    "plt.plot(np.arange(1, 20), mse_cv, '-v', label = 'Validation_MSE')\n",
    "plt.plot(np.arange(1, 20), mse_train, '-v', label = 'Train_MSE')\n",
    "plt.xlabel('Number of principal components in regression')\n",
    "plt.ylabel('MSE')\n",
    "plt.title('Salary')\n",
    "plt.xlim(left=-1)\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\"\"\"\n",
    "From the above picture, we can observe three things \n",
    "1)PC = 6, 16, 17, 18, 19 give us the small validation MSE; \n",
    "2) for all the PC dimensions,  validation MSE is higher than training MSE, is this normal? \n",
    "3) the training MSE keeps decreasing as PC number goes up, is this as expected?\n",
    "\n",
    "\n",
    "From the plot we also see that the cross-validation error is roughly the same when \n",
    "only one component is included in the model. \n",
    "This suggests that a model that uses just a small number of components might suffice. \n",
    "In the book, the authors used train/test to select the best dimension. \n",
    "I will skip that part since most of those were already covered in the previous sections.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let us move forwward with 6 PCs\n",
    "pcs = 6\n",
    "X_train, X_test , y_train, y_test = train_test_split(X, y, test_size=0.66)\n",
    "\n",
    "pca2 = PCA()\n",
    "# scale the data\n",
    "X_reduced_train = pca2.fit_transform(scale(X_train))\n",
    "X_reduced_test = pca2.transform(scale(X_test))[:,:pcs]\n",
    "# train regression model on training data \n",
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(X_reduced_train[:,:6], y_train)\n",
    "# prediction with test data\n",
    "pred = regr.predict(X_reduced_test)\n",
    "mean_squared_error(y_test, pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "We'll do a little math to get the amount of variance explained by adding each consecutive principal component. \n",
    "We can think of this as the amount of information(variance) about the data(X) or the response that is \n",
    "captured using $M$ principal components. \n",
    "For example, setting $M = 1$ only captures 38.31% of all the variance, or information, in the data.\n",
    "\"\"\"\n",
    "np.cumsum(pca.explained_variance_ratio_) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial Least Squares\n",
    "#### Scikit-learn has function PLSRegression for partial least squares regression. But, we still need to write a few line of codes to do the cross validation. The logic is same as before. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train, X_test , y_train, y_test = train_test_split(X, y, test_size=0.66)\n",
    "pls = PLSRegression(n_components=2)\n",
    "pls.fit(scale(X_train), y_train)\n",
    "\n",
    "mean_squared_error(y_test, pls.predict(scale(X_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 6"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "interpreter": {
   "hash": "4548a0e672c5b3a287feee7b2962606840aa548749d1830ef724408652b0c250"
  },
  "kernelspec": {
   "display_name": "Python 2.7.16 64-bit ('base': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}