3.Negative_Binomial_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from statsmodels.discrete.discrete_model import NegativeBinomial\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "year           2001.031593\nweekofyear       26.503434\ntotal_cases      24.675137\ndtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAELCAYAAAA1AlaNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XmUXOV55/HvU1uvarWW1o6QBGKRwGCQwU5sH9sYA7ZjxWN8EIxjZoYcxmOTZJLxcWAm5tjEmTEZT0hyTCbBwWOGmICDPTOKjY2xIYuXyAgECEmWaQRIQrtaaqnVS23P/HFvyaWmW13dXd2tvu/vc06frrr1VtX7FuJXbz/33veauyMiImFITXUHRERk8ij0RUQCotAXEQmIQl9EJCAKfRGRgCj0RUQCUlPom9m1ZrbdzDrN7PYhHm8ws0fixzeY2bJ4+zIz6zOz5+Kfv6xv90VEZDQyIzUwszRwL3A1sBt42szWu/vWqma3AEfc/VwzWwfcDdwQP/ayu19a536LiMgY1DLTvwLodPcd7p4HHgbWDmqzFnggvv0ocJWZWf26KSIi9TDiTB9YDOyqur8buHK4Nu5eNLNuYE782HIz2wQcA/7A3f/5dG82d+5cX7ZsWQ3dEhGRimeeeeaQu3eM1K6W0B9qxj547Ybh2uwFlrr7YTO7HPi/Zrba3Y+d8mSzW4FbAZYuXcrGjRtr6JaIiFSY2Wu1tKulvLMbOKvq/hJgz3BtzCwDzAS63H3A3Q8DuPszwMvAeYPfwN3vc/c17r6mo2PELyoRERmjWkL/aWClmS03sxywDlg/qM164Ob49vXAk+7uZtYR7wjGzFYAK4Ed9em6iIiM1ojlnbhGfxvwOJAGvuruW8zsLmCju68H7gceNLNOoIvoiwHgncBdZlYESsAn3L1rIgYiIiIjszNtaeU1a9a4avoiIqNjZs+4+5qR2umMXBGRgCj0RUQCotAXEQmIQl9EJCAKfRGRgNRyRm5iPLRh55Dbb7py6ST3RERkamimLyISEIW+iEhAFPoiIgFR6IuIBEShLyISEIW+iEhAFPoiIgFR6IuIBEShLyISEIW+iEhAFPoiIgFR6IuIBEShLyISEIW+iEhAFPoiIgFR6IuIBEShLyISEIW+iEhAFPoiIgFR6IuIBEShLyISEIW+iEhAFPoiIgFR6IuIBEShLyISEIW+iEhAFPoiIgFR6IuIBEShLyISkJpC38yuNbPtZtZpZrcP8XiDmT0SP77BzJYNenypmfWY2afr020RERmLEUPfzNLAvcB1wCrgRjNbNajZLcARdz8XuAe4e9Dj9wDfHX93RURkPGqZ6V8BdLr7DnfPAw8Dawe1WQs8EN9+FLjKzAzAzH4d2AFsqU+XRURkrGoJ/cXArqr7u+NtQ7Zx9yLQDcwxsxbg94HPn+4NzOxWM9toZhsPHjxYa99FRGSUagl9G2Kb19jm88A97t5zujdw9/vcfY27r+no6KihSyIiMhaZGtrsBs6qur8E2DNMm91mlgFmAl3AlcD1ZvbHQDtQNrN+d//yuHsuIiKjVkvoPw2sNLPlwOvAOuCmQW3WAzcDPwWuB550dwfeUWlgZp8DehT4IiJTZ8TQd/eimd0GPA6kga+6+xYzuwvY6O7rgfuBB82sk2iGv24iOy0iImNTy0wfd38MeGzQtjurbvcDHx3hNT43hv6JiEgd6YxcEZGAKPRFRAKi0BcRCYhCX0QkIAp9EZGAKPRFRAKi0BcRCYhCX0QkIAp9EZGAKPRFRAKi0BcRCUjQoV8qD74sgIhIsgUb+nuO9vH5v9/C3u6+qe6KiMikCTb0n9l5hGLZOXBsYKq7IiIyaYIM/VLZ2by7G4AT+eIU90ZEZPIEGfo7DvXQMxCF/YkBhb6IhKOmi6gkzfO7umnIpDCDEwOlqe6OiMikCW6mXyiV2bKnm4sWzaStMavyjogEJbjQ377vOAPFMpec1U5LQ0YzfREJSnCh/4v9x2nKplnR0UJzLq2ZvogEJbjQ7yuUaG3MkDKLZ/oKfREJR3ChXyiVyaWjYbfkMvTlSzozV0SCEVzo54tOLhOHfkMaB4725qe2UyIikyS40D9lpt8QHbHadUKhLyJhCC7088Uy2bQBUXkH4LBCX0QCEV7ol8qnlHcAjij0RSQQ4YV+sUx2UHlHM30RCUVwoV+omuk356KZvmr6IhKKoEK/7E6x7Cd35GZSKRqzKYW+iAQjqNDPF8sAJ2f6EO3MVXlHREIRVuiXotCv1PQhqut3ndCFVEQkDEGFfmHImX6awz2a6YtIGIIK/cpMP/eGmb5CX0TCEFboDzXTb8hwpDePu9bfEZHkqyn0zexaM9tuZp1mdvsQjzeY2SPx4xvMbFm8/Qozey7+ed7MPlzf7o9OoRQF+yk1/VyaQsk5rtU2RSQAI4a+maWBe4HrgFXAjWa2alCzW4Aj7n4ucA9wd7z9RWCNu18KXAv8lZlN2SUah5rpN1fW31FdX0QCUMtM/wqg0913uHseeBhYO6jNWuCB+PajwFVmZu7e6+6VKXQjMKU1lCFr+lp/R0QCUkvoLwZ2Vd3fHW8bsk0c8t3AHAAzu9LMtgCbgU9UfQlMuiGP3mnQWbkiEo5aQt+G2DZ4xj5sG3ff4O6rgbcAd5hZ4xvewOxWM9toZhsPHjxYQ5fGZrijdwAdqy8iQagl9HcDZ1XdXwLsGa5NXLOfCXRVN3D3bcAJ4KLBb+Du97n7Gndf09HRUXvvR+nkyVmZX35HqbwjIiGpJfSfBlaa2XIzywHrgPWD2qwHbo5vXw886e4ePycDYGZnA+cDr9al52OQL5ZJWbTmTkUuk6Ipm9aOXBEJwohH0rh70cxuAx4H0sBX3X2Lmd0FbHT39cD9wINm1kk0w18XP/3twO1mVgDKwCfd/dBEDKQW1WvpV5vdklNNX0SCUNPhk+7+GPDYoG13Vt3uBz46xPMeBB4cZx/rplAsn1LPr5jTmqNL18kVkQCEdUZuqXzKiVkVmumLSCjCCv3iMOWd5pwWXRORIIQV+prpi0jgggr9wnAz/dYcfYUSffnSFPRKRGTyBBX6+dIwO3JbcgAc1glaIpJwQYV+oeTDHLLZAGgpBhFJvqBCP18cvqYPOitXRJIvrNAvlcml37hMUKW8o7NyRSTpggl9dz/tjlyAIzpBS0QSLpjQHyiWcRhyR+6MhgzZtKm8IyKJF0zo98aHY2aHmOmbGbOacyrviEjiBRT60bVbhprpQ7QzVzN9EUm6YEK/cuLVUDV9iBdd03H6IpJwwYR+pbwz/Ey/Qcfpi0jiBRf6Q9X0ITpsU+UdEUm6YEK/rzByTf94f5F8fPF0EZEkCib0e0eo6VfOytWx+iKSZOGF/mlm+qD1d0Qk2YIJ/f7C6Wv6Cn0RCUEwoT/STH+OFl0TkQAEF/qZIRZcg6qZfo+O1ReR5Aom9PvyRbJpI2VDh357cw4zlXdEJNmCCf3efGnY0g5AOhWtv6PyjogkWTCh35cvDXu4ZoUukC4iSRdM6PfmS0NeNauaFl0TkaQLJ/QLI8/057TkOKLQF5EECyb0+/LF09b0AWapvCMiCRdM6PfWUNOf05LjSG+ectknqVciIpMrmNDvq6GmP7e1gbJDl9bfEZGECib0a5npL5zZCMCeo32T0SURkUkXUOiPXNNf1N4EwOtHFPoikkzBhH5fDUfvLJkVh75m+iKSUJmp7sBkKJbKFEo+bE3/oQ07AXB3cpkUT/78AM25DDdduXQyuykiMuGCmOn3Fk5/AZUKM6O9KcvR3sJkdEtEZNIFEfr9levjDrPCZrX25ixH+3T0jogkU02hb2bXmtl2M+s0s9uHeLzBzB6JH99gZsvi7Veb2TNmtjn+/Z76dr82fYXTr6Vfrb0pp5m+iCTWiCloZmngXuA6YBVwo5mtGtTsFuCIu58L3APcHW8/BPyau18M3Aw8WK+Oj0bvyZl+DaHfnKU3X9IF0kUkkWqZ6V8BdLr7DnfPAw8Dawe1WQs8EN9+FLjKzMzdN7n7nnj7FqDRzBrq0fHR6Kuxpg/RuvqASjwikki1hP5iYFfV/d3xtiHbuHsR6AbmDGrzEWCTu7/h0lRmdquZbTSzjQcPHqy17zXrH81MvykLoBKPiCRSLaE/1N7PwYvTnLaNma0mKvn8+6HewN3vc/c17r6mo6Ojhi6NzkjXx63W3qzQF5HkqiX0dwNnVd1fAuwZro2ZZYCZQFd8fwnwf4CPu/vL4+3wWFTKO7UcvTOjMUvKVN4RkWSqJfSfBlaa2XIzywHrgPWD2qwn2lELcD3wpLu7mbUD3wHucPcf16vTo3Uy9Guo6adTRpuO1ReRhBoxBeMa/W3A48A24BvuvsXM7jKzD8XN7gfmmFkn8HtA5bDO24Bzgc+a2XPxz7y6j2IEfaMo7wA6QUtEEqumZRjc/THgsUHb7qy63Q98dIjnfQH4wjj7OG6/LO/UGPrNOV49fGIiuyQiMiWCOCO3siM3U0NNH6Kducf6CpR0MRURSZggQr+/UKIpmyZlNYZ+U46yw/5j/RPcMxGRyRVE6PflSzTl0jW3rxy2qYupiEjSBBH6vflopl+ryglaWldfRJImiNDvL4x2ph8txbCrq3eiuiQiMiWCCP2+wuhm+rlMiplNWXYc0hE8IpIsQYR+b744qtAHmNua4+WDCn0RSZYgQr+vUB5VeQegY0YDOw724K7DNkUkOYII/f5R7sgFmNvawPH+Igd73rAoqIjItBVE6PcWiqOf6bdGy/7vUIlHRBIkiNDvy4+tvAPw8sGeieiSiMiUCCT0R78jt60pS2M2pZm+iCRK4kPf3Ud9yCZAyowVc1vZoZm+iCRI4kM/XypTdkZd3gFY0dGiwzZFJFESH/qVtfRHO9MHOKejld1HeumPl2YWEZnukh/6cWCPdaZfdnjtsJZjEJFkSH7oxzP95jGE/jkdrQCq64tIYiQ+9CsXUGkcQ3lnRUcLoMM2RSQ5Eh/6lXr8WGr6zbkMi2Y26rBNEUmMxId+paY/lvIOwIqOVs30RSQxEh/64ynvAFywYAY/33ecfLFcz26JiEyJxId+/ziO3gG4/OxZDBTLbN17rJ7dEhGZEokP/d5xHL0DcNnZswB45rUjdeuTiMhUSXzoj+fkLID5bY0sbm/iWYW+iCRA8kO/ML6aPkQlnmd3KvRFZPpLfujnS6QMGjJjH+plS9vZ293PnqN9deyZiMjkS37oxytsmtmYX+Pys2cDquuLyPQXRujnMuN6jQsWzqAxm1Loi8i0l/zQz5doyo1vmNl0ikuWtLNJdX0RmebCCP1x7MStuPzsWWzZc+zk0UAiItPR+Ooe08B4yjsPbdh58nZPf5Fi2fnS97fz2Q+uqlf3REQmVSAz/fEPc3lHC5mUsWWPzswVkekr+aE/huvjDqUhk+a8+TPYsqebctnr0DMRkcmX+NDvzRdpHufROxUXL57J8f4iz2iHrohMUzWFvplda2bbzazTzG4f4vEGM3skfnyDmS2Lt88xs6fMrMfMvlzfrtemv1Ae19m41S5YMINMynhs8966vJ6IyGQbMfTNLA3cC1wHrAJuNLPBezJvAY64+7nAPcDd8fZ+4LPAp+vW41GKduTW5w+ahmyalfNn8L0X96nEIyLTUi1peAXQ6e473D0PPAysHdRmLfBAfPtR4CozM3c/4e4/Igr/KVHP8g7ARYva2Nvdz6ZdR+v2miIik6WW0F8M7Kq6vzveNmQbdy8C3cCcWjthZrea2UYz23jw4MFanzaictnrWt4BuHBhG7l0iu+qxCMi01AtoT/UojWDaxu1tBmWu9/n7mvcfU1HR0etTxvRQHy1q3ocvVPRmE3zjpVz+e6L+3BXiUdEppdaQn83cFbV/SXAnuHamFkGmAl01aOD49GbLwJjv4DKcK67eCGvH+3j+d3ddX1dEZGJVkvoPw2sNLPlZpYD1gHrB7VZD9wc374eeNLPgGlwZS39es70Aa6+cD7ZtKnEIyLTzoihH9fobwMeB7YB33D3LWZ2l5l9KG52PzDHzDqB3wNOHtZpZq8CfwL8GzPbPcSRPxOmsk5OY51n+jObs/zquXP5zua9KvGIyLRS02Et7v4Y8NigbXdW3e4HPjrMc5eNo3/jUpnpN9d5pg/w/osW8plvvsCLrx/j4iUz6/76IiITIdELrp28Pm6dZ/oPbdhJb75IyuBL39/ONasXAHDTlUvr+j4iIvWW6GUYeutwfdzhNOcynNPRyubXu1XiEZFpI9Gh3x/P9Ot99E7F6kUz6TqRZ//xgQl5fRGRekt06E/U0TsV5y+YAcD2fccn5PVFROot0aHfO0E1/YqZTVkWzmxk+z6tsS8i00OiQ/9obx6IwnmiXLBgBq8d7j15IpiIyJks0aF/qCfPjIbMhOzIrTh/QRsOvLS/Z8LeQ0SkXhId+gd7Bpg7o2FC32PJrCZacml+rhKPiEwDiQ79Q8cHmNuam9D3SJlx3vwZ/GJ/DyWtsS8iZ7hEh/7hE3nmtEzsTB/ggoVt9BVKbNJlFEXkDJfo0D/UM8DcGRM70wdYOa+VlMET2/ZP+HuJiIxHYkO/UCpztLfA3NaJn+k3ZtOc09HK97fs19m5InJGS2zoH+6JDtecjNAHWLWojVcOneClAzqKR0TOXIkN/UM90dIIkxX6Fy5swwwef3HfpLyfiMhYJD70Oyahpg/Q1pjlzWe18/hWhb6InLkSHPpReWcyjt6puGb1Al58/Ri7unon7T1FREYjwaEfl3cm+OSsapV19b+/VUfxiMiZKbmhf3yAxmyKlglabG0oy+a2cP78GTyma+eKyBkquaHfM8Dc1gbMbNLe86ENOzl3XivPvHaEP/z2Vh7asHPS3ltEpBYJDv38pB25U+2K5bOZ2ZTlia06Zl9EzjwJDv2BKQn9bDrFu8+fx86uXrbv18VVROTMkuDQz0/4YmvDufzsWcxuyfHE1v2UtQibiJxBEhn6pbLTdWJqZvoA6ZTx3gvns7e7n0//3fMMFEtT0g8RkcEyU92BiXCkN0/ZmbKZPsAlS6KLpn9r0+u81tXLvTddxoKZjVPWHxERSOhMfyqO0R/MzHjPBfO496bL2LKnm1+9+0n+7f/6Gd96dje7unq1k1dEpkQiZ/qTvdja6XT3FfjUu85l42tHeHbnUZ7afhCI/gpZe+li/t3bl7O4vWmKeykioUhk6E/2YmsjmdPawDWrF3D1qqjOv6url1LZ+dpPXuVrP3mVay9awIcvXcw7z+sgl0nkH18icoZIZOgfPF4J/amr6Q8lZcbi9qaTM/vVi9r4cechnvr5Ab7zwl4aMilWzm9lQVsTqxe18c7z5nLJknYyaX0RiEh9JDL0D/XkyaaNmU3Zqe7KabU35/jAmxZxzUUL6DzQw8/3Hae7t8Dm14/yw237+bMfvsSMxgw3rDmLP/jgqqnurogkQEJDf4A5LZO7BMN4ZFIpLljQxgUL2k5u680XefngCX6wdT/3/+gVZrfm+MQ7zyGVmh5jEpEzUyLrBocn6dq4E6k5l+HixTP55LvP4eIlM/nj723nUw89S19ex/yLyNglMvR3Hemj4wzZiTteDZl0VN75wIV8b8s+brjvp+w/1j/V3RKRaSpxob9p5xE6D/Tw7gvmTXVX6sbM+M13rOArv7GGzgM9XPU//pH/9tg29nUr/EVkdBJX03/gJ6/S2pDhX122ZKq7UleVZZpvfccKntx+gPv+aQf3/+gV3rFyLh+6dBFXr1pAa0Pi/nOKSJ3VlBJmdi3wZ0Aa+Gt3/+KgxxuA/w1cDhwGbnD3V+PH7gBuAUrAb7v743Xr/SAHjvfznc17+dhbz05sAM5ra2TdW5byvlV5Nrxy+OQJX5nUC1ywYAZvWtLOnb+2isbs8BePKZWdYrlMyiz+oaad3j0DRfryJea25qbNTnIROdWIyWhmaeBe4GpgN/C0ma13961VzW4Bjrj7uWa2DrgbuMHMVgHrgNXAIuAHZnaeu0/I3siHNuykUHI+/rZlE/HyZ5TZLTmuu2gh16xewK6uXp7ffZTNrx/jxT3HWP/8Ht63aj4rOlqY2ZzjWF+BrXuO8dKB4xzqyXOkN8/gVSBSFl1PeGF7I/NmNPC+1QtozqXZ/Ho3m147yssHezh8IjrTub05y3nzZ9DWmKUhPpksXyqTNmP1ojbevHQWqxa1Mbtleu9MF0miWqbDVwCd7r4DwMweBtYC1aG/FvhcfPtR4MsWTQXXAg+7+wDwipl1xq/30/p0/5fyxTJf37CTd53fwfK5LfV++TNWyoyz57Rw9pwWPnDxInYc6qGnv8gT2/bzrU2Fk+3Omt1Ea0OWlfNaaWnIkEkZZQfHcYdiyTl4vJ/XDvfywu5ufrDtAABpMxa1N7Kio4U1y2Zz5fLZvHTgOC/t7+H1o30USmXK7vQOlCiWy3xvy76T79mcS3PhwjYWtDUyr62B+W2Np9ye3ZzDgWKpTFdvnv3HBjhyIk9focRAoYSZkU7FP1W3Uykjk4r+SklXbg9qk05B2eHIiTxHewuYQVMuTVtjlvlt0RdbQzZFJpWiN1/kUE+eQz0DHO7J03VigJaGDPPbGqM+z2wc81+O5bLTVyhxYqBIz0CREwMlTuSL0WdbNYa0GakU8XhSNGRSZNMpsmkjF9/OpVOYRV+whZKTL5YplMqkU0ZTNk0uk8IdyvE3etmdctX9xky6pjO+3Z0T+RLH+goc6y9wrK/Isb4ChVKZWS055rTkmNWSY1ZzjvSgQ4hLZafn5FiLHO8vnhx7T3+Rkju5dIrGbJpZzVnam3PMbsnR3pw95a9TdydfKtOfL9NbKJIvluPPI0UuE30+KTOK5TLFslMsOcX4s2jOZWjMRp+Fn/wcon/nlc+n7I4DXj71fnW7VApachmasulxHSpdKJXpL5ToL5QZKJbIZVK0NkSvOxV/MdfyL3kxsKvq/m7gyuHauHvRzLqBOfH2fxn03MVj7u1pbHyti0M9A9z8K8sm4uWnhXTKWDlvBgBvXjqLQqlMX6F08n+yWhVKZY73F+kvlOiY0UB20BnBFy9u5+LF7UM+ty9fYvfRXvZ393Pg+ACHT+R57XAvx/oL5IvlsQ9uijVl0ye/JDIpI5OOwroSIlFQxEHiUHKnd6BIb6H0hr+qplIm/oJozKXJxf9dPf5ycJyBYpljfQVquQyEGTRn02Qz0WdyYqBEX2Hsf8Rn03ay5DhQLNXUh8nSnEvTnMuQS78xpAcHt3v0OfYXSvQXy5SGGYhZ9KXSnEuTil/jvavm8YVfv7j+A6hSS+gP9VU0eBTDtanluZjZrcCt8d0eM9teQ7+G9J4vnvbhucChsb72NKexhyfUccM0HfsG4I/G/vSza2lUS+jvBs6qur8E2DNMm91mlgFmAl01Phd3vw+4r5YOj4eZbXT3NRP9PmcijT28sYc6bgh77COp5Tj9p4GVZrbczHJEO2bXD2qzHrg5vn098KRHC8avB9aZWYOZLQdWAj+rT9dFRGS0RpzpxzX624DHiQ7Z/Kq7bzGzu4CN7r4euB94MN5R20X0xUDc7htEO32LwKcm6sgdEREZmYV0BSczuzUuJQVHYw9v7KGOG8Ie+0iCCn0RkdAlbu0dEREZXjChb2bXmtl2M+s0s9unuj/1ZmZfNbMDZvZi1bbZZvaEmb0U/54Vbzcz+/P4s3jBzC6bup6Pj5mdZWZPmdk2M9tiZr8Tbw9h7I1m9jMzez4e++fj7cvNbEM89kfiAzCID6h4JB77BjNbNpX9Hy8zS5vZJjP7dnw/iHGPVxChX7WUxHXAKuDGeImIJPkacO2gbbcDP3T3lcAP4/sQfQ4r459bgf85SX2cCEXgP7n7hcBbgU/F/21DGPsA8B53vwS4FLjWzN5KtAzKPfHYjxAtkwJVy6UA98TtprPfAbZV3Q9l3OPj7on/Ad4GPF51/w7gjqnu1wSMcxnwYtX97cDC+PZCYHt8+6+AG4dqN91/gP9HtE5UUGMHmoFnic6WPwRk4u0n/+0THYH3tvh2Jm5nU933MY53CdGX+XuAbxOdCJr4cdfjJ4iZPkMvJTEhy0GcYea7+16A+HflIgOJ/DziP9vfTHRiYxBjj0sczwEHgCeAl4Gj7l6Mm1SP75TlUoDKcinT0Z8CnwEqa3vMIYxxj1sooV/TchABSdznYWatwDeB/+jux07XdIht03bs7l5y90uJZr5XABcO1Sz+nYixm9kHgQPu/kz15iGaJmrc9RJK6Ne0HEQC7TezhQDx7wPx9kR9HmaWJQr8r7v7t+LNQYy9wt2PAv9AtF+jPV4OBU4d38mxD1ouZbr5VeBDZvYq8DBRiedPSf646yKU0K9lKYkkql4e42aiendl+8fjI1neCnRXSiHTTbyE9/3ANnf/k6qHQhh7h5m1x7ebgPcS7dh8img5FHjj2IdaLmVacfc73H2Juy8j+n/5SXf/1yR83HUz1TsVJusHeD/wC6Ka53+Z6v5MwPj+FtgLFIhmNrcQ1S1/CLwU/54dtzWio5leBjYDa6a6/+MY99uJ/lR/AXgu/nl/IGN/E7ApHvuLwJ3x9hVEa1x1An8HNMTbG+P7nfHjK6Z6DHX4DN4FfDu0cY/nR2fkiogEJJTyjoiIoNAXEQmKQl9EJCAKfRGRgCj0RUQCotAXEQmIQl+mNTNrN7NPjtBmmZndVMNrLatemlokiRT6Mt21A6cNfaLVR0cMfZEQKPRluvsicI6ZPWdm/z3+edHMNpvZDVVt3hG3+d14Rv/PZvZs/PMrtbxRvKLll+LXfsHMfivefqeZPR2/733x0hCY2W+b2da47cPxthaLLnjzdHwBkLXx9tXxBVGei9uvrPsnJYKukSvTXLyc8rfd/SIz+wjwCaKLycwlWnPpSuB84NPu/sH4Oc1A2d3743D9W3dfU/1aw7zXfyBa3+YGdy+a2Wx376r8jts8CHzD3f/ezPYAy919wMza3f2omf1XYKu7/028bs7PiJaD/iLwL+7+9Xh9qLS7903ARyaB00xfkuTtRAFecvf9wD8CbxmiXRb4ipltJlqTpdarqL0X+EuP12yvBD3w7vgyfJuJVnxcHW9/Afi6mX2M6ApfAO8Dbo/XwP8HonVhlgI/Bf6zmf3zaXGVAAABbUlEQVQ+cLYCXyZKZuQmItPGUOumD+V3gf3AJUQTn/5RvP4pfxqbWSPwF0QLt+0ys88RBTnAB4B3Ah8CPmtmq+PX+Ii7bx/02tvMbEP8nMfN7Dfd/cka+yVSM830Zbo7DsyIb/8TcENce+8gCtyfDWoD0Xrqe929DPwGkK7xvb4PfKKyZruZzeaXAX8ovpDL9fFjKeAsd3+K6ApP7UAr0aX7fquq7v/m+PcKYIe7/znRUsBvGtWnIFIjzfRlWnP3w2b24/hQy+8SlVSeJ5qRf8bd95nZYaBoZs8TXUD+L4BvmtlHidZgP1Hj2/01cB7wgpkVgK+4+5fN7CtEyzS/SrQfAaIvkr8xs5lEs/t74pr+HxJd8OOFOPhfBT4I3AB8LH7dfcBdY/5QRE5DO3JFRAKi8o6ISEBU3hEZxMyuAe4etPkVd//wVPRHpJ5U3hERCYjKOyIiAVHoi4gERKEvIhIQhb6ISEAU+iIiAfn/UiRcwz4Ouh4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f22cb372128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = pd.read_csv(\"/media/vajira/0FD50F0F0FD50F0F/Sem7_Academic/Data Mining/dengAI/dengue_labels_train.csv\")\n",
    "sns.distplot(y.total_cases)\n",
    "#plt.show()\n",
    "print (y.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> As is evident from the distribution of our y labels in the graph, we have a highly skewed distribution, that combined with the fact our labels are count variables, makes a very convincing case for using negative binomial regression!\n",
    "You can read more about negative binomial regression [here](http://www.mathematica-journal.com/2013/06/negative-binomial-regression/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### importing the relevant packages and loading the data\n",
    "In order to implement negative binomial regression, we use the go to statistical library for python, statsmodels!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "\n",
    "from warnings import filterwarnings\n",
    "\n",
    "filterwarnings('ignore')\n",
    "\n",
    "from statsmodels.tools import eval_measures\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "# from sklearn.feature_selection import RFE\n",
    "\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.metrics import median_absolute_error\n",
    "from sklearn.model_selection import ParameterGrid\n",
    "from sklearn.model_selection import TimeSeriesSplit\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "from scipy import stats\n",
    "\n",
    "from xgboost import XGBRegressor\n",
    "\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.metrics import make_scorer\n",
    "\n",
    "from sklearn.feature_selection import f_regression\n",
    "from sklearn.feature_selection import SelectKBest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from mlxtend.feature_selection import SequentialFeatureSelector as SFS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pre Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def impute(X):\n",
    "    \n",
    "    # remove the colum that has ~20% null values, also ranks low on feature importance\n",
    "    X.drop(['ndvi_ne'], axis=1, inplace=True)\n",
    "    \n",
    "    # Filling the rest using linear interpolation\n",
    "    X.interpolate(inplace=True)\n",
    "\n",
    "def remove_outliers(df):\n",
    "    return df[(np.abs(stats.zscore(df)) < 5).all(axis=1)]\n",
    "\n",
    "def mape(Y_test, Y_pred, epsilon = 1):\n",
    "    return np.mean(np.abs((Y_test - Y_pred + epsilon) / (Y_test + epsilon))) * 100\n",
    "\n",
    "def extract_month(s):\n",
    "    return int(s[5:7])\n",
    "\n",
    "def city_indices(X):\n",
    "    # city boolean encoding\n",
    "    return X.city == 'sj'\n",
    "\n",
    "def pre_process(X, y, trees = False):\n",
    "    \"\"\"\n",
    "    Extracts the month out of date and converts it to a one hot\n",
    "    Standardizes the numerical features\n",
    "    \"\"\"\n",
    "    \n",
    "    #Extracting month from the date\n",
    "    months = X.week_start_date.apply(extract_month)\n",
    "    \n",
    "    #Response coding\n",
    "#     month = X.week_start_date.apply(extract_month)\n",
    "#     temp = pd.DataFrame(y.total_cases)\n",
    "#     temp['month'] = month\n",
    "#     for name,group in temp.groupby(month):\n",
    "#         month[group.index] = np.median(group.total_cases)\n",
    "        \n",
    "    # Removing the columns not required for classification\n",
    "    X.drop(['city', 'year', 'weekofyear', 'week_start_date'], axis=1, inplace=True)\n",
    "\n",
    "    # Standardizing the data\n",
    "    if not trees:\n",
    "        scaler = StandardScaler()\n",
    "        X[X.columns] = scaler.fit_transform(X)\n",
    "\n",
    "    # Month one hot features\n",
    "    month_features = pd.get_dummies(months, prefix='m_')\n",
    "    X = X.join(month_features)\n",
    "\n",
    "    #Alternatively use response coding \n",
    "#     X = X.join(month)\n",
    "    \n",
    "    return X\n",
    "\n",
    "def seperate_cities_data(X, is_sj):\n",
    "\n",
    "    # Seperating the cities data\n",
    "    X_sj = X.loc[is_sj]\n",
    "    X_iq = X.loc[~is_sj]\n",
    "    \n",
    "    return X_sj, X_iq\n",
    "\n",
    "def get_y_labels(X_sj, X_iq, y):    \n",
    "    \n",
    "    y = y.total_cases    \n",
    "    y_sj = y.loc[X_sj.index]\n",
    "    y_iq = y.loc[X_iq.index]\n",
    "    \n",
    "    return y_sj, y_iq\n",
    "\n",
    "def split(X_sj, X_iq, y_sj, y_iq):\n",
    "\n",
    "    # train and test split\n",
    "    sj_split_data = train_test_split(X_sj, y_sj, shuffle = False)\n",
    "    iq_split_data = train_test_split(X_iq, y_iq, shuffle = False)\n",
    "\n",
    "    return sj_split_data, iq_split_data\n",
    "\n",
    "def process(X, y = pd.Series(), train = True, trees = False, feature_selection = 0, time_shift = 0):\n",
    "    \n",
    "    is_sj = city_indices(X)\n",
    "    if not trees:\n",
    "        impute(X)\n",
    "    X = pre_process(X, y, trees)\n",
    "    \n",
    "    if feature_selection:\n",
    "        selector = SelectKBest(f_regression, k=feature_selection).fit(X,y.total_cases)\n",
    "        X = X.loc[:,selector.get_support()]\n",
    "\n",
    "    #X = remove_outliers(X)\n",
    "    X_sj, X_iq = seperate_cities_data(X, is_sj)\n",
    "    if y.empty:\n",
    "        return X_sj, X_iq\n",
    "    \n",
    "    y_sj, y_iq = get_y_labels(X_sj, X_iq, y)\n",
    "    if time_shift:\n",
    "        y_sj = y_sj.shift(time_shift).dropna()\n",
    "        y_iq = y_iq.shift(time_shift).dropna()\n",
    "        X_sj = X_sj[:-time_shift]\n",
    "        X_iq = X_iq[:-time_shift]\n",
    "        \n",
    "    if not train:\n",
    "        return X_sj, X_iq, y_sj, y_iq\n",
    "    \n",
    "    return split(X_sj, X_iq, y_sj, y_iq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = pd.read_csv(\"/media/vajira/0FD50F0F0FD50F0F/Sem7_Academic/Data Mining/dengAI/dengue_features_train.csv\")\n",
    "data = process(X,y)\n",
    "\n",
    "(X_sj_train, X_sj_test, Y_sj_train, Y_sj_test), (X_iq_train, X_iq_test, Y_iq_train, Y_iq_test) = data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Negative Binomial Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "formula = ' + '.join([str(i) for i in list(X_sj_train.columns)])\n",
    "formula = 'y ~ ' + formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We first examine results seperately on the two cities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The city of San Juan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      ndvi_nw   ndvi_se   ndvi_sw  precipitation_amt_mm  \\\n702 -1.067275 -1.470795 -1.364071              1.168702   \n703  0.298713 -0.886888 -0.924523              0.785077   \n704 -0.606574 -1.866210 -1.522569              1.181994   \n705 -0.944617 -1.016692 -0.802140              7.903910   \n706 -1.386481  0.458440  0.077812             -0.676090   \n707 -0.623109 -1.024504 -1.012159             -1.047340   \n708 -2.357911 -0.561945 -1.128550              3.142056   \n709 -1.510604  1.185482 -0.050560             -1.047340   \n710 -1.403435 -0.610019 -0.647219             -1.002423   \n711 -2.105896  0.279664 -0.304567             -1.024881   \n712 -0.926197  0.157254  0.245382             -1.047340   \n713 -1.048437  1.070689  1.053617             -0.220735   \n714 -2.565760 -2.964682 -2.296999             -1.047340   \n715 -1.974446  0.449142  0.939964             -1.047340   \n716 -1.479626  0.674015  0.179997             -1.047340   \n717 -0.876799  0.586467 -0.112693             -1.047340   \n718 -1.713221 -0.668246 -0.629779             -1.047340   \n719 -0.634622 -0.869489 -0.521432             -1.047340   \n720 -0.312486  0.236860  0.405077             -1.047340   \n721 -0.833680 -0.222763 -0.577060             -1.047340   \n722 -0.644529 -2.369800 -2.280824             -1.047340   \n723 -0.455379 -1.128451 -0.835061             -1.047340   \n724 -0.876171 -0.752501 -0.522117             -0.397652   \n725  0.112004  2.050172  0.864138             -1.047340   \n726 -0.934570 -0.086793 -0.238156             -0.915569   \n727 -1.464136  0.127812  0.839491             -0.506735   \n728 -1.821786 -1.895845 -0.558090              0.677827   \n729 -0.949815 -1.491132 -1.004400             -0.828027   \n730 -0.077844 -1.294765  0.097017             -0.685715   \n731 -0.566455 -1.186169 -0.864787              1.175348   \n..        ...       ...       ...                   ...   \n906 -0.862845 -0.241745 -1.009248              0.659952   \n907 -0.199456 -1.152468 -0.906037              1.433619   \n908  0.040558  0.580464 -0.146242             -0.644923   \n909 -1.146815 -0.161945 -0.306279             -0.669673   \n910 -0.615574 -0.627959 -0.939756              2.104848   \n911 -1.477114  0.027676 -0.348414             -0.697860   \n912 -1.593911  0.035424  0.006951             -0.084840   \n913 -1.506488  0.018379 -0.547278              0.634056   \n914 -1.419064 -1.521412 -0.873174             -0.681819   \n915 -1.424785 -1.457915 -1.353459             -0.638277   \n916 -1.451577 -0.700825 -0.704406             -0.330048   \n917 -1.396737 -0.402861 -0.620536              0.376015   \n918 -0.933733  0.056084 -0.687556             -1.047340   \n919 -3.662360  0.497371 -0.076749             -1.047340   \n920 -1.976539  1.726711  0.507091             -0.193923   \n921 -0.749536 -1.059304 -1.115712             -1.047340   \n922 -2.872825 -1.227232 -0.500208              0.813952   \n923 -2.203855  0.273467 -0.371663             -1.047340   \n924 -1.096998  0.402851  0.829221             -1.047340   \n925 -1.598516 -0.717637 -0.631320             -1.047340   \n926 -1.866160  0.582012  0.313505             -1.047340   \n927 -0.941100  0.115224 -1.013870             -1.047340   \n928 -1.824018  0.272693 -0.572439             -1.047340   \n929 -2.252416 -1.253960 -0.730252             -1.047340   \n930 -0.865496 -1.376178 -1.365441             -1.047340   \n931 -1.404272  1.455354  1.129614             -0.424235   \n932 -1.211144 -1.135618 -1.624411             -0.959798   \n933 -1.511860 -0.886341 -0.727855             -0.658673   \n934 -1.334431  0.010438  0.094416             -1.047340   \n935 -1.157002 -1.705836 -1.334460             -1.047340   \n\n     reanalysis_air_temp_k  reanalysis_avg_temp_k  \\\n702               1.658885               1.502092   \n703               1.077751               0.804786   \n704               0.997884               0.657388   \n705               0.275933              -0.056925   \n706               0.085725              -0.249676   \n707               0.301154               0.016774   \n708               0.350545              -0.028579   \n709              -0.104483              -0.487781   \n710               0.254916              -0.034249   \n711              -0.258437              -0.672028   \n712              -0.771789              -1.309808   \n713              -0.660396              -0.992336   \n714              -0.712940              -0.884622   \n715              -0.836943              -1.287132   \n716              -0.534291              -0.771239   \n717              -0.273674              -0.567149   \n718              -0.300997              -0.516127   \n719              -0.855859              -1.003674   \n720              -0.642531              -0.924306   \n721              -0.877927              -1.264455   \n722              -0.326218              -0.629510   \n723              -0.277878              -0.493450   \n724              -1.326651              -1.825701   \n725              -0.690872              -0.901629   \n726              -0.002549              -0.249676   \n727              -0.027770              -0.278022   \n728               0.269628               0.084804   \n729               0.288544              -0.011572   \n730               0.049995              -0.244007   \n731              -0.423949              -0.878953   \n..                     ...                    ...   \n906               1.710378               1.434062   \n907               1.539085               1.366033   \n908               2.056116               1.915941   \n909               1.698818               1.485085   \n910               0.559670               0.254878   \n911               1.402471               1.161943   \n912               0.909611               0.765102   \n913               0.826592               0.526998   \n914               0.288544              -0.022910   \n915               0.237051              -0.158970   \n916               0.148777              -0.198654   \n917              -0.021464              -0.368729   \n918              -0.069805              -0.374398   \n919              -0.484900              -0.873283   \n920              -0.408186              -0.788246   \n921              -0.787552              -1.071704   \n922              -0.536393              -0.867614   \n923              -0.497511              -0.827930   \n924              -1.074441              -1.542244   \n925              -0.632023              -1.043358   \n926              -0.581581              -0.924306   \n927              -0.685617              -0.912967   \n928              -0.602598              -0.963990   \n929              -0.810671              -1.253117   \n930              -0.951489              -1.417522   \n931              -1.279362              -1.797356   \n932              -0.453374              -0.788246   \n933              -0.910505              -1.315477   \n934              -0.785450              -1.145403   \n935              -0.018312              -0.419751   \n\n     reanalysis_dew_point_temp_k  reanalysis_max_air_temp_k  \\\n702                     0.922279                  -0.006237   \n703                     0.495234                  -0.222818   \n704                     0.931644                  -0.532220   \n705                     1.040279                  -0.841622   \n706                     0.337901                  -0.748801   \n707                    -0.285810                  -0.841622   \n708                    -0.139716                  -0.594100   \n709                    -0.165001                  -0.903502   \n710                    -0.143462                  -0.810682   \n711                    -1.012537                  -1.042733   \n712                    -1.881612                  -1.274785   \n713                    -0.884236                  -1.212904   \n714                    -1.035013                  -0.996323   \n715                    -1.551026                  -1.058203   \n716                    -1.119298                  -0.810682   \n717                    -1.131473                  -1.120084   \n718                    -1.294424                  -1.027263   \n719                    -1.065917                  -0.872562   \n720                    -1.048124                  -0.686921   \n721                    -2.005230                  -1.058203   \n722                    -0.651046                  -0.965383   \n723                    -0.841157                  -1.027263   \n724                    -1.286932                  -1.089143   \n725                    -1.259774                  -0.779742   \n726                    -0.597666                  -0.717861   \n727                    -1.351551                  -0.841622   \n728                    -0.128478                  -0.625041   \n729                    -0.566761                  -0.717861   \n730                     0.082236                  -0.655981   \n731                    -0.059176                  -1.058203   \n..                           ...                        ...   \n906                     0.916660                   0.024703   \n907                     1.043088                  -0.160938   \n908                     0.727487                   0.055643   \n909                     0.412821                  -0.191878   \n910                     0.508345                  -0.408459   \n911                     0.389409                  -0.253758   \n912                     0.258298                  -0.099057   \n913                     0.333219                  -0.501280   \n914                    -1.241044                  -0.841622   \n915                    -0.626697                  -0.810682   \n916                    -0.308286                  -0.779742   \n917                    -0.175303                  -0.934442   \n918                    -0.862696                  -0.810682   \n919                    -1.134282                  -1.089143   \n920                    -1.080902                  -1.058203   \n921                    -1.458312                  -1.027263   \n922                    -1.054679                  -1.212904   \n923                    -1.499519                  -0.841622   \n924                    -2.430403                  -1.274785   \n925                    -1.509820                  -1.274785   \n926                    -1.383392                  -1.181964   \n927                    -1.797327                  -1.089143   \n928                    -1.436773                  -1.089143   \n929                    -1.997738                  -1.212904   \n930                    -1.992119                  -1.120084   \n931                    -2.064230                  -1.151024   \n932                    -1.316900                  -1.120084   \n933                    -1.646550                  -1.243844   \n934                    -1.947167                  -1.151024   \n935                    -0.632316                  -0.779742   \n\n     reanalysis_min_air_temp_k  reanalysis_precip_amt_kg_per_m2  ...    m__3  \\\n702                   1.126431                        -0.532252  ...       0   \n703                   0.696430                        -0.097593  ...       0   \n704                   1.009158                         1.200381  ...       0   \n705                   0.774612                         3.270673  ...       0   \n706                   0.579157                        -0.032925  ...       0   \n707                   0.657339                        -0.689765  ...       0   \n708                   0.696430                         0.241913  ...       0   \n709                   0.461884                        -0.643111  ...       0   \n710                   0.579157                        -0.715632  ...       0   \n711                   0.364156                        -0.582139  ...       0   \n712                   0.149155                        -0.448646  ...       0   \n713                   0.305520                        -0.460194  ...       0   \n714                  -0.007209                        -0.901321  ...       0   \n715                   0.031882                        -0.880534  ...       0   \n716                  -0.046300                        -0.838962  ...       0   \n717                   0.618248                        -0.856515  ...       0   \n718                   0.579157                        -0.887463  ...       0   \n719                   0.031882                        -0.466661  ...       0   \n720                   0.110064                        -0.552577  ...       1   \n721                   0.110064                        -0.781916  ...       1   \n722                   0.110064                        -0.607775  ...       1   \n723                   0.618248                        -0.670364  ...       1   \n724                  -0.554483                         0.946330  ...       0   \n725                  -0.046300                        -0.672674  ...       0   \n726                   0.344611                        -0.866677  ...       0   \n727                   0.579157                        -0.608006  ...       0   \n728                   0.696430                        -0.652812  ...       0   \n729                   0.618248                        -0.451880  ...       0   \n730                   0.579157                        -0.074497  ...       0   \n731                   0.305520                         0.897829  ...       0   \n..                         ...                              ...  ...     ...   \n906                   1.321886                        -0.520243  ...       0   \n907                   1.400068                        -0.419084  ...       0   \n908                   1.439159                        -0.735032  ...       0   \n909                   1.282795                        -0.565048  ...       0   \n910                   0.774612                         3.106694  ...       0   \n911                   0.970067                        -0.688841  ...       0   \n912                   0.696430                        -0.432479  ...       0   \n913                   0.852794                         0.075624  ...       0   \n914                   0.618248                        -0.700851  ...       0   \n915                   0.540066                         0.643777  ...       0   \n916                   0.500975                         0.614214  ...       0   \n917                   0.188246                         0.240296  ...       0   \n918                  -0.085391                        -0.872220  ...       0   \n919                   0.344611                        -0.280279  ...       0   \n920                   0.227337                        -0.056944  ...       0   \n921                  -0.163573                        -0.832034  ...       0   \n922                   0.266428                        -0.299449  ...       0   \n923                  -0.554483                        -0.864368  ...       0   \n924                  -0.124482                        -0.471742  ...       0   \n925                   0.305520                        -0.691843  ...       0   \n926                   0.266428                        -0.735032  ...       0   \n927                   0.227337                        -0.777528  ...       0   \n928                  -0.085391                        -0.776604  ...       1   \n929                   0.149155                        -0.710089  ...       1   \n930                  -0.319937                        -0.905940  ...       1   \n931                  -0.515392                        -0.752354  ...       1   \n932                   0.305520                        -0.841965  ...       0   \n933                   0.031882                        -0.118379  ...       0   \n934                   0.070973                        -0.815405  ...       0   \n935                   0.070973                        -0.876608  ...       0   \n\n     m__4  m__5  m__6  m__7  m__8  m__9  m__10  m__11  m__12  \n702     0     0     0     0     0     0      1      0      0  \n703     0     0     0     0     0     0      0      1      0  \n704     0     0     0     0     0     0      0      1      0  \n705     0     0     0     0     0     0      0      1      0  \n706     0     0     0     0     0     0      0      1      0  \n707     0     0     0     0     0     0      0      0      1  \n708     0     0     0     0     0     0      0      0      1  \n709     0     0     0     0     0     0      0      0      1  \n710     0     0     0     0     0     0      0      0      1  \n711     0     0     0     0     0     0      0      0      0  \n712     0     0     0     0     0     0      0      0      0  \n713     0     0     0     0     0     0      0      0      0  \n714     0     0     0     0     0     0      0      0      0  \n715     0     0     0     0     0     0      0      0      0  \n716     0     0     0     0     0     0      0      0      0  \n717     0     0     0     0     0     0      0      0      0  \n718     0     0     0     0     0     0      0      0      0  \n719     0     0     0     0     0     0      0      0      0  \n720     0     0     0     0     0     0      0      0      0  \n721     0     0     0     0     0     0      0      0      0  \n722     0     0     0     0     0     0      0      0      0  \n723     0     0     0     0     0     0      0      0      0  \n724     1     0     0     0     0     0      0      0      0  \n725     1     0     0     0     0     0      0      0      0  \n726     1     0     0     0     0     0      0      0      0  \n727     1     0     0     0     0     0      0      0      0  \n728     1     0     0     0     0     0      0      0      0  \n729     0     1     0     0     0     0      0      0      0  \n730     0     1     0     0     0     0      0      0      0  \n731     0     1     0     0     0     0      0      0      0  \n..    ...   ...   ...   ...   ...   ...    ...    ...    ...  \n906     0     0     0     0     0     0      1      0      0  \n907     0     0     0     0     0     0      1      0      0  \n908     0     0     0     0     0     0      1      0      0  \n909     0     0     0     0     0     0      1      0      0  \n910     0     0     0     0     0     0      1      0      0  \n911     0     0     0     0     0     0      0      1      0  \n912     0     0     0     0     0     0      0      1      0  \n913     0     0     0     0     0     0      0      1      0  \n914     0     0     0     0     0     0      0      1      0  \n915     0     0     0     0     0     0      0      0      1  \n916     0     0     0     0     0     0      0      0      1  \n917     0     0     0     0     0     0      0      0      1  \n918     0     0     0     0     0     0      0      0      1  \n919     0     0     0     0     0     0      0      0      0  \n920     0     0     0     0     0     0      0      0      0  \n921     0     0     0     0     0     0      0      0      0  \n922     0     0     0     0     0     0      0      0      0  \n923     0     0     0     0     0     0      0      0      0  \n924     0     0     0     0     0     0      0      0      0  \n925     0     0     0     0     0     0      0      0      0  \n926     0     0     0     0     0     0      0      0      0  \n927     0     0     0     0     0     0      0      0      0  \n928     0     0     0     0     0     0      0      0      0  \n929     0     0     0     0     0     0      0      0      0  \n930     0     0     0     0     0     0      0      0      0  \n931     0     0     0     0     0     0      0      0      0  \n932     1     0     0     0     0     0      0      0      0  \n933     1     0     0     0     0     0      0      0      0  \n934     1     0     0     0     0     0      0      0      0  \n935     1     0     0     0     0     0      0      0      0  \n\n[234 rows x 31 columns]\ncv error: 22.60683760683761\ntrain error: 25.48005698005698\n"
     ]
    }
   ],
   "source": [
    "train_sj = X_sj_train.copy()\n",
    "train_sj['y'] = Y_sj_train\n",
    "test_sj = X_sj_test.copy()\n",
    "\n",
    "model = smf.glm(formula=formula,\n",
    "                data=train_sj,\n",
    "                family=sm.families.NegativeBinomial())\n",
    "model = model.fit()\n",
    "\n",
    "predictions_sj = model.predict(test_sj).astype(int)\n",
    "#print (test_sj)\n",
    "\n",
    "print (\"cv error:\", mean_absolute_error(predictions_sj, Y_sj_test))\n",
    "\n",
    "pred_train_sj = model.predict(train_sj).astype(int)\n",
    "print (\"train error:\", mean_absolute_error(pred_train_sj, Y_sj_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The city of iquitos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cv error: 8.069230769230769\n",
      "train error: 5.441025641025641\n"
     ]
    }
   ],
   "source": [
    "train_iq = X_iq_train.copy()\n",
    "train_iq['y'] = Y_iq_train\n",
    "test_iq = X_iq_test.copy()\n",
    "\n",
    "model = smf.glm(formula=formula,\n",
    "                data=train_iq,\n",
    "                family=sm.families.NegativeBinomial())\n",
    "results = model.fit()\n",
    "\n",
    "predictions_iq = results.predict(test_iq).astype(int)\n",
    "print (\"cv error:\", mean_absolute_error(predictions_iq, Y_iq_test))\n",
    "\n",
    "pred_train_iq = results.predict(train_iq).astype(int)\n",
    "print (\"train error:\", mean_absolute_error(pred_train_iq, Y_iq_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### single evaluation metric\n",
    "Now the combined score for both the cities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cv error: 17.412087912087912\n",
      "train error: 18.32234432234432\n"
     ]
    }
   ],
   "source": [
    "pred = predictions_iq.append(predictions_sj)\n",
    "true = Y_iq_test.append(Y_sj_test)\n",
    "print (\"cv error:\", mean_absolute_error(pred, true))\n",
    "\n",
    "train_pred = pred_train_iq.append(pred_train_sj)\n",
    "train_true = Y_iq_train.append(Y_sj_train)\n",
    "print (\"train error:\", mean_absolute_error(train_pred, train_true))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conclusion : \n",
    " - Negative Binomial is by far the best model we have come across\n",
    " - We had high expectations from Time shifting the data but it actually decreases the performance\n",
    " \n",
    "#### Next Steps :\n",
    "At the next stage we shall dive into some more advanced statistical models"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}