Python_ANOVA/Python_ANOVA_Factorial_Using_Statsmodels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two Way ANOVA using Statsmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                      sum_sq    df          F        PR(>F)\n",
      "C(supp)           205.350000   1.0  15.571979  2.311828e-04\n",
      "C(dose)          2426.434333   2.0  91.999965  4.046291e-18\n",
      "C(supp):C(dose)   108.319000   2.0   4.106991  2.186027e-02\n",
      "Residual          712.106000  54.0        NaN           NaN\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = 'https://vincentarelbundock.github.io/Rdatasets/csv/datasets/ToothGrowth.csv'\n",
    "\n",
    "df = pd.read_csv(data, index_col=0)\n",
    "\n",
    "# ANOVA\n",
    "formula = 'len~C(supp)+C(dose)+C(supp):C(dose)'\n",
    "model = ols(formula, df).fit()\n",
    "aov_table = anova_lm(model, typ=2)\n",
    "\n",
    "print(aov_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting Effect Sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                    sum_sq    df       F  PR(>F)  eta_sq  omega_sq\n",
      "C(supp)           205.3500   1.0  15.572  0.0002  0.0595    0.0555\n",
      "C(dose)          2426.4343   2.0  92.000  0.0000  0.7029    0.6926\n",
      "C(supp):C(dose)   108.3190   2.0   4.107  0.0219  0.0314    0.0236\n",
      "Residual          712.1060  54.0     NaN     NaN     NaN       NaN\n"
     ]
    }
   ],
   "source": [
    "def eta_squared(aov):\n",
    "    aov['eta_sq'] = 'NaN'\n",
    "    aov['eta_sq'] = aov[:-1]['sum_sq']/sum(aov['sum_sq'])\n",
    "    return aov\n",
    "\n",
    "def omega_squared(aov):\n",
    "    mse = aov['sum_sq'][-1]/aov['df'][-1]\n",
    "    aov['omega_sq'] = 'NaN'\n",
    "    aov['omega_sq'] = (aov[:-1]['sum_sq']-(aov[:-1]['df']*mse))/(sum(aov['sum_sq'])+mse)\n",
    "    return aov\n",
    "\n",
    "\n",
    "eta_squared(aov_table)\n",
    "omega_squared(aov_table)\n",
    "\n",
    "print(aov_table.round(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the Residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuclHX5//HXxUEBz4KYp10MEH9SQrqaiH1TMUuoDLXS1vKMmsdSRFpDUVABxW8CZeABlYVviVimWCpqpWa5GEcVRd1FhJBDArpy2r1+f3zugWGZ2Z1hZ3ZmZ9/Px2Mec8993zv3NaPstffndJm7IyIikqpWuQ5ARESaFyUOERFJixKHiIikRYlDRETSosQhIiJpUeIQEZG0KHGIiEhalDhERCQtShwiIpKWNrkOIBs6derkXbp0yXUYIiLNxuzZs1e5+36pnFuQiaNLly5UVFTkOgwRkWbDzKpSPVdNVSIikhYlDhERSYsSh4iIpEWJQ0RE0qLEISIiaVHiEBFpAcrLoUsXaNUqPJeX7/x7FeRwXBER2aa8HAYNgurq8LqqKrwGKC1N//10xyEiUuDKyrYljZjq6rB/ZyhxiIgUuCVL0tvfECUOEZECV1SU3v6GKHGIiBS4kSOhQ4ft93XoEPbvDCUOEZE8ksnRTzGlpTBxIhQXg1l4njhx5zrGQaOqRETyRqZHP8UrLW38e8TojkNEJE9kevRTtuQscZhZDzObE/dYZ2bX1jnnRDNbG3fOsFzFKyKSbZke/ZQtOWuqcvdFQG8AM2sNfAQ8keDUv7v7t5syNhGRXCgqCs1Tifbnk3xpquoHvOfuKRcSEREpNJke/ZQt+ZI4zgamJTnWx8zmmtkzZtYz2RuY2SAzqzCzipUrV2YnShGRLMr06KdsMXfPbQBmuwDLgJ7uvqLOsT2BWnf/1Mz6A79y9+4NvWdJSYmrdKyISJyaGnj4Ydh9d/jBD3Y4bGaz3b0klbfKhzuO04A36iYNAHdf5+6fRtszgbZm1qmpAxQRabbc4ZlnoHdvuOgi+N3vGv2W+ZA4ziFJM5WZfcHMLNo+lhDv6iaMTUSk+XrjDTjlFOjfHzZsgMceg+nTG/22OZ0AaGYdgG8Al8btuwzA3e8DzgIuN7MtwOfA2Z7rtjURkXxXWQk33RRmFHbqBOPGhZmEu+ySkbfPaeJw92qgY51998VtjwfGN3VcIiLN0n//C7ffDvfeG9Ys+cUv4IYbYK+9MnoZLTkiItLcbdwI48eHcbuffALnnw+33goHH5yVy+VDH4eIiOyM2lqYOhUOPxyuvx6++lWYMwcefDBrSQOUOEREmqcXX4Rjjw2TPPbZB557LoyeOvLIrF9aiUNEpDlZsAAGDICTT4aVK+HRR6GiIoyeaiJKHCIizcFHH8HFF0OvXvDKKzB6NCxaBOeeGzrCm5ASh4hIFmSsINO6dWFobffu8MgjcM018N57MHgwtGuXwYhTp1FVIiIZlpGCTJs3w6RJcMstoUnqnHPCqKlDD81GyGnRHYeISIY1qiCTO8yYAT17whVXwBFHwL/+FUZP5UHSACUOEZGM2+mCTK++CiecAGeeCW3bwlNPhdFTxxyT8RgbQ4lDRCTDkhVeSlqQ6Z13QrLo2xfefz+spT53bhg9FZbryytKHCIiGZZyQaaPP4YrrwzNUs8+G2Z7L14Ml1wCbfK3Czp/IxMRaaZiHeBlZaF5qqgoJI2tHePV1XDPPTBqVNgeNAhuvhn23z9nMadDiUNEJAtKSxOMoKqpgcmTYdgwWLYMvvc9uPNO6NEjFyHuNDVViYikaKfnZrjDzJmhmNLFF4dbkL//HZ54otklDVDiEBFJSWxuRlVVyAOxuRkNJo/Zs8NyIAMGbCumFBs91UwpcYiIpCDtuRmVlaGtqqQE5s0LNTIWLoSzzsrLkVLpUB+HiEgKUp6bsWZNKKY0blxo0xo6FIYMyXgxpVxS4hARSUFRUWieSrQfaPJiSrmkpioRkRQknZtxW26KKeWSEoeISApKS8OE7uLi0EVRXAx/vPZFSn+Vm2JKuaTEISKSotLS0OddO28BlT0HcMrtuSumlEs5TxxmVmlm881sjplVJDhuZnavmS02s3lmdlQu4hQR2a6Y0quvwpgxOSumlEv50jl+kruvSnLsNKB79Pgq8JvoWUSkaaxbFyrujR0bZn9fc00Yh9uxY64jy4nmkCJPBx7x4DVgbzM7INdBiUhhSTgrfPNm+PWvoVu30Dv+ve/B22+HBNJCkwbkR+Jw4Fkzm21mgxIcPwj4MO710mjfdsxskJlVmFnFypUrsxSqiBSiHWeFO09dOIN1RflbTCmX8iFx9HX3owhNUleY2f/UOZ5oiqXvsMN9oruXuHvJfvvtl404RaRAxc8K78OrvMwJTNt0JitWt4E//SkviynlUs4Th7svi54/Bp4Ajq1zylLgkLjXBwPLmiY6ESkU9S1QuGQJdONdpnMmr9KXL/I+lzCRIzbPg29/u9kvEZJpOU0cZrabme0R2wZOBRbUOe1J4CfR6KrjgLXuvryJQxWRZqzeBQo//piHdr+SNzmCU3mWYQynG4u5n0s4qDhfxg/ll1x/K/sDT1jI5m2Aqe7+ZzO7DMDd7wNmAv2BxUA1cEGOYhWRZirRAoVeXc2yK++BmlH8+LNq7m8ziF9uuZmPCcWUElbsEyDHdxzu/r6794oePd19ZLT/vihpEI2musLdu7r7l919h7keIiIxiZqk4hcibEUNF/Ag79KdwZ/cBP360erNhew2+de0L95/66zwiRMTFGISAMx9h37mZq+kpMQrKpRfRFqaWJNU/N1Fhw7Qvj2sXu2cxjOMYghfZgH/4Dju3n8M0//TfOtiZJKZzXb3klTOzXnnuIhIpiSrmdFry2xeatWPmQygPZ9zFo9xSvtXGXi3ksbOUOIQkYJRtzZGMZVMoZRZa0v46m7zGbbPOHryJhXFZzFxkqkpaicpcYhIwYjVxtiHNYzhehbRgzOYwfg9f0G7Dxdz65or2ei7bC3OJztHiUNECsYdN29gaNu7eI+u/JyxlFPKke3eZZ9fjyyoCny5luvhuCIijVdbC9Omcc7wMs7ZXMVL7b7F1RtGs674y4wcqbuLTFPiEJHm7YUXYPBgeOMN+MpX4IEHOLFfP+blOq4CpqYqEWmWnrpzAS+27w/9+rF07ipeuSwqptSvX65DK3hKHCLSvHz0EYtPvIjThvbiKxteZTCj6VaziFMfOZfyafqV1hT0LYtIXovNBN/L1jFur5vY8sXuHPLXKfwv19KV97iLwWykHdXVYR6HZJ/6OEQkb5WXw08v2cy5n0/kZobTed1Kft/6HIYwkkp2rItRdx6HZEeDdxxm1tXMdo22TzSzq81s7+yHJiKFLNGaUtvtK3aevWwGr3/ekwlcyUJ6UsLr/LBmKh+2TlxMKTaPQ7Irlaaqx4EaM+sGPAAcCkzNalQiUlDqJomf/nTHZc4vuAAuvDBsH+evMnVJXx7+9Ew205YBPMXJvMBswlJKNTVhDap4Ws226aSSOGrdfQswEPhfd/8ZoJrfIrJVosQQe92p07aEEEsS992345pSmzdD8aZ3thZT6kIlFzOJXsxlJgOILwYaW722uBitZpsDqfRxbDazc4DzgO9E+9pmLyQRyXfl5aEjeskS2HdfWL8eNm0Kx6qq4De/2Xbu6tU7/nzdRbn342NuZjiX8ls+pz2/5FbG8nOq2W2Hn43dWZSWKlHkSiqJ4wLgMmCku39gZocCU7Iblojkq7pLlydKDKlqTzU/4x6GMIoOVPNbLuVWhm0tpgTQsSPsvntIUkVFaCZ4HkipHoeZtQeK3H1R9kNqPNXjEMmeLl3CXUVjtKKG83mYW/klB7GMxzmDm9vczrutemy9c4Fwd6EmqKaR0XocZvYdYA7w5+h1bzN7snEhikhz1bghr85pzGS+9eIBLuI/uxZzAi9zXfHjDJ3cgwcfVL9Fc5BKU9UtwLHASwDuPidqrhKRFqioKP07jrZt4WsdZlO2djAn8yLrOneDCdM5+owzeNlsu3OVKPJfKqOqtrj72jr7Cq/erIikZOTIHYfCtm0b+iJidwqXX77tzqHvQZUsOiYUUzq503wYN449P3wTzjwznCDNTip3HAvM7EdAazPrDlwNvJrdsEQkX8XuCGKjqpJ2WP/3v+HAuHGwuhX84hcwZAjsuWeTxyyZlcodx1VAT2AjMA1YB1zb2Aub2SFm9qKZvWVmC83smgTnnGhma81sTvQY1tjrikj66s7TAKisDGUwdqimt2ED3H03dO0KY8eGg+++G5KIkkZBaPCOw92rgbLokUlbgOvc/Q0z2wOYbWbPufubdc77u7t/O8PXFpEU1R1+W1UVXkOdhBEVU6KsLJx02mlw551w5JFNHrNkV9LEYWZ/op6+DHf/bmMu7O7LgeXR9nozews4CKibOEQkh8rKdpzlHVuJdmviSFBMSXUxCld9dxx3NVUQZtYF+ArwzwSH+5jZXGAZcL27L0zyHoOAQQBFWulMJGOSDb9dsgRYsABuuAGeeSZ0djz6KPzoR6FNSwpW0sTh7n9tigDMbHfCQorXuvu6OoffAIrd/VMz6w/8Aeie6H3cfSIwEcIEwCyGLNKiJBp+eyAfMbbDMOg1OfRbjBkDV14J7drlJEZpWkn/LDCz30fP881sXt1HJi5uZm0JSaPc3WfUPe7u69z902h7JtDWzDpl4toikpr44bd7sI7buIl36c5ZG6fANdfA4sVw/fVKGi1IfU1VsVFOWemYNjMjLNP+lruPTXLOF4AV7u5mdiwh0TViZRwRSabuwoUAa9aEO44Lzt3MPtMnctWa4XRmJZV9zqFL+Ug4VHOBW6L6mqqWR5s/dfch8cfMbBQwZMefSktf4MfAfDObE+37BVAUXf8+4CzgcjPbAnwOnO2pLK4lIg2qb4XbbQsXOkdVPcE1k26ku78LJ54IY8bQpSSlJY2kQDW4yKGZveHuR9XZN8/d83aMnRY5FKlf3SG2ifThVe7ieo7nHyzkCO7abzQPreiv2d4FKiOLHJrZ5WY2H+hRp3/jAyAjfRwikhuJhtjGdCdxMaWHVw1Q0hCg/j6OqcAzwB3AjXH717v7mqxGJSJZlWiIbUPFlIo1yl0i9fVxrAXWAueYWWtg/+j83c1sd3dv1OLKIpI78UNsO/DZ1mJK7fk8YTEl1fOWeKnU47gSWAE8BzwdPZ7KclwikkUjR8Lu7Wu4kAd4h8MYwS+ZZd/ghL0XcpVNoKbj/tutdqu6GBIvldVxrwV6uLuGwYoUAndK93mGAR1vYO+lC/kHx3HN/r9j4N0n8JqSg6QglcTxIaHJSkSau9mzwxIhL7zA3t26wfTp9DnjDKar01vSkErieB94ycyeJiytDkCySXsikocqK8NQqqlTWd2qE7cwjj9vGsQtG3ahVDlD0pRK4lgSPXaJHiLSTDz22zX894bbOW/dOGpoza9alXFH7Q2sZ09YkmR5dJEGNDgBsDnSBEBp8TZs4I0Lx3PotJHsxVomcz7DuJWPOHiHU4uLww2JtGzpTABs8I7DzPYDbiBUAdy6ipm7n7zTEYpIdsQVUzqqqoqZnMYQRrGALyf9kWTLposkk8qi+eXA28ChwHCgEng9izGJyM544QU45hg491zYd1++wXMMYGa9SQPCnA6RdKSSODq6+wPAZnf/q7tfCByX5bhEJFULFkD//qHi3qpVvHL5FA5dXcHznNLgj2pin+yMVDrHN0fPy81sAKES344NpSLStD76CIYNg8mTQzGl0aOZ1ukqLr6yXdJ1qNq2DafGlksfOVId45K+VBLHCDPbC7gOGAfsCfwsq1GJSHLr1sHo0TB2LGzZEooplZVBx44M7ZJ88cLiYiUKyYwGE4e7x5YXWQuclN1wRCSpzZvD2h/Dh8PKlaG294gRcOihW2tr1C3xGmOmkVOSOamMqnoI2GHMbtTXISJZVF4OZb9wSpbMYFSroXStfZcXOZGRe41hzl9KWDNtxyJMiagDXDIplaaq+AUN2wEDCf0cIpJF5eXw4EWvMHXj4FBMqfYIBvAUM+kPa7dN917dwCpy6gCXTGtwVJW7Px73KAd+AHwp+6GJtBzl5dClC7RqBZ06wVf3eYf2557BrI0nbFdMaSYDgNTXCNHKtpINqdxx1NWdqC64iOy8+H4JM3CHzqxg2OpbkxZTSodmhEu2pNLHsZ7Qx2HR83+AIVmOS6Sg1a353d4bLqaUDjVPSTalMqpqj6YIRKQlidX8bkUNF/AQtzKMA1nODAYylDt4hx5pvZ/mZ0hTqjdxmFl7oBQ4ItpVAUx393rGb6TOzL4F/ApoDdzv7nfWOb4r8AhwNLAa+KG7V2bi2iK5tKTKOY1nGM0NfImFvEofvs9jvErfpD/TsWN4XrMmjKSKbStRSFNL2jluZl8G3gK+Rlifqgr4JvCKme1tZiMac+GojvkE4DRCYjrHzI6oc9pFwH/dvRtwDzCqMdcUyQsVFby8az9mMoBd2ciZTKcvryRNGh06wJQpsGpVeNTWbr9dWamkIU2rvjuOe4FL3P25+J1mdgqwAFjYyGsfCyx29/ej9/0/4HTgzbhzTgduibanA+PNzLwQ14KXwvfBB6GNato0jtqjEz+vHce4zZeyhbbAtg7y+DsL3U1IPqpvOO4BdZMGgLs/T1i/amAjr30QoSxtzNJoX8Jz3H0LYfZ6x0ZeV6RprVkD110Hhx8Of/gDDB1Kuw8Xc/RDV3JQcVvMwgioRx8NiUN3E5Lv6rvjaGVmu7r7xvidZtaOsFJukhVxUpZoMHrdO4lUzonFNQgYBFCkabKSD6JiSl1/N5I9atfy2G4XMPeM4UydejBL7tTdhDRf9d1xPAI8bmZdYjui7d8Dj2bg2kuBQ+JeH8yOM9K3nmNmbYC9gDWJ3szdJ7p7ibuX7LfffhkIT2Qn1dbC1Kl8esjhHDVtMK/U9qEXczn7swe449GDqaoKdxZVVWFIbnl5rgMWSU/SxOHuI4A/A38zs1Vmtgr4K/Ccu9+WgWu/DnQ3s0PNbBfgbODJOuc8CZwXbZ8FvKD+DclrsWJKpaVUrtuXfjxfbzGl6urQ7SHSnNQ7HNfdxxM6pPeIXq/P1IXdfYuZXQn8hTAc90F3X2hmtwIV7v4k8ADwqJktJtxpnJ2p64tk1Pz5MGQIPPMMn3Us4hcdpzBu9Tl4CrXSVLpVmpuUlhzJZMKo874zgZl19g2L294AfD8b1xbJiKiYUu1Dk1lvezKC0YxffRUbaJfyW6hLTpqbnVmrSkTWrmXBT0bT9U/30MprmMC1jPAy/su+ab2NlgaR5iiVmuMiErNpE4wfz4ZDuvGlJ29nhg/kcN7mOu5uMGkUF8Pll4fn2BBcrVwrzVEqixx2IJSNLXL3S8ysO9AjrjKgSOFzhxkz4MYbYfFi/r3riVzFGGZTktKPa6VaKSSp3HE8BGwE+kSvlwKNWm5EpFl55RXo2xfOOotPNuzK+Z2f5viNL6ScNNQcJYUmlcTR1d1HE2aL4+6fk04lGZHmatEiOOMMOOEEqKzktYsmUbR6Dg9/3J+G/glYdFjNUVKIUkkcm6JVch3AzLoS7kBECtOKFXDFFdCzJzz3HNx2G7z7Lmc/fzHrP0/euhufLGLLh2jJEClEqYyqupkwEfAQMysH+gLnZzMokZz47DMYOxZGj4YNG+DSS+Hmm6FzZ6D++RbFxVo+RFqOVAo5PWdmbwDHEe7Pr3H3VVmPTKSpbNkCkyfDsGGwfDmccQZP9rmDq8cfxpLfbKt9kWzNAnV8S0uTNHGY2VF1di2PnovMrMjd38heWCJNwB1mzgwzvhcuhD594LHHKK/su11Z19Wrk7+FOr6lJarvjuPueo45cHKGYxFpOhUVMHgwvPQSdOsG06eHjnAzykq3JY36qHlKWqqkicPdT2rKQESyobw8LCK4ZElY2uPHJ3xAyRNlnF49jZW2H2N2G8//Lh7Enpe2hUtD6YxUltE0U/OUtFypTABsB/wUOIFwp/F34L5oHSmRvFVeztYmp31Yw1VVI7myajw1tGYEZYz2G1j/2Z5A/c1RiWh9KWnJUhlV9QiwHhgXvT6HUI9Diw9KXisrg5rqDfycCdzECPZkHQ9xATcznGU7FJtMnfo1pKVLJXH0cPdeca9fNLO52QpIJCNqazmhahojKKMLVczkNIYwKmldjGTMto2qUg1wkSCVxPFvMzvO3V8DMLOvAq9kNyyRRpg1CwYPZgr/ZjZHcSEP8uJOjOXQMFuRxFKZOf5V4FUzqzSzSuAfwNfNbL6ZzctqdCLpmD8fTjsNTjkF1qzhlcun8PX2r+9U0lBzlEhyqdxxfCvrUYg0RlRMicmTYc894a674Ior6NuuHb/tu/2oqv79w9SNJUu2b4JSc5RI6lKZOV5lZvsAh8SfrwmAknNr14blQe65B2pq4NprQ5bYd1tdjNJSJQCRTGuwqcrMbgPmAfcSJgXeDdyV5bhEkouKKdGtG9x+OwwcGFayvftu2HdfysuhSxdo1So8l5fnOmCRwpJKU9UPCEurb8p2MCL1qlNMiRNPhDFjKF9UQtmJ25qf1q8PuQWgqirM5QDdeYhkSiqd4wuAvbMdiEi94oopseuu8PTT8MILlC8qYdCgkCDcw0S+TXX+xKmuDi1YIpIZqdxx3EEYkruAuDoc7v7dnb2omY0BvgNsAt4DLnD3TxKcV0mYfFgDbHH31EquSeFYtAiGDoUnnoADDoBJk+D886FN+F+3rCy1daXqWxJdRNKTSuJ4GBgFzAdqM3Td54Ch7r7FzEYBQ4EhSc49Scu4t0ArVsDw4aF8Xvv2oZjSz34Gu+0GbFuDqqoqtbfTEiEimZNK4ljl7vdm8qLu/mzcy9eAszL5/tKMNVBMCbZfgyoVmpMhklmp9HHMNrM7zKyPmR0Ve2QwhguBZ5Icc+BZM5ttZoMyeE3JN1u2wP33Q/fuYU7GqaeGGhkTJkDnztuNlDrvvPqTRtu20LFjWC5ENb9FMi+VO46vRM/Hxe1rsB6HmT0PfCHBoTJ3/2N0ThmwBUg2YLKvuy8zs87Ac2b2trv/Lcn1BgGDAIrULtF8uIeO7iFD4M03tyumVHZq4pFSNTXJ3041MkSyzzyV4gPZuLDZecBlQD93b7DRwcxuAT519wbnkJSUlHhFRUXjg5Tsmj0brr8eXnqJdV/ozuCaO5m0ciD7drTtEkWqtLaUyM4zs9mpDkBK5Y4DMxsA9ATaxfa5+607Fx6Y2bcIneFfT5Y0zGw3oJW7r4+2TwV2+pqSRz74IPRsT5sG++3H6+eN55TfD2Ld522B9GtjgPoxRJpSKjPH7wN+CFwFGKEOR3Ejrzse2IPQ/DQnugZmdqCZzYzO2R94OVrC/V/A0+7+50ZeV3KkvBx6HbKGsXYdG794OJse+wP37nUTe61cTJ8pV2xNGulo3Vr9GCK5kModx/HufqSZzXP34WZ2NzCjMRd1925J9i8D+kfb7wO9Ep0nzcu0hzaw8NLxvLR55LZiSluGs2xtVEypnj6LZDp0ULIQyZVURlV9Hj1Xm9mBwGbg0OyFJM1JonWhYvtaWy3XdprCCRf34PbNg/kHfejNHC7h/rQr8GmklEj+SOWO4ykz2xsYA7xBGFE1KatRSbNQdz5FVRVccEH45d530wvM4HqOWh2KKZ3HQ2nVxWjbNqyQrmXORfJPKsuq3xZtPm5mTwHt3H1tdsOS5iDRch89Ns9nFEPozzNUUkwpU5jGOXgKN7etW0NtrRKFSL5L+q/ZzI4xsy/Evf4J8HvgNjPbN9nPScsRv/7TQSzlAS5kLr3owz+4jrs4nLeZSmlKSaNDB3j44ZA4KiuVNETyWX3/on9LWIQQM/sf4E7gEWAtMDH7oUmuJKtnUXf/vvvCnqxlBGW8w2GUUs5Yfk5X3mMs17Fx2+htOnYMfROxPorLL9/+tfosRJqP+pqqWrv7mmj7h8BEd3+c0GQ1J/uhSbbFFgqMlVWNzYOo228xaFBY1fzhh7ftX1a1ictbTeQmhrMfqyjnR5QxkmVtu2BG9CdH0KED/OpXSgwihaLexGFmbdx9C9CPaDmPFH5OmoFEHduDBoWFaOv2W1RXhzuCsNSHcwYzuJMb6V67mL+3OYkLOo1h5oqjKSqCh6LkUzchKWmIFI76mqqmAX81sz8ShuT+HcDMuhGaqyRPNFQqNdHxRB3b1dXJZ23X1MDxvMIr9OVxzmIju9Kfp/n6llk8tfzo7fomSkvDtvorRApT0jsHdx9pZrOAA4BnfduiVq0Is8glDyS7c4DwCzvZ8VSXJAc4jEXcyVAG8gQfcSAXcT+TOZ9aWlPc2DUERKTZydkih9nUkhY57NIlcTGj2IJ/yY63bp14ldmOHeHzz0Ni6cwKbmY4g5jI5rYdGMUQxmy+lmpCMSXN3hYpHOkscpjKzHHJY8lKosb2JzteUxN+8ceLdWI/cO9n3LXXbSymG5cwife+cRntly6m+0Nl7Fe8m0ZCibRwShzNXLLSI7H9yY7HfvHHD4md9JstlH5+P2f/sjvXrR3GHgO/QdtFC+nx7Hjo3Fl9FyICKHE0eyNHJr5ziA2tre/41kRQ41SOf4ofjeoFl1wS2rdeeQVmzIDDDmuKjyEizYgSRzNXWrrjnUN8E1JDx6mogJNPhu98BzZvhscfD0nj+ONz9plEJL+pc7ylii+m1KkT3HwzXHppWF1QRFqcjFcAlAKyZg2MGAETJoShVWVlcMMNYSlaEZEUKHG0FBs2wLhxcPvtsG5dWP98+HA4KL26GCIi6uPIcw3NCm9QbS1MmQI9eoQ7i+OPhzlz4P77lTREZKcoceSx2Kzvqipw3zbrO+Xk8fzzUFICP/5x6MeYNQuefhq+/OWsxi0ihU2JI48lW0+qrKyBH5w3D047Db7xjdCnUV4Or78eRk+JiDSSEkcea2hW+A6WLoULL4TeveG11+Cuu+Dtt+FHPwptXSIiGZCT3yZmdouZfWRmc6JH/yTnfcvMFpnZYjO7sanjzLWGZoVvtW5duA057LBwd3HddfDee+G5XbuE7yEisrNy+WfoPe7eO3rMrHvQzFoDE4DTgCOAc8zsiKYOMhtS7fBuaFazf5ICAAAOk0lEQVQ4mzbB+PHQtWsYLTVwICxaBGPGhPJ8IiJZkM/tF8cCi939fXffBPwfcHqOY2q0dDq8k876/pGHGd49e8JVV4XO7oqKbRlJRCSLcpk4rjSzeWb2oJntk+D4QcCHca+XRvuatXQ7vHdYWLBLtBzIWWfBrruGUVKzZsHRR2c7dBERIIuJw8yeN7MFCR6nA78BugK9geXA3YneIsG+pOujmNkgM6sws4qVK1dm5DNkQ9od3jGLFsEZZ8AJJ4STH3gA5s6F/v3D7YiISBPJ2sxxdz8llfPMbBLwVIJDS4FD4l4fDCyr53oTgYkQ1qpKPdKmVVSUuLBSso5wVqwIM7wnTgwdHCNGwLXXwm67ZTVOEZFkcjWq6oC4lwOBBQlOex3obmaHmtkuwNnAk00RXzY12OEd89lncNtt0K0bTJoEl10GixeHNi0lDRHJoVz1cYw2s/lmNg84CfgZgJkdaGYzAdx9C3Al8BfgLeD37r4wR/FmTIPLnG/ZEpYD6d4dhg2Db34TFi4Mo6c6d85p7CIioGXV84d76OgeMgTefBP69AkT+FQXQ0SagGqONzcqpiQizYgSRy598AGccw4cc8y25qiFC8PoKY2UEpE8pXocubB6degNHz8e2rRRMSURaVaUOJpSrJjSyJGwfj2cfz7ceqvqYohIs6KmqqZQt5hS376hmNIDDyhpiEizo8SRbSqmJCIFRokjW1RMSUQKlBJHpi1dChdcoGJKIlKw1DmeKWvXwqhRcM89oU/juutg6FDVxRCRgqPE0VibNsFvfxtGR61aFdYOGTFCdTFEpGCp7WRnucP06aGY0tVXw5FHhhngU6YoaYhIQVPi2BkvvxyWA/n+90NN76efDqOnVExJRFoAJY50LFoU6np/7WvbiinNmaNiSiLSoihxpGLFCvjpT0Oz1KxZoQ/j3XfhwguhdetcRyci0qTUOV6fzz6Du++GMWPCciGXXRZqZKguhoi0YEociWzZApMnhySxfDmceSbcfjscdliuIxMRyTkljnh1iykdf3wYOaW6GCIiW6mPI6ZuMaUZM7aNnhIRka2UOOoWU5owITwPHKiRUiIiCbTcpqo1a8LoqAkTwsiom26CwYNVTElEpAEtL3HEiindfjusWxeG1A4fDgcemOvIRESahZwkDjP7HdAjerk38Im7905wXiWwHqgBtrh7yU5ftLYWpk4NZVqXLIEBA+DOO+FLX9rptxQRaYlykjjc/YexbTO7G1hbz+knufuqRl4QvvnNbcuCTJ4MJ53UqLcUEWmpcto5bmYG/ACYluULwbnnhjuOf/0rY0mjvDysZ9iqVXguL8/I24qI5LVc93F8DVjh7u8mOe7As2bmwG/dfeJOX+m883b6RxMpL4dBg6C6OryuqgqvIaysLiJSqMzds/PGZs8DX0hwqMzd/xid8xtgsbvfneQ9DnT3ZWbWGXgOuMrd/5bk3EHAIICioqKjq6qqMvExkurSJSSLuoqLobIyq5cWEck4M5udaj9y1hJHgxc2awN8BBzt7ktTOP8W4FN3v6uhc0tKSryioqLxQdajVavQdVKXWeiHFxFpTtJJHLns4zgFeDtZ0jCz3cxsj9g2cCqwoAnjq1dRUXr7RUQKRS4Tx9nU6RQ3swPNbGb0cn/gZTObC/wLeNrd/9zEMSY1ciR06LD9vg4dwn4RkUKWs85xdz8/wb5lQP9o+32gVxOHlbJYB3hsWkhRUUga6hgXkULXItaqytaw2dLS0BFeWxuelTREpCXI9XDcrNOwWRGRzCr4O46ysm1JI6a6OuwXEZH0FXziWLIkvf0iIlK/gk8cGjYrIpJZBZ84NGxWRCSzCj5xlJbCxIlhKRCz8DxxojrGRUR2VsGPqoKQJJQoREQyo+DvOEREJLOUOEREJC1KHCIikhYlDhERSYsSh4iIpCVnhZyyycxWAo0pAdgJWJWhcDJJcaUuH2MCxZWOfIwJCjeuYnffL5UTCzJxNJaZVaRaCaspKa7U5WNMoLjSkY8xgeICNVWJiEialDhERCQtShyJTcx1AEkortTlY0yguNKRjzGB4lIfh4iIpEd3HCIikhYlDsDMxpjZ22Y2z8yeMLO9k5z3LTNbZGaLzezGJojr+2a20MxqzSzpaAkzqzSz+WY2x8wq8iiuJvu+zGxfM3vOzN6NnvdJcl5N9D3NMbMnsxhPvZ/dzHY1s99Fx/9pZl2yFUsaMZ1vZivjvp+Lsx1TdN0HzexjM1uQ5LiZ2b1R3PPM7Kg8iOlEM1sb910Ny3ZM0XUPMbMXzeyt6N/gNQnOyf735e4t/gGcCrSJtkcBoxKc0xp4D/gisAswFzgiy3H9P6AH8BJQUs95lUCnJvy+Goyrqb8vYDRwY7R9Y6L/htGxT5vg+2nwswM/Be6Lts8GfpcHMZ0PjG+q/4/irvs/wFHAgiTH+wPPAAYcB/wzD2I6EXgqB9/VAcBR0fYewDsJ/jtm/fvSHQfg7s+6+5bo5WvAwQlOOxZY7O7vu/sm4P+A07Mc11vuviib19gZKcbV1N/X6cDD0fbDwPeyeK2GpPLZ4+OdDvQzM8txTDnh7n8D1tRzyunAIx68BuxtZgfkOKaccPfl7v5GtL0eeAs4qM5pWf++lDh2dCEhW9d1EPBh3Oul7PgfLFcceNbMZpvZoFwHE2nq72t/d18O4R8X0DnJee3MrMLMXjOzbCWXVD771nOiP1rWAh2zFE+qMQGcGTVvTDezQ7IYTzry9d9eHzOba2bPmFnPpr541Lz5FeCfdQ5l/ftqEYWcAMzseeALCQ6Vufsfo3PKgC1AeaK3SLCv0UPSUokrBX3dfZmZdQaeM7O3o7+YchlXxr+v+mJK422Kou/qi8ALZjbf3d9rTFwJpPLZs/L/Uz1Sud6fgGnuvtHMLiPcEZ2cxZhS1dTfVSreICzR8amZ9Qf+AHRvqoub2e7A48C17r6u7uEEP5LR76vFJA53P6W+42Z2HvBtoJ9HDYV1LAXi/wI7GFiW7bhSfI9l0fPHZvYEoVmiUYkjA3Fl/PuqLyYzW2FmB7j78ui2/OMk7xH7rt43s5cIf7FlOnGk8tlj5yw1szbAXmS3aaTBmNx9ddzLSYT+vnyQlX97jRH/y9rdZ5rZr82sk7tnfQ0rM2tLSBrl7j4jwSlZ/77UVEUYbQIMAb7r7tVJTnsd6G5mh5rZLoQOzayNykmVme1mZnvEtgkd/QlHgjSxpv6+ngTOi7bPA3a4KzKzfcxs12i7E9AXeDMLsaTy2ePjPQt4IckfLE0WU5128O8S2s/zwZPAT6LRQscBa2PNkrliZl+I9UmZ2bGE36Wr6/+pjFzXgAeAt9x9bJLTsv99NfWogHx8AIsJbYJzokdstMuBwMy48/oTRjG8R2iyyXZcAwl/PWwEVgB/qRsXYZTM3OixMF/iaurvi9A/MAt4N3reN9pfAtwfbR8PzI++q/nARVmMZ4fPDtxK+OMEoB3wWPT/3r+ALzbBf7eGYroj+n9oLvAicHi2Y4quOw1YDmyO/r+6CLgMuCw6bsCEKO751DPCsAljujLuu3oNOL6JvqsTCM1O8+J+X/Vv6u9LM8dFRCQtaqoSEZG0KHGIiEhalDhERCQtShwiIpIWJQ4REUmLEofkLTPrGLf66H/M7KNo+xMzy8b8i/pi6R3NEI69/m6iFWZTfK/KaB5J3f17mdkjZvZe9Ci3JKv8NkZ9n8XMbjGz6zN9TSksShySt9x9tbv3dvfewH3APdF2b6A209eLZnAn05swXj4W25PufmeGQ3gAeN/du7p7V8Icj8kZvgY0zWeRAqbEIc1VazObFNUkeNbM2gOYWVcz+3O04OPfzezwaH+xmc2KFvCbZWZF0f7JZjbWzF4ERkUz8R80s9fN7N9mdno00/pW4IfRHc8PLdSuGB+9x/4W6rjMjR7HR/v/EMWx0BpYfNLMugFHA7fF7b4V6GVmPSzUf3gq7vzxZnZ+tD0sineBmU2Mm9H8kpmNMrN/mdk7Zva1hj5LnZiSfZffj64118watbSNNE9KHNJcdQcmuHtP4BPgzGj/ROAqdz8auB74dbR/PGGp6SMJi1jeG/dehwGnuPt1hAUTX3D3Y4CTgDFAW2AYoWZGb3f/XZ1Y7gX+6u69CDUcFkb7L4ziKAGuNrP6Vr89Apjj7jWxHdH2vwn1T+oz3t2PcfcvAe0Ja67FtHH3Y4FrgZs9LKle32eJl+y7HAZ8M/q8320gNilALWaRQyk4H7j7nGh7NtDFwoqhxwOP2bbSFrtGz32AM6LtRwmFn2Iei/uFfSrw3bh2/nZAUQOxnAz8BLb+sl8b7b/azAZG24cQkl2y9YyMxCuYplKj4yQzuwHoAOxLSFx/io7FFsGbDXRJ4b3CRev/Ll8BJpvZ7+PeX1oQJQ5prjbGbdcQ/tJuBXwS9YM0JP6X9Gdx2wac6XUKVZnZV9MJzsxOBE4B+rh7tYWVeNvV8yMLga+YWSt3r43eoxVwJGEJ7yK2byFoF53TjnAnUOLuH5rZLXWuE/ueakjv33vS79LdL4u+jwHAHDPr7duvrCsFTk1VUjA8LHX9gZl9H7bWXu4VHX6VsCIsQCnwcpK3+QtwVVw/wVei/esJpToTmQVcHp3f2sz2JCyT/t8oaRxOKOFZX+yLCc1SN8XtvgmY5e5LgCrgCAu1yvcC+kXnxJLEqugu4az6rpPCZ4nFk/S7NLOu7v5Pdx8GrGL7JbylBVDikEJTClxkZrHVgmPlUa8GLjCzecCPgWuS/PxthD6NeWa2gG2d1S8SfnHPMbMf1vmZawjNRfMJTUI9gT8DbaLr3UZYQbUhFxKWPl9sZisJyeYyAHf/EPg9YVXUckKSwd0/IdTOmE8oJvR6Ctep77PES/ZdjjGz+dH38zfCCrHSgmh1XJE8ZGY9gJmEzumZuY5HJJ4Sh4iIpEVNVSIikhYlDhERSYsSh4iIpEWJQ0RE0qLEISIiaVHiEBGRtChxiIhIWv4/RuCfzPfBAosAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = model.resid \n",
    "fig = sm.qqplot(res, line='s')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}