08_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Linear Regression\n",
    "\n",
    "*Adapted from Chapter 3 of [An Introduction to Statistical Learning](http://www-bcf.usc.edu/~gareth/ISL/)*\n",
    "\n",
    "||continuous|categorical|\n",
    "|---|---|---|\n",
    "|**supervised**|**regression**|classification|\n",
    "|**unsupervised**|dimension reduction|clustering|\n",
    "\n",
    "## Motivation\n",
    "\n",
    "Why are we learning linear regression?\n",
    "- widely used\n",
    "- runs fast\n",
    "- easy to use (not a lot of tuning required)\n",
    "- highly interpretable\n",
    "- basis for many other methods\n",
    "\n",
    "## Libraries\n",
    "\n",
    "Will be using [Statsmodels](http://statsmodels.sourceforge.net/) for **teaching purposes** since it has some nice characteristics for linear modeling. However, we recommend that you spend most of your energy on [scikit-learn](http://scikit-learn.org/stable/) since it provides significantly more useful functionality for machine learning in general."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# this allows plots to appear directly in the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Advertising Data\n",
    "\n",
    "Let's take a look at some data, ask some questions about that data, and then use linear regression to answer those questions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>radio</th>\n",
       "      <th>newspaper</th>\n",
       "      <th>sales</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>230.1</td>\n",
       "      <td>37.8</td>\n",
       "      <td>69.2</td>\n",
       "      <td>22.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>44.5</td>\n",
       "      <td>39.3</td>\n",
       "      <td>45.1</td>\n",
       "      <td>10.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>17.2</td>\n",
       "      <td>45.9</td>\n",
       "      <td>69.3</td>\n",
       "      <td>9.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>151.5</td>\n",
       "      <td>41.3</td>\n",
       "      <td>58.5</td>\n",
       "      <td>18.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>180.8</td>\n",
       "      <td>10.8</td>\n",
       "      <td>58.4</td>\n",
       "      <td>12.9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  radio  newspaper  sales\n",
       "1  230.1   37.8       69.2   22.1\n",
       "2   44.5   39.3       45.1   10.4\n",
       "3   17.2   45.9       69.3    9.3\n",
       "4  151.5   41.3       58.5   18.5\n",
       "5  180.8   10.8       58.4   12.9"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read data into a DataFrame\n",
    "data = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the **features**?\n",
    "- TV: advertising dollars spent on TV for a single product in a given market (in thousands of dollars)\n",
    "- Radio: advertising dollars spent on Radio\n",
    "- Newspaper: advertising dollars spent on Newspaper\n",
    "\n",
    "What is the **response**?\n",
    "- Sales: sales of a single product in a given market (in thousands of widgets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(200, 4)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the shape of the DataFrame\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 200 **observations**, and thus 200 markets in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'Sales'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/core/indexes/base.py\u001b[0m in \u001b[0;36mget_loc\u001b[0;34m(self, key, method, tolerance)\u001b[0m\n\u001b[1;32m   2441\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2442\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2443\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'Sales'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-1feceee91561>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# visualize the relationship between the features and the response using scatterplots\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msharey\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkind\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'scatter'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'TV'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Sales'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m16\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkind\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'scatter'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Radio'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Sales'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkind\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'scatter'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Newspaper'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Sales'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/plotting/_core.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, x, y, kind, ax, subplots, sharex, sharey, layout, figsize, use_index, title, grid, legend, style, logx, logy, loglog, xticks, yticks, xlim, ylim, rot, fontsize, colormap, table, yerr, xerr, secondary_y, sort_columns, **kwds)\u001b[0m\n\u001b[1;32m   2625\u001b[0m                           \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfontsize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolormap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcolormap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtable\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtable\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2626\u001b[0m                           \u001b[0myerr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0myerr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxerr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mxerr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msecondary_y\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msecondary_y\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2627\u001b[0;31m                           sort_columns=sort_columns, **kwds)\n\u001b[0m\u001b[1;32m   2628\u001b[0m     \u001b[0m__call__\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__doc__\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplot_frame\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__doc__\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2629\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/plotting/_core.py\u001b[0m in \u001b[0;36mplot_frame\u001b[0;34m(data, x, y, kind, ax, subplots, sharex, sharey, layout, figsize, use_index, title, grid, legend, style, logx, logy, loglog, xticks, yticks, xlim, ylim, rot, fontsize, colormap, table, yerr, xerr, secondary_y, sort_columns, **kwds)\u001b[0m\n\u001b[1;32m   1867\u001b[0m                  \u001b[0myerr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0myerr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxerr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mxerr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1868\u001b[0m                  \u001b[0msecondary_y\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msecondary_y\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msort_columns\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msort_columns\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1869\u001b[0;31m                  **kwds)\n\u001b[0m\u001b[1;32m   1870\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1871\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/plotting/_core.py\u001b[0m in \u001b[0;36m_plot\u001b[0;34m(data, x, y, subplots, ax, kind, **kwds)\u001b[0m\n\u001b[1;32m   1650\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDataFrame\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1651\u001b[0m             plot_obj = klass(data, x=x, y=y, subplots=subplots, ax=ax,\n\u001b[0;32m-> 1652\u001b[0;31m                              kind=kind, **kwds)\n\u001b[0m\u001b[1;32m   1653\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1654\u001b[0m             raise ValueError(\"plot kind %r can only be used for data frames\"\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/plotting/_core.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, data, x, y, s, c, **kwargs)\u001b[0m\n\u001b[1;32m    808\u001b[0m             \u001b[0;31m# the handling of this argument later\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    809\u001b[0m             \u001b[0ms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 810\u001b[0;31m         \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mScatterPlot\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    811\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mis_integer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mholds_integer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    812\u001b[0m             \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/plotting/_core.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, data, x, y, **kwargs)\u001b[0m\n\u001b[1;32m    784\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_numeric_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    785\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_kind\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m' requires x column to be numeric'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 786\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_numeric_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    787\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_kind\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m' requires y column to be numeric'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    788\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   1962\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_multilevel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1963\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1964\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1965\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1966\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_getitem_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_getitem_column\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   1969\u001b[0m         \u001b[0;31m# get column\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1970\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_unique\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1971\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_item_cache\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1972\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1973\u001b[0m         \u001b[0;31m# duplicate columns & possible reduce dimensionality\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36m_get_item_cache\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   1643\u001b[0m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1644\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mres\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1645\u001b[0;31m             \u001b[0mvalues\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_data\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1646\u001b[0m             \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_box_item_values\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1647\u001b[0m             \u001b[0mcache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/core/internals.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self, item, fastpath)\u001b[0m\n\u001b[1;32m   3588\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3589\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misnull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3590\u001b[0;31m                 \u001b[0mloc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3591\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3592\u001b[0m                 \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0misnull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/pandas/core/indexes/base.py\u001b[0m in \u001b[0;36mget_loc\u001b[0;34m(self, key, method, tolerance)\u001b[0m\n\u001b[1;32m   2442\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2443\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2444\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_cast_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2445\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2446\u001b[0m         \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtolerance\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtolerance\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'Sales'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADIlJREFUeJzt3V+Ipfddx/H3p7uuhVgbMCOU3dWkuDUuUkg6xEBBI62w\nycXuhUV2odRI7FI0etEiRCpR1guxvSgUV+uqobZg0jQXusqWBTVSEBN3QtuYTViZrn922ECmachN\nMenC14s5SQ+zMztn5zyTPTnf9wsG5nnOb878nn7hnbPnzDlNVSFJmn/vuNEbkCS9NQy+JDVh8CWp\nCYMvSU0YfElqwuBLUhNbBj/JI0leSvLcJrcnyeeTLCd5Nsmdw29TkjStSR7hfxE4dI3b7wUOjL6O\nA382/bYkSUPbMvhV9XXgu9dYcgT4Uq15Crg5yXuG2qAkaRi7B7iPvcClseOV0bkX1y9Mcpy1fwVw\n0003feD2228f4NdrWs8888x3qmphuz/vXGeTc51P08x1iOBng3Mbfl5DVZ0CTgEsLi7W0tLSAL9e\n00ryP9P8vHOdTc51Pk0z1yH+SmcF2D92vA+4PMD9SpIGNETwTwMfG/21zt3Aq1V11dM5kqQba8un\ndJI8CtwD3JJkBfh94IcAquoLwBngPmAZ+B7wazu1WUnS9m0Z/Ko6tsXtBfzmYDuSJO0I32krSU0Y\nfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYM\nviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMG\nX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSExMFP8mhJBeSLCd5aIPbfyLJk0m+keTZ\nJPcNv1VJ0jS2DH6SXcBJ4F7gIHAsycF1y34PeLyq7gCOAn869EYlSdOZ5BH+XcByVV2sqteBx4Aj\n69YU8KOj798NXB5ui5KkIUwS/L3ApbHjldG5cX8AfDTJCnAG+K2N7ijJ8SRLSZZWV1e3sV3NIuc6\nn5zr/Jkk+NngXK07PgZ8sar2AfcBX05y1X1X1amqWqyqxYWFhevfrWaSc51PznX+TBL8FWD/2PE+\nrn7K5gHgcYCq+jfgncAtQ2xQkjSMSYJ/DjiQ5LYke1h7Ufb0ujX/C3wIIMnPsBZ8/w0oSTNky+BX\n1RXgQeAs8AJrf41zPsmJJIdHyz4FfDzJt4BHgfurav3TPpKkG2j3JIuq6gxrL8aOn3t47PvngQ8O\nuzVJ0pB8p60kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJ\nasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4k\nNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqYmJgp/kUJILSZaTPLTJml9J\n8nyS80n+ZthtSpKmtXurBUl2ASeBXwJWgHNJTlfV82NrDgC/C3ywql5J8uM7tWFJ0vZM8gj/LmC5\nqi5W1evAY8CRdWs+DpysqlcAquqlYbcpSZrWJMHfC1waO14ZnRv3PuB9Sf41yVNJDm10R0mOJ1lK\nsrS6urq9HWvmONf55FznzyTBzwbnat3xbuAAcA9wDPjLJDdf9UNVp6pqsaoWFxYWrnevmlHOdT45\n1/kzSfBXgP1jx/uAyxus+buq+n5V/RdwgbX/AEiSZsQkwT8HHEhyW5I9wFHg9Lo1fwv8IkCSW1h7\niufikBuVJE1ny+BX1RXgQeAs8ALweFWdT3IiyeHRsrPAy0meB54EfqeqXt6pTUuSrt+Wf5YJUFVn\ngDPrzj089n0Bnxx9SZJmkO+0laQmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYM\nviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMG\nX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqYmJ\ngp/kUJILSZaTPHSNdR9JUkkWh9uiJGkIWwY/yS7gJHAvcBA4luTgBuveBfw28PTQm5QkTW+SR/h3\nActVdbGqXgceA45ssO4Pgc8A/zfg/iRJA5kk+HuBS2PHK6Nzb0pyB7C/qv7hWneU5HiSpSRLq6ur\n171ZzSbnOp+c6/yZJPjZ4Fy9eWPyDuBzwKe2uqOqOlVVi1W1uLCwMPkuNdOc63xyrvNnkuCvAPvH\njvcBl8eO3wX8LPAvSf4buBs47Qu3kjRbJgn+OeBAktuS7AGOAqffuLGqXq2qW6rq1qq6FXgKOFxV\nSzuyY0nStmwZ/Kq6AjwInAVeAB6vqvNJTiQ5vNMblCQNY/cki6rqDHBm3bmHN1l7z/TbkiQNzXfa\nSlITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHw\nJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4\nktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KamCj4SQ4luZBkOclDG9z+ySTPJ3k2yT8l\n+cnhtypJmsaWwU+yCzgJ3AscBI4lObhu2TeAxap6P/AE8JmhNypJms4kj/DvApar6mJVvQ48BhwZ\nX1BVT1bV90aHTwH7ht2mJGlakwR/L3Bp7HhldG4zDwBf2+iGJMeTLCVZWl1dnXyXmmnOdT451/kz\nSfCzwbnacGHyUWAR+OxGt1fVqaparKrFhYWFyXepmeZc55NznT+7J1izAuwfO94HXF6/KMmHgU8D\nv1BVrw2zPUnSUCZ5hH8OOJDktiR7gKPA6fEFSe4A/hw4XFUvDb9NSdK0tgx+VV0BHgTOAi8Aj1fV\n+SQnkhweLfss8CPAV5N8M8npTe5OknSDTPKUDlV1Bjiz7tzDY99/eOB9SZIG5jttJakJgy9JTRh8\nSWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+\nJDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZf\nkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJamKi4Cc5lORCkuUkD21w+w8n+cro9qeT3Dr0RiVJ\n09ky+El2ASeBe4GDwLEkB9ctewB4pap+Cvgc8MdDb1SSNJ1JHuHfBSxX1cWqeh14DDiybs0R4K9H\n3z8BfChJhtumJGlauydYsxe4NHa8AvzcZmuq6kqSV4EfA74zvijJceD46PC1JM9tZ9Mz5BbWXePb\n1E9P88POdWY516vNw2y3PddJgr/RI/Xaxhqq6hRwCiDJUlUtTvD7Z9Y8XAOsXcc0P+9cZ5Nzvdo8\nXMc0c53kKZ0VYP/Y8T7g8mZrkuwG3g18d7ubkiQNb5LgnwMOJLktyR7gKHB63ZrTwK+Ovv8I8M9V\nddUjfEnSjbPlUzqj5+QfBM4Cu4BHqup8khPAUlWdBv4K+HKSZdYe2R+d4HefmmLfs2IergGGvY55\n+N9kHq4BnOtG5uE6tn0N8YG4JPXgO20lqQmDL0lN7Hjw5+FjGSa4hvuTrCb55ujr12/EPq8lySNJ\nXtrsb6mz5vOja3w2yZ1b3J9znQHO9WrO9Rqqase+WHuR99vAe4E9wLeAg+vW/AbwhdH3R4Gv7OSe\nduga7gf+5EbvdYvr+HngTuC5TW6/D/gaa++puBt42rk6V+f69p/r+NdOP8Kfh49lmOQaZl5VfZ1r\nvzfiCPClWvMUcHOS92yy1rnOCOd6Fed6DTsd/I0+lmHvZmuq6grwxscyzIpJrgHgl0f/tHoiyf4N\nbp91k17npGud62xwrs71TTsd/ME+luEGmmR/fw/cWlXvB/6RHzwCeju5njk417cP5+pc37TTwZ+H\nj2XY8hqq6uWqem10+BfAB96ivQ1pklldz1rnOhucq3N9004Hfx4+lmHLa1j33Nlh4IW3cH9DOQ18\nbPTq/93Aq1X14iZrnevbh3N1rj/wFrzafB/wn6y9cv7p0bkTwOHR9+8EvgosA/8OvPdGv0K+jWv4\nI+A8a38R8CRw+43e8wbX8CjwIvB91h4dPAB8AvjE6Paw9n90823gP4BF5+pcnet8zPWNLz9aQZKa\n8J22ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhP/D1waPo5c5UySAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113f74a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize the relationship between the features and the response using scatterplots\n",
    "fig, axs = plt.subplots(1, 3, sharey=True)\n",
    "data.plot(kind='scatter', x='TV', y='Sales', ax=axs[0], figsize=(16, 8))\n",
    "data.plot(kind='scatter', x='Radio', y='Sales', ax=axs[1])\n",
    "data.plot(kind='scatter', x='Newspaper', y='Sales', ax=axs[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions About the Advertising Data\n",
    "\n",
    "Let's pretend you work for the company that manufactures and markets this widget. The company might ask you the following: On the basis of this data, how should we spend our advertising money in the future?\n",
    "\n",
    "This general question might lead you to more specific questions:\n",
    "1. Is there a relationship between ads and sales?\n",
    "2. How strong is that relationship?\n",
    "3. Which ad types contribute to sales?\n",
    "4. What is the effect of each ad type of sales?\n",
    "5. Given ad spending in a particular market, can sales be predicted?\n",
    "\n",
    "We will explore these questions below!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Linear Regression\n",
    "\n",
    "Simple linear regression is an approach for predicting a **quantitative response** using a **single feature** (or \"predictor\" or \"input variable\"). It takes the following form:\n",
    "\n",
    "$y = \\beta_0 + \\beta_1x$\n",
    "\n",
    "What does each term represent?\n",
    "- $y$ is the response\n",
    "- $x$ is the feature\n",
    "- $\\beta_0$ is the intercept\n",
    "- $\\beta_1$ is the coefficient for x\n",
    "\n",
    "Together, $\\beta_0$ and $\\beta_1$ are called the **model coefficients**. To create your model, you must \"learn\" the values of these coefficients. And once we've learned these coefficients, we can use the model to predict Sales!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating (\"Learning\") Model Coefficients\n",
    "\n",
    "Generally speaking, coefficients are estimated using the **least squares criterion**, which means we are find the line (mathematically) which minimizes the **sum of squared residuals** (or \"sum of squared errors\"):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"08_estimating_coefficients.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What elements are present in the diagram?\n",
    "- The black dots are the **observed values** of x and y.\n",
    "- The blue line is our **least squares line**.\n",
    "- The red lines are the **residuals**, which are the distances between the observed values and the least squares line.\n",
    "\n",
    "How do the model coefficients relate to the least squares line?\n",
    "- $\\beta_0$ is the **intercept** (the value of $y$ when $x$=0)\n",
    "- $\\beta_1$ is the **slope** (the change in $y$ divided by change in $x$)\n",
    "\n",
    "Here is a graphical depiction of those calculations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"08_slope_intercept.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use **Statsmodels** to estimate the model coefficients for the advertising data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "PatsyError",
     "evalue": "Error evaluating factor: NameError: name 'Sales' is not defined\n    Sales ~ TV\n    ^^^^^",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/compat.py\u001b[0m in \u001b[0;36mcall_and_wrap_exc\u001b[0;34m(msg, origin, f, *args, **kwargs)\u001b[0m\n\u001b[1;32m    116\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 117\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    118\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/eval.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(self, expr, source_name, inner_namespace)\u001b[0m\n\u001b[1;32m    165\u001b[0m         return eval(code, {}, VarLookupDict([inner_namespace]\n\u001b[0;32m--> 166\u001b[0;31m                                             + self._namespaces))\n\u001b[0m\u001b[1;32m    167\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<string>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'Sales' is not defined",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mPatsyError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-f9c0a330dd35>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m# create a fitted model in one line\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mlm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msmf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mols\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mformula\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Sales ~ TV'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;31m# print the coefficients\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/statsmodels/base/model.py\u001b[0m in \u001b[0;36mfrom_formula\u001b[0;34m(cls, formula, data, subset, drop_cols, *args, **kwargs)\u001b[0m\n\u001b[1;32m    153\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    154\u001b[0m         tmp = handle_formula_data(data, None, formula, depth=eval_env,\n\u001b[0;32m--> 155\u001b[0;31m                                   missing=missing)\n\u001b[0m\u001b[1;32m    156\u001b[0m         \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mendog\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexog\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmissing_idx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdesign_info\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtmp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    157\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/statsmodels/formula/formulatools.py\u001b[0m in \u001b[0;36mhandle_formula_data\u001b[0;34m(Y, X, formula, depth, missing)\u001b[0m\n\u001b[1;32m     63\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdata_util\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_is_using_pandas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     64\u001b[0m             result = dmatrices(formula, Y, depth, return_type='dataframe',\n\u001b[0;32m---> 65\u001b[0;31m                                NA_action=na_action)\n\u001b[0m\u001b[1;32m     66\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m             result = dmatrices(formula, Y, depth, return_type='dataframe',\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/highlevel.py\u001b[0m in \u001b[0;36mdmatrices\u001b[0;34m(formula_like, data, eval_env, NA_action, return_type)\u001b[0m\n\u001b[1;32m    308\u001b[0m     \u001b[0meval_env\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mEvalEnvironment\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcapture\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0meval_env\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreference\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    309\u001b[0m     (lhs, rhs) = _do_highlevel_design(formula_like, data, eval_env,\n\u001b[0;32m--> 310\u001b[0;31m                                       NA_action, return_type)\n\u001b[0m\u001b[1;32m    311\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mlhs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    312\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mPatsyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"model is missing required outcome variables\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/highlevel.py\u001b[0m in \u001b[0;36m_do_highlevel_design\u001b[0;34m(formula_like, data, eval_env, NA_action, return_type)\u001b[0m\n\u001b[1;32m    163\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0miter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    164\u001b[0m     design_infos = _try_incr_builders(formula_like, data_iter_maker, eval_env,\n\u001b[0;32m--> 165\u001b[0;31m                                       NA_action)\n\u001b[0m\u001b[1;32m    166\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mdesign_infos\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    167\u001b[0m         return build_design_matrices(design_infos, data,\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/highlevel.py\u001b[0m in \u001b[0;36m_try_incr_builders\u001b[0;34m(formula_like, data_iter_maker, eval_env, NA_action)\u001b[0m\n\u001b[1;32m     68\u001b[0m                                       \u001b[0mdata_iter_maker\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     69\u001b[0m                                       \u001b[0meval_env\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 70\u001b[0;31m                                       NA_action)\n\u001b[0m\u001b[1;32m     71\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     72\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/build.py\u001b[0m in \u001b[0;36mdesign_matrix_builders\u001b[0;34m(termlists, data_iter_maker, eval_env, NA_action)\u001b[0m\n\u001b[1;32m    694\u001b[0m                                                    \u001b[0mfactor_states\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    695\u001b[0m                                                    \u001b[0mdata_iter_maker\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 696\u001b[0;31m                                                    NA_action)\n\u001b[0m\u001b[1;32m    697\u001b[0m     \u001b[0;31m# Now we need the factor infos, which encapsulate the knowledge of\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    698\u001b[0m     \u001b[0;31m# how to turn any given factor into a chunk of data:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/build.py\u001b[0m in \u001b[0;36m_examine_factor_types\u001b[0;34m(factors, factor_states, data_iter_maker, NA_action)\u001b[0m\n\u001b[1;32m    441\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdata_iter_maker\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    442\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mfactor\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexamine_needed\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 443\u001b[0;31m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfactor\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfactor_states\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mfactor\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    444\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mfactor\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcat_sniffers\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mguess_categorical\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    445\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mfactor\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcat_sniffers\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/eval.py\u001b[0m in \u001b[0;36meval\u001b[0;34m(self, memorize_state, data)\u001b[0m\n\u001b[1;32m    564\u001b[0m         return self._eval(memorize_state[\"eval_code\"],\n\u001b[1;32m    565\u001b[0m                           \u001b[0mmemorize_state\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 566\u001b[0;31m                           data)\n\u001b[0m\u001b[1;32m    567\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    568\u001b[0m     \u001b[0m__getstate__\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mno_pickling\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/eval.py\u001b[0m in \u001b[0;36m_eval\u001b[0;34m(self, code, memorize_state, data)\u001b[0m\n\u001b[1;32m    549\u001b[0m                                  \u001b[0mmemorize_state\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"eval_env\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    550\u001b[0m                                  \u001b[0mcode\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 551\u001b[0;31m                                  inner_namespace=inner_namespace)\n\u001b[0m\u001b[1;32m    552\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    553\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmemorize_chunk\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwhich_pass\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/compat.py\u001b[0m in \u001b[0;36mcall_and_wrap_exc\u001b[0;34m(msg, origin, f, *args, **kwargs)\u001b[0m\n\u001b[1;32m    122\u001b[0m                                  origin)\n\u001b[1;32m    123\u001b[0m             \u001b[0;31m# Use 'exec' to hide this syntax from the Python 2 parser:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 124\u001b[0;31m             \u001b[0mexec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"raise new_exc from e\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    125\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    126\u001b[0m             \u001b[0;31m# In python 2, we just let the original exception escape -- better\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/patsy/compat.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mPatsyError\u001b[0m: Error evaluating factor: NameError: name 'Sales' is not defined\n    Sales ~ TV\n    ^^^^^"
     ]
    }
   ],
   "source": [
    "# this is the standard import if you're using \"formula notation\" (similar to R)\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# create a fitted model in one line\n",
    "lm = smf.ols(formula='Sales ~ TV', data=data).fit()\n",
    "\n",
    "# print the coefficients\n",
    "lm.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting Model Coefficients\n",
    "\n",
    "How do we interpret the TV coefficient ($\\beta_1$)?\n",
    "- A \"unit\" increase in TV ad spending is **associated with** a 0.047537 \"unit\" increase in Sales.\n",
    "- Or more clearly: An additional $1,000 spent on TV ads is **associated with** an increase in sales of 47.537 widgets.\n",
    "\n",
    "Note that if an increase in TV ad spending was associated with a **decrease** in sales, $\\beta_1$ would be **negative**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the Model for Prediction\n",
    "\n",
    "Let's say that there was a new market where the TV advertising spend was **$50,000**. What would we predict for the Sales in that market?\n",
    "\n",
    "$$y = \\beta_0 + \\beta_1x$$\n",
    "$$y = 7.032594 + 0.047537 \\times 50$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.409444"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# manually calculate the prediction\n",
    "7.032594 + 0.047537*50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we would predict Sales of **9,409 widgets** in that market.\n",
    "\n",
    "Of course, we can also use Statsmodels to make the prediction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td> 50</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   TV\n",
       "0  50"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# you have to create a DataFrame since the Statsmodels formula interface expects it\n",
    "X_new = pd.DataFrame({'TV': [50]})\n",
    "X_new.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.40942557])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# use the model to make predictions on a new value\n",
    "lm.predict(X_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the Least Squares Line\n",
    "\n",
    "Let's make predictions for the **smallest and largest observed values of x**, and then use the predicted values to plot the least squares line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>   0.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td> 296.4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV\n",
       "0    0.7\n",
       "1  296.4"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a DataFrame with the minimum and maximum values of TV\n",
    "X_new = pd.DataFrame({'TV': [data.TV.min(), data.TV.max()]})\n",
    "X_new.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  7.0658692 ,  21.12245377])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make predictions for those x values and store them\n",
    "preds = lm.predict(X_new)\n",
    "preds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x14625128>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt4VOW1/z9vwMTIxRiigDfAW1VEGWkVy7HEKkRrDwo8\n2urRH9gqeurdIFGsFmtQQaOeaqsHDwq2Xh8pHmxrLh6JLXrUygFFxFsL1EulXJSWGg2Q9ftj78lM\nkpnJntnvntk7rM/z7Cez98zs/Z03M+/a71rrXa8RERRFURSlqNACFEVRlHCgBkFRFEUB1CAoiqIo\nLmoQFEVRFEANgqIoiuKiBkFRFEUBAjQIxpjdjTGvGmNWGmPeNsbc5h4vN8Y0GWPeM8Y0GmPKgtKg\nKIqieMcEOQ/BGLOHiHxhjOkNLAOmAxOATSIy1xhTA+wlItcFJkJRFEXxRKAuIxH5wn1YDPQCPsMx\nCAvd4wuBM4PUoCiKongjUINgjCkyxqwENgBLRWQ1MFBENrgv2QAMDFKDoiiK4o3eQZ5cRNqAkcaY\nPYEGY8xJnZ4XY4zWzlAURQkBgRqEOCKy1RjzW2AUsMEYM0hEPjXGDAb+1vn1aiQURVFyQ0RMru8N\nMsuoIp5BZIwpBcYBK4AlwBT3ZVOAZ1K9X0RCv/3kJz8puAbVqTpVp2qMb34JcoQwGFhojCnCMTy/\nFJH/McasAJ4yxvwQWAecHaCGQFm3bl2hJXhCddpFddolCjqjoNEGgRkEEVkFHJvi+BbglKCuqyiK\nouSGzlT2wdSpUwstwROq0y6q0y5R0BkFjTYIdGJarhhjJIy6FEVRwowxBgljUHlXoLm5udASPKE6\n7aI67RIFnVHQaAM1CIqiKAqgLiNFUZQeg7qMFEVRFCuoQfBBVPyKqtMuqtMuUdAZBY02UIOgKIqi\nABpDUBRF6TFoDEFRFEWxghoEH0TFr6g67aI67RIFnVHQaAM1CIqiKAqgMQRFUZQeg8YQFEVRFCuo\nQfBBVPyKqtMuqtMuUdAZBY02UIOgKEqPoqGhgfHjJzN+/GQaGhoKLSdSaAxBUZQeQ0NDAxMnTqGl\nZQ4ApaU1LF68kKqqqgIryw9+YwhqEBRF6TGMHz+ZpqYJJJZtX8i4cUtobFxUSFl5Q4PKBSQqfkXV\naRfVaZco6IyCRhsEtqayoihKvqmunsayZVNoaXH2S0trqK5eWFhREUJdRoqi9CgaGhqoq5sHOAZi\nV4kfgMYQFEVRsqanGg2NIRSQqPgVVaddVKdd8q0znonU1DSBpqYJTJw4pdv01Ki0pV80hqAoyi5F\nXd08Ny3VyURqaXGO9ZRRgh/UZaQoyi5FT05N9esy0hGCoii7FJqJlB6NIfggKn5F1WkX1WmXfOmM\nl7Soq5vHDTdczrhxSxg3bomnmcxRaUu/qEFQFCVQbNcWyuV8nQPJs2ffS3X1NKqrp1FXN0/rHsUR\nkUA24ABgKbAaeAu4wj0+C/gIWOFup6Z4ryiKEn3q6+ultHSgwAKBBVJaOlDq6+vzfr5x4ya57xF3\nWyCx2Fir2sKA23fm3G8HGUPYDlwtIiuNMX2B5caYJkCAu0TkrgCvrShKCLCd0WPzfOvXf6TZRp0I\nzGUkIp+KyEr38TZgDbCf+3TOUfAwERW/ouq0i+q0Sz50VldPo7S0BlgILKS0tIYhQwZ5fn9U2tIv\neckyMsYMBWLAK8AY4HJjzP8DXgeqReTzfOhQFCW/2M7oyfV8VVVVLF68MGl2svMep1S2HW09gcDn\nIbjuomagVkSeMcbsA2x0n74FGCwiP+z0Hglal6Io6bFZ2sF2mYgways0oa5lZIzZDfgN8JyI3JPi\n+aHAsyIyotNxmTJlCkOHDgWgrKyMkSNHUllZCSSGb7qv+7pvf3/u3LnceOPttLbeDUBx8dXccst1\nzJgxIxT6dD+x39zczIIFCwAYOnQoN998sy+DEGSWkQEeAe7udHxw0uOrgcdSvNdCvD14li5dWmgJ\nnlCddunpOlNl5IwbN8muuCSCbs/6+noZN26SjBs3Kecsoqj8zwlxltEY4DzgTWPMCvfYTOAcY8xI\nnGyjtcDFAWpQFGUXpvOSmsuWTdmlltTMFq1lpChKB8K0LrFfH39PrluUCq1lpCiKVVJl5BTKGOjd\nfZ7x428KakNjCFZRnXZRnXZJp9NGLMPWTOmotCU+Ywhay0hRehDJdX5ee+01T6/ryTV84qOdbArZ\n7SptkxI/1iSojYiMEBQlTHi9G7ZdXygoCqEzKm2TDnyOEAre+acUpQZBUbLGq4sl32mlfrCRMpoN\nUWqbVPg1COoy8kF8gkjYUZ12iYpOp3xY+MnUnlVVVTQ2LqKxcVFBg8nR+Z/7Q7OMFKWH0LnOT3Hx\nPKqrH+/2dUHU8IlqSYhdfTU1nYegKD0Irx1xkB12mOYx5EJUjRmEvJZRrqhBUJTosqtNBgsTfg2C\nxhB8EBW/ouq0i+q0SxR0RkGjDTSGoCiKVXZ1P3yUUZeRovRwvPjEw7xmgeIdjSEoyi5Ops7XS4A3\n6kFgJYFfg1DwSWipNiIyMS0q9U1Up13CpDPTzNqlS5d6mmiVeE29wCSB0RKLjcnbZwhTe6YjChpF\ndGKaouzS1NXNc+/spwDOXX58tJAdq9xzTAAu4Y033s5Yx8dGvZ/4OaZPv2nXqxkUVvxYk6A2IjJC\nUJRC090IwEttnvr6eikqGuC5ZIONej9RrxkUVtBaRoqy6+K1w++uHlAsNtazQbBR7yfqNYPCil+D\noC4jH0QlN1l12iVMOjOVd47rzFQPKO62ASguvgpYCCx0U0Wn+dbnzbXU7Ps6QROm/3mQ6DwERYk4\nVVVVOWUEdc4uKi6+lljsQSoqBmZcJc3rPINMK551PMcaSksX6FyFEKBpp4qyi5JLiYl4iuumTZuB\nHa7xSD3PoLvz61wF++iayoqi5IXU8xVuzLkjz3VkEwRqnFz8BCCC2ohIUDkqucmq0y750GljYZju\ndGab6ZNtINjr+Qv9f/eis9AavYLPoLKOEBQlZGTyvdskHpBO3BnbvUbQ57dFx7kc0NLiHAuj1qDR\nGIKihIywlo/uqSUuUrX3WWOf4qmqf4FLL4X+/QspLyu0/LWiKHkhU4profEzc7q6ehqlpTXAQr7B\nT/hVr0t44uUmmDkTHnkkGMFhxY+/KagNjSFYRXXaJWidHX3a1VJUNEBisbFZxxIK1Z7Zxj/86vQ9\n67mlRd6cPl1W7dE3HhyRNmNEJkwQeeklKxrzBToxTVF6FvE78VjsQYqKHqKtrY4VKy5g4sQpoa/5\nE3crNTVNoKlpQl4051zP6cMPnVHAAQcw4s47OeqLbWymD3P4DkcUD6DhRz+Cb34zUO2hw481CWoj\nIiMERQmSKJZ3KITmrK7Z1ibywgsikyaJFBW1jwje67enTOWHsjtfRKatU4GOEBRF2ZVJjgGkLbux\nbRvcfz8cdRR8+9vw619DURGccw68/DKXHv9tFnAiX1JaiI8QHvxYk0wbcACwFFgNvAVc4R4vB5qA\n94BGoCzFewOyn3aJil9RddolXzr9+sb96Mx1HkQuxfa8zJfoTkva17z7rsgVV4j0798+GpDBg0Vu\nvlnkk088647Kd5OwVjsFBgEj3cd9gXeBI4C5wAz3eA1we4r3BtNalonKl0R12iWfOv1MUMtVp19D\nlElzqnPPmTPHrpYdO0SefVakqkqS/EgiJ54o8uSTIq2tWeuOyncztAahy4XgGeAU4B1goCSMxjsp\nXhtAUymK4oUg4wCZzp2qQ85Ky+bNInfeKXLQQQkjUFoqcuGFIitXWtEfdvwahLzMVDbGDAViwKuu\nMdjgPrUBGJgPDYqi2CGIuj/pZmd74o034L774NFHaS/BOmyYM6nsggugvNy3vl2FwA2CMaYvsAi4\nUkT+YUxiEp2IiDEm5ZTkqVOnMnToUADKysoYOXIklZWVQKI2eaH348fCoifd/j333BPK9ttV2nPu\n3Lk89dRvKC/fm+rqaZSUlORFb/xYtu8/4IAyjLkMcX+ZxcVXc8op1wFOxz1hwjm0tk4DjmDZsinM\nmnUNxx13nKfzV1dP48UXz6G1dQ1wBKWlNQwbdgY33DA7KXW0mZaWqdTVzUv5+lNOuYbm55+ncvNm\nuO8+mpctc84PUFVF89ixcNxxVJ58srX2XLlyJVdddZW189nab25uZsGCBQDt/aUv/AwvutuA3YAG\n4KqkY+8Ag9zHg4mwyygqfkXVaZdsdBZyqchc2rO2tlaKivYSqBYYLUVFA6S2trb9eRvupFRBZa+u\npBcee8wJCA8enHAL9e8vcuWVTgA5IKLy3SSsMQTAAI8Ad3c6PheocR9fR4SDyorSHVGaS+BlbeWg\nPk9Gw9nW5swYPucckd12SxiC4cNF7r9f5B//8H39noJfgxDkPIQxwHnAScaYFe52KnA7MM4Y8x7w\nbXdfUQqGnzo4hcS27rq6ebS1HZrxNZ5y/nMgZZ2kb30LHnoIRo2CMWPg8cdh506YNAleeAFWrYJL\nLoG+fX1fX3HxY02C2ojICCEqw0jVmZ5cXDp+XUa1tbW+1zrwWsM/m7RV5+6/WiBx3qKivbqdR+CX\nLu25dq3IjBki5eWJ0UBFhcjMmSLr1/u+Xi5E5TdEWF1GvkSpQbCK6kxPLi6QbHUmd6C1tbVWYgpe\ndM+ZMyerayWMTOr4QbaftbtrxV83Z84cxy3U2OgUlDMmYQi+/nWRhQtFWlqy0mGbqPyG1CAoig/y\n7eO3dT0v58nlWn7u/r2OtpJf14/75erd+sk/9t8/YQSKi0XOP1/k1Vezur7i3yDoimnKLk119TSW\nLZvSnr7u+MQ95r8XkKB0+1nn2OvKY3V18xjSciWX8kemsJB+27fBR/+A/faDf/93uOgi2Gcfvx9F\nyQU/1iSojYiMEKIyjFSdmcnmrri+vl5GjToxZ/+5zTTU7nRn6zLyS7cjkh07RBYvluXle0vSi+Qu\nDpCbjz5eZPv2wLT5JSq/IdRlVDii8iVRnXZIdOY1vjpY20HZdGQKKgehIa2x27hR5PbbRQ48sN0I\nbAP5BSfJcGqluLhr4DpshP27GcevQdA1lRXFI2Fd6zgZL2UlglwbOfn6N084iROWL3fSRb/6ynnB\nIYfApZfyPwccwJz/fCyjTiV7/K6prDEERekhpKsHlMqH78XXnwtVJ51ElVtSgssvdw4aA6efDpdd\nBuPHQ1ERJwMnT57s+3qKXXSBHB8k14wJM6rTwe9ErsSkrOvoblJWISa7dV1Kcqq3pSSTyFn3xx/D\nTTfBgQfCv/0b/O//QlkZXHMNvPce/OY3cOqpzqI0nYjC9zMKGq3gx98U1IbGEKyiOu0Fc70ElQtV\nv6hrULcmZZppOn1Z625rE/n970XOOkukV69EoHjECJF580S2bfOkOwrfzyhoFNGgsqJ4Ip/zDQpV\nvyibDt3X2gPbtjkd/tFHJ4xAr14iZ5/tGIi2tiA/ppIBvwZBYwiK0kOI1wNKBJXTB4pzmm/wpz/B\nL37h1Bf6/HPn2D77wMUXO9t+++WsPYg1FpQc8GNNgtqIyAghKsNI1WnXjeNlDeBClbxOxmt7xkcL\nsdgYKS7eu6Pu3/1O5LnnRL7znY4lJU44QeTRR0W+/NK3zmzmS+QrZbczUfkNoS6jwhGVL4nqdPDb\nmcTfP2rUiZ4msAXVcXk9t5f27Gy8iovLJBYbK2dWflfWXHyxyCGHJIxASYnI1Kkir79u8dOIjBp1\noidXVdTWligEahAUJQey7bDDctdvW0fnuMFwamXJ/sNE+vRJGIIDD3Qmlm3caPGTpNeQziBEaW2J\nQqEGQVGyJJdONcjOKPsS1fZ0xGJjpBfHyyS+Ji+wpySdWOTkk0WeecYpOREgXv8fahC6x69B0HkI\nPohKbrLq7EjXfP05WebrN3f7Cq/5/PHJZE1NE2hqmsDEiVOszVvorj2XPvkk331zBWtZzSLe5SS2\n8g925/6i3Vk2bx48/zyccQb06mVFTzpKSkq6Lo6TIqgc1OI8XojKb8g3fqxJUBsRGSFExa+oOjuS\na1noxF1sTbeBT68jkGy1ZHPutO356qsi558vX5mi9tHAGgbJpdwr/dia9zvvbBcc0qByelCXkdLT\nsd0J5OqH96ojm07e5poFnY8n7zcuWeIsNPONb7QbgZ0gzxCTk/kXgYfVFdMDUIOg9GiCCuYGeaeZ\nTSdvcwZ1x2yhvaW4uEz2p05m8135G0kpo+XlIjNmSPOCBUkrpFUUPGCu+EcNQgGJyjAyyjrDGEi0\nPQ/BhnHq2E5tUkmNzKKf7CDhFnqv354iDz0k8sUXXa4di42RWGxs3l0xItH4fkZBo4h/g6AzlRXF\nMtnMGI6/3sbM3D58yXk8wGXcx1GsphloozdP8j3u4yD6Hv82jRdcEMi1lZ5BVushGGN6AX1E5O/B\nSdL1EJQEQdbu7zG89x7ra2ooe+a/2RPnd/MJRTzYq5gHdt7Op5Ttsu22q5XECHw9BGPM48DFwE7g\nj8Cexpj/EJG5uV5UUbwSv9u+/vrbWL/+I4YMOaTQksLBzp3w3HPOugMNDQxxD68qG8CSAw7mG7fd\nxOjevXmpbh4j6H6U0hPxuj6EkkR3PiXgDffvvwF1wG7AKj9+Kg/XtOlWC4yo+BWjrjPorCBbOvPC\n5s0id94pctBBiSBxaanIhReKrFjR4aWZ2jPo1M1srhFUe6aLP+Xy+aPyGyIPMYTexpjdgDOBn4vI\ndmOM+nOUvJHLCl897u7wjTec0cCjjzoNADBsGFx6KVxwAZSXezpNPtqloaGBCRO+T2vr4QC8+OL3\nWbLkiVC0/aZNG3rW98I23VkM4ArgY+A5nBXWhgJ/8GOFPFzTvulUIksumUa5ZicVauJTSlpbZeX1\n18ubZQMk6YPIawP2keU335xTSQm/WVte2icWG9MhjRUqJBYbk7VWv6QaWcZiY0OXtWYTgh4hiMjP\ngJ/F940x64GTrFsmRUlDdfU0li2b0n5j7JQsWGj9OqEZVXz6Kcybx5f/8R8cs2ULAFsp5mHgF/yU\n9zcPovT2GhYff3xetXltn/XrPwXuJD6ic47dkjedcVJle2W7pOguR3cWAxgEzAfq3f0jgR/6sUIe\nrmndcgZBVPyKYdKZ6Q4zk858VCf1evccSHu2tYm89JLIOeeI7LZb+2jgLfaVS/iF9OVfs76zTaXT\nz0Q4r+2T6i48Fhublc6gyPXzh+k3lAnyEENYADwM3ODuvw885RqJjBhjHgJOB/4mIiPcY7OAC4GN\n7suuF5F6DzqUiOPnDtxrvnxymuENN1zOiy8uAUKcZdPSAk884cQH/u//nGNFRTBpEteu38idy38A\nTAWe7/LWTZs2M378ZMB7SmW2cyRy4bbbrmfChPNpbXX2i4uv5bbbfmn1GrmSj88fabqzGMDr7t8V\nScdWerE2wIlAjKSsJOAnwDXdvC8A26kUmqBnHfstA9H5/UVFe0ltba01fR1Yu1ZkxgynjES8QSoq\nRGbOFFm/PoWeaoH+XUpT5LPchN81mwtNGDXZhqBLV+DU+h0QNwjAaOBFzxdwgtCdDUJ1N+8JprWU\nghK0QbBx/traWikqGiAwWqDabkfb1ibS2CgyYUKH5Sjf6V8mc4aPcgrQdSK5E5syZYqUlx8s5eUH\ny8EHH1mQ4GhUO9WwLHAUNPkwCKOAl4Gt7t/3gWM8XyC1QVgHvIHjdipL8Z7gWswiUfEr2tTpp0Po\n7kfpV6cNg+DlHFnr3LpV5N57Rb72tcRooLhYPv72t+XEknLPd9zJ6x0bU+aOGizq7IYozuuIay4v\nP9jXdyMqv3W/BsFLltFyY8xY4GvuoXdFZLtHj1Qq7gd+6j6+BWey2w87v2jq1KkMHToUgLKyMkaO\nHEllZSWQWKyi0PtxwqIn3f7KlSutnO+rr75yYwBTASe/fPjwYygqauPss7/LjBkzMr4/7r+94YbZ\nAMye7fhvbbVnIhtpDQClpQuorl6Y9flgDc7A2Nlfu/ZPfP3r36K8fG+qq6exZs0ab+cbOBB+/nOa\n58+HL790zrbffjRXVcHpp3PrA4/yh6/uAneecXyhnpKSki7nu+yyGlpb78DJ3GlG5GKMmYcTmltD\ncfE8qqsf73B9v+2ZvP/aa68xa9ZdbvxnDS++eA5Lljze4f9X6O9n5u/rmuQW6bDv5XwrV64Mze85\neb+5uZkFCxYAtPeXvkhnKYDJwKROf+OPJ3m1OHQaIXh5joiMEHY1Ot4910sYSyb7vYtNV0ba8+fc\nscNZdvKUUyTpdlQ2H320yNNPi2zf3v7SbEY0qe5w+/U7MG/umzBWne2OKHxfbUOAI4R/BbdSVmp+\nnYsBMsYMFpG/ursTgVW5nEcpNPNIzjX3Mns4H/it3tk5C2XTpsNYseIiuv2cmzbB/Plw//2wfj0A\n/wR+SSU/52T+9P59LO7bl6rezk+uoaGBTZs2U1RUTVvbKmBExvkVQ4YMYsuW6UlHpnPIIV+jsXFR\nzp9116IKmEJ5+S2MGnWMZhelw4816W4DHgc+AVqBD4EfAI8Ab+LEEJ4BBqZ4XzDm0zJR8Sva0tnx\n7nm09TtGPzqD8m+nujMeNerExAuWLxeZOlWkpCQxIjjkEPnFYUfLnvw8ZfukymaKxcZ0u2aCM1IZ\nLTBaiovL8lojKMigbFC/I5uao/JbJx8L5ADfBWYAN8U3Pxf1cD37LRUAUfmSBBFUjsXGdAhy2ugg\nctUZZGeV6txzb71V5NFHRU44IWEEjBE5/XSR554T2bkzo4sl1XPl5Qd3qzlbo6dBZXuao/JbD9wg\nAP/p3tV/hJMh9BYw389FPVwzgKZSbBPEWse5nC8f8xvGjZsk3z/xNPng3HNFBg5MGIKyMpFrrhF5\n//0u70lnpFLphdE91q+t5I98GIRV7t833b99gWV+LurhmtYbSgk3+SipkBNtbSK//73IWWeJ9OqV\nMAQjRojMmyeybVvGz5TKwHX+rDDQDXqGP1CrhJt8GITX3L+vAPsBuwMf+Lmoh2vab6kAiMowMgo6\nnU69JqdOPRCX0bZtTod/9NEJI9Crl8hZZ8nSe+5xDIUP6uvr3cyh0a4xsD+yicL/XSQaOqOgUcS/\nQSjyEHd+1hizF3AHsBxY6waLlR5OQ0MD48dPZvz4yTQ0NBRaTlrimUHjxi1h3Lgl/iqU/ulPMH06\n7L8/TJsGb74J++wDN97oZA899RQccwyYnFcpbNf82GM/p7R0LfApsNDNMprm67yK4ot0lgI4Dhic\ntD8FaALuBQb4sULdbURkhNCTyfdU/4KWFti5U16vrZX/rRgkO5PmDsjo0U7w+MsvA7t0VEpBZHJ/\nRUH/rgJBuYyAFUC5+/hbwF9xJqbVAk/7uWi3otQgFBw/fvlcO4m8dy6ffSZy992ybd99241AC73l\nkV67y8s/+1nw148I6Yz1rlIfKEoEaRDeSHr8c2BWqueC2KJiEKLiV8xFp58Vx3LtJPI2D2HVKpFL\nLpHtu+/ebgjWMUBmcLsMYGO3n7Un/99Tke67kOk7EoY1lW0SBY0i/g1CppnKvYwxu4lTt+gUINm5\n6WUdBSXC5LpKWS7rH/vF0zoLO3bAf/+3s+6AWwumN/A8R3IfO3iW62jjgsA07kqEZuU5JXvSWQqc\nBXFeBpbguI+K3OOHAi/5sULdbURkhNDTycWFk+ukKz9kHM1s2CBSWyuy//7to4F/9uotD5b0lyOY\nKmGscxOE6yyIKrXpjkex7lFPgSDTToETcOoN9Uk6dhhwrJ+LditKDUJk6ZpjXyHW1xXoRKoO6LLj\nKkXOP1+kuLjdEGzbbz+5erd+0o/7k7TFUz6rpbz84IIER5M769raWut+eRu+/myCykEZhKBiTD0p\nMB6oQSjUFhWDEBW/Yr7XQ8g1x95v6YoS5sn5XCR/NL3bjYAY4yxI09go40+Z2KWjitdkCirW0V17\npapr1N06B16J68z3HXu2BshLewYVwPZ63qj81v0aBI0FKClJXps4vl6vV99wVVUVo0YdQ1PTBJwq\nk8FSdeSRrDrjJAYsupyy7V85NXrLy+HCC+GSS2DYMADkjge6vLe8fCOjRi0JpPqll/bqHHNpawPo\nqjMqxL83hx9+CPAwFRUDrLRtULGpQsS8Qo0faxLURkRGCD0VG75hm3d0Ke+y29pEXnhBZNIkkaKi\nxIhg5EiR+fNFvvgiUE1e8NJeqV7jLOEZLpdR9teplqKiARKLjbVyraBGOT0t3oG6jBTb5JJmmAob\nvtnOnVnF7nvLW5dfLjJ8eMII9O4tcs45Ii+91G1JiXxOsPLSXqk669raWutaamtr29djrq2tTfs6\nP+2Q+Lz14tRnCr9R62lzKdQgFJCo+BWz1ZmuIwv6x5NKZ1zLobwrd3OlfE5pwhAMGiQya5bIJ5/4\num4QPu9szht0WelsdPj5/ya+N9ndOGTTnoUKKkflt64GoYBE5UuSrc5MHUN9fb3EYmOlvPzgbhd1\n8a1zxw758cgT5DmOkqTeRV7uvbvUjjhObp01y0oHke3Ix2ZQOUiyDSr7nZ0ei4113V2HW23PMGQB\nReW3rgZByRqvmULpXCvJC+MUF+9t/0e6ebPInXeKHHRQuxH4J8Uyj2/JMfQRJwunWqC/ldFKT/Mj\ndyZog9D5BsKYMjGmr5X/TS6jljAYkEKhBkFJSaYO3Y9bIBYb26XTiMXG2hG9cqXIhReKlCa5hYYN\nk3cuvFAmVX7XTWWNp2R27bxisbE511AK2o8c/3/EYmPadQZ1jVz/551fV1y8t6s3s+ZUhiTX/4WX\nc3cXt+pJMYFsUYNQQMI6jOz6w97L2qQhp1PuOhM5Z1pbRZ58UuTEE2VpkltIqqpk+c03S9UpE9s7\nlY7a7WbnZHNXmbsLrlqCmhGdqiOcM2dOh+e9fL6E4RorvXv3EThKoNzVnlqz3+9UpvbM9txBjfbC\n+lvvjBqEAhLWL0nXH0VN+4/C7w8mFhvToVODConFxmSt8YXHHpMFBx0hm4oTBeaWlpaKXHGFyLvv\nps2+SU5rTHYZ2ZzQ1R25B+mDc02l+r/267dvzgan6/85/apufu/KM7VntudWg6AGQelEph9FNj+w\nVHeVtbUysm0cAAAV+UlEQVS1Ykw/gf0FjpLi4jLvP/62NpGXXpJPKivlq6TRwGrTS1ZfdpnI3//e\n7WdI1pScnul0YPlNifVKoQyCn3WaU40E4/pTaQ6yPbM5t7qM1CAonejuR+E1qJz5Dt25K8+U097O\nF1+IPPSQyLHHthuBHRhZxESp5AWBhz1N2LLlO853p1EIl5HfdZpTxYr8GJh8okFlNQgFIczDyOQf\nRbIv2SupOuRUd40ZO5u1a0VmzBApL283BFJRIY8NPUwOoK7LeZLbMxejlqkjSH4uVWdny+edjnwF\nlTvWkFqas0HonE1mTJn1NOM4Yf4dxYmCRhE1CAWlEF8SJ997jPTtO1j69TvQU2kAWwvk9O072O1s\nJqW/+2xrE2lqEjnjjI4lJb7+dZGFC0VaWtJ29p112sqU6vz6RPA5fwYhX3T8rDX2S4YEQJjbM04U\nNIqoQdilcO7aysQJpiZcD0HMBeiaqVQmvXsnsnigomP8YOtWkXvvFTn88IQRKC4WOe88kVdfTXn+\nVLEAL5/Df+ZJtRuE7pl+5l3ZZbKr49cgaLXTCFFXN4/W1sPdvUuIV2hsbbVfobGqqorFixe2Vzzd\ntOkYVqy4oP2aAMOHP0zVkCFw2WWwcCFs2+Y8sd9+vH/KKdywdiMfrP4QLrmGioqB7VVT4+fPpoKq\nXUZwzDFHUlGxBCCQSqeKEkn8WJOgNiIyQsj3MNK50x0t8Rr+mTJwku8MbaxVnBw/KGKHTOAKWV6+\njySJEBk7VuTpp6Xht7/1vEhOx7v3pQILJBYbk/EO16/LyO+IoHOsI0x34zZdRvkiCu6YKGgUUZdR\nQcn3lySTyyjTSlt+F56J5/0PoK/M4CxZy4CEEdhjD5GLLxZ5883296VOgUydstjVIHhz52TbEXf3\neic2461GU7ZF4/JJKgMb9jIcUehso6BRJOQGAXgI2ACsSjpWDjQB7wGNQFmK9wXSWD2BdEHlICbk\nxM8ZY7nM5wJpoVe7Idi2774id98t8tlnad/nxSDYDvjmQueMmi7xkTQENQnKD2HUpOQPvwYh6BjC\nw8C9wCNJx64DmkRkrjGmxt2/LmAdPYa47z2ZhoYGli9/A3gbp8kHAMP8Xai1lZP++iGzqOWbfABA\nG4ZXKwZy/CMP06eqCoqKUr61unoay5ZNoaUlfmQ6MIXS0hqqqxd2+TwdYxVHsWKFP+nZ4sRm7iA5\nPtLa+kCXuEznVeTCSOe2T9XmipIWP9bEywYMpeMI4R1goPt4EPBOivfYN50BEIZhZLpJT9C/fdJY\nVjo//ljkpptEBg5sHw1sYQ+5kyoZXjJA6uvrPU9si+f8x2MCXtJjC+GGSTfLN/nOOpVvPpObrpDE\n237UqBNDoac7wvA76o4oaBQJuctIUhuEz5Iem+T9pOP2WyoAwvAl8VIWoVudbW0iv/+9yNlnO6uP\nxU8yYoS8deWV8q/fntDeoQfZYSf75vMZqPXiMkrnmw9bUDmZMHw/vRAFnVHQKBJ+l1FGRESMMZLq\nualTpzJ06FAAysrKGDlyJJWVlQA0NzcD6L67D78DXgEmuPvNwBqSaW5u7vL+HVu3sur6m6n4YA0H\nbP+SSqCtqIhnRozgia8Mnw88mOrTTuOaM0varzd+/GRaWqYCQ4BKWlrghhtmU1JS4klvQ0MDN9ww\nG4CjjjqIZ59dxvbtLZx11qnMnz8fgJKSEmbOvDxv7VdSUsItt0zniSceZv36jygvH8RFF53f7i5q\nbm5my5aNdMRp36qqKkpKSjqcr9Dfh6jtx4+FRU+6/WStYdBTWVlJc3MzCxYsAGjvL33hx5p42Ujt\nMhrkPh5MhF1Gtsm2HIOIU2zOyTrKok7OBx/InydPli2Y9tHApxj5KaNl2G79O9wtdz6Pn6Bl13o7\n/QUmd3FxhZEwZhQpSmeIoMtoLlDjPr4OuD3FewJoKvvYHEZ6KSbXfR5/vcBoKS8/uOs8hJ07RZ57\nTuT000VMwhC8zGg5l19JMQ+6bqfUcxySdea6YlpqX/0AiZfB6NdvXyttGRTqmw+GKOiMgkaRkLuM\njDGPA2OBCmPMh8BNwO3AU8aYHwLrgLOD1BA2OmeqxGfrnnvupe5sXSfTpaUF7rrrli7HMs9IrgI+\nZdSoJYnXbN0KTz8NF10EHzjZQpSU0DBgIDd8cgHLmeW+N1UmSgPwAMuXb6ShoSHputuBB5Ie++FQ\nYB4Jd1d4iWd4Jbs3FKVH4ceaBLURkRFCtmQeBXS9M/dSXTStK2PVKpFLLhHp06d9NCAHHihy220i\nGzemcN84s4mLi8vcEUBqF5Rfl1HypDOnRHO1+9nD7TJSlChA2F1GOYnqoQYhc0npereDzM5lJJJw\nZVSdfKas+PGPRSorJekiIiefLLJ4scj27Snf1zk1NFFGuWvH73fiU21trTv5bLRrDPaUvn0HqzFQ\nFAuoQSgguS+lmG4UUC9wlPTuvU97+QRPaY0bNojMni2y//4JI9Cnj8iPfiSyerW18tdxDZ2roGZb\n3z+Imkv5RHXaJQo6o6BRRA1CQcn2S9L57ri7dYK7zWR59VWR8893ykzHe+7DDhP52c9EPv88pU6v\nefOZsmqSF3vJlJGU7TWj8qNTnXaJgs4oaBRRgxAZutbsSSw/maqiaFp3TEuLs9DMN76RMALGiEyY\nINLY6GQTedTQXQfe3epeXtxHmq6pKPlDDUJE6K7zzOS3FxGRv/xFZOZMkb33ThiCvfYSufZakT//\nOeU1O9+Z2+7AvZxPi60pSv7waxBSVydTPNF5BmOuxBeJ2bLlTJxCcAuBhZTuPoPZpxwPkyfD0KFw\n662wcSOMHAnz58NHH8HcuTCsayG7+DmbmibQ1HQoEydOYdOmzd1qqaubl5Tq6ixcE0+T7Ux19TRK\nS2sSektrfBV9s9WeQaM67RIFnVHQaANdMS1PZKpC2bETHkcfbuTiPdYya+8+9Kupcd7Quzd873vO\n6mQnnADGZLxex3M209JyBPAgpaU11iphdq5UmmrlMa2+qSgRws/wIqiNHugyEkmsZeAswjK23RUT\nd6scyrtyN1fK55Qm3EKDBonMmiXyySdZXStdRlN3cYEgfP5hLgCnKD0JfLqMjHOOcGGMkTDq8kvn\n9YNLS2tY/PRDVLz+Optvns34ttb213525JHsddNNMHEiFBf7vhZcAVwEjHCum2Hd4lSzqRVFCT/G\nGEQks/sgE36sSVAbERkh+JmHsBeb5Rq+Jx+XJmYSf1lUJL/bd4i8dN99VvTF78z79dvXTWkNd2A3\nKql9qtMuUdAZBY0iGlSOHEfzF+ZxER+xL3U8yb4t/+SLgQPhjjso2biR0z5exzcvvTTjORoaGhg/\nfjLjx0+moaEh7euqqqpobFzEYYcdDIyw/EkURelx+LEmQW1EZITgmdZWkSeflM1HHSVJt+lSz1Fy\nOldJn933yWqWb7Y+fp0LoCi7Bug8hPDywmOPyYKDjpCNJbu3G4Hte+whD5TsKYdxW8Y5CemCsLnm\n9WtgV1F6Pn4NgrqMfJAyN1mEV+6+m9+V782Yc89lyp/XUPHVl7xtevH2ZZfR+9NPWfStk3mPwR3e\ntmnTZsaPn8yxx1YyYcL33fkDE5g4cUpGt5BXnXH3UWPjotAGiaOS66067RIFnVHQaAOdh2CLlhZ4\n4gm23noro911B3Zi+DVnch+XsVT+wrh3n6WxX78uufnFxdeyevV2WlvvcU82HRgEVHVZA0Hz+hVF\nCQpNO/XLunVw//3wX/8FW7YAsJF+PMhgHuBiPuQa94ULGTduCY2Ni2hoaOD6629h/fpPGTJkf2AH\nK1ZcRHwhHGfm7xJgUYf3xdG0UEVRUuE37VRHCLnywQcwfTo8+yy0tTnHRo3iji/hxtUX8xUH4nTw\nA4DEnXzn+QEtLTUcfvghKS7wCbCQ4uJrqa7+ZYdn4it3KYqi2ERjCLnSpw/Nzz4LvXrBeefBK6/A\nH//I0XWzKSq9EfgUOI+iompisYfbJ4I5JSXOwxkBLHEf9+5QEwiuBXoBDzB8+GG+O/+o+D9Vp11U\npz2ioNEGOkLIlcGDYdYsmDYNBg7s8NThhx/O+vW3MGTIIG677dEOHfqmTRuA3wN3ukem8/e/D2bx\n4oWce+6lbNmyN/BLnPWRF1JRsSQvH0dRFEVjCBZJWZqiU4mIY4+tZMWKC0iOFxhzNc899zhAt+9X\nFEVJh8YQQkTHCqN0yRACqKgY0OV9Il+jrm4ejY2Luq0eqiiKEhQaQ/BBLn7F6uppFBVdTSJeUAOM\naX8+iPkCUfF/qk67qE57REGjDdQgWMTLgjFVVVX89KfVFBVVAw8A51Fa+itfC8soiqLYQGMIlulu\njkD8eWf1sh1UVAzUuQSKolhBYwgRInXQ+UY1BoqihAJ1Gfmgs1+x4zrGTh2i2bNnt5eqvv762zyv\nVxykzrCiOu2iOu0RBY020BGCRTp2+E6W0U03VdPWVgfgxg0URVHCicYQLDF79mx+/OO7gLvoWJPo\nTuAwd99QVPQCbW13AzrPQFEUu/iNIRTMIBhj1gF/B3YC20XkuKTnImUQGhoa+M53zqGt7QckjADA\nlYAAP3P3p3Pwwftw0EGHA1qYTlEUu/g1CIWMIQhQKSKxZGMQJeJ+xbq6ee5d/53Ar3DSSX9MaWkx\njjGY4m530r//3nlflyAq/k/VaRfVaY8oaLRBoYPKOVuy8FIFXEJRUQuHH35kl2dTzVRWFEUJA4V0\nGf0Z2IrjMvpPEXkw6bnIuYyS00mNuYqDDtqX/v33YvXq92htvQPQmIGiKMES5RjCYBH5qzFmb6AJ\nuFxE/uA+FymDAB0nnK1e/Ub76mfFxVcxfPgxVFQM0JiBoiiBEtmJaSLyV/fvRmPMYuA44A/x56dO\nncrQoUMBKCsrY+TIkVRWVgIJf16h9+PHmpubKSkpobFxEePHT6a19WJgCFBJaysUFc1n5szLC6b3\nnnvuCWX7ZWrPMOjR9szPfhTac+XKlVx11VWh0RPfb25uZsGCBQDt/aUvRCTvG7AH0M993Ad4CRif\n9LxEgaVLl3Y5Nm7cJIEFAuJuC2TcuEn5F5dEKp1hRHXaRXXaIwoaRUTcvjPnvrkgLiNjzDBgsbvb\nG3hURG5Lel4KocsGXtZEUBRFCYLIxhAyEWWDAN0XuFMURQmCKM9DiDzJvtpkgljTwA/pdIYN1WkX\n1WmPKGi0gRoERVEUBVCXkaIoSo9BXUaKoiiKFdQg+CAqfkXVaRfVaZco6IyCRhuoQVAURVEAjSFY\nQdNMFUUJAzoPocDoRDRFUcKCBpULSHNzM3V18wqyTnI2RMX/qTrtojrtEQWNNlCDoCiKogDqMvKN\nuowURQkLGkMIARpUVhQlDGgMoYDE/Yphq13Umaj4P1WnXVSnPaKg0QZqEBRFURRAXUaKoig9BnUZ\nKYqiKFZQg+CDqPgVVaddVKddoqAzChptoAZBURRFATSGkDOaaqooStjQeQgFQCejKYoSRjSoXAAS\n9YuGENb6RclExf+pOu2iOu0RBY02UIOgKIqiAOoyygl1GSmKEkY0hlAgNKisKErY0BhCgaiqqmLm\nzMtDW78omaj4P1WnXVSnPaKg0QZqEBRFURRAXUaKoig9BnUZKYqiKFYoiEEwxpxqjHnHGPO+Maam\nEBpsEBW/ouq0i+q0SxR0RkGjDfJuEIwxvYD7gFOBI4FzjDFH5FuHDVauXFloCZ5QnXZRnXaJgs4o\naLRBIUYIxwEfiMg6EdkOPAGcUQAdvvn8888LLcETqtMuqtMuUdAZBY02KIRB2A/4MGn/I/eYoiiK\nUkAKYRB6TPrQunXrCi3BE6rTLqrTLlHQGQWNNsh72qkxZjQwS0ROdfevB9pEZE7Sa3qM0VAURckn\nkSpdYYzpDbwLnAx8ArwGnCMia/IqRFEURelA73xfUER2GGMuAxqAXsB8NQaKoiiFJ5QzlRVFUZT8\nE6qZysaYWcaYj4wxK9zttKTnrncnsr1jjBlfSJ2untBOrjPGrDPGvOm24WvusXJjTJMx5j1jTKMx\npizPmh4yxmwwxqxKOpZWU6H+32l0hu57aYw5wBiz1Biz2hjzljHmCvd4qNo0g85QtakxZndjzKvG\nmJXGmLeNMbe5x0PTnhk02mtLEQnNBvwEuCbF8SOBlcBuwFDgA6CogDp7uRqGuppWAkcUuv2S9K0F\nyjsdmwvMcB/XALfnWdOJQAxY1Z2mQv6/0+gM3fcSGASMdB/3xYnLHRG2Ns2gM4xtuof7tzfwCvAv\nIWzPVBqttWWoRgguqSLkZwCPi8h2EVmH88GOy6uqjkRhcl3ndpwALHQfLwTOzKcYEfkD8JlHTQX7\nf6fRCSH7XorIpyKy0n28DViDM58nVG2aQSeEr02/cB8W49z0fUb42jOVRrDUlmE0CJcbY94wxsxP\nGp7tizOBLU6hJ7OFfXKdAM8bY143xlzkHhsoIhvcxxuAgYWR1oF0msL2/4YQfy+NMUNxRjWvEuI2\nTdL5insoVG1qjCkyxqzEabelIrKakLVnGo1gqS0LUcuoyRizKsU2AbgfGAaMBP4K1GU4VSGj4WGP\nxI8RkRhwGnCpMebE5CfFGU+G6jN40FRIvaH9Xhpj+gKLgCtF5B8dhISoTV2dT+Po3EYI21RE2kRk\nJLA/8C1jzEmdni94e6bQWInFtixE2uk4L68zxvwX8Ky7+zFwQNLT+7vHCkVnPQfQ0RIXFBH5q/t3\nozFmMc4wcYMxZpCIfGqMGQz8raAiHdJpCtX/W0Ta2ypM30tjzG44xuCXIvKMezh0bZqk81dxnWFt\nUwAR2WqM+S0wihC2ZyeNXxeR5vhxv20ZKpeR2+BxJgLxTI8lwPeNMcXGmGHAoTgT2grF68Chxpih\nxphi4HuuxoJjjNnDGNPPfdwHGI/TjkuAKe7LpgDPpD5DXkmnKVT/7zB+L40xBpgPvC0i9yQ9Fao2\nTaczbG1qjKmIu1qMMaXAOGAFIWrPdBqNMYOSXuavLYOOimcZQX8EeBN4A6fhByY9NxMnKPIOUBUC\nrafhZEx8AFxfaD1JuobhZBasBN6KawPKgeeB94BGoCzPuh7HmZneihN/uSCTpkL9v1Po/EEYv5c4\n2SVt7v95hbudGrY2TaPztLC1KTAC+D9X55vAte7x0LRnBo3W2lInpimKoihAyFxGiqIoSuFQg6Ao\niqIAahAURVEUFzUIiqIoCqAGQVEURXFRg6AoiqIAahAUJSXGmAFJ5YT/mlReuK1zGWFjzFXGmF8U\nSqui2EINgqKkQEQ2i0hMnJpQDwB3uY8vBr7f6eXfAx7Lt0ZFsY0aBEXxRry88CLgdOOsDR6v4Lmv\niCwrkC5FsYYaBEXJAhHZglMP5jvuoe8DTxZOkaLYQw2ComTP4yTcRt9z9xUl8qhBUJTsWQKcbIyJ\n4SxpuKLQghTFBmoQFCVLxFngZSnwMBpMVnoQahAUxRudywI/jlOOWN1FSo9By18riqIogI4QFEVR\nFBc1CIqiKAqgBkFRFEVxUYOgKIqiAGoQFEVRFBc1CIqiKAqgBkFRFEVxUYOgKIqiAPD/AUd7480b\nVee2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbf92b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# first, plot the observed data\n",
    "data.plot(kind='scatter', x='TV', y='Sales')\n",
    "\n",
    "# then, plot the least squares line\n",
    "plt.plot(X_new, preds, c='red', linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence in our Model\n",
    "\n",
    "**Question:** Is linear regression a high bias/low variance model, or a low bias/high variance model?\n",
    "\n",
    "**Answer:** High bias/low variance. Under repeated sampling, the line will stay roughly in the same place (low variance), but the average of those models won't do a great job capturing the true relationship (high bias). Note that low variance is a useful characteristic when you don't have a lot of training data!\n",
    "\n",
    "A closely related concept is **confidence intervals**. Statsmodels calculates 95% confidence intervals for our model coefficients, which are interpreted as follows: If the population from which this sample was drawn was **sampled 100 times**, approximately **95 of those confidence intervals** would contain the \"true\" coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td> 6.129719</td>\n",
       "      <td> 7.935468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TV</th>\n",
       "      <td> 0.042231</td>\n",
       "      <td> 0.052843</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  0         1\n",
       "Intercept  6.129719  7.935468\n",
       "TV         0.042231  0.052843"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the confidence intervals for the model coefficients\n",
    "lm.conf_int()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind that we only have a **single sample of data**, and not the **entire population of data**. The \"true\" coefficient is either within this interval or it isn't, but there's no way to actually know. We estimate the coefficient with the data we do have, and we show uncertainty about that estimate by giving a range that the coefficient is **probably** within.\n",
    "\n",
    "Note that using 95% confidence intervals is just a convention. You can create 90% confidence intervals (which will be more narrow), 99% confidence intervals (which will be wider), or whatever intervals you like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis Testing and p-values\n",
    "\n",
    "Closely related to confidence intervals is **hypothesis testing**. Generally speaking, you start with a **null hypothesis** and an **alternative hypothesis** (that is opposite the null). Then, you check whether the data supports **rejecting the null hypothesis** or **failing to reject the null hypothesis**.\n",
    "\n",
    "(Note that \"failing to reject\" the null is not the same as \"accepting\" the null hypothesis. The alternative hypothesis may indeed be true, except that you just don't have enough data to show that.)\n",
    "\n",
    "As it relates to model coefficients, here is the conventional hypothesis test:\n",
    "- **null hypothesis:** There is no relationship between TV ads and Sales (and thus $\\beta_1$ equals zero)\n",
    "- **alternative hypothesis:** There is a relationship between TV ads and Sales (and thus $\\beta_1$ is not equal to zero)\n",
    "\n",
    "How do we test this hypothesis? Intuitively, we reject the null (and thus believe the alternative) if the 95% confidence interval **does not include zero**. Conversely, the **p-value** represents the probability that the coefficient is actually zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    1.406300e-35\n",
       "TV           1.467390e-42\n",
       "dtype: float64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the p-values for the model coefficients\n",
    "lm.pvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the 95% confidence interval **includes zero**, the p-value for that coefficient will be **greater than 0.05**. If the 95% confidence interval **does not include zero**, the p-value will be **less than 0.05**. Thus, a p-value less than 0.05 is one way to decide whether there is likely a relationship between the feature and the response. (Again, using 0.05 as the cutoff is just a convention.)\n",
    "\n",
    "In this case, the p-value for TV is far less than 0.05, and so we **believe** that there is a relationship between TV ads and Sales.\n",
    "\n",
    "Note that we generally ignore the p-value for the intercept."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How Well Does the Model Fit the data?\n",
    "\n",
    "The most common way to evaluate the overall fit of a linear model is by the **R-squared** value. R-squared is the **proportion of variance explained**, meaning the proportion of variance in the observed data that is explained by the model, or the reduction in error over the **null model**. (The null model just predicts the mean of the observed response, and thus it has an intercept and no slope.)\n",
    "\n",
    "R-squared is between 0 and 1, and higher is better because it means that more variance is explained by the model. Here's an example of what R-squared \"looks like\":"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"08_r_squared.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the **blue line** explains some of the variance in the data (R-squared=0.54), the **green line** explains more of the variance (R-squared=0.64), and the **red line** fits the training data even further (R-squared=0.66). (Does the red line look like it's overfitting?)\n",
    "\n",
    "Let's calculate the R-squared value for our simple linear model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.61187505085007099"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the R-squared value for the model\n",
    "lm.rsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is that a \"good\" R-squared value? It's hard to say. The threshold for a good R-squared value depends widely on the domain. Therefore, it's most useful as a tool for **comparing different models**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Linear Regression\n",
    "\n",
    "Simple linear regression can easily be extended to include multiple features. This is called **multiple linear regression**:\n",
    "\n",
    "$y = \\beta_0 + \\beta_1x_1 + ... + \\beta_nx_n$\n",
    "\n",
    "Each $x$ represents a different feature, and each feature has its own coefficient. In this case:\n",
    "\n",
    "$y = \\beta_0 + \\beta_1 \\times TV + \\beta_2 \\times Radio + \\beta_3 \\times Newspaper$\n",
    "\n",
    "Let's use Statsmodels to estimate these coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    2.938889\n",
       "TV           0.045765\n",
       "Radio        0.188530\n",
       "Newspaper   -0.001037\n",
       "dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a fitted model with all three features\n",
    "lm = smf.ols(formula='Sales ~ TV + Radio + Newspaper', data=data).fit()\n",
    "\n",
    "# print the coefficients\n",
    "lm.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we interpret these coefficients? For a given amount of Radio and Newspaper ad spending, an **increase of $1000 in TV ad spending** is associated with an **increase in Sales of 45.765 widgets**.\n",
    "\n",
    "A lot of the information we have been reviewing piece-by-piece is available in the model summary output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>Sales</td>      <th>  R-squared:         </th> <td>   0.897</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.896</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   570.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 21 Jan 2015</td> <th>  Prob (F-statistic):</th> <td>1.58e-96</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>18:22:21</td>     <th>  Log-Likelihood:    </th> <td> -386.18</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   200</td>      <th>  AIC:               </th> <td>   780.4</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   196</td>      <th>  BIC:               </th> <td>   793.6</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    2.9389</td> <td>    0.312</td> <td>    9.422</td> <td> 0.000</td> <td>    2.324     3.554</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>TV</th>        <td>    0.0458</td> <td>    0.001</td> <td>   32.809</td> <td> 0.000</td> <td>    0.043     0.049</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Radio</th>     <td>    0.1885</td> <td>    0.009</td> <td>   21.893</td> <td> 0.000</td> <td>    0.172     0.206</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Newspaper</th> <td>   -0.0010</td> <td>    0.006</td> <td>   -0.177</td> <td> 0.860</td> <td>   -0.013     0.011</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>60.414</td> <th>  Durbin-Watson:     </th> <td>   2.084</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td> 151.241</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-1.327</td> <th>  Prob(JB):          </th> <td>1.44e-33</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 6.332</td> <th>  Cond. No.          </th> <td>    454.</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                  Sales   R-squared:                       0.897\n",
       "Model:                            OLS   Adj. R-squared:                  0.896\n",
       "Method:                 Least Squares   F-statistic:                     570.3\n",
       "Date:                Wed, 21 Jan 2015   Prob (F-statistic):           1.58e-96\n",
       "Time:                        18:22:21   Log-Likelihood:                -386.18\n",
       "No. Observations:                 200   AIC:                             780.4\n",
       "Df Residuals:                     196   BIC:                             793.6\n",
       "Df Model:                           3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      2.9389      0.312      9.422      0.000         2.324     3.554\n",
       "TV             0.0458      0.001     32.809      0.000         0.043     0.049\n",
       "Radio          0.1885      0.009     21.893      0.000         0.172     0.206\n",
       "Newspaper     -0.0010      0.006     -0.177      0.860        -0.013     0.011\n",
       "==============================================================================\n",
       "Omnibus:                       60.414   Durbin-Watson:                   2.084\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              151.241\n",
       "Skew:                          -1.327   Prob(JB):                     1.44e-33\n",
       "Kurtosis:                       6.332   Cond. No.                         454.\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print a summary of the fitted model\n",
    "lm.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are a few key things we learn from this output?\n",
    "\n",
    "- TV and Radio have significant **p-values**, whereas Newspaper does not. Thus we reject the null hypothesis for TV and Radio (that there is no association between those features and Sales), and fail to reject the null hypothesis for Newspaper.\n",
    "- TV and Radio ad spending are both **positively associated** with Sales, whereas Newspaper ad spending is **slightly negatively associated** with Sales. (However, this is irrelevant since we have failed to reject the null hypothesis for Newspaper.)\n",
    "- This model has a higher **R-squared** (0.897) than the previous model, which means that this model provides a better fit to the data than a model that only includes TV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Selection\n",
    "\n",
    "How do I decide **which features to include** in a linear model? Here's one idea:\n",
    "- Try different models, and only keep predictors in the model if they have small p-values.\n",
    "- Check whether the R-squared value goes up when you add new predictors.\n",
    "\n",
    "What are the **drawbacks** to this approach?\n",
    "- Linear models rely upon a lot of **assumptions** (such as the features being independent), and if those assumptions are violated (which they usually are), R-squared and p-values are less reliable.\n",
    "- Using a p-value cutoff of 0.05 means that if you add 100 predictors to a model that are **pure noise**, 5 of them (on average) will still be counted as significant.\n",
    "- R-squared is susceptible to **overfitting**, and thus there is no guarantee that a model with a high R-squared value will generalize. Below is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.89719426108289557"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# only include TV and Radio in the model\n",
    "lm = smf.ols(formula='Sales ~ TV + Radio', data=data).fit()\n",
    "lm.rsquared"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.89721063817895219"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add Newspaper to the model (which we believe has no association with Sales)\n",
    "lm = smf.ols(formula='Sales ~ TV + Radio + Newspaper', data=data).fit()\n",
    "lm.rsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**R-squared will always increase as you add more features to the model**, even if they are unrelated to the response. Thus, selecting the model with the highest R-squared is not a reliable approach for choosing the best linear model.\n",
    "\n",
    "There is alternative to R-squared called **adjusted R-squared** that penalizes model complexity (to control for overfitting), but it generally [under-penalizes complexity](http://scott.fortmann-roe.com/docs/MeasuringError.html).\n",
    "\n",
    "So is there a better approach to feature selection? **Cross-validation.** It provides a more reliable estimate of out-of-sample error, and thus is a better way to choose which of your models will best **generalize** to out-of-sample data. There is extensive functionality for cross-validation in scikit-learn, including automated methods for searching different sets of parameters and different models. Importantly, cross-validation can be applied to any model, whereas the methods described above only apply to linear models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression in scikit-learn\n",
    "\n",
    "Let's redo some of the Statsmodels code above in scikit-learn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.93888936946\n",
      "[ 0.04576465  0.18853002 -0.00103749]\n"
     ]
    }
   ],
   "source": [
    "# create X and y\n",
    "feature_cols = ['TV', 'Radio', 'Newspaper']\n",
    "X = data[feature_cols]\n",
    "y = data.Sales\n",
    "\n",
    "# follow the usual sklearn pattern: import, instantiate, fit\n",
    "from sklearn.linear_model import LinearRegression\n",
    "lm = LinearRegression()\n",
    "lm.fit(X, y)\n",
    "\n",
    "# print intercept and coefficients\n",
    "print lm.intercept_\n",
    "print lm.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('TV', 0.04576464545539765),\n",
       " ('Radio', 0.18853001691820462),\n",
       " ('Newspaper', -0.0010374930424762799)]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pair the feature names with the coefficients\n",
    "zip(feature_cols, lm.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.202667011892373"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# predict for a new observation\n",
    "lm.predict([100, 25, 25])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.89721063817895208"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the R-squared\n",
    "lm.score(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that **p-values** and **confidence intervals** are not (easily) accessible through scikit-learn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling Categorical Predictors with Two Categories\n",
    "\n",
    "Up to now, all of our predictors have been numeric. What if one of our predictors was categorical?\n",
    "\n",
    "Let's create a new feature called **Size**, and randomly assign observations to be **small or large**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "      <th>Size</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td> 230.1</td>\n",
       "      <td> 37.8</td>\n",
       "      <td> 69.2</td>\n",
       "      <td> 22.1</td>\n",
       "      <td> large</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>  44.5</td>\n",
       "      <td> 39.3</td>\n",
       "      <td> 45.1</td>\n",
       "      <td> 10.4</td>\n",
       "      <td> small</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>  17.2</td>\n",
       "      <td> 45.9</td>\n",
       "      <td> 69.3</td>\n",
       "      <td>  9.3</td>\n",
       "      <td> small</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td> 151.5</td>\n",
       "      <td> 41.3</td>\n",
       "      <td> 58.5</td>\n",
       "      <td> 18.5</td>\n",
       "      <td> small</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td> 180.8</td>\n",
       "      <td> 10.8</td>\n",
       "      <td> 58.4</td>\n",
       "      <td> 12.9</td>\n",
       "      <td> large</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  Radio  Newspaper  Sales   Size\n",
       "1  230.1   37.8       69.2   22.1  large\n",
       "2   44.5   39.3       45.1   10.4  small\n",
       "3   17.2   45.9       69.3    9.3  small\n",
       "4  151.5   41.3       58.5   18.5  small\n",
       "5  180.8   10.8       58.4   12.9  large"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# set a seed for reproducibility\n",
    "np.random.seed(12345)\n",
    "\n",
    "# create a Series of booleans in which roughly half are True\n",
    "nums = np.random.rand(len(data))\n",
    "mask_large = nums > 0.5\n",
    "\n",
    "# initially set Size to small, then change roughly half to be large\n",
    "data['Size'] = 'small'\n",
    "data.loc[mask_large, 'Size'] = 'large'\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For scikit-learn, we need to represent all data **numerically**. If the feature only has two categories, we can simply create a **dummy variable** that represents the categories as a binary value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "      <th>Size</th>\n",
       "      <th>IsLarge</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td> 230.1</td>\n",
       "      <td> 37.8</td>\n",
       "      <td> 69.2</td>\n",
       "      <td> 22.1</td>\n",
       "      <td> large</td>\n",
       "      <td> 1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>  44.5</td>\n",
       "      <td> 39.3</td>\n",
       "      <td> 45.1</td>\n",
       "      <td> 10.4</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>  17.2</td>\n",
       "      <td> 45.9</td>\n",
       "      <td> 69.3</td>\n",
       "      <td>  9.3</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td> 151.5</td>\n",
       "      <td> 41.3</td>\n",
       "      <td> 58.5</td>\n",
       "      <td> 18.5</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td> 180.8</td>\n",
       "      <td> 10.8</td>\n",
       "      <td> 58.4</td>\n",
       "      <td> 12.9</td>\n",
       "      <td> large</td>\n",
       "      <td> 1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  Radio  Newspaper  Sales   Size  IsLarge\n",
       "1  230.1   37.8       69.2   22.1  large        1\n",
       "2   44.5   39.3       45.1   10.4  small        0\n",
       "3   17.2   45.9       69.3    9.3  small        0\n",
       "4  151.5   41.3       58.5   18.5  small        0\n",
       "5  180.8   10.8       58.4   12.9  large        1"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a new Series called IsLarge\n",
    "data['IsLarge'] = data.Size.map({'small':0, 'large':1})\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's redo the multiple linear regression and include the **IsLarge** predictor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('TV', 0.045719820924362803),\n",
       " ('Radio', 0.18872814313427874),\n",
       " ('Newspaper', -0.0010976794483515545),\n",
       " ('IsLarge', 0.057423850854828061)]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create X and y\n",
    "feature_cols = ['TV', 'Radio', 'Newspaper', 'IsLarge']\n",
    "X = data[feature_cols]\n",
    "y = data.Sales\n",
    "\n",
    "# instantiate, fit\n",
    "lm = LinearRegression()\n",
    "lm.fit(X, y)\n",
    "\n",
    "# print coefficients\n",
    "zip(feature_cols, lm.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we interpret the **IsLarge coefficient**? For a given amount of TV/Radio/Newspaper ad spending, being a large market is associated with an average **increase** in Sales of 57.42 widgets (as compared to a Small market, which is called the **baseline level**).\n",
    "\n",
    "What if we had reversed the 0/1 coding and created the feature 'IsSmall' instead? The coefficient would be the same, except it would be **negative instead of positive**. As such, your choice of category for the baseline does not matter, all that changes is your **interpretation** of the coefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling Categorical Predictors with More than Two Categories\n",
    "\n",
    "Let's create a new feature called **Area**, and randomly assign observations to be **rural, suburban, or urban**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "      <th>Size</th>\n",
       "      <th>IsLarge</th>\n",
       "      <th>Area</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td> 230.1</td>\n",
       "      <td> 37.8</td>\n",
       "      <td> 69.2</td>\n",
       "      <td> 22.1</td>\n",
       "      <td> large</td>\n",
       "      <td> 1</td>\n",
       "      <td>    rural</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>  44.5</td>\n",
       "      <td> 39.3</td>\n",
       "      <td> 45.1</td>\n",
       "      <td> 10.4</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "      <td>    urban</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>  17.2</td>\n",
       "      <td> 45.9</td>\n",
       "      <td> 69.3</td>\n",
       "      <td>  9.3</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "      <td>    rural</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td> 151.5</td>\n",
       "      <td> 41.3</td>\n",
       "      <td> 58.5</td>\n",
       "      <td> 18.5</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "      <td>    urban</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td> 180.8</td>\n",
       "      <td> 10.8</td>\n",
       "      <td> 58.4</td>\n",
       "      <td> 12.9</td>\n",
       "      <td> large</td>\n",
       "      <td> 1</td>\n",
       "      <td> suburban</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  Radio  Newspaper  Sales   Size  IsLarge      Area\n",
       "1  230.1   37.8       69.2   22.1  large        1     rural\n",
       "2   44.5   39.3       45.1   10.4  small        0     urban\n",
       "3   17.2   45.9       69.3    9.3  small        0     rural\n",
       "4  151.5   41.3       58.5   18.5  small        0     urban\n",
       "5  180.8   10.8       58.4   12.9  large        1  suburban"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set a seed for reproducibility\n",
    "np.random.seed(123456)\n",
    "\n",
    "# assign roughly one third of observations to each group\n",
    "nums = np.random.rand(len(data))\n",
    "mask_suburban = (nums > 0.33) & (nums < 0.66)\n",
    "mask_urban = nums > 0.66\n",
    "data['Area'] = 'rural'\n",
    "data.loc[mask_suburban, 'Area'] = 'suburban'\n",
    "data.loc[mask_urban, 'Area'] = 'urban'\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have to represent Area numerically, but we can't simply code it as 0=rural, 1=suburban, 2=urban because that would imply an **ordered relationship** between suburban and urban (and thus urban is somehow \"twice\" the suburban category).\n",
    "\n",
    "Instead, we create **another dummy variable**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TV</th>\n",
       "      <th>Radio</th>\n",
       "      <th>Newspaper</th>\n",
       "      <th>Sales</th>\n",
       "      <th>Size</th>\n",
       "      <th>IsLarge</th>\n",
       "      <th>Area</th>\n",
       "      <th>Area_suburban</th>\n",
       "      <th>Area_urban</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td> 230.1</td>\n",
       "      <td> 37.8</td>\n",
       "      <td> 69.2</td>\n",
       "      <td> 22.1</td>\n",
       "      <td> large</td>\n",
       "      <td> 1</td>\n",
       "      <td>    rural</td>\n",
       "      <td> 0</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>  44.5</td>\n",
       "      <td> 39.3</td>\n",
       "      <td> 45.1</td>\n",
       "      <td> 10.4</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "      <td>    urban</td>\n",
       "      <td> 0</td>\n",
       "      <td> 1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>  17.2</td>\n",
       "      <td> 45.9</td>\n",
       "      <td> 69.3</td>\n",
       "      <td>  9.3</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "      <td>    rural</td>\n",
       "      <td> 0</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td> 151.5</td>\n",
       "      <td> 41.3</td>\n",
       "      <td> 58.5</td>\n",
       "      <td> 18.5</td>\n",
       "      <td> small</td>\n",
       "      <td> 0</td>\n",
       "      <td>    urban</td>\n",
       "      <td> 0</td>\n",
       "      <td> 1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td> 180.8</td>\n",
       "      <td> 10.8</td>\n",
       "      <td> 58.4</td>\n",
       "      <td> 12.9</td>\n",
       "      <td> large</td>\n",
       "      <td> 1</td>\n",
       "      <td> suburban</td>\n",
       "      <td> 1</td>\n",
       "      <td> 0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      TV  Radio  Newspaper  Sales   Size  IsLarge      Area  Area_suburban  \\\n",
       "1  230.1   37.8       69.2   22.1  large        1     rural              0   \n",
       "2   44.5   39.3       45.1   10.4  small        0     urban              0   \n",
       "3   17.2   45.9       69.3    9.3  small        0     rural              0   \n",
       "4  151.5   41.3       58.5   18.5  small        0     urban              0   \n",
       "5  180.8   10.8       58.4   12.9  large        1  suburban              1   \n",
       "\n",
       "   Area_urban  \n",
       "1           0  \n",
       "2           1  \n",
       "3           0  \n",
       "4           1  \n",
       "5           0  "
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create three dummy variables using get_dummies, then exclude the first dummy column\n",
    "area_dummies = pd.get_dummies(data.Area, prefix='Area').iloc[:, 1:]\n",
    "\n",
    "# concatenate the dummy variable columns onto the original DataFrame (axis=0 means rows, axis=1 means columns)\n",
    "data = pd.concat([data, area_dummies], axis=1)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how we interpret the coding:\n",
    "- **rural** is coded as Area_suburban=0 and Area_urban=0\n",
    "- **suburban** is coded as Area_suburban=1 and Area_urban=0\n",
    "- **urban** is coded as Area_suburban=0 and Area_urban=1\n",
    "\n",
    "Why do we only need **two dummy variables, not three?** Because two dummies captures all of the information about the Area feature, and implicitly defines rural as the baseline level. (In general, if you have a categorical feature with k levels, you create k-1 dummy variables.)\n",
    "\n",
    "If this is confusing, think about why we only needed one dummy variable for Size (IsLarge), not two dummy variables (IsSmall and IsLarge).\n",
    "\n",
    "Let's include the two new dummy variables in the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('TV', 0.04574401036331379),\n",
       " ('Radio', 0.18786669552525814),\n",
       " ('Newspaper', -0.0010876977267108706),\n",
       " ('IsLarge', 0.077396607497479411),\n",
       " ('Area_suburban', -0.10656299015958708),\n",
       " ('Area_urban', 0.26813802165220019)]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create X and y\n",
    "feature_cols = ['TV', 'Radio', 'Newspaper', 'IsLarge', 'Area_suburban', 'Area_urban']\n",
    "X = data[feature_cols]\n",
    "y = data.Sales\n",
    "\n",
    "# instantiate, fit\n",
    "lm = LinearRegression()\n",
    "lm.fit(X, y)\n",
    "\n",
    "# print coefficients\n",
    "zip(feature_cols, lm.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we interpret the coefficients?\n",
    "- Holding all other variables fixed, being a **suburban** area is associated with an average **decrease** in Sales of 106.56 widgets (as compared to the baseline level, which is rural).\n",
    "- Being an **urban** area is associated with an average **increase** in Sales of 268.13 widgets (as compared to rural).\n",
    "\n",
    "**A final note about dummy encoding:** If you have categories that can be ranked (i.e., strongly disagree, disagree, neutral, agree, strongly agree), you can potentially use a single dummy variable and represent the categories numerically (such as 1, 2, 3, 4, 5)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What Didn't We Cover?\n",
    "\n",
    "- Detecting collinearity\n",
    "- Diagnosing model fit\n",
    "- Transforming predictors to fit non-linear relationships\n",
    "- Interaction terms\n",
    "- Assumptions of linear regression\n",
    "- And so much more!\n",
    "\n",
    "You could certainly go very deep into linear regression, and learn how to apply it really, really well. It's an excellent way to **start your modeling process** when working a regression problem. However, it is limited by the fact that it can only make good predictions if there is a **linear relationship** between the features and the response, which is why more complex methods (with higher variance and lower bias) will often outperform linear regression.\n",
    "\n",
    "Therefore, we want you to understand linear regression conceptually, understand its strengths and weaknesses, be familiar with the terminology, and know how to apply it. However, we also want to spend time on many other machine learning models, which is why we aren't going deeper here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources\n",
    "\n",
    "- To go much more in-depth on linear regression, read Chapter 3 of [An Introduction to Statistical Learning](http://www-bcf.usc.edu/~gareth/ISL/), from which this lesson was adapted. Alternatively, watch the [related videos](http://www.dataschool.io/15-hours-of-expert-machine-learning-videos/) or read my [quick reference guide](http://www.dataschool.io/applying-and-interpreting-linear-regression/) to the key points in that chapter.\n",
    "- To learn more about Statsmodels and how to interpret the output, DataRobot has some decent posts on [simple linear regression](http://www.datarobot.com/blog/ordinary-least-squares-in-python/) and [multiple linear regression](http://www.datarobot.com/blog/multiple-regression-using-statsmodels/).\n",
    "- This [introduction to linear regression](http://people.duke.edu/~rnau/regintro.htm) is much more detailed and mathematically thorough, and includes lots of good advice.\n",
    "- This is a relatively quick post on the [assumptions of linear regression](http://pareonline.net/getvn.asp?n=2&v=8)."
   ]
  }
 ],
 "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
}