curriculum/ml_intro/linear-regression-code-intro/intro-to-regression-starter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intro to Regression\n",
    "\n",
    "We're going through a few regression exercises that will help you see how regression works in Python.  \n",
    "\n",
    "## Objectives: \n",
    "At the end of this notebook the students should:\n",
    "- Be able to visualize data\n",
    "- Look for correlations and multicollinearity\n",
    "- Understand how linear regression models work\n",
    "- Interpret basic regression statistics like R^2\n",
    "- Do basic feature engineering and selection to improve models\n",
    "\n",
    "\n",
    "Be able to create linear regression in:\n",
    "- [***statsmodels***](http://statsmodels.sourceforge.net/): a package mainly best at doing regressions with traditional R formula syntax\n",
    "- [***scikit-learn***](http://scikit-learn.org/dev/index.html): This is the main machine learning package we'll be using throughout the course.  It has a multitude of machine learning algorithms and helpful machine learning pipeline tools.  sklearn has a tremendous amount of functionality, to get the most out of this course it will help to really explore the depth of the documentation on your own and watch as you understand more and more of the functionality as the course progresses.\n",
    "\n",
    "\n",
    "Gain familiarity with the following:\n",
    "- [***R formulas***](http://science.nature.nps.gov/im/datamgmt/statistics/r/formulas/): R formulas are a convenient way for encapsulating functional relationships for regressions\n",
    "- [***seaborn***](http://stanford.edu/~mwaskom/software/seaborn/): We'll use seaborn for **visualization** as we go along\n",
    "- [***Variable Preprocessing and Polynomial Regression***](http://scikit-learn.org/dev/modules/preprocessing.html#preprocessing) with scikit-learn:  We'll be **\"standardizing\"** or **\"normalizing\"** many of our variables to yield better model data.  We'll show how the \"linear\" models can be extended to basically any type of function by using functions of the different fields as the inputs to the linear model.\n",
    "\n",
    "## Datasets\n",
    "We'll take a look at a few different datasets:\n",
    "1. [Survey Responses](http://www.ats.ucla.edu/stat/examples/chp/p054.txt) (done together)\n",
    "2. Manufactured random dataset  (done together)\n",
    "3. [Car price predictor dataset](https://archive.ics.uci.edu/ml/datasets/Automobile) (exercise done in pairs)\n",
    "\n",
    "# Setup\n",
    "Try running the following imports and check that they all import successfully.  If they do not, install commands are given below. If necessary, at a command line window use `pip` to install the ones that are failing for you and then retry the imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:10:30.961980Z",
     "start_time": "2019-01-29T16:10:30.931958Z"
    }
   },
   "outputs": [],
   "source": [
    "# Necessary imports\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import patsy\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.linear_model import RidgeCV\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installations (if necessary)\n",
    "\n",
    "```` bash\n",
    "conda install pandas numpy statsmodels seaborn scikit-learn\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survey Data\n",
    "We will use this [simple survey data](http://www.ats.ucla.edu/stat/examples/chp/p054.txt) to demonstrate a few basic features of ***statsmodels*** and ***seaborn*** and how they might be used in a data science workflow for regression.\n",
    "\n",
    "The dataset is simply the results of a survey where the question responses are all numeric.  This leads to 6 numeric independent variable (predictor) fields and 1 numeric dependent variable (response) field.  The predictors are labeled ***X<sub>i</sub>*** and the response is labeled ***Y***.\n",
    "\n",
    "Let's load the dataset in using ***pandas*** and take a look at it.  Here we use [***pandas.read_table***](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_table.html) to load the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:11:12.287313Z",
     "start_time": "2019-01-29T16:11:11.922764Z"
    }
   },
   "outputs": [],
   "source": [
    "# Load data\n",
    "df = pd.read_table(\n",
    "    'https://stats.idre.ucla.edu/wp-content/uploads/2016/02/p054.txt')\n",
    "# Take a look at the datatypes\n",
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:11:13.382808Z",
     "start_time": "2019-01-29T16:11:13.341662Z"
    }
   },
   "outputs": [],
   "source": [
    "# Take a look at the first few rows\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look at the column names, we'll notice we have the trailing whitespace problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:11:20.316221Z",
     "start_time": "2019-01-29T16:11:20.305682Z"
    }
   },
   "outputs": [],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can remove this by calling map on the columns list and stripping the whitespace with strip.  The ***map*** function is applied to Series objects, whereas the ***apply*** and ***applymap*** functions are called on Dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:12:21.274828Z",
     "start_time": "2019-01-29T16:12:21.258283Z"
    }
   },
   "outputs": [],
   "source": [
    "df.columns = df.columns.map(str.strip)\n",
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:12:22.728544Z",
     "start_time": "2019-01-29T16:12:22.713083Z"
    }
   },
   "outputs": [],
   "source": [
    "# How many rows and columns does the dataset have?\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing with Seaborn\n",
    "We see that the data has 30 responses with 7 fields (6 independent, 1 dependent) each.  Let's use pandas to check out the correlations between the different variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:16:14.945316Z",
     "start_time": "2019-01-29T16:16:14.903888Z"
    }
   },
   "outputs": [],
   "source": [
    "# View the correlations\n",
    "df.corr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:16:32.651443Z",
     "start_time": "2019-01-29T16:16:31.098847Z"
    }
   },
   "outputs": [],
   "source": [
    "# example of a better corr matrix\n",
    "sns.heatmap(df.corr(), cmap=\"seismic\", annot=True, vmin=-1, vmax=1);\n",
    "# more cmaps: https://matplotlib.org/examples/color/colormaps_reference.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlation and Multicollinearity\n",
    "We notice that some of the variables are highly correlated.  This is something we might want to take into consideration later.  When 2 predictor variables are highly correlated this is called [multicollinearity](https://en.wikipedia.org/wiki/Multicollinearity) and it is something we want to watch out for as it can destabilize our model.  In the extreme case, when 2 predictors are perfectly correlated then there is absolutely nothing gained by making both variables part of our regression.\n",
    "\n",
    "The other takeaway from this table is that some of our predictors are highly correlated with our ***target variable Y***.  This is a good thing, it means that these are the variables that we most likely want to include as part of our model as they explain a large amount of the variance in the target variable (correlation=R, variance_explained=R<sup>2</sup>).\n",
    "\n",
    "Let's try to visualize these correlations all together by using the [***seaborn pairplot***](http://stanford.edu/~mwaskom/software/seaborn/generated/seaborn.pairplot.html) function.\n",
    "\n",
    "> What do you notice?\n",
    "\n",
    "> Almost all correlations are positive, somewhat normal distributions, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:19:20.146145Z",
     "start_time": "2019-01-29T16:19:04.991374Z"
    }
   },
   "outputs": [],
   "source": [
    "# Plot all of the variable-to-variable relations as scatterplots\n",
    "sns.pairplot(df, height=1.2, aspect=1.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ordinary Least Squares Regression with Statsmodels\n",
    "Now that we have a feel for our data, let's jump right in and try a basic regression model.  \n",
    "\n",
    "#### Statsmodels\n",
    "We are going to use the [**`statsmodels`**](http://statsmodels.sourceforge.net/) library first.  `statsmodels` is a Python package for implementing [**linear models**](https://en.wikipedia.org/wiki/Linear_model), of which **Linear Regression** is one.  It has a bunch of nice features for evaluating and executing such models.  Essentially, a linear model is one that is a **linear function of the parameters**.  For Linear Regression, this means the parameter $\\beta$ (note that here $\\beta$ is a **vector** of parameters, which includes $\\beta_0$, $\\beta_1$, $\\beta_2$, etc).  We'll discuss linear models generally later but for  now just accept that Linear Regression is one of these.\n",
    "\n",
    "#### Modeling with Statsmodels\n",
    "There are 2 main ways you can generate models with stats models:\n",
    "- Via the `statsmodels.api` package\n",
    "- Via the `statsmodels.formula.api` package\n",
    "\n",
    "For both approaches, you'll need somewhere to use the [R formula](http://science.nature.nps.gov/im/datamgmt/statistics/r/formulas/) styles formulas for defining the relationship between target variable and feature variables in your model.  ***Statsmodels*** uses [***patsy***](http://patsy.readthedocs.org/en/latest/) to convert this syntax into the proper data matrices for input into its linear models under the covers.  There are a variety of interactions and functions of variables that you can incorporate with this syntax, so feel free to check out the docs.\n",
    "\n",
    "Here we'll just start by defining a regression model that takes as its inputs each of the 6 predictor variables.  The other parameter of course is the data that the model is to be built from, our pandas dataframe.\n",
    "\n",
    "This first model fitting is done for you, it fits a multiple linear regression model of the following form (notice the use of [MathJax](https://www.mathjax.org/) for rendering such lovely math equations in markdown):\n",
    "\n",
    "$$\n",
    "\\widehat{Y} = \\beta_0 + \\beta_1X_1 + \\beta_2X_2 + \\beta_3X_3 + \\beta_4X_4 + \\beta_5X_5 + \\beta_6X_6\n",
    "$$\n",
    "\n",
    "##### `statsmodels.api`\n",
    "To use this method, you need to generate a **matrix** of **features**, **`X`** and a **vector** of **targets**, **`y`** where each row represents a single **observation**.  In statsmodels, you can do this with a call to **`patsy.dmatrices`**:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:28:09.616670Z",
     "start_time": "2019-01-29T16:28:09.592646Z"
    }
   },
   "outputs": [],
   "source": [
    "y, X = patsy.dmatrices('Y ~ X1 + X2 + X3 + X4 + X5 + X6', data=df, return_type=\"dataframe\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:28:52.012588Z",
     "start_time": "2019-01-29T16:28:51.948265Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create your feature matrix (X) and target vector (y)\n",
    "y, X = patsy.dmatrices('Y ~ X1 + X2 + X3 + X4 + X5 + X6', data=df, return_type=\"dataframe\")\n",
    "\n",
    "# Create your model\n",
    "model = sm.OLS(y, X)\n",
    "\n",
    "# Fit your model to your training set\n",
    "fit = model.fit()\n",
    "\n",
    "# Print summary statistics of the model's performance\n",
    "fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `statsmodels.formula.api`\n",
    "The formula approach handles the creation of the `X` and `y` matrices internally, so all you have to do is supply the R formula for your model when you create your `ols` object.  \n",
    "\n",
    "**NOTE:** We'll use this for the remainder of these exercises, but many people like to stick with the `dmatrices` and `X`, `y` matrix creation approach because that is the way `sklearn` works (see later).  Notice that here we've built the same model and it's deterministic, so the results are the same.\n",
    "\n",
    "[Nonrobust vs. robust regression](https://en.wikipedia.org/wiki/Robust_regression#Methods_for_robust_regression)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:38:58.049795Z",
     "start_time": "2019-01-29T16:38:57.959068Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define the model\n",
    "lm1 = smf.ols('Y ~ X1 + X2 + X3 + X4 + X5', data=df)\n",
    "\n",
    "# Fit the model\n",
    "fit1 = lm1.fit()\n",
    "\n",
    "# Print summary statistics of the model's performance\n",
    "fit1.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Statistics\n",
    "From this we get a handful of useful statistics describing our regression. A google search on statistics should yield details for those interested, and we will cover them in depth soon. We will focus on the ***R<sup>2</sup>*** values and the middle table here.\n",
    "\n",
    "***R<sup>2</sup>*** is the square of the correlation coefficient and represents the estimated percentage of the variance in our target variable ***Y*** that can be explained by our regression model.  ***Adjusted R<sup>2</sup>*** also penalizes for things such as large coefficients and extra variables to try and limit ***overfitting*** so it is often a better measure of model efficacy. We will use that here.\n",
    "\n",
    "The middle table provides the **coefficients** that are regression has found, along with the **standard error** for each coefficient. This defines our model, aka these are the model parameters that our algorithm was seeking to determine.  \n",
    "\n",
    "The **t-scores** are values that the coefficients score in the [Student's T Distribution](https://en.wikipedia.org/wiki/Student's_t-distribution) and the **P>|t|** field represents the probability of finding such a t-score if the actual value of the coefficient were 0. In other words, if we had a coefficient whose true value should be 0 (aka the predictor has no impact on the response) then this P-value is the probability of finding such a coefficient value in our regression by random chance. In essence, it measures our degree of belief that the coefficient for each variable should be zero. Thus, the lowest P-values represent the most likely predictors to be impacting the response."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Putting it all together, the final column returns a **95% Confidence Interval** for the value of each coefficient.\n",
    "\n",
    "Given these stats, lets remove the highest 3 P-values from our regression model, from ***X<sub>2</sub>***, ***X<sub>4</sub>***, and ***X<sub>5</sub>*** and see how our model performs now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:54:17.230096Z",
     "start_time": "2019-01-29T16:54:17.167701Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Define the model removing X2, X4, and X5\n",
    "lm2 = smf.ols('Y ~ X1 + X3 + X6', data=df)\n",
    "\n",
    "# Fit the model.  Let's call this fit2.\n",
    "fit2 = lm2.fit()\n",
    "\n",
    "# Check out the results\n",
    "fit2.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see our **Adjusted R<sup>2</sup>** has increased, and our P-values are lower so we likely have a better model.  Let's just try removing ***X<sub>6</sub>*** to see if that might improve our model a little bit more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:54:18.904112Z",
     "start_time": "2019-01-29T16:54:18.890436Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Define the model removing X6 this time\n",
    "\n",
    "\n",
    "# Fit the model.  Let's call this fit3.\n",
    "\n",
    "\n",
    "# Check out the results\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nope, both **R<sup>2</sup>** coefficients decreased so lets stick with the 2nd model.\n",
    "\n",
    "### Plotting Residuals\n",
    "Before we call it a day with this model and dataset, let's take a quick look at a plot of our residuals *(actual value - predicted value)* with this model.  We do this because in a good model we essentially want our errors to be random.  If our residuals look systematic (e.g. missing high for one range and low for another) then we probably are missing the actual functional dependency underlying the data (perhaps it's not really linear).\n",
    "\n",
    "Take a look [here](http://blog.minitab.com/blog/adventures-in-statistics/regression-analysis-how-do-i-interpret-r-squared-and-assess-the-goodness-of-fit) for an example of a bad residual plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:49:55.959948Z",
     "start_time": "2019-01-29T16:49:55.948885Z"
    }
   },
   "source": [
    "residuals $\\epsilon = y - \\hat{y}$\n",
    "\n",
    "linear regression equation $ \\hat{y} = \\beta_0 + \\beta_1 x_1 \\dots + \\epsilon$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:49:17.304975Z",
     "start_time": "2019-01-29T16:49:16.097702Z"
    }
   },
   "outputs": [],
   "source": [
    "# Use statsmodels to plot the residuals vs the fitted values\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.scatter(fit2.predict(), fit2.resid);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks pretty random! Let's move on to modeling functions with sklearn!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression with sklearn\n",
    "Statsmodels has decent functionality for linear models, and is great for statistical summaries. But, scikit-learn has more modeling options for all sorts of algorithms as well as data preparation and is growing every day, so we will generally be working with that from here on.\n",
    "\n",
    "### Regression with sklearn\n",
    "Before we jump into some of the additional features of sklearn, let's try to repeat our basic survey example using sklearn's built in **LinearRegression**.\n",
    "\n",
    "You should still have your Dataframe loaded from earlier.  Let's try repeating some of the different models we tried earlier with sklearn.  Here's the first model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T16:56:10.470147Z",
     "start_time": "2019-01-29T16:56:10.435527Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create an empty model\n",
    "lr = LinearRegression()\n",
    "\n",
    "# Choose the predictor variables, here all but the first which is the response variable\n",
    "# This model is analogous to the Y ~ X1 + X2 + X3 + X4 + X5 + X6 model\n",
    "X = df.iloc[:, 1:]\n",
    "\n",
    "# Choose the response variable(s)\n",
    "y = df.iloc[:, 0]\n",
    "\n",
    "# Fit the model to the full dataset\n",
    "lr.fit(X, y)\n",
    "\n",
    "# Print out the R^2 for the model against the full dataset\n",
    "lr.score(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that this is the same **R<sup>2</sup>** value that was reported for the first model above. \n",
    "\n",
    "## Exercise\n",
    "Let's quickly run the best model from earlier (***X1***, ***X3***, and ***X6***) and see how it performs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:00:27.476541Z",
     "start_time": "2019-01-29T17:00:27.422374Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Create an empty model\n",
    "\n",
    "\n",
    "# Choose the predictor variables, here all but the first which is the response variable\n",
    "# This model is analogous to the Y ~ X1 + X3 + X6 model\n",
    "\n",
    "\n",
    "# Choose the response variable(s)\n",
    "\n",
    "\n",
    "# Fit the model to the full dataset\n",
    "\n",
    "\n",
    "# Print out the R^2 for the model against the full dataset\n",
    "\n",
    "\n",
    "# Note:  adjusted_r_squared = 1 - (1-r_squared)*(len(y)-1)/(len(y)-X.shape[1]-1)\n",
    "\n",
    "# Create an empty model\n",
    "lr = LinearRegression()\n",
    "\n",
    "y, X = patsy.dmatrices(\"Y ~ X1 + X2 + X3 + X4 + X5 + X6\", data=df, return_type=\"dataframe\")\n",
    "\n",
    "# Fit the model to the full dataset\n",
    "lr.fit(X, y)\n",
    "\n",
    "# Print out the R^2 for the model against the full dataset\n",
    "lr.score(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the **R<sup>2</sup>** value is the same again.  It's slightly lower, the same as we saw earlier, but the **Adjusted R<sup>2</sup>** value that we saw earlier was higher for this one.  \n",
    "\n",
    "So we've seen how can do simple models with statsmodels and sklearn, but does anyone see an issue with what we've done here?  What data have we fit and scored our models on?  What claims can we make about the probable performance of our model going forward?\n",
    "\n",
    "## sklearn: What's in a model?\n",
    "Essentially all models in `sklearn` inherit from the same type of \"estimator\" interface.  This means they'll share common methods that we'll see over and over again like:\n",
    "- `fit()`: Fit a model to a set of training data\n",
    "- `score()`: Score the performance of a model on a given sample of data with known _ground truth_ dependent variables\n",
    "- `predict()`: Predict target/response variables based on a sample of independent variables (features, predictors, etc)\n",
    "\n",
    "Additionally, models are usually loaded with other goodies once they've been fit, which can provide information about the resulting fitted model.  For instance, the following might be of interest in regression:\n",
    "- `intercept_`: our $\\beta_0$ intercept in our regression model\n",
    "- `coef_`: the other $\\beta$s in our model\n",
    "\n",
    "## Exercise\n",
    "Let's print those out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:02:38.783935Z",
     "start_time": "2019-01-29T17:02:38.771010Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# print out intercept\n",
    "print(lr.intercept_)\n",
    "\n",
    "# print out other coefficients\n",
    "lr.coef_\n",
    "\n",
    "# note how much more elegant statsmodel renders results; yet, sklearn's raw values are easier to capture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Pickling for Later\n",
    "We're going to come back to this dataset!  So let's pickle it away for later...\n",
    "\n",
    "`sklearn`, `pandas` and `statsmodels` have their own respective methods for \"pickling\" their objects.\n",
    "\n",
    "**Pandas:**\n",
    "\n",
    "To pickle a `pandas.DataFrame` use the [`to_pickle()`](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.to_pickle.html) method.  Use this now to pickle your training data `df` now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:03:13.362958Z",
     "start_time": "2019-01-29T17:03:13.351290Z"
    }
   },
   "outputs": [],
   "source": [
    "# Pickle your pandas dataframe\n",
    "df.to_pickle('data/survey_data.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Statsmodels**\n",
    "\n",
    "To pickle a statsmodels object us the [`save()`](http://statsmodels.sourceforge.net/devel/generated/statsmodels.regression.linear_model.OLSResults.save.html#statsmodels.regression.linear_model.OLSResults.save) method.  Use this to pickle your best fit model `fit2` to `survey_sm_model.pkl` now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:03:48.309173Z",
     "start_time": "2019-01-29T17:03:48.300817Z"
    }
   },
   "outputs": [],
   "source": [
    "# Pickle fit2 to a file \n",
    "fit2.save('data/survey_sm_model.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**sklearn**\n",
    "\n",
    "Pickle an `sklearn` model to the file `survey_sk_model.pkl` using the `sklearn` replacement for pickle `joblib` as seen [here](http://scikit-learn.org/stable/modules/model_persistence.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:03:59.316430Z",
     "start_time": "2019-01-29T17:03:59.299785Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.externals import joblib\n",
    "joblib.dump(lr, 'data/survey_sk_model.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial Regression\n",
    "So far we've only tried to create regression models that are linear functions of the predictor variables.  However, there's no reason we can't transform the predictor variables by any type of function we want before inputting them to linear regression.  This is the idea behind [**Polynomial Regression**](https://en.wikipedia.org/wiki/Polynomial_regression) and it allows us (along with similar functional regressions) to essentially model our response variables as any function of our predictor variables that we like.  Viewed in this way, Linear Regression is just a special instance of Polynomial Regression with a polynomial of degree 1.\n",
    "\n",
    "## Polynomial Regression with sklearn\n",
    "sklearn has built-in options for converting your predictor variables to polynomial functions of them.  In this exercise we'll use the [**PolynomialFeatures**](http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html) class of sklearn to manipulate incoming predictors into nth-order polynomials of those features.  We'll combine this with the [***make_pipeline***](http://scikit-learn.org/stable/modules/generated/sklearn.pipeline.make_pipeline.html) function to string together a pipeline of operations that is able to first transform our linear features into polynomial features and then run a linear regression against the resulting polynomial features. \n",
    "\n",
    "### Generating Random Data\n",
    "The first thing we're going to do is manufacture some data from a known distribution with a little additive noise.  This allows us to compare our results to the known ground truth.  Let's create that data from a sine curve as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:05:34.603572Z",
     "start_time": "2019-01-29T17:05:33.887940Z"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.pylabtools import figsize\n",
    "figsize(5,5)\n",
    "plt.style.use('fivethirtyeight')\n",
    "\n",
    "# We start by seeding the random number generator so that everyone will have the same \"random\" results\n",
    "np.random.seed(9)\n",
    "\n",
    "# Function that returns the sin(2*pi*x)\n",
    "def f(x):\n",
    "    return np.sin(2 * np.pi * x)\n",
    "\n",
    "# generate points used to plot\n",
    "# This returns 100 evenly spaced numbers from 0 to 1\n",
    "x_plot = np.linspace(0, 1, 100)\n",
    "\n",
    "# generate points and keep a subset of them\n",
    "n_samples = 100\n",
    "# Generate the x values from the random uniform distribution between 0 and 1\n",
    "X = np.random.uniform(0, 1, size=n_samples)[:, np.newaxis]\n",
    "# Generate the y values by taking the sin and adding a random Gaussian (normal) noise term\n",
    "y = f(X) + np.random.normal(scale=0.3, size=n_samples)[:, np.newaxis]\n",
    "\n",
    "# Plot the training data against what we know to be the ground truth sin function\n",
    "fig,ax = plt.subplots(1,1);\n",
    "ax.plot(x_plot, f(x_plot), label='ground truth', color='green')\n",
    "ax.scatter(X, y, label='data', s=100)\n",
    "ax.set_ylim((-2, 2))\n",
    "ax.set_xlim((0, 1))\n",
    "ax.set_ylabel('y')\n",
    "ax.set_xlabel('x')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting an nth-degree Polynomial\n",
    "Now that we have our data and know the ground truth, let's try fitting a 3rd degree polynomial to our training data and see how it looks.  3rd degree makes sense for this interval because the sin function has 2 turning points over the interval [0,1] and a 3rd degree polynomial will general have 2 (or less) turning points.\n",
    "\n",
    "We first define a function `plot_approximation` that takes a pipeline of steps from make_pipeline and some plotting info and will plot the results of the sklearn pipeline on the specified plot with the ground truth and data in the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:05:36.430577Z",
     "start_time": "2019-01-29T17:05:36.387692Z"
    }
   },
   "outputs": [],
   "source": [
    "# import PolynomialFeatures and make_pipeline for Polynomial Regression\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "# Plot the results of a pipeline against ground truth and actual data\n",
    "def plot_approximation(est, ax, label=None):\n",
    "    \"\"\"Plot the approximation of ``est`` on axis ``ax``. \"\"\"\n",
    "    ax.plot(x_plot, f(x_plot), label='ground truth', color='green')\n",
    "    ax.scatter(X, y, s=100)\n",
    "    ax.plot(x_plot, est.predict(x_plot[:, np.newaxis]), color='red', label=label)\n",
    "    ax.set_ylim((-2, 2))\n",
    "    ax.set_xlim((0, 1))\n",
    "    ax.set_ylabel('y')\n",
    "    ax.set_xlabel('x')\n",
    "    ax.legend(loc='upper right',frameon=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's generate our pipeline for a 3rd degree polynomial and try it out in our plotting function.  Note that the steps are:\n",
    "- Use PolynomialFeatures(3) to create a generator of 3rd degree polynomials\n",
    "- Feed this generator to make_pipeline along with a LinearRegression object to tell it to string together these operations when given a new set of input predictor variables.  This results in a new model object that has the same `fit()`, `score()`, `predict()`, etc functions\n",
    "- Call `fit()` on our new object to fit a 3rd degree polynomial regression\n",
    "- Send the result to our plotting function to view the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:05:38.328895Z",
     "start_time": "2019-01-29T17:05:37.909638Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set up the plot\n",
    "fig,ax = plt.subplots(1,1)\n",
    "# Set the degree of our polynomial\n",
    "degree = 3\n",
    "\n",
    "# Generate the model type with make_pipeline\n",
    "# This tells it the first step is to generate 3rd degree polynomial features in the input features and then run\n",
    "# a linear regression on the resulting features\n",
    "est = make_pipeline(PolynomialFeatures(degree), LinearRegression())\n",
    "\n",
    "# Fit our model to the training data\n",
    "est.fit(X, y)\n",
    "# Plot the results\n",
    "plot_approximation(est, ax, label='degree=%d' % degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> How would you characterize this fit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the fit of a polynomial of degree 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:05:40.612733Z",
     "start_time": "2019-01-29T17:05:40.209182Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set up the plot\n",
    "fig,ax = plt.subplots(1,1)\n",
    "# Set the degree of our polynomial\n",
    "degree = 2\n",
    "# Generate the model type with make_pipeline\n",
    "# This tells it the first step is to generate 3rd degree polynomial features in the input features and then run\n",
    "# a linear regression on the resulting features\n",
    "est = make_pipeline(PolynomialFeatures(degree), LinearRegression())\n",
    "# Fit our model to the training data\n",
    "est.fit(X, y)\n",
    "# Plot the results\n",
    "plot_approximation(est, ax, label='degree=%d' % degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the fit of a polynomial of degree 9, or heck, how bout 27!?:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:05:42.912055Z",
     "start_time": "2019-01-29T17:05:42.461780Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set up the plot\n",
    "fig,ax = plt.subplots(1,1)\n",
    "# Set the degree of our polynomial\n",
    "degree = 27\n",
    "# Generate the model type with make_pipeline\n",
    "# This tells it the first step is to generate 3rd degree polynomial features in the input features and then run\n",
    "# a linear regression on the resulting features\n",
    "est = make_pipeline(PolynomialFeatures(degree), LinearRegression())\n",
    "# Fit our model to the training data\n",
    "est.fit(X, y)\n",
    "# Plot the results\n",
    "plot_approximation(est, ax, label='degree=%d' % degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> What happens as we increase the degree of polynomial?\n",
    "\n",
    "> Which polynomial should we choose?\n",
    "\n",
    "To gain some insight into this, let's plot polynomials from degree 1 to 9 and examine how the errors in our predictions change vs. the training and test sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:05:47.688273Z",
     "start_time": "2019-01-29T17:05:44.364748Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step through degrees from 0 to 9 and store the training and test (generalization) error.\n",
    "# This sets up 5 rows of 2 plots each (KEEP)\n",
    "fig, ax_rows = plt.subplots(5, 2, figsize=(15, 20))\n",
    "for degree in range(10):\n",
    "    est = make_pipeline(PolynomialFeatures(degree), LinearRegression())\n",
    "    est.fit(X, y)\n",
    "    # This sets the appropriate axis for each degree (KEEP)\n",
    "    ax_row_left, ax_row_right = ax_rows[degree//2]\n",
    "    if degree%2 == 0:\n",
    "        ax = ax_row_left\n",
    "    else:\n",
    "        ax = ax_row_right\n",
    "    plot_approximation(est, ax, label='degree=%d' % degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> What do you notice?\n",
    "\n",
    "### Pickling for later\n",
    "We may return to this fake data for further exploration later, so let's pickle our `X` and `y` so we can do just that.  To save `numpy` arrays, you use the [`numpy.save()`](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.save.html) method.  We can save multiple arrays to one file with the [`numpy.savez()`](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.savez.html) method.  Use this now to save `X` and `y` to the file `poly_data.npz`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pickle the numpy arrays X and y\n",
    "np.savez('data/poly_data.npz', X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises: Car Price Predictor Dataset\n",
    "For these exercises we'll be exploring the auto data available [here](https://archive.ics.uci.edu/ml/datasets/Automobile).  The goal is to be able to predict auto price from the \n",
    "\n",
    "## Data Exploration\n",
    "Use pandas `read_csv()` to load the data into a dataframe and then call `head()` to make sure everything looks good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:21:04.300159Z",
     "start_time": "2019-01-29T17:21:04.000907Z"
    }
   },
   "outputs": [],
   "source": [
    "# read in the car dataset\n",
    "df=pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/autos/imports-85.data',header=None)\n",
    "\n",
    "columns= ['symboling','normalized-losses','make','fuel-type','aspiration','num-of-doors','body-style','drive-wheels','engine-location','wheel-base','length','width','height','curb-weight','engine-type','num-of-cylinders','engine-size','fuel-system','bore','stroke','compression-ratio','horsepower','peak-rpm','city-mpg','highway-mpg','price']\n",
    "df.columns=columns\n",
    "# Use head to view the first few rows\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Use [`shape`](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.shape.html#pandas.DataFrame.shape) to check out how many rows and columns the dataframe has."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:24:02.314847Z",
     "start_time": "2019-01-29T17:24:02.302661Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# How many rows and columns do we have? \n",
    "\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use [`info()`](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.info.html#pandas.DataFrame.info) to get a summary of the dataframe and its datatypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:24:23.154665Z",
     "start_time": "2019-01-29T17:24:23.132400Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Let's examine the datatypes\n",
    "\n",
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Lets keep numeric data ** \n",
    "- Create a list of columns to keep\n",
    "- Select out only those columns from the dataframe and reassign the dataframe to that selection\n",
    "- Use `head()` & `info()` to make sure everything worked as expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:34:05.441449Z",
     "start_time": "2019-01-29T17:34:05.415446Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a list of columns to keep.  These are the numeric columns.\n",
    "subset=['wheel-base','length','width','height','curb-weight',\n",
    "        'engine-size','bore','stroke','compression-ratio',\n",
    "        'horsepower','peak-rpm','city-mpg','highway-mpg','price']\n",
    "\n",
    "# Now create a new dataframe named cars that only contains your column subset.\n",
    "small_df = df.loc[:, subset]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like some of our features (even our target feature :'price') is listed as an object.    \n",
    "Run cars.head(10), and we will see why!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:34:06.540706Z",
     "start_time": "2019-01-29T17:34:06.511444Z"
    }
   },
   "outputs": [],
   "source": [
    "# check data head\n",
    "small_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:35:26.972025Z",
     "start_time": "2019-01-29T17:35:26.927547Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION : replace \"?\" data in order to turn our features into numerics\n",
    "\n",
    "objects=['bore','stroke','horsepower','peak-rpm','price']\n",
    "\n",
    "for o in objects:\n",
    "    small_df[o]=small_df[o].replace('?', np.nan)\n",
    "    small_df[o]=small_df[o].astype(float)\n",
    "\n",
    "    \n",
    "small_df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the above output, that only a few entries were unknown in the first place.   \n",
    "To keep things simple for now, let's just go ahead and drop the entries that were unknown: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:36:47.238815Z",
     "start_time": "2019-01-29T17:36:47.219826Z"
    }
   },
   "outputs": [],
   "source": [
    "# drop nulls\n",
    "\n",
    "small_df.dropna(inplace=True)\n",
    "small_df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we begin modeling, use the [`corr()`](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.corr.html#pandas.DataFrame.corr) function to get a feel for the correlations among the different variables, especially with regard to 'price'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:42:11.444121Z",
     "start_time": "2019-01-29T17:42:11.438771Z"
    }
   },
   "outputs": [],
   "source": [
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:42:13.035774Z",
     "start_time": "2019-01-29T17:42:11.926399Z"
    }
   },
   "outputs": [],
   "source": [
    "# build correlations\n",
    "sns.heatmap(np.abs(small_df.corr()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at only the 'price' column of the correlations and order it in descending order wih [`sort_values()`](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.sort_values.html#pandas.DataFrame.sort_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-29T17:43:21.612208Z",
     "start_time": "2019-01-29T17:43:21.527478Z"
    }
   },
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Get the correlations with 'price' sorted in descending order\n",
    "small_df.corr()['price'].sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now have a better feel for which variables might be most valuable for your model.\n",
    "Q :  Do correlations provide the 'entire picture' of what is happening with our model? \n",
    "A:  Nope. It can give us an idea : but corrs will only provide the relationship with the response variable (all other factors being held constant) \n",
    "\n",
    "Now use ***seaborn's*** `pairplot()` function to visualize these correlations for the variables.  What do you think of our distributions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "#Let's try visualizing some of these pairwise correlations with seaborn.  Hint: Look at a subset of the columns.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Modeling with statsmodels\n",
    "Let's try some exploration with statsmodels.  As a first model, try creating an ordinary least squares model with statsmodels by incorporating all of the variables that had at least a .10 absolute value of correlation with cnt above \n",
    "- Create your model with the `ols()` function with the appropriate **R Formula** syntax and your dataframe\n",
    "- Fit the model\n",
    "- Print the fit summary to check out the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Let's jump right in and try a model with statsmodels using all variables above .10 correlation\n",
    "# lsm = smf.ols('price~ ....', data = cars)\n",
    "\n",
    "# You might have issues with some of the features .. Fix that\n",
    "## ex:\n",
    "cars.rename(\n",
    "    inplace=True,\n",
    "    columns={\n",
    "        \"curb-weight\": \"curb_weight\",\n",
    "        'engine-size': 'engine_size',\n",
    "        'wheel-base': 'wheel_base'\n",
    "    })\n",
    "\n",
    "# this also works\n",
    "#cars.rename(columns=lambda name: name.replace(\"-\", \"_\"), inplace=True)\n",
    "\n",
    "cars.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit and summarize in statsmodels\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Seaborn for Exploring Distributions\n",
    "Your **R<sup>2</sup>** should be .834, not bad.  That means we believe we can explain about 83.4% of the variance in price with this model.  \n",
    "\n",
    "One thing we will be discussing next week is Linear Regression assumptions and one being: normal distribution of the predictor variable.  Perhaps you noticed from our pairplot above that our 'price' variable is skewed.  Transform the y variable and rerun your OLS model.  Are there any other variables we should transform, why? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take log of price as target and graph.  Call this log_price.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# refit and summarize\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# log-transform the engine size.  Call this new column log_engine_size.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we see that the transformed engize-size values have a stronger correlation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adjust model to improve R^2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "#Let's try visualizing some of these pairwise correlations with seaborn\n",
    "sns.pairplot(cars[['log_price','log_engine_size', 'curb_weight', 'horsepower', 'width', 'length', 'wheel_base','bore']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What can we do with this observation?\n",
    "##### Indicator (Dummy) Variables\n",
    "As a first attempt, let's try to create [***indicator variables***](https://en.wikipedia.org/wiki/Dummy_variable_(statistics)) Indicator variables are 0/1 binary variables that indicate whether a condition is met or not, and they are quite useful in regression as they have an easy conceptual understanding: a value of 1 for a given field increases the target variable by the amount of its model coefficient.  These are often useful when provided with categorical attributes.  Any field with n unique categorical values can be reformulated into n indicator variable fields where each represents whether or not that attribute value is present.\n",
    "\n",
    "Let's see how to transform a categorical variable into a dummy variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's add 'make'\n",
    "print(df.make.value_counts())\n",
    "print('\\n')\n",
    "print(len(df.make.value_counts()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use patsy to convert the feature 'make' to categorical variables\n",
    "\n",
    "X=patsy.dmatrix('make',data=df,return_type='dataframe')\n",
    "X.head()\n",
    "\n",
    "# Q for Students: What do you notice about the returned matrix ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why are there three alfa romeos but more intercept 1 values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A:  A couple of things:  1) First of all, we've added an intercept (comes for free w/ patsy!)  \n",
    "#  2) 2nd: our favorite 'make' ('alfo-romero' of course) is missing from the above matrix\n",
    "\n",
    "## This is because patsy knows about (The Dummy Variable Trap)[http://www.algosome.com/articles/dummy-variable-trap-regression.html]\n",
    "# The idea behind DVT: Take a categorical variable that has two outcomes (example: boy & girl for the feature 'sex) - \n",
    "# We only need one column: \"Girl\", we can obviously solve for 'Boy' given this column -- if we had \n",
    "# both columns, we would automatically introduce multicolinearity  -- this idea can be extended to categorical features with \n",
    "# 20+ categories ! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STUDENT SECTION\n",
    "# Append the X 'make' matrix to your cars matrix.  Call your result cars2.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Cool thing to note: even though cars & X are different lengths (due to our previous .dropna()) ~ \n",
    "# joins allows us to 'merge' on their common index "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## For fun, let's see a quick way to throw cars2 into a OLS model,\n",
    "\n",
    "# sm (vs smf) allows for the following shortcut!\n",
    "import statsmodels.api as sm\n",
    "\n",
    "y = cars2.log_price\n",
    "x = cars2.drop(['log_price', 'price'], 1)\n",
    "\n",
    "lsm = sm.OLS(y,x)\n",
    "\n",
    "fit = lsm.fit()\n",
    "\n",
    "fit.summary()\n",
    "\n",
    "#  Looks like 'make' gives our R-squared a boost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  And this is just the beginning !\n",
    "## We'll pick up where we left off here soon :)\n",
    "\n",
    "# we can pickle our dataframe:\n",
    "cars2.to_pickle('data/cars_dataframe.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## What have we done?\n",
    "- Learned how to use statsmodels with R formula syntax for creating linear models\n",
    "- Learned how to evaluate models using adjusted R^2\n",
    "- Understood linear and polynomial regression in sklearn\n",
    "- Used seaborn for visualizing relationships in data\n",
    "- Used pandas for manipulating data as we move through our workflow\n",
    "- Gotten a peak into a genuine data science workflow\n",
    "- Seen how curiosity and creativity can yield big gains in a data science modeling pipeline\n",
    "\n",
    "## Play Time\n",
    "See if you can improve the model by trying out whatever methods you like!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}