code/08.06-statsmodels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Statistical Modeling with Python\n",
    "\n",
    "`statsmodels` is better suited for traditional stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# the statsmodels.api uses numpy array notation\n",
    "# statsmodels.formula.api use formula notation (similar to R's formula notation)\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A minimal OLS example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Four pairs of points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([1,2,3,4])\n",
    "y = np.array([2,6,4,8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAADEtJREFUeJzt3F1oZPUZx/HfL5kN6qqYpqlYlWhAVsSLagY77YKIL8WqaC96oajY4rI3fdHaIvZKelHoRRF7IYWw2lIMSlFLxYJV7EoRnNUZXetLtEpqdOu2G9O0ainNpvP0IiORdV+SOWf3bJ75fiAkk5yZ83DQL4f/nnMcEQIArH8DVQ8AACgHQQeAJAg6ACRB0AEgCYIOAEkQdABI4pBBt32f7T22X/nE7z5j+0nbb3a/Dx/eMQEAh7KaM/RfSrp8n9/dIempiDhL0lPd1wCACnk1NxbZPkPSYxFxbvf1G5Iuiojdtk+R9HREbDqcgwIADq7W4/tOjojdktSN+ucOtKHtrZK2StLGjRsnzj777B53CQD9qd1uvx8Ro4fartegr1pETEqalKR6vR6tVutw7xIAUrE9u5rter3K5e/dpRZ1v+/p8XMAACXpNeiPSrqp+/NNkn5bzjgAgF6t5rLFByQ9K2mT7V22b5b0E0mX2X5T0mXd1wCACh1yDT0irjvAny4peRYAQAHcKQoASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0AShYJu+3u2X7X9iu0HbB9T1mAAgLXpOei2T5X0XUn1iDhX0qCka8saDED/as8u6J7tb6k9u1D1KOtKrYT3H2t7r6TjJL1XfCQA/aw9u6DrtzW1uNTRUG1AU1samhgbrnqsdaHnM/SI+Kukn0p6R9JuSf+KiCf23c72Vtst2625ubneJwXQF5oz81pc6qgT0t6ljpoz81WPtG4UWXIZlnSNpDMlfV7SRts37LtdRExGRD0i6qOjo71PCqAvNMZHNFQb0KClDbUBNcZHqh5p3Siy5HKppL9ExJwk2X5E0pcl3V/GYAD608TYsKa2NNScmVdjfITlljUoEvR3JDVsHyfpP5IukdQqZSoAfW1ibJiQ96DIGvoOSQ9JekHSy93PmixpLgDAGhW6yiUi7pR0Z0mzAAAK4E5RAEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0AShYJu+yTbD9l+3fa07S+VNRgAYG1qBd//M0mPR8TXbQ9JOq6EmdBH2rMLas7MqzE+oomx4arHAda1noNu+0RJF0r6hiRFxKKkxXLGQj9ozy7o+m1NLS51NFQb0NSWBlEHCiiy5DIuaU7SL2y/aHub7Y37bmR7q+2W7dbc3FyB3SGb5sy8Fpc66oS0d6mj5sx81SMB61qRoNcknS/p5xFxnqR/S7pj340iYjIi6hFRHx0dLbA7ZNMYH9FQbUCDljbUBtQYH6l6JGBdK7KGvkvSrojY0X39kPYTdOBAJsaGNbWlwRo6UJKegx4Rf7P9ru1NEfGGpEskvVbeaOgHE2PDhBwoSdGrXL4jaap7hcuMpG8WHwkA0ItCQY+InZLqJc0CACiAO0UBIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIonDQbQ/aftH2Y2UMBADoTRln6LdImi7hc4C+1p5d0D3b31J7dqHqUbBO1Yq82fZpkq6U9GNJt5UyEdCH2rMLun5bU4tLHQ3VBjS1paGJseGqx8I6U/QM/W5Jt0vqHGgD21ttt2y35ubmCu4OyKk5M6/FpY46Ie1d6qg5M1/1SFiHeg667ask7YmI9sG2i4jJiKhHRH10dLTX3QGpNcZHNFQb0KClDbUBNcZHqh4J61CRJZfNkq62fYWkYySdaPv+iLihnNGA/jExNqypLQ01Z+bVGB9huQU9cUQU/xD7Ikk/iIirDrZdvV6PVqtVeH8A0E9styOifqjtuA4dAJIodJXLxyLiaUlPl/FZAIDecIYOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEii56DbPt32dtvTtl+1fUuZgwEA1qZW4L1Lkr4fES/YPkFS2/aTEfFaSbOl1Z5dUHNmXo3xEU2MDVc9DoAkeg56ROyWtLv784e2pyWdKomgH0R7dkHXb2tqcamjodqAprY0iDqAUpSyhm77DEnnSdqxn79ttd2y3Zqbmytjd+tac2Zei0sddULau9RRc2a+6pEAJFE46LaPl/SwpFsj4oN9/x4RkxFRj4j66Oho0d2te43xEQ3VBjRoaUNtQI3xkapHApBEkTV02d6g5ZhPRcQj5YyU28TYsKa2NFhDB1C6noNu25LulTQdEXeVN1J+E2PDhBxA6YosuWyWdKOki23v7H5dUdJcAIA1KnKVyzOSXOIsAIACuFMUAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAg6ACRB0AEgCYIOAEkQdABIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASAJgg4ASRB0AEiCoANAEgQdAJIg6ACQBEEHgCQIOgAkQdABIAmCDgBJEHQASIKgA0ASBB0AkiDoAJAEQQeAJAoF3fbltt+w/ZbtO8oaCgCwdj0H3fagpHskfVXSOZKus31OWYMBANamyBn6BZLeioiZiFiU9KCka8oZCwCwVrUC7z1V0rufeL1L0hf33cj2Vklbuy//a/uVAvvM5LOS3q96iKMEx2IFx2IFx2LFptVsVCTo3s/v4lO/iJiUNClJtlsRUS+wzzQ4Fis4Fis4Fis4Fitst1azXZEll12STv/E69MkvVfg8wAABRQJ+vOSzrJ9pu0hSddKerScsQAAa9XzkktELNn+tqTfSxqUdF9EvHqIt032ur+EOBYrOBYrOBYrOBYrVnUsHPGpZW8AwDrEnaIAkARBB4AkjkjQeUTACtv32d7T79fj2z7d9nbb07ZftX1L1TNVxfYxtp+z/VL3WPyo6pmqZnvQ9ou2H6t6lirZftv2y7Z3rubSxcO+ht59RMCfJV2m5Usdn5d0XUS8dlh3fJSyfaGkjyT9KiLOrXqeqtg+RdIpEfGC7RMktSV9rR//u7BtSRsj4iPbGyQ9I+mWiGhWPFplbN8mqS7pxIi4qup5qmL7bUn1iFjVDVZH4gydRwR8QkT8UdI/qp6jahGxOyJe6P78oaRpLd993Hdi2Ufdlxu6X317tYLt0yRdKWlb1bOsN0ci6Pt7REBf/o+L/bN9hqTzJO2odpLqdJcYdkraI+nJiOjbYyHpbkm3S+pUPchRICQ9YbvdfYzKQR2JoK/qEQHoT7aPl/SwpFsj4oOq56lKRPwvIr6g5TuuL7Ddl8txtq+StCci2lXPcpTYHBHna/mptt/qLtke0JEIOo8IwH5114sfljQVEY9UPc/RICL+KelpSZdXPEpVNku6urt2/KCki23fX+1I1YmI97rf90j6jZaXsA/oSASdRwTgU7r/EHivpOmIuKvqeapke9T2Sd2fj5V0qaTXq52qGhHxw4g4LSLO0HIr/hARN1Q8ViVsb+xeMCDbGyV9RdJBr4477EGPiCVJHz8iYFrSr1fxiIC0bD8g6VlJm2zvsn1z1TNVZLOkG7V8Braz+3VF1UNV5BRJ223/ScsnQE9GRF9frgdJ0smSnrH9kqTnJP0uIh4/2Bu49R8AkuBOUQBIgqADQBIEHQCSIOgAkARBB4AkCDoAJEHQASCJ/wPk3NKhIeUPUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x,y, marker = '.')\n",
    "plt.xlim(0,5)\n",
    "plt.ylim(0,10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   x  y\n",
      "0  1  2\n",
      "1  2  6\n",
      "2  3  4\n",
      "3  4  8\n"
     ]
    }
   ],
   "source": [
    "# make a dataframe of our data\n",
    "d = pd.DataFrame({'x':x, 'y':y})\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Seaborn lmplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x17581f122e8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XtwXPd53vHvu1fsYnHlRaRIgDZjxrZsWRRJMUk9k8p27MqqYzWtp0NPm8hOOnTS0VidNk2cpI1TN51e0lzUOrWi2JrYqWM649qR4sgXJWqiZBpbpKiLZUu2ZEYkIVIicSNui72+/WMX0BJcAAtgd89ens8MB4vdA+DlIfHgxe+c8x5zd0REpPlCQRcgItKtFMAiIgFRAIuIBEQBLCISEAWwiEhAFMAiIgFRAIuIBEQBLCISEAWwiEhAIkEX0Ai33Xabf/WrXw26DBHpXlbLRh3ZAY+PjwddgojIujoygEVE2oECWEQkIApgEZGAKIBFRAKiABYRCUjDA9jM7jezS2b2TMVzv2ZmL5nZk+U/t6/ysbeZ2XfN7AUz+0ijaxURaaZmdMB/ANxW5fnfdveD5T8PrXzRzMLA7wLvBm4A3m9mNzS0UhGRJmp4ALv7o8DkJj70KPCCu59x9yxwArijrsWJiAQoyDXgu8zs6fISxVCV1/cA5yveHys/V5WZHTezU2Z26vLly/WuVUSk7oIK4E8APwAcBC4Cv1llm2qX8q16B1F3v8/dj7j7kR07dtSnShGRBgokgN39FXcvuHsR+H1Kyw0rjQEjFe/vBS40oz4RkWYIJIDNbHfFuz8BPFNls5PAATN7rZnFgGPAg82oT0Rks9xX/UX9Gg2fhmZmnwNuBbab2RjwUeBWMztIaUnhReBD5W2vBz7p7re7e97M7gK+BoSB+939242uV0RksxZzBS7PZhgZTta0vW0krdvFkSNH/NSpU0GXISJdolh0JheyzKRzAOzfkappHGVHzgMWEWmWdLbA+FyGXKG44Y9VAIuIbMLKrnczFMAiIhu0la63kgJYRKRGxaIzMZ9ldnHzXW8lBbCISA0WsnnGZ7Pki1vreispgEVE1lAoOhPzGeYW83X/3ApgEZFVzGfyTMzVt+utpAAWEVmhUHQm5jLMZerf9VZSAIuIVJjP5Bmfy1AoNv4iNQWwiAjN63orKYBFpOvNZfJMNKnrraQAFpGulS8UmZjPMt/ErreSAlhEutLsYo7J+WzTu95KCmAR6Sr5QpHxuSwL2WC63koKYBHpGjOLOSbnshRbZAyvAlhEOl4rdb2VFMAi0tFareutpAAWkY6UKxQZn8uQzhaCLmVVCmAR6ThX0jmm5luz662kABaRjpErFLk8m2Ex17pdbyUFsIh0hCsLOSYXshu6LXzQFMAi0tay+dJab7t0vZUUwCLSttqx662kABaRtpPNF7k8lyHThl1vJQWwiLQNdy+d4bCQa9uut5ICWETaQiZfYHwu2/Zdb6WGB7CZ3Q+8B7jk7m8uP/cbwI8DWeD7wAfdfbrKx74IzAIFIO/uRxpdr4i0FndneiHHdLozut5KoSZ8jT8Ablvx3MPAm939LcD3gF9a4+Pf5u4HFb4i3SeTL/DSdJqpNj7QtpaGB7C7PwpMrnju6+6+NBXjG8DeRtchIu3D3Zmcz3JhepFsvjF3JG4FzeiA1/PTwFdWec2Br5vZ42Z2vIk1iUhAFnOlrne6Q7veSoEehDOzXwHywGdX2eSt7n7BzHYCD5vZc+WOutrnOg4cBxgdHW1IvSLSOO7O1EKO6YVs0KU0TWAdsJndSeng3D/zVX7MufuF8ttLwJeAo6t9Pne/z92PuPuRHTt2NKJkEWmQxVyBsal0V4UvBBTAZnYb8IvAe919YZVtes2sb+kx8C7gmeZVKSKN5l66FfyF6TS5Queu9a6m4QFsZp8D/hZ4vZmNmdnPAB8H+igtKzxpZveWt73ezB4qf+h1wN+Y2VPAY8CfuftXG12viDTHUtd7JZ0LupTANHwN2N3fX+XpT62y7QXg9vLjM8BNDSxNRALg7kzMZ5np4uBdoivhRKRp0tkC43OZrlxuqEYBLCINVyyWut7ZRXW9lRTAItJQ6npXpwAWkYZQ17s+BbCI1N1CNs/4bJZ8UV3vWhTAIlI3haIzMZ9hbjG//saiABaR+pjP5JmYU9e7EQpgEdmSQrF0NdtcRl3vRimARWTT5jN5xucyFIqdPbWsURTAIrJh6nrrQwEsIhsyl8kzoa63LhTAIlKTQtEZn8swr663bhTAIrKu2cUck/NZdb11pgAWkVXlC0XG57IsZNX1NoICWESqml3MMTGXpdjh92ULkgJYRK6irrd5FMAismxmMcekut6mUQCLCLlCkfG5DOlsIehSuooCWKTLXUnnmJpX1xsEBbBIl8oVilyezbCYU9cbFAWwSBe6spBjciGLq+sNlAJYpItk86W1XnW9rUEBLNIl1PW2HgWwSIfL5otcnsuQUdfbchTAIh3K3UtnOCzk1PW2KAWwSAfK5Atcns2Qzev2QK1MASzSQdyd6YUc02l1ve0g1IwvYmb3m9klM3um4rlhM3vYzJ4vvx1a5WPvLG/zvJnd2Yx6RdpRJl/gpek0UzrQ1jaaEsDAHwC3rXjuI8BfuPsB4C/K71/FzIaBjwI/BBwFPrpaUIt0K3dncj7LhelFLTm0maYEsLs/CkyuePoO4NPlx58G/lGVD/0HwMPuPunuU8DDXBvkIl1rMVfqeqfV9balINeAr3P3iwDuftHMdlbZZg9wvuL9sfJz1zCz48BxgNHR0TqXKtJa3J2phRzTC9mgS5EtaNYSxGZZleeq/ph39/vc/Yi7H9mxY0eDyxIJzmKuwNhUWuHbAYIM4FfMbDdA+e2lKtuMASMV7+8FLjShNpGW4166FfyF6TS5gtZ6O0GQAfwgsHRWw53AA1W2+RrwLjMbKh98e1f5OZGustT1Xknngi5F6qhZp6F9Dvhb4PVmNmZmPwP8F+CdZvY88M7y+5jZETP7JIC7TwL/EThZ/vOx8nMiXcG9dCt4db2dyTrxyOmRI0f81KlTQZchsiXpbIHxuYyCtw3t35GqdvzqGroSTqTFFIvO5EKWGS03dDwFsEgLUdfbXRTAIi2gWHQm5rPMLqrr7SYKYJGALWTzjM9myRfV9XYbBbBIQIpFZ3w+w9xiPuhSJCAKYJEAzGfyTMyp6+12CmCRJioUS1ezzWXU9YoCWKRp5jN5xucyFIqdd+69bI4CWKTB1PXKahTAIg00l8kzoa5XVqEAFmmAQrE0w2FeXa+sQQEsUmezizkm57PqemVdCmCROskXiozPZVnIquuV2iiARepgZjHH5FyWYgdOF5TGUQCLbIG6XtkKBbDIJqnrla1SAItsUK5QZHwuQzpbCLoUaXMdGcD5gjOfyZOMhTGraTC9SE2upHNMzavrlfroyAAuuPPKzCIhM5KxML3xiMJYtiRXKHJ5NsNiTl2v1E9HBvCSojtzmTxzmfxyGCfjEZLRMKGQwlhqc2Uhx+RClk68f6Jc7bEzk5w4eZ6LM2l29yc4dssIR/cPN+zrdXQAV6oMY6vsjBXGsopsvrTWq663Ozx2ZpJ7HnmeSMjo74kwMZ/hnkee524ONCyEuyaAK7mX1ojnFcayCnW93efEyfNEQkYiGgYgEQ2TzhU4cfK8ArhRqoVxMhamNxZRGHehbL7I5bkMGXW9XefiTJr+nqsjsSca4uWZdMO+ZtcHcKXKMB63LIlomN64wrhbTC9kmVrIqevtUrv7E0zMZ5Y7YIDFXJFd/YmGfc1Qwz5zm3N3FrJ5Ls9mODu5wMtXFplZzGnASgfK5Au8NJ1mcl5LDt3s2C0j5ItOOlfAKb3NF51jt4w07GsqgGuwFMbjsxnOKYw7hrszNZ/lwvSilhyEo/uHufvtB9jWG2d2Mc+23jh3v732A3DuzkvTab789IWav2ZgSxBm9nrg8xVP7Qd+1d1/p2KbW4EHgL8rP/VFd/9Y04qsYimMF7J5JixLTzREbzxCbyxCWMsUbSOTL3B5NkM2r5tiyquO7h/e0AG3qYUsT5yb5vS5KU6fneblmUUAPvyOH6zp4wMLYHf/LnAQwMzCwEvAl6ps+tfu/p5m1lYrdyedLZDOFhgnQ6J8NoXCuHW5O1MLOa6ktdYrG5fOFnj6pWlOny2F7vcvz1+zzcoDeWtplYNw7wC+7+5ngy5kKyrDuCe6FMZhImGt9LSCxVyB8Tl1vVK7fKHIcy/Pljrcc9N858IM+RVLj/FIiBv3DHBo3xCHRgd53c5UzZ+/VQL4GPC5VV77ETN7CrgA/Ly7f7vaRmZ2HDgOcP3exi2a12oxV2AxV2ACFMYBW+p6pxeyQZciLc7deXFigcfPTnH63BRPj11hYcXQpZDBG3b1cfNoKXDfdP0Ascjmvq8t6F/DzCxGKVzf5O6vrHitHyi6+5yZ3Q7c4+4H1vucNx485A88/GhjCt6inmjptLbeuMK4GRZzpbXeXEFdr1T3yswip89Nc7oculMLuWu22TecXO5wb9o7SGqdZYb9O1I1rUG2Qgf8buD0yvAFcPeZiscPmdn/MrPt7j7e1ArraLkznod4NExKYdwQ7s7kfJYr6Wu/maS7zS7meOL8q+u4Y1PXXmixPRXj0OjQcuhuT8UbUksrBPD7WWX5wcx2Aa+4u5vZUUqnzU00s7hGyuQKZFaEcTIeJqow3hJ1vVIpkyvwzIWZ5TMVvvfKLCt/7++Nhzk4Msih0SEOjw4xMpxoyvTEQAPYzJLAO4EPVTz3swDufi/wPuDnzCwPpIFjHvSaSYOsDOPe8hkVCuPauTsT81lm1PV2tULR+d4rszxxbprHz03xzEtXyBWujo1o2HjT9QMc3lcK3R+8ri+QM5cCDWB3XwC2rXju3orHHwc+3uy6grYUxpPzWWKREKl4RGG8jnS2dIaDut7u4+6cn0qX13CnefL8NHOZq+/RZ8CB61KlZYXRQd68Z4CeikuOg9IKSxCyhmy+yGQ+e1UYJ2ORTR917TTFojO5oK6320zMZUoHzsrLCpfnMtdss2cwwaFyh3vzyCD9iWgAla5NAdxGVoZx6WyK7g1jdb3dYy6T56nz08uhe3Zi4ZptBhNRbh4d5PC+0sGzXf09AVS6MQrgNpXNF8nms0wtZImGy51xPEw8EvyvVY1WLJbWemcX1fV2qmy+yLMXZ3j83BSnz07x3MuzrBy90hMN8Za9pcA9PDrIa7f3tt1txxTAHSBXKDK1kGVqAaLh8myKDg3j0lCkLPmiut5OUnTn+5fmljvcp8eukFlxxWI4ZNywu3QBxOHRId6wu6/tj4sogDtMrlBkeiHLdIeFcbHojM9nmFvMr7+xtDx358KVRZ44N8XjZ6d54twUM1X+bffv6OXQaGkd96a9gyRi7f3/eCUFcAerFsbJWLgljv5uxHwmz8Scut52N12eHPb4islhlXb2xUtruKND3Dw6yHBvLIBKm0cB3CVWhvHSffBaOYwLRWdiLnPNKUXSHtK5Ak+PrT857ODoIIdHS6F7/WBP263jboUCuAvlCkWupItcSedaNoznM3nG5zIaet9GNjQ5bHSQQ/uGeN3OFKEuCtyVFMBdrjKMI6FQ6R54AYaxut72sXJy2FPnr5DOXTs57PW7+pYvgNjK5LBOpACWZfni1WGcjIdJNTGM5zJ5JtT1trRLM4s8Xp4c9sT5aSbnrx3xuW84uXw+bi2Tw7qZ9oxUlS8WmUkXmWlCGBeKzvhchnl1vS1nJp3jybG1J4dtS8XKa7ilZYVGTQ7rRApgWVdlGIdDRjIWIRWP1OWUoNnFHJPzWXW9LWJpctjSssLzr8xVnxy2d3B5VOPocLKrDpzVkwJYNqRQdGYXc8wuXh3GPdHQhr4J84Ui43NZFrLqeoNUKDrPX5rl9Nn1Jof1l0Y17gtuclgnUgDLptUSxn/53CV+79EznJ9aYGQoyYd+dD+HXzPExFyWYmdOFm1pS5PDli6AaKfJYZ1IASx1US2MT784ya8/9CyxSIjBRJRXZtL88p98iw+/7cCGbv0tWzNenhxWCt0pxueuPXBWOTns4MggAy04OawTKYCl7pbC+BN/dQaASChEruBEwiHCBefEyfMK4AZamhy2dNVZtclhQ8loeabCIDe3yeSwTqQAloa5OJOmvydC5U1MeqIhXp659ki6bN7S5LDT5WWF516eqTo57Kby5LBDbTo5rBMpgKVhdvcnmJjPkKhYP1zMFdnVnwiwqvZXdOfM5XkePzvFE+XJYYurTA47VL7EtxMmh3UiBbA0zLFbRrjnkedJ5wr0REMs5orki86xW0aCLq3tXJhOL9/94Ynz01Xv9rx/e+/yOu5b9g6QjOnbu9XpX0ga5uj+Ye7mACdOnuflmTS7+hMcu2VE6781qJwc9sS5aS5euXZy2HX98eUOtxsmh3UiBbA01NH9wwrcGqSzBZ5+SZPDuo0CWCQAlZPDHj87zbMXr50cFouEeIsmh3U0BbBIEyxNDjtdPhdXk8MEFMAiDfPKzOLyBRCnz60+OWxppoImh3Uf/WuL1MlMOseTFbdOX3Ny2L4hbh4ZZEefJod1MwWwyCZVTg574tw033tltvrksJHSqWGHR4cYGU7owJksUwCL1Khyctjpc1N8a43JYUs3ltTkMFlL4AFsZi8Cs0AByLv7kRWvG3APcDuwAHzA3U83u07pPrVODnvdzhSH95XOxb1Rk8NkAwIP4LK3ufv4Kq+9GzhQ/vNDwCfKb0XqbqI8OWzpqrPLc5lrtrl+sGd5HVeTw2Qr1g1gM7sL+Ky7TzWhnmruAD7jpYku3zCzQTPb7e4XA6pHOsh8Js9T5VvurDU57ODI4PKywq4BTQ6T+qilA94FnDSz08D9wNfc6zpJ24Gvm5kDv+fu9614fQ9wvuL9sfJzVwWwmR0HjgNcv1ezBqS6WieHHRwZXB7XqMlh0ijrBrC7/zsz+/fAu4APAh83sz8GPuXu369DDW919wtmthN42Myec/dHK16v9j//mh8A5eC+D+DGg4d0qwUBNjc57I27+4hocpg0QU1rwO7uZvYy8DKQB4aAL5jZw+7+C1spwN0vlN9eMrMvAUeBygAeAypb2r3Aha18Telsmhwm7aKWNeAPA3cC48AngX/r7jkzCwHPA5sOYDPrBULuPlt+/C7gYys2exC4y8xOUDr4dkXrv1KpcnLY6bPTvDyz+uSwpbMVhpKaHCbBq+XH/nbgH7v72con3b1oZu/Z4te/DvhSeX0tAvyRu3/VzH62/DXuBR6idAraC5ROQ/vgFr+mtLl0rsC3xq4s3zpdk8OkXVl9j6e1hhsPHvIHHn50/Q2lLVRODjt9bprvXKg+OezGPQMc1uQwaQH7d6Rq+s+nhS9pOSsnhz09doWFrCaHSedRAEtLuFSeHHZ6nclhN4+WzsfV5DDpBPofLIGYXczxxPlX7wCx2uSwpQ730OiQJodJx1EAS1MsTQ5bOj1Mk8NEFMDSILVPDivdcufwPk0Ok+6jAJa6cHfGptLLl/hqclhjPXZmkhMnz3NxJs1u3W26bSmAZdNqnhxWHmKjyWH18diZSe555HkiIaO/J8LEfIZ7HnmeuzmgEG4zCmCp2Vwmz1Pnp5evOlttctjNFQfONDms/k6cPE8kZCTKvz0komHSuQInTp5XALcZBbCsamly2NIlvqtNDrtp7+DyjSX3a3JYw12cSdPfE1leL3dKIfzyzLVnkkhrUwDLsqI73780t7ys8K1VJoe9cVffcuC+cXc/UU0Oa6rrBxJMp7P0RF/d7wuFPK/Z1su2VJz5TJ7FFbe8l9akAO5ytU4OW7oAQpPDgmNmDCSifPjtr+Ojf/odFrL55eWHXMH52b//AwwkogwkouQLReazBRayedJZhXGr0ndSl5kqTw47vcbksJ198eUDZzePDjLcq8lhQYtFQuzoixOPhHnbG6/DzPi9R88wNrXA3qEkH/rR/dz6hp3L20fCIQYSIQYSUQpFZz6bZz6jMG41CuAOl84WePqlV684W3Vy2EhpHfewJoe1FDNjMBFlMBm96t/k1jfsvCpw1xIOGf09Ufp7rg7jxVyRThzG1U4UwB2m1slhb9lTugBCk8NaV2XXWy/VwnghUyCdKyiMA6AAbnOVk8NOn53mqbHpNSeHHd43xA27+zU5rIWZGUPJ0lpuI38TWRnGC9k88wrjplIAt6FXypPDnlhjctjocHL5XNyDI5oc1i7i0TA7UvGm/4AMh4y+nih9PVGKy8sUCuNG03dlG9jI5LDDo6W7+WpyWHtZ6noHW+BWSaEqYbyQLbCQVRjXmwK4BWVyBb594dULIFadHLb31QNnmhzWvoLqemuxMowXcgXmM3mFcZ0ogFtAoeh875XZ5Ut8n9HksK5gZgwnYwwk22M+RihkpOIRUvHIchgvZPLMK4w3TQEcAHfn/FSaJzQ5rGv1RMNsb9GutxaVYezupYs+yp1xUWFcMwVwk2hymACEzBjqjXXUv63Z1WG8kH11mUJhvDYFcIPUOjns4MggR/YNcfO+IXb1a3JYJ0vESl1vJ8/OMDN64xF6K8O4fK6xwvhaCuA62dDksPI67ms1OawrhMwYTsXo7+mcrrcWV4VxyknnCsxlFMaVFMCbVHTnzOV5Hj87tebksBt29y3Px9XksO6TiJXOcIh0+b+7mZGMRUjGrg7jdLZAYWWn0kUUwBugyWFSq27temtxVRh7KYznM6XJbd0WxkqHNUyXJ4c9XtPksNIFEJocJslYhO2pWNd3vbW4OoxjLOaKpWWKLgnjwALYzEaAzwC7gCJwn7vfs2KbW4EHgL8rP/VFd/9Yo2pK5wp8a+zK8rLCqpPDypf4HhodZM+gLoCQknDIGO6N0aeud1PMjEQsTCIWBuKks4WOD+MgO+A88G/c/bSZ9QGPm9nD7v6dFdv9tbu/pxEFFIrOsxdnlrvc1SaH3bhngMOaHCZr6I1H2NarrreeVobx0tkU+WJx3Y9tF4EFsLtfBC6WH8+a2bPAHmBlANfzay5PDnv87BRPj11Zc3LYodFB3nT9QNueLC+NFw4Z21JxUnGt5jXSchinYLHibIp2D+OW+F9jZq8Bbga+WeXlHzGzp4ALwM+7+7dX+RzHgeMA1+8dWX6+lslh+4aT3KzJYbJBvfEI21NxXRLeZD3RcOmq0HIYz2dKk9vaMYwt6Gu4zSwF/BXwn9z9iyte6weK7j5nZrcD97j7gfU+52vfcKP/+Ec/s+bksMOjQxzaN8TNI4OaHCYboq63NbVSGO/fkarpp3KgAWxmUeDLwNfc/bdq2P5F4Ii7j6+1XXz3Ad995+8sv780Oezm0SEO7xtkdDipA2eyKal4hG3qelveYsXUtlyh+WFcawAHeRaEAZ8Cnl0tfM1sF/CKu7uZHQVCwMS6nxs4ODLI4X2lZQVNDpOtioRCbEvF6FXX2xaWlim2UQrjpfkUQYTxWoL83/RW4CeBb5nZk+XnfhkYBXD3e4H3AT9nZnkgDRzzGlr21+1M8Vv/9KYNF/TYmUlOnDzPxZk0u/sTHLtlhKP7hzf8eaSzpHoibOtV19uulsJ4uDdGJl+66KNVwjjwNeBGuPHgIX/g4Uc39DGPnZnknkeeJxIyeqIhFnNF8kXn7rcfUAh3qUgoxPa+mK5m7FCNDONalyB0flXZiZPniYSMRDSMUXobCRknTp4PujQJQKonwt6hhMK3g8Ujpa54ZDjJnqEEQ8lY02e16H9X2cWZNP0rTj/riYZ4eebasyikc6nr7U7xSJh4JMxQb4xsvlg6myKbJ5tv7DKF/peV7e5PMDGfIVFx14nFXJFd/YkAq5Jm6uuJsq03RkhrvV0tFgkRi8SaEsZagig7dssI+WJpMpNTepsvOsduGVn/g6WtRcMhdg8k2NEXV/jKVWKREEO9MfYOJdk7lGS4N1bXK2PVAZcd3T/M3RzgxMnzvDyTZpfOgugK/Ykow0l1vbK+pc54MBkjV1jqjAtkcoX1P3gVCuAKR/cPK3C7RDQcYkdfXDc6lU2JhkMMJmMMJiFXKLKQKTCXzW84jBXA0nUGElGGe2O6GlLqIhoOMZAMMZCMLodxrRTA0jXU9UqjLYVxrRTA0hXU9UorUgBLR1PXK61MASwdycwYSEQZSkbV9UrLUgBLx4lFSl1vPKKuV1qbAlg6hpkxmIgyqK5X2oQCWDqCul5pRwpgaWtmxlAyykBCXa+0HwWwtK14NMz2VExdr7QtBbC0naWudzAZC7oUkS1RAEtbiUfD7EjF6zqRSiQoCmBpC2bGcDLGQDIadCkidaMAlpbXEw2zXV2vdCAFsLSskBlDvTEGEup6pTMpgKUlJWKlrrfZN0kUaSYFsLSUkBnDqRj9Pep6pfMpgKVlqOuVbqMAlsCp65VupQCWQCVjEbanYkTU9UoXUgBLIMIhY7g3Rp+6XuligbYdZnabmX3XzF4ws49UeT1uZp8vv/5NM3tN86uUeuuNR9gzmFD4StcLLIDNLAz8LvBu4Abg/WZ2w4rNfgaYcvfXAb8N/NfmVin1FA4ZO/t7uK6/R0sOIgTbAR8FXnD3M+6eBU4Ad6zY5g7g0+XHXwDeYZo52JZ64xH2DiVJxbXqJbIkyADeA5yveH+s/FzVbdw9D1wBtlX7ZGZ23MxOmdmpyYnxBpQrm1HZ9YZD+tkpUinIAK723eib2Kb0pPt97n7E3Y8Mb9u+5eJk61LqekXWFOR3xhgwUvH+XuDCKtuMmVkEGAAmm1OebFYkFGJbKkavgldkTUF2wCeBA2b2WjOLAceAB1ds8yBwZ/nx+4BH3L1qByytIdUTYc9QQuErUoPAvkvcPW9mdwFfA8LA/e7+bTP7GHDK3R8EPgX8oZm9QKnzPRZUvbK2SCjE9r4YyZiCV6RW1okN5Y0HD/kDDz8adBldI9UTYXtvnJAOsoksqembQe2KbJq6XpGt0XeObEpfT5RtvTF1vSJboACWDYmGQ2xPxUnEdCt4ka1SAEvN+hNRhpPqekXqRQEs64qGQ+zoi9MTVdcrUk8KYFnTQCLKcG8MjeAQqT8FsFSlrlek8RTAcg11vSLNoQCWZep6RZpLASyYGQOJKEPJqLpekSZo6lu7AAAJlElEQVRSAHe5WKR0Xq+6XpHmUwB3KTNjMBFlUF2vSGAUwF0oFimt9cYj6npFgqQA7iJmxlAyykBCXa9IK1AAd4l4NMyOVJxYRHcjFmkVCuAOt9T1DiZjQZciIisogDuYul6R1qYA7kBmxnAyxkAyGnQpIrIGBXCH6YmG2a6uV6QtKIA7RMiMod4YAwl1vSLtQgHcARKxUtcbDavrFWknCuA2FjJjOBWjv0ddr0g7UgC3KXW9Iu1PAdxm1PWKdA4FcBtJxiJsT8WIqOsV6QgK4DYQDhnDvTH61PWKdJRAAtjMfgP4cSALfB/4oLtPV9nuRWAWKAB5dz/SzDpbQW88wrZedb0inSio7+qHgTe7+1uA7wG/tMa2b3P3g90WvuGQsbO/h+v6exS+Ih0qkO9sd/+6u+fL734D2BtEHa2qNx5h71CSVFwrRCKdrBVaq58GvrLKaw583cweN7Pja30SMztuZqfM7NTkxHjdi2yGyq43HNK8XpFO17AWy8z+HNhV5aVfcfcHytv8CpAHPrvKp3mru18ws53Aw2b2nLs/Wm1Dd78PuA/gxoOHfMt/gSZLxSNsS8UVvCJdpGEB7O4/ttbrZnYn8B7gHe5eNTDd/UL57SUz+xJwFKgawO0qEgqxLRWjV8sNIl0nkCUIM7sN+EXgve6+sMo2vWbWt/QYeBfwTPOqbLxUT4Q9QwmFr0iXCmoN+ONAH6VlhSfN7F4AM7vezB4qb3Md8Ddm9hTwGPBn7v7VYMqtr0goxK6BHnb2aa1XpJsF0nq5++tWef4CcHv58RngpmbW1Qypngjbe+OEFLwiXU+/+zZJJBRie1+MZEy7XERKlAZN0NcTZVtvTF2viFxFAdxA0XCI7ak4iVg46FJEpAUpgBukPxFlOKmuV0RWpwCus2g4xI6+OD1Rdb0isjYFcB0NJKIM98YwU9crIutTANeBul4R2QwF8Bap6xWRzVIAb5K6XhHZKgXwBpkZA4koQ8moul4R2RIF8AbEIqXzetX1ikg9KIBrYGYMJqIMqusVkTpSAK8jFimt9cYj6npFpL4UwKswM4aSUQYS6npFpDEUwFXEo2G2p2LqekWkoRTAFZa63sFkLOhSRKQLKIDL4tEwO1JxYpFWuFG0iHSDrg9gM2M4GWMgGQ26FBHpMl0dwD3RMNvV9YpIQLoygENmDPXGGEio6xWR4HRdACdipa43GlbXKyLB6poAVtcrIq2mKwJYXa+ItKKODuCQGcOpGP096npFpPV0bAAnYxG2p2JE1PWKSIvqyACOho1dAz1BlyEisqZA2kMz+zUze8nMniz/uX2V7W4zs++a2Qtm9pFaP39Iw3NEpA0E2QH/trv/99VeNLMw8LvAO4Ex4KSZPeju32lWgSIijdTKC6RHgRfc/Yy7Z4ETwB0B1yQiUjdBBvBdZva0md1vZkNVXt8DnK94f6z8nIhIR2hYAJvZn5vZM1X+3AF8AvgB4CBwEfjNap+iynO+xtc7bmanzOzU5cuX6/J3EBFppIatAbv7j9WynZn9PvDlKi+NASMV7+8FLqzx9e4D7gM4cuTIqkEtItIqgjoLYnfFuz8BPFNls5PAATN7rZnFgGPAg82oT0SkGYI6C+K/mdlBSksKLwIfAjCz64FPuvvt7p43s7uArwFh4H53/3ZA9YqI1J25d95v60eOHPFTp04FXYaIdK+aLkZo5dPQREQ6mgJYRCQgCmARkYAogEVEAqIAFhEJSEeeBWFml4GzW/gU24HxOpXTSO1SJ6jWRlGtjbHVWsfd/bb1NurIAN4qMzvl7keCrmM97VInqNZGUa2N0axatQQhIhIQBbCISEAUwNXdF3QBNWqXOkG1NopqbYym1Ko1YBGRgKgDFhEJiAJYRCQgXRvA5VshXTKzarOIsZL/Ub4j89NmdqjZNVbUsl6tt5rZlYq7TP9qs2ss1zFiZv/XzJ41s2+b2d1VtmmJ/Vpjra2yX3vM7DEze6pc63+osk3czD5f3q/fNLPXNL/Smmv9gJldrtiv/yKIWivqCZvZE2Z2zY0hGr5f3b0r/wA/ChwCnlnl9duBr1AaK/fDwDdbuNZbgS+3wD7dDRwqP+4Dvgfc0Ir7tcZaW2W/GpAqP44C3wR+eMU2/xK4t/z4GPD5Fq71A8DHg96vFfX8a+CPqv1bN3q/dm0H7O6PApNrbHIH8Bkv+QYwuOJOHk1TQ60twd0vuvvp8uNZ4FmuvZFqS+zXGmttCeV9NVd+N1r+s/Lo+R3Ap8uPvwC8w8xqmklbTzXW2jLMbC/wD4FPrrJJQ/dr1wZwDdrtrsw/Uv617ytm9qagiyn/qnYzpQ6oUsvt1zVqhRbZr+Vfk58ELgEPu/uq+9Xd88AVYFtzqyypoVaAf1JegvqCmY1Ueb1Zfgf4BaC4yusN3a8K4NVt6K7MATsN7HP3m4D/CfxJkMWYWQr4P8C/cveZlS9X+ZDA9us6tbbMfnX3grsfpHRz2qNm9uYVm7TMfq2h1j8FXuPubwH+nFc7zKYys/cAl9z98bU2q/Jc3farAnh1G7orc5DcfWbp1z53fwiImtn2IGoxsyilQPusu3+xyiYts1/Xq7WV9mtFTdPAXwIrB70s71cziwADBLxstVqt7j7h7pnyu78PHG5yaUveCrzXzF4ETgBvN7P/vWKbhu5XBfDqHgR+qnzU/oeBK+5+MeiiqjGzXUvrUmZ2lNK/60QAdRjwKeBZd/+tVTZrif1aS60ttF93mNlg+XEC+DHguRWbPQjcWX78PuARLx85aqZaal2x5v9eSuvvTefuv+Tue939NZQOsD3i7v98xWYN3a9B3RU5cGb2OUpHubeb2RjwUUoHDHD3e4GHKB2xfwFYAD4YTKU11fo+4OfMLA+kgWNBfPNR6ih+EvhWeQ0Q4JeB0YpaW2W/1lJrq+zX3cCnzSxM6YfAH7v7l83sY8Apd3+Q0g+TPzSzFyh1aMcCqLPWWj9sZu8F8uVaPxBQrVU1c7/qUmQRkYBoCUJEJCAKYBGRgCiARUQCogAWEQmIAlhEJCAKYBGRgCiARUQCogAWKTOzW8oDYnrMrLc8z3blHAORutGFGCIVzOzXgR4gAYy5+38OuCTpYApgkQpmFgNOAovA33P3QsAlSQfTEoTI1YaBFKW7ZPQEXIt0OHXAIhXM7EFKowlfC+x297sCLkk6WNdOQxNZycx+Csi7+x+Vp3n9PzN7u7s/EnRt0pnUAYuIBERrwCIiAVEAi4gERAEsIhIQBbCISEAUwCIiAVEAi4gERAEsIhKQ/w+2X1lUZVvYJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.lmplot(x = 'x', y = 'y', data = d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## formula notation with statsmodels\n",
    "use statsmodels.formula.api (often imported as smf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# data is in a dataframe\n",
    "model = smf.ols('y ~ x', data = d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# estimation of coefficients is not done until you call fit() on the model\n",
    "results = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.640\n",
      "Model:                            OLS   Adj. R-squared:                  0.460\n",
      "Method:                 Least Squares   F-statistic:                     3.556\n",
      "Date:                Fri, 09 Nov 2018   Prob (F-statistic):              0.200\n",
      "Time:                        15:32:30   Log-Likelihood:                -6.8513\n",
      "No. Observations:                   4   AIC:                             17.70\n",
      "Df Residuals:                       2   BIC:                             16.48\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      1.0000      2.324      0.430      0.709      -8.998      10.998\n",
      "x              1.6000      0.849      1.886      0.200      -2.051       5.251\n",
      "==============================================================================\n",
      "Omnibus:                          nan   Durbin-Watson:                   3.400\n",
      "Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.308\n",
      "Skew:                          -0.000   Prob(JB):                        0.857\n",
      "Kurtosis:                       1.640   Cond. No.                         7.47\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\miles\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\stattools.py:72: ValueWarning: omni_normtest is not valid with less than 8 observations; 4 samples were given.\n",
      "  \"samples were given.\" % int(n), ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "print(results.summary()) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the abline_plot function for plotting the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAH91JREFUeJzt3Xl81PW97/HX1yxkgwQIWza2LIhsgbgVQdmtS8V9r7X20tYe6kpvPfec28fpuT3tORAEtC7U/eipx1Zq19OERTZBlEVFxJkkJJAFCAFC9m3me/8giK2gJDOZ32Tm/Xw8eECGyfw+juTF8OM3H4y1FhER6f3Oc3oAERHxDwVdRCREKOgiIiFCQRcRCREKuohIiFDQRURCxFcG3RjzgjGm2hjz8eduG2CMWW2MKer8vn/PjikiIl/lXF6hvwRc+Xe3/RhYa63NAtZ2fiwiIg4y5/LGImPMCOBP1tpxnR+7gCustQeNMcOA9dbanJ4cVEREvlxkNz9viLX2IEBn1Aef7Y7GmAXAAoD4+PgpY8aM6eYhRUTC044dO2qstYO+6n7dDfo5s9auBFYC5OXl2e3bt/f0IUVEQooxZv+53K+7V7kc7jzVQuf31d18HBER8ZPuBv0PwD2dP74H+L1/xhERke46l8sWfw1sBXKMMRXGmPuAXwBzjDFFwJzOj0VExEFfeQ7dWnv7WX5qlp9nERERH+idoiIiIUJBFxEJEQq6iEiIUNBFREKEgi4iEiIUdBGREKGgi4iECAVdRCREKOgiIiFCQRcRCREKuohIiFDQRURChIIuIhIiFHQRkRChoIuIhIge/zdFRUSke9yH61la6D7n+yvoIiJBpqymkWVr3Pz+wyoSos890wq6iEiQqKpt5ol1RbyxvYKoCMN3p4/mu9NHMeCn5/b5CrqIiMOO1Lfy1PpiXtt2ACzcfclw7p8xmsF9Y7r0OAq6iIhDTjS18+zGEl58p4w2j5ebJqexcFYmaf3juvV4CrqISIA1tHbw4uZSVm7aR0NrB9dOSOHB2VmMGpTg0+Mq6CIiAdLS7uHVd/fz1PoSjjW2MWfsEB6Zm82Yof388vgKuohID2vr8PLG9nKeXFfMoboWpmUl88jcHCalJ/n1OAq6iEgP8Xgtb+2qZNlaN+XHmskb3p9lt03iklEDe+R4CrqIiJ95vZa/7jnE0tVuiqsbGJfaj5/eO44rsgdhjOmx4yroIiJ+Yq1lvesISwpd7KmqI3NwAk/fOZkrxw3t0ZCfoqCLiPjB1pKj5Be62L7/OBkD4lh6y0Sum5RKxHk9H/JTFHQRER98UF7LkgIXm4trGNovhp9dP45b8tKJigj87kMFXUSkG/YerCO/0M2avYcZEB/NP119PnddMpyYqAjHZlLQRUS6YN+RBh5fU8SfPqoioU8kj8zJ5t7LRpLQx/mcOj+BiEgvUHG8iRVri3hzZyV9Is/j/itGs2DaaBLjopwe7TMKuojIl6iub+GX64r59XvlYOCeS0dw/4zRJCf0cXq0L1DQRUTO4HhjG89u3MdLW0pp91huyUtn4cxMUpJinR7trBR0EZHPqW9p5/nNpTy/qZSGtg6um5jCg7OzGZEc7/RoX0lBFxEBmts8vLK1jGc2lHC8qZ0rLxjKw3OzyR7S1+nRzplPQTfGPAR8B7DAbuBea22LPwYTkfD01q5KFhe4qKptJiUplkXzcpifm9pjx2vr8PL6+wd4cl0x1fWtTM8exKNzs5mQ5t/FWYHQ7aAbY1KBHwJjrbXNxpg3gNuAl/w0m4iEmbd2VfLYqt00t3sAqKxt5rFVuwH8HvUOj5dVuypZvqaIytpmLhoxgCfvmMxFIwf49TiB5Ospl0gg1hjTDsQBVb6PJCLhanGB67OYn9Lc7mFxgctvQfd6LX/efZDH17jZd6SRCWmJ/NsN45melRyQfSs9qdtBt9ZWGmOWAAeAZqDQWlv49/czxiwAFgBkZGR093AiEgaqapu7dHtXWGtZ92k1Swrd7D1YR86Qvjx79xTmjh3S60N+ii+nXPoD1wEjgVrgN8aYu6y1r37+ftbalcBKgLy8POvDrCIS4lKSYqk8Q7x9vVRwS3ENiwtd7DpQy/CBcSy/bRLXTEgJ6OKsQPDllMtsoNRaewTAGLMK+Brw6pd+lojIWSyal/M359ABYqMiWDQvp1uPt2P/cfILXWwpOcqwxBh+fsN4bpqS5sjirEDwJegHgEuMMXGcPOUyC9jul6lEJCydOk/u61Uue6pOkF/oZt2n1SQnRPN/rxnLHRdnOLo4KxB8OYe+zRjzW2An0AHsovPUiohId83PTe32X4AWVzfw+Bo3f/7oIP1iIlk0L4d7p44gLjo83nLj03+ltfYnwE/8NIuISLeUH2ti+doiVu2sIDYqgoUzM/nOtFEkxgbP4qxACI/ftkQkJB2ua+HJdcW8/v4BjDF8e+pIvn/FaAYG4eKsQFDQRaTXOdbYxjMbSnh5Sxker+XWC9NZODOLoYkxTo/mKAVdRHqNupZ2nttUygubS2lq62B+bioPzsomY2Cc06MFBQVdRIJeU1sHL2/ZzzMbSjjR3M5V44fy8JxsMgf3nsVZgaCgi0jQau3w8OttB3jy7RJqGlqZkTOIR+bmMC410enRgpKCLiJBp8Pj5c2dFaxYW0xlbTOXjBrAs3dPZsrw3rs4KxAUdBEJGl6v5Y8fVbFsTRGlNY1MTE/i32+cwNTMgSGzb6UnKegi4jhrLas/OczS1W4+PVTPmKF9+dU385h9/mCFvAsUdBFxjLWWTUU15Be6+LDiBKOS41lxey7XjB/GeSG2OCsQFHQRccT2smMsLnCxrfQYqUmx/MeNE7hhciqRIbo4KxAUdBEJqI8rT7Ck0MV61xEG9e3Dv3zjAm67KJ0+kaG9OCsQFHQRCYiiw/UsXe3mfz4+RFJcFD/++hjuuXQEsdEKub8o6CLSow4cbWLZGje/+6CS+OhIHpiVxX3TRtIvJrwWZwWCgi4iPeLgiWaeWFfMG++XExlhWDBtFN+7fDT946OdHi1kKegi4lc1Da08vb6E/3x3P9Za7rg4g3+YkcngfuG9OCsQFHQR8YsTze38auM+XninlJZ2DzdOTuOHs7JIH6DFWYGioIuITxpbO3hpSxnPbiihrqWDayYM46E52YwelOD0aGFHQReRbmlp9/DatgM8vb6YmoY2Zp8/mIfn5DA2pZ/To4UtBV1EuqTd4+U32yt4Yl0RB0+0MDVzICvn5jA5o7/To4U9BV1EzonHa/nDh5UsW1PE/qNNTM5IIv+WiXxtdLLTo0knBV1EvpS1loI9h1i62o37cANjh/XjhW/lMSNHi7OCjYIuImdkrWWD+wj5hW52V55g9KB4fnnHZL4+bqgWZwUpBV1EvmDbvqMsKXTxftlx0vrHsuTmicyflKLFWUFOQReRz3xYXsuSQhebimoY3LcP/zp/HLfmpRMdqZD3Bgq6iOA6VE9+oYvCTw7TPy6K/3PV+dx96XBiorQ4qzdR0EXCWFlNI4+vcfOHD6tIiI7k4TnZfPuykST0URp6I/1fEwlDVbXNrFhbxG92VBAdcR7fu3w0350+iqQ4Lc7qzRR0kTBypL6VX75dzH9tOwDA3ZcM5/4ZoxncV4uzQoGCLhIGapvaeHbjPl56p4w2j5ebp6SxcFYWqUmxTo8mfqSgi4SwhtYOXthcyq827qOhrYNvTEzhwdnZjEyOd3o06QEKukgIamn38J9b9/P0hhKONbYxd+wQHp6bzZihWpwVyhR0kRDS1uHlv7eX8+S6Ig7XtTItK5lH5+YwMT3J6dEkABR0kRDg8Vp+t6uS5WvdlB9r5sIR/VlxWy4Xjxro9GgSQAq6SC/m9Vr+5+NDLF3touRII+NTE/nXe8dxefYgLc4KQz4F3RiTBDwHjAMs8G1r7VZ/DCah761dlSwucFFV20xKUiyL5uUwPzfV6bF6BWstb7uqyS90s6eqjqzBCTxz12TmXTBUIQ9jvr5CXw781Vp7kzEmGtA/Hijn5K1dlTy2ajfN7R4AKmubeWzVbgBF/StsKakhv9DNjv3HyRgQx+O3TuQbE1OJ0AbEsNftoBtj+gHTgW8BWGvbgDb/jCWhbnGB67OYn9Lc7mFxgUtBP4tdB46zpNDFO8VHGdovhn+7fjw356URpQ2I0smXV+ijgCPAi8aYicAO4AFrbePn72SMWQAsAMjIyPDhcBJKqmqbu3R7ONt7sI78Qhdr9lYzMD6af75mLHdenKHFWfIFvgQ9EpgMLLTWbjPGLAd+DPzz5+9krV0JrATIy8uzPhxPQkhKUiyVZ4h3it65+JmSIw08vtrNnz46SN+YSB6dm829U0cSr8VZcha+/MqoACqstds6P/4tJ4Mu8pUWzcv5m3PoALFRESyal+PgVMGh4ngTy9cU8ebOCmKiIviHGZn8r2mjSIyLcno0CXLdDrq19pAxptwYk2OtdQGzgE/8N5qEslPnyXWVy2nVdS08+XYxv37vAMYY7p06ku9fMZrkhD5Ojya9hK9/dlsIvNZ5hcs+4F7fR5JwMT83NawDfsrxxjae2VDCy1vL6PBYbrkwnYUzMxmWqNNP0jU+Bd1a+wGQ56dZRMJKfUs7z20q5fnNpTS2dXD9pFQemJ3F8IFanCXdo79dEQmw5jYPL28t45kNJdQ2tfP1cUN5eE42WUP6Oj2a9HIKukiAtHZ4eP29cp58u5gj9a1ckTOIR+bkMD4t0enRJEQo6CI9rMPjZdXOSpavLaKytpmLRg7gqTsnc+GIAU6PJiFGQRfpIV6v5U+7D7JstZt9NY1MTEvk5zeMZ1pWsvatSI9Q0EX8zFrL2r3VLCl08emhenKG9GXl3VOYM3aIQi49SkEX8aN3imtYXODig/JaRgyMY/ltk7h2QgrnaXGWBICCLuIHO/YfY0mBm637jpKSGMMvbhjPjVO0OEsCS0EX8cHHlSfIL3TxtusIyQl9+Mm1Y7nj4gz6RGpxlgSegi7SDcXV9Sxd7eYvuw+RGBvF/75yDPd8bThx0fqSEufoV59IF5Qfa2LZmiJ+t6uC2KgIfjgri+9MG0m/GC3OEucp6CLn4NCJFp5YV8R/v19OxHmG+y4byfcuH81ALc6SIKKgi3yJY41tPL2+mFe27sdrLbddlM7CmVkM6Rfj9GgiX6Cgi5xBXUs7z23cx/ObS2lu93B9bhoPzs4ifYD+2VwJXgq6yOc0tXXw0pYynt2wjxPN7Vw9fhgPzckic7AWZ0nwU9BFgJZ2D/+17QBPrS+mpqGNmWMG8/CcbMalanGW9B4KuoS1do+XN3dUsGJtEVUnWrh01ECevTuHKcP7Oz2aSJcp6BKWvF7LHz+q4vHVbsqONjEpPYnFN09kamay06OJdJuCLmHFWkvhJ4dZWujGdbieMUP78tw385h1/mAtzpJeT0GXsGCtZVNRDfmFLj6sOMGo5HieuD2Xq8cP0+IsCRkKuoS898uOsbjAxXulx0hNiuU/bprADbmpRGpxloQYBV1C1u6KEywpdLHBfYRBffvw0+su4NYL07U4S0KWgi4hx324nqWFbv665xBJcVE89vUxfPPSEcRGK+QS2hR0CRn7jzaybE0Rb31QSXx0JA/OzuK+y0bSV4uzJEwo6NLrHTzRzIq1xfxmezmREYYF00fxvemj6R8f7fRoIgGloEuvVdPQylNvl/Dqtv1Ya7nz4gx+MCOTwVqcJWFKQZde50RTOys3lfDiO2W0dni5cXIqP5yVRVp/Lc6S8KagS6/R2NrBi++UsnLjPupaOrh2YgoPzc5i1KAEp0cTCQoKugS9lnYPr767n6fXl3C0sY3Z5w/hkbnZnD+sn9OjiQQVBV2CVrvHyxvby3libTGH6lq4LDOZR+Zmk5uhxVkiZ6KgS9DxeC2//6CSZWuKOHCsiSnD+/P4rZO4dPRAp0cTCWoKugQNay1//fgQS1e7Kapu4IKUfrz4rQu5ImeQFmeJnAMFXRxnrWW9+wj5hS4+rqwjc3ACT905mSsvGKrFWSJdoKCLo97dd5QlBS627z9O+oBY8m+eyPzcVCIUcpEuU9DFER+W17Kk0MWmohqG9OvD/5s/jlvy0omO1AZEke5S0CWgPj1UR36hm9WfHGZAfDT/dPX53HXJcGKitDhLxFc+B90YEwFsByqttdf4PpKEotKaRh5f7eaPH1WR0CeSR+Zkc+9lI0noo9cUAG/tqmRxgYuq2mZSkmJZNC+H+bmpTo8lvYw/vpoeAPYCepeHfEFlbTMr1hTx250VREecx/cvH82C6aNIitPirFPe2lXJY6t209zuAU4+Z4+t2g2gqEuX+BR0Y0wacDXwM+Bhv0wkIaG6voWn3i7hv7YdAOCblw7n/isyGdS3j8OTBZ/FBa7PYn5Kc7uHxQUuBV26xNdX6MuAHwF9z3YHY8wCYAFARkaGj4eTYFfb1MYzG/bx8pYy2jxebslLY+HMLFKSYp0eLWhV1TZ36XaRs+l20I0x1wDV1todxpgrznY/a+1KYCVAXl6e7e7xJLg1tHbw/KZSntu0j4a2Dq6bmMKDs7MZkRzv9GhBLyUplsozxFu/CUpX+fIKfSrwDWPMVUAM0M8Y86q19i7/jCa9QUu7h1e2lvH0+hKON7Uz74IhPDwnh5yhZ/1Dm/ydRfNy/uYcOkBsVASL5uU4OJX0Rt0OurX2MeAxgM5X6I8q5uGjrcPLf79/gCfWFVNd38r07EE8OjebCWlJTo/W65w6T66rXMRXumZMuqTD4+V3uypZvraIiuPNXDRiAE/cnsvFo7Q4yxfzc1MVcPGZX4JurV0PrPfHY0lw8notf/n4IEtXu9l3pJHxqYn87PrxTM9K1uIskSChV+jypay1rPu0mvxCN58crCN7SALP3DWFeRcMUchFgoyCLme1pbiGJYUudh6oZfjAOJbdOolrJ6ZocZZIkFLQ5Qt2HjjOkgIXW0qOMiwxhp/fMJ6bpqQRFaHFWSLBTEGXz3xSVUd+oYu1n1aTnBDN/71mLHdcnKHFWSK9hIIulBxpYOlqN3/+6CD9YiJZNC+Hb31tBPFanCXSq+grNoyVH2ti+doiVu2sICYqgoUzM/nOtFEkxkY5PZqIdIOCHoaq61p4Yl0xr79/AGMM3546ku9fMZqBCVqcJdKbKehh5FhjG89sKOHlLWV4vJZbL0xn4cwshibGOD2aiPiBgh4G6lraeW5TKS9sLqWprYP5uak8OCubjIFxTo8mIn6koIew5jYPL28t45kNJdQ2tXPV+KE8NDubrCFanCUSihT0ENTa4eH198p58u1ijtS3MiNnEI/MzWFcaqLTo4lID1LQQ0iHx8ubOytYsbaYytpmLh45gKfvnEzeiAFOjyYiAaCghwCv1/LHj6pYtqaI0ppGJqYn8Ysbx3NZphZniYQTBb0Xs9ayZm81+YUuPj1Uz5ihffnVN/OYff5ghVwkDCnovZC1lneKj7K40MWH5bWMTI5nxe25XDN+GOdpcZZI2FLQe5ntZcdYXOBiW+kxUpNi+fcbx3Pj5DQitThLJOwp6L3Ex5UnyC908bbrCMkJffiXb1zAbRel0ydSi7NE5CQFPcgVV9ezdLWbv+w+RFJcFD/++hjuuXQEsdEKuYj8LQU9SB042sSytW7e2lVJbFQED8zK4r5pI+kXo8VZInJmCnqQOXSihRXrinjj/XIizjN8Z9oovnf5aAbERzs9mogEOQU9SBxtaOXp9SW88u5+rLXccXEGP5iRyZB+WpwlIudGQXfYieZ2ntu0jxc2l9Lc7uGGyWk8MCuL9AFanCUiXaOgO6SxtYOXtpTx7IYS6lo6uHrCMB6anU3m4ASnRxORXkpBD7CWdg+vbTvA0+uLqWloY9aYwTw8N5sLUrQ4S0R8o6AHSLvHy293VLBibREHT7QwNXMgK+fmMDmjv9OjiUiIUNB7mMdr+eOHVTy+xs3+o03kZiSRf/NEvpaZ7PRoIhJiFPQeYq2lYM9hlq524T7cwPnD+vH8PXnMHKPFWSLSMxR0P7PWsrGohvxCFx9VnGDUoHievCOXq8ZpcZaI9CwF3Y/eKz3GkgIX75UdI61/LItvmsD1ualanCUiAaGg+8FHFbUsKXSz0X2EwX378K/XXcCtF2YQHamQi0jgKOg+cB2qZ+lqFwV7DtM/Lop/vGoMd1+ixVki4gwFvRvKahpZtsbN7z+sIiE6kodmZ/Pty0bQV4uzRMRBCnoXVNU288S6It7YXkFUhOG700fzvctHkRSnxVki4jwF/RwcqW/lqfXFvPbuAQDuvmQ4988YzeC+WpwlIsGj20E3xqQDrwBDAS+w0lq73F+DBYMTTe08u7GEF98po83j5abJafxwdhapSbHdfsy3dlWyuMBFVW0zKUmxLJqXw/zcVD9OLSLhypdX6B3AI9bancaYvsAOY8xqa+0nfprNMQ2tHby4uZSVm/bR0NrBtRNSeGhONiOT43163Ld2VfLYqt00t3sAqKxt5rFVuwEUdRHxWbeDbq09CBzs/HG9MWYvkAr02qC3tHt49d39PLW+hGONbcwZO4RH5mYzZmg/vzz+4gLXZzE/pbndw+ICl4IuIj7zyzl0Y8wIIBfYdoafWwAsAMjIyPDH4fyurcPLG9vLeWJdEYfrWpmWlcwjc3OYlJ7k1+NU1TZ36XYRka7wOejGmATgTeBBa23d3/+8tXYlsBIgLy/P+no8f/J4LW/tqmTZWjflx5rJG96f5bflcsmogT1yvJSkWCrPEO8UH87Ji4ic4lPQjTFRnIz5a9baVf4Zqed5vZa/7jnE0tVuiqsbGJfaj5/eO44rsgf16OKsRfNy/uYcOkBsVASL5uX02DFFJHz4cpWLAZ4H9lprl/pvpJ5jrWW96whLCl3sqaojc3ACT985mSvHDQ3IBsRT58l1lYuI9ARfXqFPBe4GdhtjPui87R+ttX/xfSz/21pylCWFLnbsP07GgDiW3jKR6yalEhHgDYjzc1MVcBHpEb5c5bIZCPp9sB+U17KkwMXm4hqG9ovhZ9eP45a8dKK0AVFEQkzIvlN078E68gvdrNl7mIHx0fzT1edz1yXDiYnS4iwRCU0hF/R9Rxp4fE0Rf/qoioQ+kTw6N5t7p44kvk/I/aeKiPyNkKlcxfEmVqwt4s2dlfSJPI/7rxjNgmmjSYzTBkQRCQ+9PujV9S38cl0xv36vHAzcc+kI7p8xmuSEPk6PJiISUL026Mcb23hmYwkvbymjw2O5OS+dhTMz9SYdEQlbvS7o9S3tPL+5lOc3ldLQ1sH8Sak8ODuL4QN9W5wlItLb9ZqgN7d5eGVrGc9sKOF4UztXXjCUh+dmkz2kr9OjiYgEhaAPeluHl9ffP8CT64qprm/l8uxBPDo3h/FpiU6PJiISVII26B0eL6t2VbJ8TRGVtc1cNHIAT94xmYtGDnB6NBGRoBR0Qfd6LX/efZDH17jZd6SRCWmJ/PyG8UzLSg7IvhURkd4qaIJurWXt3mryV7vZe7COnCF9efbuKcwdO0QhFxE5B0ER9HeKa1hS6GLXgVpGDIxj+W2TuGZCSsAXZ4mI9GaOBn3H/uMsKXCxdd9RUhJj+MUN47lxSpoWZ4mIdIMjQd9TdYL8QjfrPq0mOSGan1w7ltsvytDiLBERHwQ06K0dXn7w2k7+vPsgibFR/OjKHL71tRHERQfFmR8RkV4toCV1H67H66rmhzMzuW/aKBJjtThLRMRfAhr05IQ+bPzRDAZqcZaIiN8F9G8fhyXGKOYiIj1El5OIiIQIBV1EJEQo6CIiIUJBFxEJEQq6iEiIUNBFREKEgi4iEiIUdBGREKGgi4iECAVdRCREKOgiIiFCQRcRCREKuohIiFDQRURChIIuIhIiFHQRkRChoIuIhAifgm6MudIY4zLGFBtjfuyvoUREpOu6HXRjTATwS+DrwFjgdmPMWH8NJiIiXePLK/SLgGJr7T5rbRvwOnCdf8YSEZGuivThc1OB8s99XAFc/Pd3MsYsABZ0fthqjPnYh2OGkmSgxukhgoSei9P0XJym5+K0nHO5ky9BN2e4zX7hBmtXAisBjDHbrbV5PhwzZOi5OE3PxWl6Lk7Tc3GaMWb7udzPl1MuFUD65z5OA6p8eDwREfGBL0F/H8gyxow0xkQDtwF/8M9YIiLSVd0+5WKt7TDG/ANQAEQAL1hr93zFp63s7vFCkJ6L0/RcnKbn4jQ9F6ed03NhrP3CaW8REemF9E5REZEQoaCLiISIgARdKwJOM8a8YIypDvfr8Y0x6caYt40xe40xe4wxDzg9k1OMMTHGmPeMMR92Phf/4vRMTjPGRBhjdhlj/uT0LE4yxpQZY3YbYz44l0sXe/wceueKADcwh5OXOr4P3G6t/aRHDxykjDHTgQbgFWvtOKfncYoxZhgwzFq70xjTF9gBzA/HXxfGGAPEW2sbjDFRwGbgAWvtuw6P5hhjzMNAHtDPWnuN0/M4xRhTBuRZa8/pDVaBeIWuFQGfY63dCBxzeg6nWWsPWmt3dv64HtjLyXcfhx17UkPnh1Gd38L2agVjTBpwNfCc07P0NoEI+plWBITlF66cmTFmBJALbHN2Eud0nmL4AKgGVltrw/a5AJYBPwK8Tg8SBCxQaIzZ0blG5UsFIujntCJAwpMxJgF4E3jQWlvn9DxOsdZ6rLWTOPmO64uMMWF5Os4Ycw1Qba3d4fQsQWKqtXYyJ7fa/qDzlO1ZBSLoWhEgZ9R5vvhN4DVr7Sqn5wkG1tpaYD1wpcOjOGUq8I3Oc8evAzONMa86O5JzrLVVnd9XA7/j5CnsswpE0LUiQL6g8y8Cnwf2WmuXOj2Pk4wxg4wxSZ0/jgVmA586O5UzrLWPWWvTrLUjONmKddbauxweyxHGmPjOCwYwxsQDc4EvvTqux4Nure0ATq0I2Au8cQ4rAkKWMebXwFYgxxhTYYy5z+mZHDIVuJuTr8A+6Px2ldNDOWQY8LYx5iNOvgBaba0N68v1BIAhwGZjzIfAe8CfrbV//bJP0Fv/RURChN4pKiISIhR0EZEQoaCLiIQIBV1EJEQo6CIiIUJBFxEJEQq6iEiI+P9k1ngXAmQo/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sm.graphics.abline_plot(model_results = results)\n",
    "plt.scatter(d.x, d.y)\n",
    "\n",
    "plt.xlim(0,5)\n",
    "plt.ylim(0,10)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Generating an anova table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           df  sum_sq  mean_sq         F  PR(>F)\n",
      "x         1.0    12.8     12.8  3.555556     0.2\n",
      "Residual  2.0     7.2      3.6       NaN     NaN\n"
     ]
    }
   ],
   "source": [
    "print(sm.stats.anova_lm(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Making predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    4.2\n",
       "dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.predict({'x' : 2})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## numpy array notation\n",
    "similar to sklearn's notation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4]\n"
     ]
    }
   ],
   "source": [
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X = sm.add_constant(x)  \n",
    "# need to add a constant for the intercept term.\n",
    "# because we are using the numpy notation, we use sm rather than smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1.]\n",
      " [1. 2.]\n",
      " [1. 3.]\n",
      " [1. 4.]]\n"
     ]
    }
   ],
   "source": [
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$y_i = \\beta_0 + \\beta_1 x_i + \\epsilon_i$$\n",
    "\n",
    "$$\\mathbf{\\hat{Y}} = \\boldsymbol{\\beta} \\mathbf{X}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# OLS is capitalized in the numpy notation\n",
    "model2 = sm.OLS(y, X)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "results2 = model2.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.640\n",
      "Model:                            OLS   Adj. R-squared:                  0.460\n",
      "Method:                 Least Squares   F-statistic:                     3.556\n",
      "Date:                Fri, 09 Nov 2018   Prob (F-statistic):              0.200\n",
      "Time:                        15:41:47   Log-Likelihood:                -6.8513\n",
      "No. Observations:                   4   AIC:                             17.70\n",
      "Df Residuals:                       2   BIC:                             16.48\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.0000      2.324      0.430      0.709      -8.998      10.998\n",
      "x1             1.6000      0.849      1.886      0.200      -2.051       5.251\n",
      "==============================================================================\n",
      "Omnibus:                          nan   Durbin-Watson:                   3.400\n",
      "Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.308\n",
      "Skew:                          -0.000   Prob(JB):                        0.857\n",
      "Kurtosis:                       1.640   Cond. No.                         7.47\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\miles\\Anaconda3\\lib\\site-packages\\statsmodels\\stats\\stattools.py:72: ValueWarning: omni_normtest is not valid with less than 8 observations; 4 samples were given.\n",
      "  \"samples were given.\" % int(n), ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "print(results2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## OLS solution:\n",
    "\n",
    "$$(X^TX)^{-1}X^TY$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1.],\n",
       "       [1., 2.],\n",
       "       [1., 3.],\n",
       "       [1., 4.]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1. , 1.6])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(X.T @ X) @ (X.T @ y)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}