LinearRegression_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear Regression Tutorial\n",
    "===\n",
    "\n",
    "Some problems don't have discrete (categorical) labels (e.g. color, plant species), but rather a continuous range of numbers (e.g. length, price). For these types of problems, regression is usually a good choice. Rather than predicting a categorical label for each example, it fits a continuous line (or plane, or curve) to the data in order to give a predicition as a number. \n",
    "\n",
    "If you've ever found a \"line of best fit\" using Excel, you've already used regression!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setup\n",
    "===\n",
    "Tell matplotlib to print figures in the notebook. Then import numpy (for numerical data), matplotlib.pyplot (for plotting figures), linear_model (for the scikit-learn linear regression algorithm), datasets (to download the Boston housing prices dataset from scikit-learn), and cross_validation (to create training and testing sets)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import linear_model, datasets # Import the linear regression function and dataset from scikit-learn\n",
    "from sklearn.model_selection import train_test_split, KFold\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "# Print figures in the notebook\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Import the dataset\n",
    "===\n",
    "Import the dataset and store it to a variable called diabetes. This dataset is similar to a python dictionary, with the keys: ['DESCR', 'target', 'data', 'feature_names']\n",
    "\n",
    "The data features are stored in diabetes.data, where each row is an example from a single patient, and each column is a single feature. The feature names are stored in diabetes.feature_names. Target values are stored in diabetes.target.\n",
    "\n",
    "Below, we load the labels into y, the data into X, and the names of the features into featureNames. We also print the description of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import some data to play with\n",
    "diabetes = datasets.load_diabetes()\n",
    "\n",
    "# Store the labels (y), features (X), and feature names\n",
    "y = diabetes.target       # Labels are stored in y as numbers\n",
    "X = diabetes.data\n",
    "featureNames = diabetes.feature_names\n",
    "\n",
    "print(diabetes.DESCR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create Training and Testing Sets\n",
    "---\n",
    "\n",
    "In order to see how well our classifier works, we need to divide our data into training and testing sets"
   ]
  },
  {
   "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.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize The Data\n",
    "===\n",
    "\n",
    "There are too many features to visualize the whole training set, but we can plot a single feature (e.g. body mass index) against the quantified disease progression.\n",
    "\n",
    "Remember that the features have been normalized to center it around 0 and scaled by the standard deviation. So the values shown aren't the actual BMI levels, but a standardized version of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bmi_ind = 2\n",
    "\n",
    "plt.scatter(X_train[:,bmi_ind], y_train)\n",
    "plt.ylabel('Disease Progression')\n",
    "plt.xlabel('Body Mass Index')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train A Toy Model\n",
    "===\n",
    "Next, we train a linear regression algorithm on our data. \n",
    "\n",
    "In scikit learn, we use [linear_model.LinearRegression()](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html) to create our model. Models are trained using a `fit()` method. For linear regression, `fit()` expects two arguments: the training examples `X_train` and the corresponding labels `y_train`.\n",
    "\n",
    "To start, we will train a toy model using only the body mass index feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "x_train = X_train[:,bmi_ind].reshape(-1, 1) # regression expects a (#examples,#features) array shape\n",
    "regr.fit(x_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Making Predictions\n",
    "---\n",
    "\n",
    "To make predictions on new data we use the `predict()` method. It expects a single input: an array-like object containing the features for the examples we want to predict.\n",
    "\n",
    "Here, we get the predictions for two body mass index values: -.05 and 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bmi = [[-0.05],[0.1]]\n",
    "predictions = regr.predict(bmi)\n",
    "print(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Visualize the Toy Model\n",
    "---\n",
    "\n",
    "Here we plot the linear regression line  of our trained model on top of the data. We do this by predicting the output of the model on a range of values from -0.1 to 0.15. These predictions are plotted as a line on top of the training data. \n",
    "\n",
    "This can't tell us how well it will perform on new, unseen, data, but it can show us how well the line fits the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = np.arange(-0.1, 0.15, 0.01).reshape(-1, 1) # Reshape so each feature is a separate row\n",
    "\n",
    "plt.scatter(x_train, y_train)\n",
    "plt.plot(values, regr.predict(values), c='r')\n",
    "plt.ylabel('Disease Progression')\n",
    "plt.xlabel('Body Mass Index')\n",
    "plt.title('Regression Line on Training Data')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the Toy Model\n",
    "===\n",
    "\n",
    "Next, we will test the ability of our model to predict the disease progression in our test set, using only the body mass index.\n",
    "\n",
    "First, we get our predictions for the training data, and plot the predicted model on top of the test data. Since we trained our model using only one feature, we need to get only that feature from the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_test = X_test[:,bmi_ind][np.newaxis].T # regression expects a (#examples,#features) array shape\n",
    "predictions = regr.predict(x_test)\n",
    "\n",
    "plt.scatter(x_test, y_test)\n",
    "plt.plot(x_test, predictions, c='r')\n",
    "plt.ylabel('Disease Progression')\n",
    "plt.xlabel('Body Mass Index')\n",
    "plt.title('Regression Line on Test Data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the Toy Model\n",
    "---\n",
    "\n",
    "Next, we evaluate how well our model worked on the training dataset. Unlike with discrete classifiers (e.g. KNN, SVM), the number of examples it got \"correct\" isn't meaningful here. We certainly care if the predicted value is off by 100, but do we care as much if it is off by 1? What about 0.01?\n",
    "\n",
    "There are many ways to evaluate a linear classifier, but one popular way is the mean-squared error, or MSE. As the name implies, you find the error for each example (the distance between the point and the predicted line), square each of them, and then add them all together. \n",
    "\n",
    "Scikit-learn has a function that does this for you easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mse = mean_squared_error(y_test, predictions)\n",
    "\n",
    "print('The MSE is ' + '%.2f' % mse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MSE isn't as intuitive as the accuracy of a discrete classifier, but it is highly useful for comparing the effectiveness of different models. Another option is to look at the $R^2$ score, which you may already be familiar with if you've ever fit a line to data in Excel. A value of 1.0 is a perfect predictor, while 0.0 means there is no correlation between the input and output of the regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2score = r2_score(y_test, predictions)\n",
    "\n",
    "print('The R^2 score is ' + '%.2f' % r2score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train A Model on All Features\n",
    "===\n",
    "\n",
    "Next we will train a model on all of the available features and use it to predict the progression of diabetes after a year. We can then see how this compares to using only a single feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(X_train, y_train)\n",
    "\n",
    "predictions = regr.predict(X_test)\n",
    "\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "print('The MSE is ' + '%.2f' % mse)\n",
    "\n",
    "r2score = r2_score(y_test, predictions)\n",
    "print('The R^2 score is ' + '%.2f' % r2score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Using Crossvalidation\n",
    "===\n",
    "\n",
    "To properly compare these two models, we need to use crossvalidation to select the best model. We can then get our final result using our test data.\n",
    "\n",
    "First we create our two linear regression models, then we divide our training data into folds. We loop through the sets of training and validation folds. Each time, we train each model on the training data and evaluate on the validation data. We store the MSE and $R^2$ score of each classifier on each fold so we can look at them later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create our regression models\n",
    "regr_1 = linear_model.LinearRegression()\n",
    "regr_all = linear_model.LinearRegression()\n",
    "\n",
    "# Create arrays to store the MSE and R^2 score\n",
    "mse_1 = []\n",
    "r2score_1 = []\n",
    "mse_all = []\n",
    "r2score_all = []\n",
    "\n",
    "# Loop through 5 folds\n",
    "kf = KFold(n_splits=5)\n",
    "for trainInd, valInd in kf.split(X_train):\n",
    "    X_tr = X_train[trainInd,:]\n",
    "    y_tr = y_train[trainInd]\n",
    "    X_val = X_train[valInd,:]\n",
    "    y_val = y_train[valInd]\n",
    "    \n",
    "    # Train our models\n",
    "    regr_1.fit(X_tr[:,bmi_ind].reshape(-1, 1), y_tr) # Train on only one feature\n",
    "    regr_all.fit(X_tr, y_tr) # Train on all features\n",
    "    \n",
    "    # Make our predictions\n",
    "    pred_1 = regr_1.predict(X_val[:,bmi_ind].reshape(-1, 1))\n",
    "    pred_all = regr_all.predict(X_val)\n",
    "    \n",
    "    # Calculate the MSE\n",
    "    mse_1.append(mean_squared_error(y_val, pred_1))\n",
    "    mse_all.append(mean_squared_error(y_val, pred_all))\n",
    "    \n",
    "    # Calculate the R^2 score\n",
    "    r2score_1.append(r2_score(y_val, pred_1))\n",
    "    r2score_all.append(r2_score(y_val, pred_all))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a Model\n",
    "---\n",
    "\n",
    "To select a model, we look at the average $R^2$ score and MSE across all folds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('One Feature:\\nMSE: ' + '%.2f' % np.mean(mse_1))\n",
    "print('R^2: ' + '%.2f' % np.mean(r2score_1))\n",
    "\n",
    "print('\\nAll Features:\\nMSE: ' + '%.2f' % np.mean(mse_all))\n",
    "print('R^2: ' + '%.2f' % np.mean(r2score_all))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Final Evaluation\n",
    "---\n",
    "\n",
    "The model using all features has a higher $R^2$ score and lower MSE, so we select it as our best model. Now we can evaluate it on our training set and get our final MSE an $R^2$ score values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regr_all.fit(X_train, y_train)\n",
    "predictions = regr_all.predict(X_test)\n",
    "\n",
    "mse = mean_squared_error(y_test, predictions)\n",
    "r2score = r2_score(y_test, predictions)\n",
    "\n",
    "print('The final MSE is ' + '%.2f' % mse)\n",
    "print('The final R^2 score is ' + '%.2f' % r2score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sidenote: Randomness and Results\n",
    "===\n",
    "\n",
    "Every time you run this notebook, you will get slightly different results. Why? Because data is randomly divided among the training/testing/validation data sets. Running the code again will create a different division of the data, and will make the results slightly different. However, the overall outcome should remain consistent and have approximately the same values. If you have drastically different results when running an analysis multiple times, it suggests a problem with your model or that you need more data.\n",
    "\n",
    "If it's important that you get the exact same results every time you run the code, you can specify a random state in the `random_state` argument of `train_test_split()` and `KFold`."
   ]
  },
  {
   "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.3, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "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
}