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": 122,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-13T15:18:36.966142Z",
     "start_time": "2020-05-13T15:18:36.963232Z"
    },
    "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)\n"
   ]
  },
  {
   "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)\n"
   ]
  },
  {
   "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()\n",
    "\n"
   ]
  },
  {
   "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": [
    "\n",
    "print(results.summary()) \n",
    "\n"
   ]
  },
  {
   "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))\n"
   ]
  },
  {
   "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)  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "results2 = model2.fit()\n"
   ]
  },
  {
   "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())\n"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Interaction of Categorical Factors\n",
    "\n",
    "https://www.statsmodels.org/dev/examples/notebooks/generated/categorical_interaction_plot.html\n",
    "\n",
    "In this example, we will visualize the interaction between categorical factors. First, we will create some categorical data. Then, we will plot it using the interaction_plot function, which internally re-codes the x-factor categories to integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Consumption  Gender  Income\n",
      "0            51    Male      80\n",
      "1            52    Male      80\n",
      "2            53    Male      90\n",
      "3            54    Male      90\n",
      "4            56    Male     100\n",
      "5            57    Male     100\n",
      "6            55  Female      80\n",
      "7            56  Female      80\n",
      "8            58  Female      90\n",
      "9            59  Female      90\n",
      "10           62  Female     100\n",
      "11           63  Female     100\n"
     ]
    }
   ],
   "source": [
    "# https://stackoverflow.com/questions/55663474/interaction-plot-from-statsmodels-formula-api-using-python\n",
    "\n",
    "import pandas as pd\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "Consumption = [51, 52, 53, 54, 56, 57, 55, 56, 58, 59, 62, 63]\n",
    "Gender  = [\"Male\", \"Male\", \"Male\", \"Male\", \"Male\", \"Male\", \"Female\",\n",
    "           \"Female\", \"Female\", \"Female\", \"Female\", \"Female\"]\n",
    "Income = [80, 80, 90, 90, 100, 100, 80, 80, 90, 90, 100, 100]\n",
    "\n",
    "df = pd.DataFrame( {\"Consumption\": Consumption, \"Gender\": Gender, \"Income\": Income})\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/lib/python3.9/site-packages/scipy/stats/stats.py:1541: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=12\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>Consumption</td>   <th>  R-squared:         </th> <td>   0.976</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.967</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   108.4</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Thu, 13 Jan 2022</td> <th>  Prob (F-statistic):</th> <td>8.11e-07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>14:41:06</td>     <th>  Log-Likelihood:    </th> <td> -9.9135</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    12</td>      <th>  AIC:               </th> <td>   27.83</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>     8</td>      <th>  BIC:               </th> <td>   29.77</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</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>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>             <td>   27.3333</td> <td>    3.059</td> <td>    8.935</td> <td> 0.000</td> <td>   20.279</td> <td>   34.387</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Gender[T.Male]</th>        <td>    4.0000</td> <td>    4.326</td> <td>    0.925</td> <td> 0.382</td> <td>   -5.976</td> <td>   13.976</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Income</th>                <td>    0.3500</td> <td>    0.034</td> <td>   10.340</td> <td> 0.000</td> <td>    0.272</td> <td>    0.428</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Gender[T.Male]:Income</th> <td>   -0.1000</td> <td>    0.048</td> <td>   -2.089</td> <td> 0.070</td> <td>   -0.210</td> <td>    0.010</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 2.522</td> <th>  Durbin-Watson:     </th> <td>   3.273</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.283</td> <th>  Jarque-Bera (JB):  </th> <td>   0.970</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.055</td> <th>  Prob(JB):          </th> <td>   0.616</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 1.612</td> <th>  Cond. No.          </th> <td>2.62e+03</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 2.62e+03. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            Consumption   R-squared:                       0.976\n",
       "Model:                            OLS   Adj. R-squared:                  0.967\n",
       "Method:                 Least Squares   F-statistic:                     108.4\n",
       "Date:                Thu, 13 Jan 2022   Prob (F-statistic):           8.11e-07\n",
       "Time:                        14:41:06   Log-Likelihood:                -9.9135\n",
       "No. Observations:                  12   AIC:                             27.83\n",
       "Df Residuals:                       8   BIC:                             29.77\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "=========================================================================================\n",
       "                            coef    std err          t      P>|t|      [0.025      0.975]\n",
       "-----------------------------------------------------------------------------------------\n",
       "Intercept                27.3333      3.059      8.935      0.000      20.279      34.387\n",
       "Gender[T.Male]            4.0000      4.326      0.925      0.382      -5.976      13.976\n",
       "Income                    0.3500      0.034     10.340      0.000       0.272       0.428\n",
       "Gender[T.Male]:Income    -0.1000      0.048     -2.089      0.070      -0.210       0.010\n",
       "==============================================================================\n",
       "Omnibus:                        2.522   Durbin-Watson:                   3.273\n",
       "Prob(Omnibus):                  0.283   Jarque-Bera (JB):                0.970\n",
       "Skew:                          -0.055   Prob(JB):                        0.616\n",
       "Kurtosis:                       1.612   Cond. No.                     2.62e+03\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 2.62e+03. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Reg = ols(formula = \"Consumption ~ Gender + Income + Gender*Income\", data = df)\n",
    "Fit = Reg.fit()\n",
    "Fit.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEJCAYAAAByupuRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6dElEQVR4nO3dd3hUZdr48e/0ySSEkAIYYEEENBTpNZGWARELvhBFFAmw6i6wqCC236uyuqiwoqigYkHeBFZEl2rBklAFkSZNkCZRIBhKGkmmz/n9cWAQksgAyUySuT/X5XWRmczMnceTc+c8z3PuW6MoioIQQoiQpA12AEIIIYJHkoAQQoQwSQJCCBHCJAkIIUQIkyQghBAhTJKAEEKEMH2wA7gS2dnZAfus2NhYTp06FbDPqy5kXEqTMSmbjEvZAj0u8fHxZT4uVwJCCBHCJAkIIUQIkyQghBAhrFquCVxMURTsdjterxeNRlOh752Tk4PD4ajQ9wwURVHQarWYzeYKHxchRM1QI5KA3W7HYDCg11f8j6PX69HpdBX+voHidrux2+2EhYUFOxQhRBVUI6aDvF5vpSSAmkCv1+P1eoMdhhCiiqoRSUCmOv6cjI8Qojw1IgkIIYS4MiGbBIzr11O3SxeM69dXyPs1atSIfv36+f47cuRIhbxvWbp27Upubm6lvb8QInSE5ES6cf16olNT0dpsRKemkpuWhjMx8are02w28+2331ZQhEIIERghdyXwxwQA+BJBRV0R/NHOnTsZMmQIAwYM4N577yUnJweAlJQUJk+ezODBg+nVqxfbt2/ngQceIDExkWnTpvleP3r0aAYMGECfPn2YP39+mZ+xaNEibr31Vvr168cTTzyBx+Op8J9DCFFz1bgrgcjnnsOwZ0+Zz2kKCjD8/DOai3bLaG02Yu65B9cNN6DUrn3hazQanAkJFL7wwp9+rt1up1+/fgD85S9/Yfbs2TzzzDPMnTuXmJgYli1bxrRp03jttdcAMBqNLF68mA8++IDRo0ezYsUKoqKi6NGjBw8++CDR0dG8+uqr1KlTB5vNxq233srAgQOJjo72feaBAwdYvnw5S5cuxWAw8PTTT7N48WLuuuuuyx43IURoqnFJ4M8YDh4slQDO0Xi9GA4exNmx4xW998XTQT///DP79u3jnnvuAdRtrHXr1vU9379/fwBuuOEGWrRoQb169QBo3Lgx2dnZREdH8+GHH7JixQpALZp3+PDhC5LAd999x65duxg4cCCgJqLY2Ngril8IETjG9esxTJqEcfr0q56Kvlo1Lgn82V/sF08F/ZE3LKzMtQG9Xo/b7b7sOBRFoUWLFnz22Wdlx2I0AqDVan3/Pve1x+Nhw4YNrFu3js8++4ywsDBSUlJK3bmsKAp33XUXTz/99GXHJ4QIjnPnIU0FrklejZBaE3AmJpKblob3ortny0sAV+O6664jNzeXLVu2AOByudi3b5/frz9z5gy1a9cmLCyMgwcPsm3btlLfk5SUxOeff+4rR5uXl8fRo0cr5gcQQlS4QK5J+iukkgCUTgSVkQBA/Uv/3Xff5aWXXsJqtdK/f39fQvBH79698Xg8WK1W/v3vf9OhQ4dS39OiRQueeOIJhg0bhtVqZdiwYb7FZyFE1VLeTESwE4FGURQlKJ98FS5uKlNSUoLFYrms9zCuX0/UhAnkz5jxpwngSqeDqpIrGZ9LkUYhpcmYlE3GRVW3c2f0f9IQy92gASc2baq0zy+vqUyNWxPwlzMxsVIHXAghAHS//opl/nw0hYXlfo83LIz8GTMCGNV5IZsEhBCi0ng8mFauJDw9HdOqVaDVYr/5ZpwdO1Jr+vQLpoQqa0raXwFLAsXFxcyePZsjR46g0WgYM2YMP/zwA1u3bkWv11OvXj3Gjh1LeHh4oEISQogKpT11CsvHH2OZNw/90aN46tWjaMIEiu+9F+811wDgatPGtzYQ7AQAAVwTmDVrFgkJCSQnJ+N2u3E4HBw8eJDWrVuj0+l8d8QOHz78ku9VEWsC/pI1gbLJPG9pMiZlq/HjoigYt2zBkpZG2Oefo3G5cCQmUjxiBPabbwaDodRLjOvXEzNpEqcDeJ9AUNcESkpK2Lt3L+PGjVM/VK9Hr9fTtm1b3/e0aNGCjRs3BiIcIYS4apriYsIWLyY8LQ3D3r14a9WieMQISu6/H3fz5n/6WmdiIq4DB3BWgeQYkCRw4sQJIiMjefvtt/n1119p2rQpI0eOxGw2+75n5cqV9OjRo8zXZ2RkkJGRAcDUqVNL3RWbk5NTqU1lqnvDGpPJVOF3Euv1erk7+SIyJmWraeOi2bMH7bvvov3Pf9CcOYO3bVvc77yDd+hQjOHhGC/9FkDVGZeAnN08Hg+HDx9m9OjRNG/enLlz57J06VJfSYXFixej0+m46aabyny91WrFarX6vr740tLhcFRaC0h/p4MaNGjA4MGDmTlzJqC2dWzfvj3t27cnPT293Ndt2LCB2bNn/+n3XC2Hw1Hhl+M1/hL/CsiYlK1GjIvTiXnFCsLnzcP4/fcoRiO222+nODUVV4cOoNGAzab+56dAj0t500EBuVksJiaGmJgYmp+9ROrWrRuHDx8GYPXq1WzdupWHH344oB2wcnK0DBkSw4kTFTMEFouFffv2YTt7EKxdu5b69etXyHsLIYJDe+wYtf79b+p16UL02LHojh2j8H//l5wtW8h/801cHTuqCaAaC0gSiIqKIiYmxregu2vXLho2bMj27dtZtmwZTz75JCaTKRCh+Lz+egQ//GDk9dcjKuw9+/TpQ2ZmJgBLly7lzjvv9D33448/cscdd9C/f3/uuOMODh48WOr1JSUlTJw4kYEDB9K/f3++/vrrCotNCOEnrxfT2rXU+etfqdetGxFvvomrbVtOz5vHifXrKRo7Fm9MTLCjrDABm+wePXo0b775Jm63m7p16zJ27Fiefvpp3G43//rXvwBo3rw5Dz300FV9znPPRbJnT+nV+D9yOmHbNiOKomHevHB27zZgLGciT6PRkJDg5IUXyr/R45xBgwYxY8YMrFYre/fu5Z577uGHH34AoFmzZixevBi9Xs/atWuZNm0a77///gWvf+ONN0hMTOS1116joKCAW2+9lZtuuqnSdj4JIc7T5OVh+eQTwufNQ3/4MJ6YGIrGjqVk+HA8jRoFO7xKE7Ak0KRJE6ZOnXrBY+fmzwPt6NHz6weKon7dtOnVN2Np2bIlR48eZdmyZfTt2/eC5woLC3n00Uc5fPgwGo0Gl8tV6vVr167l22+/Zfbs2YA6l3/s2DHfNJoQouIZduwgPC2NsGXL0NjtODp35sxjj2EbOBACPEMRDNV720sZLvUXe06Olu7d66Eo6jyeomgoKNDx9tunqFu3dK+By71PoH///rzwwgv897//JS8vz/f4K6+8Qo8ePZgzZw5HjhwhJSWl1GsVReG9996jWbNmfn+eEOLyaWw2zMuXE56WhnHHDrwWCyUpKRSPGIG7VatghxdQIVdF9PXXI7j49jivlwpbGxg6dCgTJkwgISHhgsfPnDnjWyj+5JNPynxtr169mDt3Lufu39u9e3eFxCSEUOl++YXI55+nXqdO1Jk4EY3NRv6LL5KzdSsF06aFXAKAGnglcClbtxpxOi9czXc6NWzZ4u/u3j8XHx/PAw88UOrxMWPG8Oijj/Lee++RWM4dgo8++iiTJ0/GarWiKAoNGzas1K2jQoQEtxtzRgaW9HTMa9ag6PXYb7mF4tRUnN26VfvdPVcrZEtJ+0vKRpStRuz9rmAyJmUL1rhoT5zA8tFHhM+fj+74cTzXXEPx8OGUDBuG92w712CqKvcJhNyVgBCiBlMUjBs3Ep6ejvnLL9G43dh79qRgyhTsVitU87v/K4OMiBCi2tOcOUPYokVqHZ/9+/HWrk3x6NEU338/nqZNgx1elVYjkkA1nNEKKBkfUVPp9+xRt3cuXoy2pARn27bkvfYa9jvuQLmol7goW41IAlqtFrfbXe0LvVUGt9uNVhtym8BETeZwEPbll1jS0jBt3oxiNmMbNIjiESNwtWsX7OiqnRpx1jSbzdjtdhwOR4XXHzKZTDgcjgp9z0BRFAWtVntBtVYhqivdkSNY5s/HsmAButOncTdpQsFzz1Fy990odeoEO7xqq0YkAY1GQ1glXfrJjg8hgsjrxbR6NeFpaZgyM0Gjwd6/PyWpqTiSkkCucq9ajUgCQoiaRZubS9jChWodn19/xRMXR9HDD1N83314GzQIdng1iiQBIUTVoCgYtm1TF3o//xyNw4Gje3cKn3oK+4ABlFvlUVwVSQJCiKDSlJQQtmQJlvR0jLt3442IoGTYMLWOz/XXBzu8Gk+SgBAiKPQHD2JJT8fy6adoCwtxJSSQ//LL2AYPRomouD4f4s9JEhBCBI7LhfnrrwlPT8e0fj2KwYDtttsoSU3F2alTyNfxCQZJAkKISqc9fpzwjz7C8p//oMvJwd2wIYVPPaXW8akCzdZDmSQBIUTlUBSM69ej//hj6i1fDl4vjj59yJ82DUffvqDTXfo9RKWTJCCEqFCaggIsn36KJT0dw6FDKNHRFD/0EMXDh+Np0iTY4YmLSBIQQlQIw65dWNLT1To+djvODh3Ie+MNwkeOpLCoKNjhiXJIEhBCXDm7nbDPPlPbNP74I16zGdvgwZSMGIGrTRsAws1mkCRQZUkSEEJcNl1WFuHz5xP28cfo8vJwXXcdBS+8QElKCkrt2sEOT1wGSQJCCP94PJgyM9XtnatXg1aL/eab1TaNiYmyvbOakiQghPhT2lOnsCxYgGX+fPRHj+KpX58zEyeq2zuvuSbY4YmrJElACFGaomDcvBlLWhphX3yBxuXCkZhI4XPPYe/fHwyGYEcoKkjAkkBxcTGzZ8/myJEjaDQaxowZQ3x8PDNmzODkyZPExcUxYcIEIuR2cSGCRlNURNjixYSnp2PYuxdvZCTFI0ZQMmIE7mbNgh2eqAQBSwJz586lXbt2PPbYY7jdbhwOB0uWLKFNmzbceeedLF26lKVLlzJ8+PBAhSSEOEv/88+Ep6cTtmgR2qIinK1bk//KK9juvBPFYgl2eKISBaQjQ0lJCXv37qVv374A6PV6wsPD2bx5M7169QKgV69ebN68ORDhCCEAnE7My5YRM2QIdZOTsXz8MfYBAzj52Wec+uorSu69VxJACAjIlcCJEyeIjIzk7bff5tdff6Vp06aMHDmSgoIC6pxtC1enTh0KCwvLfH1GRgYZGRkATJ06ldgA1hrR6/UB/bzqQsaltGozJr/9hm7OHLRz56LJyUFp0gT3Sy/hTU1FHxtLRW/wrDbjEmBVZVwCkgQ8Hg+HDx9m9OjRNG/enLlz57J06VK/X2+1WrFarb6vA9nuUdpLlk3GpbQqPSZeL6Z167CkpWH+9ltQFBxWK8UjRuDo3ft8m8ZKiL9Kj0sQBXpc4uPjy3w8IEkgJiaGmJgYmjdvDkC3bt1YunQptWvXJi8vjzp16pCXl0dkZGQgwhEiZGjy8rCca9OYlYUnJoaisWMpGT4cT6NGwQ5PVAEBSQJRUVHExMSQnZ1NfHw8u3btomHDhjRs2JA1a9Zw5513smbNGjp37hyIcISo8Qzbt6ttGpcvR2O34+jShTOTJmEbOBBMpmCHJ6qQgO0OGj16NG+++SZut5u6desyduxYFEVhxowZrFy5ktjYWCZOnBiocISocTQ2G+ZlywhPT8e4Ywdei4WSu+5S2zS2bBns8EQVpVEURQl2EJcrOzs7YJ8l85llk3EpLVhjojt0iPB587B88gnaggJcLVpQnJqKbcgQlFq1Ah7PxeRYKVtIrQkIISqY2405I0Nd6F27FkWvxz5woFrHp2tXqeMj/CZJQIhqRJuTg+Wjjwj/z3/QHT+O55prKHz8cUruvRdv3brBDk9UQ5IEhKjqFAXjxo2Ep6VhXrECjduNvVcvCl58EXtyMujl11hcOTl6hKiiNIWFhC1apNbx2b8fb1QUxaNHU3z//XiaNg12eKKGkCQgRBWj/+kntY7P4sVoS0pwtmtH3muvYbvjDggLC3Z4ooaRJCBEVeBwEPbFF2qbxi1bUMxmbIMGUTxiBK527YIdnajBJAkIEUS6I0ewzJ+P5aOP0OXm4m7ShILJkym56y6Us3W1hKhMkgSECDSPB9Pq1YSnpWFauRI0Guz9+1OSmoojKel8HR8hAkCSgBABos3NxfLxx1jmzUP/22944uIoeuQRiu+9F2+DBsEOT4QoSQJCVADj+vUYJk3COH262nT9HEXBsHWrWsfn88/ROJ04unen8OmnsQ8YAEZj8IIWgstIAm63m9WrV5OVlYXdbr/guX/84x8VHpgQ1YVx/XqiU1PR2GxEp6aSm5aGq317wpYsITwtDcNPP+GNiKD4vvsouf9+3NdfH+yQhfDxOwnMmjWLX3/9lY4dO1K7dkW3nRCiejqXALQ2GwBam42YYcNQjEa0NhuuhATyp07FNngwSnh4kKMVojS/k8COHTuYNWsW4XIgCwGUTgDnaDwecDop+Ne/KB41Sur4iCrN720IsbGxuFyuyoxFiGqlzsMPl0oA52g8HsJnz5YEIKo8v68EevbsySuvvMItt9xCVFTUBc+1bt26ouMSompSFIzffUd4ejraEydQgLJO896wMPJnzAh0dEJcNr+TwFdffQXAggULLnhco9Ewa9asio1KiCpGk5+P5dNPscybh+HQITx16lD097/jSkgg6oknLrgi8IaFkZuWduEuISGqKL+TwFtvvVWZcQhRJRl27cKSlkbYkiVo7XacHTuS98Yb2G67DcxmAHLr1fOtDUgCENXNZd0n4PF42LdvH7m5ucTExNCiRQt0Ol1lxSZEcNhshH32mdqm8ccf8YaFYRsyRG3TWMbUpzMxkdy0NGImTSL34vsEhKji/E4Cx44dY9q0aTidTmJiYjh9+jQGg4Enn3yShg0bVmaMQgSELitLbdP48cdo8/NxNWtGwQsvUJKSgnKJbdHOxERcBw7glDaKoprxOwl88MEHWK1Wbr/9djRndzwsX76cOXPmMHny5EoLUIhK5fFgyswkPD0d86pVKDod9gED1DaNPXrI7h5R4/mdBLKysnj22Wd9CQDg1ltvZcmSJZUSmBCVSXvyJJYFC7DMn4/+2DE89etT+NhjlAwbhveaa4IdnhAB43cSiI6OZs+ePRdsB927dy91pNytqC4UBeOmTVjS0wn74gs0LheOpCQK//lP7P36gcEQ7AiFCDi/k8CwYcOYNm0aHTt2JDY2llOnTrFt2zbGjx9fmfEJcdU0RUXn2zT+/DPeyEiKR4ygZMQI3M2aBTs8IYLK7yTQqVMnpk2bxvfff09eXh6NGjXi7rvvJj4+vjLjE+KK6X/+WW3T+N//oi0uxtm6NfnTp2MbNAjFYgl2eEJUCZe1RTQ+Pp4hQ4Zc0QeNGzcOs9mMVqtFp9MxdepUsrKyeP/993E6neh0Oh544AGayV9m4mo4nZhXrFAbtvzwA4rJhO322ylOTcXVvr0s9ApxkT9NAu+++y5/+9vfAJg5c+YFi8J/5G8p6cmTJxMZGen7ev78+aSkpNC+fXu2bdvG/Pnz+ec//+ln6EKcpzt27HybxlOncDduTMGzz1Jy990o0dHBDk+IKutPk0DdunV9/65fv36Ff7hGo8F29nb7kpISWWQWl8frxbR2LZa0NMwZGaAoOKxW8lNTcfTqJW0ahfCDRlEUxZ9vzM/PL1U47s8ev9i4ceOIiIgAoF+/flitVo4ePcqLL74IgNfrZcqUKcTFxZV6bUZGBhkZGQBMnToVp9PpT8gVQq/X43a7A/Z51UVQx+X0abTp6ejeew/NL7+gxMXhHT0az1//Co0bBycm5Fgpj4xL2QI9LsZyutj5nQRSU1NJS0sr9fioUaOYO3fuJV+fm5tLdHQ0BQUFTJkyhVGjRrFx40ZatmxJt27d2LBhA5mZmTz77LOXfK/s7Gx/Qq4Q53ZCiQsFfFwUBcP27WqbxuXL0TgcOLp0oSQ1Fdstt4DJFLhYyiHHStlkXMoW6HEpbxOP3wvDZeWKkpIStH5eckefnZetXbs2nTt35uDBg6xZs4ZRo0YB0L17d959911/wxEhQmOzYV62jPC0NIw7d+IND6dk6FC1jk9CQrDDE6Lau2QSGDNmDABOp9P373OKiopI9KNYlt1uR1EUwsLCsNvt7Ny5k5SUFN8NaK1atWL37t2Vsu4gqifdoUOEp6dj+fRTtAUFuK6/nvwXX8Q2ZAhKrVrBDk+IGuOSSWD8+PEoisLLL79c6sawqKgov+4TKCgoYPr06YBaiTQpKYl27dphNpuZO3cuXq8Xg8Hg24kkQpTbjfnbb9XtnevWoRgM2AYOpGTECJxdu8r2TiEqgd9rAg6HA1MVmHcFWROoCipyXLQ5OVg++ojw+fPR/f477vh4SoYPV+v4/GGHWlUnx0rZZFzKVu3WBHQ6HQsXLmT9+vXk5eVRp04devToweDBg8tddRaiXIqC8fvvCU9Lw/zVV2jcbuy9e1Pw0kvYk5NBf1n3MQohrpDfv2nvv/8+2dnZjBo1iri4OE6ePMnSpUv54IMPGDt2bGXGKGoQTWEhlv/+F0t6OoYDB/BGRVH8179SPHw4nqZNgx2eECHH7ySwefNmZs6cSXh4OAANGzakefPmUkBO+EW/e7dax2fJErQlJTjbtyfvtdew3XEHhIUFOzwhQpbfSSAqKgqHw+FLAqDuGJK7fEW57HbCvvhCbdO4ZQuK2UzJnXdSMmIErrZtgx2dEILLSAI9e/bkpZdeYsCAAb72kl9//TU9e/Zk9+7dvu9rXUYPVhFadL/9ptbxWbAAXW4u7muvpeCf/6TkrrtQ/Li7XAgROH4ngW+//RagVCexb7/91vecRqNh1qxZFRieqDY8HkyrVhGeno5p5UrQaLDffDPFI0bgTEqSOj5CVFF+J4G33nqrMuMQ1ZT29GksH3+MZd489EeO4Klbl6JHHqH43nvxNmgQ7PCEEJcg+/DE5VMUNN9/T9SbbxL2+edonE4c3btT+L//i33AAGnTKEQ1clmN5tPS0sjKysJut1/w3IIFCyo8MFH1aIqLCVuyhPC0NAx79qCrVYvi++5T2zS2aBHs8IQQV8DvJPDGG2/QtWtXRo0aJTeHhRj9gQNYztXxOXMGV0IC7rfe4mS/fih/2C0mhKh+/E4C+fn5DB06tNzuYqKGcbkwf/WVWsfn++9RjEZst91G8YgRuDp1IjYuDkVKAQhR7fmdBHr16sV3333HTTfdVJnxiCDTZmcT/tFHapvGnBzcjRpR+P/+HyX33IM3JibY4QkhKpjfSeDOO+/kmWeeYcmSJdSuXfuC5yZPnlzhgYkA8noxfvcd4enpmL/5BrxeHH36kP/vf+Po0wd0umBHKISoJH4ngddee426devSpUsXWROoITT5+Vg+/ZTw9HT0v/yCJzqaor//nZLhw/H85S/BDk8IEQCXtTvoww8/RC/VHas9w86dWNLSCFu6FK3djrNjR/LefBPbrbeC2Rzs8IQQAeT3GT0hIYGjR4/SpEmTSgxHVBqbjbDPPlPr+Pz4I96wMGxDhqhtGqXUhxAhy+8kEBcXx5QpU+jSpUupNYGhQ4dWeGCiYugOHyZ83jwsCxeizc/H1awZBf/6FyUpKSiRkcEOTwgRZH4nAafTSYcOHXC73Zw+fboyYxJXy+PBlJmpNmxZvRpFr8c+YADFqak4u3eXNo1CCB+/k4A0jqn6tCdPYvnoIyz/+Q/6Y8fw1K9P4aRJapvG+vWDHZ4Q4qycHC1Dh+qZOVNL3breoMbidxLIyckp97l69epVSDDiCigKxk2b1IXeL79E43LhuOkmCp9/Hnu/ftKmUYgq6PXXI9iwQcPrr0fw0kuFQY3F7zPEww8/XO5zCxcurJBghP80Z84QtmgR4fPmYfj5Z7yRkRSnplJ8//14mjULdnhCiHL8/LOejz4Kx+vVsHBhOI8+WhTUqwG/k8DFJ/r8/Hw+/fRTEhISKjwoUT793r1qm8ZFi9AWF+Ns04a8V1/FPmgQirRpFKLKURT46Sc9GRlmMjPNbNtmANR1Oa+XoF8NXPFcQVRUFCNHjuSRRx4hKSmpImMSF3M6CfvySyzp6Zh++AHFZMJ2xx0Up6biatdOFnqFqGJsNg3r1hl9J/7ff1fvum/Z0olOBx6P+n1OZ/CvBq5qwjg7OxuHw1FRsYiL6I4dwzJvntqm8dQp3E2aUPDss5TcfTdKdHSwwxNC/MHRozoyMkxkZppZv96Ew6EhPNxLr14OkpPt9O3rYMaMCA4eNPiSAAT/asDvJPDcc89dUEHU4XBw5MgRUlJS/Hr9uHHjMJvNaLVadDodU6dOBWDFihV89dVX6HQ6OnTowPDhwy/zR6hhvF5Ma9ZgSU/HnJEBgN1qpSQ1FUfPntKmUYgqwu2GbduMvhP/zz+rzZSaNHEzfHgxVquDrl0dmEznX7N1qxGn88Ird6dTw5YtwSvF43cS6Nu37wVfm81mGjduzDXXXOP3h02ePJnIP9ygtHv3brZs2cL06dMxGAwUFBT4/V41jSY3F8snnxA+bx76rCw8sbEUjRun1vFp2DDY4QkhgLw8DatXm8nMNLFqlZn8fC16vUKXLk6efbYAq9XOddd5yp2h/eab8+XXY2NjOVUFyrH7nQR69+5d4R/+zTffMGjQIAxn2xFefCdyjacoGLZvJzwtjbDly9E4HDi6dqXwiSew33ILSKE+IYJKUWD//nOLuiY2bzbi9WqIjvZgtdqxWu306uUgMlIJdqhXTKMoil/Rf/7557Ru3ZomTZqwf/9+ZsyYgU6nY/z48Vx//fWXfP24ceOIiIgAoF+/flitVh5//HE6d+7M9u3bMRgM3H///TQrY3tjRkYGGWenRqZOnYrT6bycn/Gq6PV63G53xb5pSQnahQvRvvsu2h9/RImIwHvvvXj/9jeUalLHp1LGpZqTMSlbdRsXux3WrNHw5ZdaVqzQ8uuv6p/1bdt6ueUWhYEDvXTqpFx1hfVAj0t51Z/9TgJjxozh1VdfxWKx8Pzzz9OpUyfCwsLIyMjgpZdeuuTrc3NziY6OpqCggClTpjBq1CjmzJlDq1atGDVqFIcOHWLGjBnMmjXrkt3LsrOz/Qm5QlTkJZvu4EG1js+nn6ItKMB1ww0UjxiBbcgQlLMJsrqoKpeyVYmMSdmqw7gcP64lM1P9a3/dOhM2mxaz2ctNNzmxWu307WsnPr5id+8Eelzi4+PLfNzv6aCSkhIsFgs2m42srCyeffZZtFot6enpfr0++uxultq1a9O5c2cOHjxIdHQ0Xbt2RaPR0KxZM7RaLWfOnLlg3aDac7sxf/ON2qbxu+9QDAZsAwdSkpqKs0sX2d4pRBB4vbB9u4GMDDMZGWZ++kmdkm7QwM3dd9uwWu107+4gFG698TsJxMTEsG/fPo4cOUJCQgJarZaSkhK0fuxWsdvtKIpCWFgYdrudnTt3kpKSgtlsZvfu3bRq1Yrs7Gzcbje1atW6qh+oqtD+/juWBQsInz8f3e+/446Pp/DJJ9U6PnFxwQ5PiJBTWKhhzRp1J8/KlSZOn9ah1Sp06uTk6acLsVrtXH+9O+T+LvM7CQwfPpzXXnsNvV7PY489BsC2bdvKnMO/WEFBAdOnTwfA4/GQlJREu3btcLvdvP322zz22GPo9XrGjRtXvRvZKwrGDRvU6p1ff43G7cbeuzf5L7+MIzlZ2jQKEWCHDul8f+1v2mTE7dYQFeWld287VquDXr3sREdX30XdiuD3mkBZzi1qBLrbWFVbE9AUFmL573+xpKdjOHAAb1QUJffcQ/Hw4XiuvTZAkQZWdZjnDTQZk7IFclycTti40Uhmpnriz8pSz03XX+/CarWTnOygY0dnlairWO3WBEBdF8jOzsZut1/weOtqsqOloul371br+CxejNZmw9m+PXkzZmC7/XZCYjJRiCrg5EktK1eayMgws3atiaIiLSaTQo8eDh58sIjkZAeNGnku/UYhyu8ksHr1aubMmYPZbL5gq5FGo2HWrFmVElywGdevxzBpEsbp03EmJqoP2u2EffEF4WlpGLduxWs2Y/uf/6FkxAhcN94Y3ICFCAFeL+zebSAzUz3xb9+uno/q1/cwaJC6qJuU5MRiCe1pHn/5nQQWLFjAxIkTad++fWXGU2UY168nOjUVjc1GdGoq+a+8gmHPHrWOT14e7qZNKfjnPym56y6UqKhghytEjVZcrGHdOhOZmerCbk6ODo1GoV07F48/ri7qtmoVeou6FcHvJOD1emnbtm1lxlJlnEsAWpsNAK3NRp1//AO0WrVN44gROJOSZHunEJXo1191Z+f2TXz/vQmnU0OtWl569nSc3bvvIDY2uF25agK/k8CgQYNYtGgRQ4YM8WtbaHV1cQI4RwMoRiPFI0eenxoSQlQYlwu2bDH6TvwHDqh795s2dTNyZDHJyXa6dHFKNZUK5ncS+OKLL8jPz2f58uW+8g/nvPPOOxUeWLBETZhQKgGco7HbiZowgRObNgU4KiFqptxcLatWqXP7a9aYKCjQYjAodOvm5L77CkhOttO0qSzqVia/k8D48eMrM44qI3/GjDKvBAC8YWHkz5gRhKiEqBkUBfbu1fu2cG7bZsDr1RAb62HAADvJyXZ69nRQq5Ys6gaK30mgZcuWlRlHleFMTCQ3La1UIvCGhZGbliZTQUJcJpsNMjJMvkqc2dnqaadNGyePPFKE1Wrnxhtd0iojSPxOAm63m8WLF7N27Vry8vKoU6cOPXv2ZPDgwQG/WayyXZwIJAEIcXmOHdP6/tpfv96A3R6DxaIu6k6YUETfvnbq15dF3arA77P3/PnzOXToEA8++CBxcXGcPHmSRYsWUVJSwsiRIysxxOA4lwhiJk0i94/3CQghSvF4YNs2g6+n7t696qLuX/7iZvRoL4mJ+XTr5sBsDnKgohS/k8DGjRt55ZVXfAXe4uPjufbaa3n88cdrZBIANRG4DhzAKaUAhCiloEDD6tXqNM+qVSby8nTodGqXrWeeKcBqddCsmZu4uFhOnZJe5FWV30ngKkoMCSFqAEWBgwf1vp66mzYZ8XjUgmx9+57vshUVJeeK6sTvJNC9e3emTZtGSkqKr/DRokWL6NatW2XGJ4QIIocDNm40+RZ2f/tNPWUkJLgYM0Zd1O3QwSUFcquxyyolvWjRIubMmUNeXh7R0dEkJiYyZMiQyoxPCBFgv/+uZeVKdSfP2rUmSkq0mM0KiYkO/v73IqxWBw0ayN79muKSSeDnn39my5YtDB8+nKFDhzJ06FDfc/Pnz+eXX36hRYsWlRqkEKLyeL2wc6fBt4Vz5071ltz4eDdDhqgF2RITnYSFyTRPTXTJJLBkyRJuvvnmMp9r3bo1ixcv5qmnnqrwwIQQlaeoSMPateoUz8qVJk6eVAuydezo4skn1YJsCQlSkC0UXDIJZGVl0a5duzKfa9OmTY0qGSFETXb4sM63hXPjRiMul4bISC+9e6sF2fr0cRAdLXv3Q80lk4DNZsPtdl/QQ+Acj8eDrZw6O0KI4HI6YdOm8122fvlF/XVv3tzFAw+oBdk6dXJiMAQ5UBFUl0wCDRo0YMeOHXTu3LnUczt27KBBgwaVEpgQ4vKdOqV22crMVAuynTmjxWhU6N7dwahR6om/cWNZ1BXnXTIJ3Hrrrbz33nt4vV46d+6MVqvF6/WyefNm5syZw4gRIwIRpxCiDIoCP/2k9zVT377dgKJoqFvXw2232bBaHdx0k4PwcFnUFWW7ZBJISkoiPz+ft956C5fLRWRkJIWFhRiNRu666y6SkpICEacQ4qySEg3ffWf0ze///ru6Sb9dOyePPXaG5GQHrVtLQTbhH7/uE7jtttvo27cv+/fvp6ioiIiICFq0aIHFYqns+IQQwJEjOl9P3Q0bTDgcGsLDvfTqdX5Rt25dWdQVl8/vm8UsFku5u4SEEBXL7YatW42+E/++ferqbZMmbu6/X53b79rVickU5EBFtVezakALUY3l5WlYvVptrbh6tZn8fC16vVqQ7bnnCrBa7Vx3nSzqiooVsCQwbtw4zGYzWq0WnU7H1KlTfc8tX76c+fPn88EHHxAZGRmokIQIKkWBffv0vp66W7YY8Xo1xMR46NdP7bLVq5eDyEhZ1BWVJ6BXApMnTy51kj916hS7du0iNjY2kKEIERR2O2zYcL7L1tGj6q9gq1Yuxo8vIjnZTrt2UpBNBE7Qp4PS0tK47777eOWVV4IdihBXLCdHy9ChembO1JZaoD1+/HyXre++M2KzaTGb1S5b48erXbbi42VRVwRHQJPAiy++CEC/fv2wWq1s2bKF6OhomjRp8qevy8jIICMjA4CpU6cG9KpBr9fLVUoZZFwu9PzzOjZs0DB7dhwzZnjYskXDl19qWbFCw44d6l7Nxo0VUlO93HKLi169FMLCdIDl7H81lxwrZasq46JRAtQtJjc3l+joaAoKCpgyZQqjRo1i/vz5PPPMM1gsFsaNG8fLL7/s15pAdnZ2ACJWneudIC4k43JeTo6W7t3r4XBo0OkUIiO95OXp0GoVOnd2kpysbuNs0SI0C7LJsVK2QI9LfHx8mY8H7EogOjoagNq1a9O5c2f27NnDiRMnePzxxwE4ffo0Tz75JC+//DJRUVGBCkuIK6IocOiQWpDtvffCcZztnujxQFSUlylTCunVy06dOrKoK6q2gCQBu92OoiiEhYVht9vZuXMnKSkpfPDBB77vuZwrASGCwemEjRvP36mblaX++mg0CnDuT3wNx4/r6dHDIQlAVAsBSQIFBQVMnz4dUCuPJiUlyY1nolo4ceLCgmzFxVpMJrXL1oMPFrF1q5HPPw/D6Tz/Gq8XXn89gpdeKgxe4EL4KSBJoF69epfc/fPWW28FIhQh/pTXC7t3G3zN1LdvV0uo16/v4c471S5bSUlOLBb1r/yPPrLgdF440e90atiypXTpdSGqoqBvERUi2IqLNaxbZ/Kd+E+cULtstW/v4vHH1S5brVqVvaj7zTfnF/ZkAVRUR5IEREjKytKRmanesPX99yacTg21al1YkC02Vvbui5pPkoAICS4XbN5s9JVoOHhQLch23XUuRo4sxmq106WLdNkSoUeSgKixcnPPL+quXm2isFCLwaDQrZuT++8vIDnZzrXXSkE2EdokCYgaQ1Fg797zXba2bVO7bMXFebjlFjtWq52ePR1ERMjWTSHOkSQgqjWbTe2yda42z/HjauW1G290MmGCWpDtxhuly5YQ5ZEkIKqdY8d0ZGSc77Jlt2uwWNSCbI895qBvXzv16smirhD+kCQgqjyPB7ZtM/q2cO7dq67eNm7s5r77iklOdtCtm0O6bAlxBSQJiCopP1/DmjXqX/urVpnIy9Oh06ldtp59tgCr1cF114VmQTYhKpIkAVElKAocOKD39dTdvNmIx6OhTh0Pffs6SE6207u3g9q1ZVFXiIokSUAEjd0OGzeev1P3t9/UwzEhwcXYseqibocO0mVLiMokSUAE1O+/a1m5Ur1ha906EyUlWsxmtSDbmDFFJCc7aNBA9u4LESiSBESl8nphxw6D707dXbvUwmrx8W5SUmwkJ9tJTHQSFibTPEIEgyQBUeHOnNGwdq06t79ypYlTp9QuWx07OnnqKbUg2w03yKKuEFWBJAFRIX75RedrtvLDD0ZcLg21a3vp3dtOcrKDPn0cREfL3n0hqhpJAuKKOJ3w3Xfnu2z98ot6KLVo4eKBB9SCbJ06OdHLESZElSa/osJvp06pBdkyMsysXWvgzJlYjEaFHj0cjBpVTHKyncaNZVFXiOpEkoAol6LATz/p+fZb89kuW2pBtnr1PNx1l5ekpAKSkhyEh8uirhDVlSQBcYGSEs0F0zy//65u0m/f3sljj53BanXQqpWLunVjOXXKHuRohRBXS5KA4LffdGRmqjdsbdhgwuHQEBGhFmSzWu307esgLk4WdYWoiSQJhCC3G7ZsMfpKNOzfrxZka9LEzf33q4u6Xbs6MUqvdCFqPEkCISI3V8Pq1WpP3dWrzeTna9HrFbp2dTJsmNpl67rrZFFXiFAjSaCGUhTYt+9cly0TW7ca8Xo1xMR46NfvfJetyEhZ1BUilEkSqEFsNtiwweQr0XDsmPq/t3VrJ+PHF2G12mnXTrpsCSHOC1gSGDduHGazGa1Wi06nY+rUqcybN4+tW7ei1+upV68eY8eOJTw8PFAh1QjZ2VoyM9WdPOvWGbHbtYSFebnpJgePPFJE3752rrlGFnWFEGUL6JXA5MmTiYyM9H194403cu+996LT6Zg/fz5Llixh+PDhgQyp2vF44McfDb6eunv2qIu6jRq5GTashORkB927OzCbgxyoEKJaCOp0UNu2bX3/btGiBRs3bgxiNFVXQcGFXbZyc9UuW506Ofnf/y0kOdlOixZSkE0IcfkCmgRefPFFAPr164fVar3guZUrV9KjR48yX5eRkUFGRgYAU6dOJTY2tnID/QO9Xh/Qz4Nzi7qwYoWWFSu0rF+vwe3WEB2tcPPNXm65xU2/fl6iozWA+ex/gRWMcanqZEzKJuNStqoyLhpFUQKyPSQ3N5fo6GgKCgqYMmUKo0aNomXLlgAsXryYQ4cOMWnSJDR+/DmbnZ1d2eH6xMbGcurUqUr/HIcDfvjhfJetrCw1P99wgwurVa3E2aFD1SnIFqhxqU5kTMom41K2QI9LfHx8mY8H7JQSHR0NQO3atencuTMHDx6kZcuWrF69mq1bt/Lcc8/5lQBqkhMn/liQzURxsRaTSe2y9eCDRVitDho2lL37QojKE5AkYLfbURSFsLAw7HY7O3fuJCUlhe3bt7Ns2TKef/55TCZTIEIJKq8Xdu0y+O7U3bFDvSW3fn0P//M/apetpCQnFovs3RdCBEZAkkBBQQHTp08HwOPxkJSURLt27Rg/fjxut5t//etfADRv3pyHHnooECEFTFGRhnXr1GmelSvNnDihQ6NR6NDBxRNPqF22WraURV0hRHAEJAnUq1ePV155pdTjM2fODMTHB1xW1rkuWya+/96Ey6WhVi0vvXs7SE5WC7LFxMjefSFE8FWRZcbqzeWCzZuNvhP/wYPq3v3rrnMxerRakK1zZycGQ5ADFUKIi0gSuEK5uecXddesMVFYqMVgUOje3cH995eQnGzn2mtlUVcIUbVJEvgTOTlahg7VM3Omlrg4L3v26H3NVrZtU7tsxcV5GDjQhtXq4KabHEREyKKuEKL6kCTwJ6ZPr8X69RqGDo3hzBktx4+rXbbatnUyYYJakK1NGynIJoSoviQJlOPYMS0ffWQBNOzfr6dvXzuTJtnp08dBvXqyqCuEqBkkCZRj1qwIdDq1YJvBAI0aebjnHluwwxJCiAolExllyMnRsnBhOB6Punnf5dKwcGE4J07IcAkhahY5q5Xh9dcjuLiikterPi6EEDWJJIEybN1qxOm88BZep1PDli3SeV0IUbPImkAZvvnmfGU/qYAohKjJ5EpACCFCmCQBIYQIYZIEhBAihEkSEEKIECZJQAghQpgkASGECGGSBIQQIoRJEhBCiBAmSUAIIUKYJAEhhAhhkgSEECKESRIQQogQJklACCFCmEZRLq6cL4QQIlTIlcAlPPXUU8EOoUqScSlNxqRsMi5lqyrjIklACCFCmCQBIYQIYZIELsFqtQY7hCpJxqU0GZOyybiUraqMiywMCyFECJMrASGECGGSBIQQIoTpgx1AVfL555+zcuVKNBoNjRo1YuzYsTidTmbMmMHJkyeJi4tjwoQJREREBDvUgCprXJYuXUpmZiaRkZEADBs2jA4dOgQ50sD68ssvyczMRFEUkpOTufXWWykqKgrp46WsMfnkk09C7lh5++232bZtG7Vr1+bVV18F+NNjY8mSJaxcuRKtVsuoUaNo165d4IJVhKIoinL69Gll7NixisPhUBRFUV599VVl1apVyrx585QlS5YoiqIoS5YsUebNmxfEKAOvvHFZuHChsmzZsiBHFzy//vqrMnHiRMVutytut1t54YUXlOzs7JA+Xsobk1A8Vn766Sfl0KFDysSJE32PlXdsHDlyRJk0aZLidDqVnJwc5R//+Ifi8XgCFqtMB/2B1+vF6XTi8XhwOp3UqVOHzZs306tXLwB69erF5s2bgxxl4JU1LqHu2LFjNG/eHJPJhE6nIyEhgU2bNoX08VLemISili1blroCLO/Y2Lx5Mz169MBgMFC3bl3q16/PwYMHAxarJIGzoqOjuf322xkzZgwPPfQQFouFtm3bUlBQ4Dvp1alTh8LCwiBHGljljQvA119/zaRJk3j77bcpKioKcqSB1ahRI/bu3cuZM2dwOBz8+OOPnD59OqSPl/LGBEL7WDmnvGMjNzeXmJgY3/dFR0eTm5sbsLhkTeCsoqIiNm/ezFtvvYXFYuG1115j7dq1wQ4r6Mobl/79+5OSkgLAwoULSU9PZ+zYsUGONnAaNmzIoEGDmDJlCmazmcaNG6PVhvbfVOWNSagfK5eiBHmXfmgftX+wa9cu6tatS2RkJHq9nq5du7J//35q165NXl4eAHl5eb7FrVBR3rhERUWh1WrRarUkJydz6NChYIcacH379mXatGk8//zzREREcM0114T88VLWmMixoirv2IiJifFdMYF6ZRAdHR2wuCQJnBUbG8uBAwdwOBwoisKuXbto0KABnTp1Ys2aNQCsWbOGzp07BznSwCpvXM4dzACbNm2iUaNGQYwyOAoKCgA4deoUmzZtIjExMeSPl7LGRI4VVXnHRqdOndiwYQMul4sTJ05w/PhxmjVrFrC45I7hP/jkk0/YsGEDOp2OJk2a8Pe//x273c6MGTM4deoUsbGxTJw4MaS2/EHZ4zJ79myysrLQaDTExcXx0EMPhdyC8XPPPceZM2fQ6/WMGDGCNm3acObMmZA+Xsoak5kzZ4bcsfL666+zZ88ezpw5Q+3atbn77rvp3LlzucfG4sWLWbVqFVqtlpEjR9K+ffuAxSpJQAghQphMBwkhRAiTJCCEECFMkoAQQoQwSQJCCBHCJAkIIUQIkyQghBAhTJKACEnjxo1j586dwQ5DiKCTJCCEECFMCsiJkLZ69WoyMzNp3rw5q1atwmKx8MADD/ju2CwqKiI9PZ0dO3bgdDpJSEjgiSeeACAjI4Nly5ZRVFTEDTfcwIMPPuir+XL33Xfz17/+lS+++IL8/HwGDhxI7969mTlzJkePHqVt27Y8/PDD6PXqr+DWrVv5+OOPOXnyJA0bNuTBBx+kcePGwRkUEVLkSkCEvIMHDxIfH8+cOXMYNGgQs2fP9lV2nDlzJg6Hg1dffZX333+f2267DYDdu3ezYMECJkyYwHvvvUdcXBxvvPHGBe+7fft2pk6dyosvvsjy5ct57733ePjhh3nnnXc4cuQI3333HQC//PIL77zzDg899BAffvghVquVf//737hcrsAOhAhJkgREyIuNjcVqtaLVaunVqxd5eXkUFBSQl5fH9u3befDBB4mIiECv19OyZUsA1q1bR58+fWjatCkGg4F7772X/fv3c+LECd/7Dho0CIvFQqNGjWjUqBE33ngj9erVw2Kx0L59e7KysgDIzMzEarXSvHlztFotvXv3Rq/Xc+DAgWAMhwgxMh0kQl5UVJTv3yaTCQC73U5RURERERFlFoDLy8vj2muv9X1tNpuJiIggNzeXunXrlnpfo9FY6uv8/HxArbi5Zs0avvrqK9/zbrc7oI1FROiSJCBEOWJiYigqKqK4uJjw8PALnqtTpw6nTp3yfX0uaVxJHfiYmBgGDx7M4MGDrzpmIS6XTAcJUY46derQrl07PvjgA4qKinC73ezZsweApKQkVq1aRVZWFi6XiwULFtCsWTPfVcDlSE5O5ttvv+XAgQMoioLdbmfbtm3YbLaK/pGEKEWuBIT4E+PHj+f//u//mDBhAm63m1atWtGyZUvatGnD0KFDefXVVykqKuL666/n0UcfvaLPuO666/jb3/7Ghx9+yPHjxzEajdxwww0kJCRU7A8jRBmkn4AQQoQwmQ4SQogQJklACCFCmCQBIYQIYZIEhBAihEkSEEKIECZJQAghQpgkASGECGGSBIQQIoT9fyuSUiWkpTLGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.graphics.factorplots import interaction_plot\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "fig = interaction_plot(\n",
    "   x        = Income,\n",
    "   trace    = Gender,\n",
    "   response = Fit.fittedvalues,\n",
    "   colors   = ['red','blue'],\n",
    "   markers  = ['D','^'])\n",
    "plt.xlabel('Income')\n",
    "plt.ylabel('Consumption')\n",
    "plt.legend().set_title(None)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  },
  "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
}