_2_poisson_statsmodels_overdispersion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Another simple example with `statsmodels`, plus showing how to adjust for over-dispersion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example is taken from here: https://online.stat.psu.edu/stat504/node/169/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as tm\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "plt.style.use('seaborn')\n",
    "mpl.rcParams['font.family'] = 'serif'\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('./crab.txt', names=[\"Obs\",\"C\",\"S\",\"W\",\"Wt\",\"Sa\"], sep='\\s+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "endog_var = endog = df.Sa.values\n",
    "exog_vars = np.concatenate([np.ones(df.shape[0]).reshape(-1,1),\n",
    "                            df.W.values.reshape(-1,1)],\n",
    "                           axis=1)\n",
    "m = sm.GLM(endog = endog_var,\n",
    "           exog = exog_vars,\n",
    "           family = sm.families.Poisson()).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>y</td>        <th>  No. Observations:  </th>  <td>   173</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>   171</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>         <td>Poisson</td>     <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>          <td>log</td>       <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -461.59</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Sun, 12 Apr 2020</td> <th>  Deviance:          </th> <td>  567.88</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>12:11:58</td>     <th>  Pearson chi2:      </th>  <td>  544.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>5</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>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -3.3048</td> <td>    0.542</td> <td>   -6.095</td> <td> 0.000</td> <td>   -4.368</td> <td>   -2.242</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.1640</td> <td>    0.020</td> <td>    8.216</td> <td> 0.000</td> <td>    0.125</td> <td>    0.203</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   No. Observations:                  173\n",
       "Model:                            GLM   Df Residuals:                      171\n",
       "Model Family:                 Poisson   Df Model:                            1\n",
       "Link Function:                    log   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -461.59\n",
       "Date:                Sun, 12 Apr 2020   Deviance:                       567.88\n",
       "Time:                        12:11:58   Pearson chi2:                     544.\n",
       "No. Iterations:                     5                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -3.3048      0.542     -6.095      0.000      -4.368      -2.242\n",
       "x1             0.1640      0.020      8.216      0.000       0.125       0.203\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['pred'] = m.predict(exog_vars)\n",
    "df['lindpred'] = m.predict(exog_vars, linear=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note that `m.predict(exog_vars)` is the same as `np.exp(m.predict(exog_vars, linear=True))`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Obs</th>\n",
       "      <th>C</th>\n",
       "      <th>S</th>\n",
       "      <th>W</th>\n",
       "      <th>Wt</th>\n",
       "      <th>Sa</th>\n",
       "      <th>pred</th>\n",
       "      <th>lindpred</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>28.3</td>\n",
       "      <td>3.05</td>\n",
       "      <td>8</td>\n",
       "      <td>3.810341</td>\n",
       "      <td>1.337719</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>26.0</td>\n",
       "      <td>2.60</td>\n",
       "      <td>4</td>\n",
       "      <td>2.612781</td>\n",
       "      <td>0.960415</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>25.6</td>\n",
       "      <td>2.15</td>\n",
       "      <td>0</td>\n",
       "      <td>2.446839</td>\n",
       "      <td>0.894797</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>21.0</td>\n",
       "      <td>1.85</td>\n",
       "      <td>0</td>\n",
       "      <td>1.150492</td>\n",
       "      <td>0.140190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>29.0</td>\n",
       "      <td>3.00</td>\n",
       "      <td>1</td>\n",
       "      <td>4.274001</td>\n",
       "      <td>1.452550</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>168</th>\n",
       "      <td>169</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>28.3</td>\n",
       "      <td>3.20</td>\n",
       "      <td>0</td>\n",
       "      <td>3.810341</td>\n",
       "      <td>1.337719</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>169</th>\n",
       "      <td>170</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>26.5</td>\n",
       "      <td>2.35</td>\n",
       "      <td>4</td>\n",
       "      <td>2.836122</td>\n",
       "      <td>1.042438</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>170</th>\n",
       "      <td>171</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>26.5</td>\n",
       "      <td>2.75</td>\n",
       "      <td>7</td>\n",
       "      <td>2.836122</td>\n",
       "      <td>1.042438</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>171</th>\n",
       "      <td>172</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>26.1</td>\n",
       "      <td>2.75</td>\n",
       "      <td>3</td>\n",
       "      <td>2.655995</td>\n",
       "      <td>0.976820</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>172</th>\n",
       "      <td>173</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>24.5</td>\n",
       "      <td>2.00</td>\n",
       "      <td>0</td>\n",
       "      <td>2.042853</td>\n",
       "      <td>0.714347</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>173 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Obs  C  S     W    Wt  Sa      pred  lindpred\n",
       "0      1  2  3  28.3  3.05   8  3.810341  1.337719\n",
       "1      2  3  3  26.0  2.60   4  2.612781  0.960415\n",
       "2      3  3  3  25.6  2.15   0  2.446839  0.894797\n",
       "3      4  4  2  21.0  1.85   0  1.150492  0.140190\n",
       "4      5  2  3  29.0  3.00   1  4.274001  1.452550\n",
       "..   ... .. ..   ...   ...  ..       ...       ...\n",
       "168  169  2  3  28.3  3.20   0  3.810341  1.337719\n",
       "169  170  2  3  26.5  2.35   4  2.836122  1.042438\n",
       "170  171  2  3  26.5  2.75   7  2.836122  1.042438\n",
       "171  172  3  3  26.1  2.75   3  2.655995  0.976820\n",
       "172  173  2  2  24.5  2.00   0  2.042853  0.714347\n",
       "\n",
       "[173 rows x 8 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAFyCAYAAADLZb9oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxcddn//9fs2SYlaaZbCgUacgTKkrKJLIVawX2BGyt4gwveLrciKLjd4i0gqAjIov7Ur+jNjcptEQRXREsR2aEQSssyaRHapgudLE0mmeTMzJn5/ZE2NMkkk8nMnDOTvJ+PB4+S+ZxzPtfnOudMrpzVlU6nERERERF7uJ0OQERERGQmUfElIiIiYiMVXyIiIiI2UvElIiIiYiMVXyIiIiI2UvElIiIiYiOv0wFMRiQS1fMwgLq6Krq7Y06HUdaUw8JQHvOnHOZPOcyfcpi/TDkMhYKuiebRka8y4vV6nA6h7CmHhaE85k85zJ9ymD/lMH9TyaGKLxEREREbqfgSERERsZGKLxEREREbqfgSERERsZGKLxEREREbqfgSERERsZGKLxEREREbqfgSERGRnP3tb3/lxhu/x7XXXsO9997NN7/5teG2O++8Y8S0o3/OZt265/j4xz/Ms8+uzSvGxx57hHPOeS87dmzPazmFpuJLREREcvaXv/yBD37wPL70pa/xnve8nyuu+PZw2513/t+IaUf/nM1RRx3N4sWH5B3jW95yMnPnzst7OYVWFq8XEhGZycyERU+fyayaAAGfnkguI1VfcTmBP947tZndLupTY9/gZ77n/fRfcfW4sz3yyEO0t2/lrrtWcfzxJ9De3s6qVb/mrrv+yAMP/J2+vig///lPWbToQFwu94ifV6w4k1tv/QmWZeF2u6mqquLDH/4IADfddB3JZJIFCxrZtWvXmH6ffvpJrrzycs466xw+/vFPctttt/L88+u48spv8+1vX8mhhx7Grl27OPLIozjjjHeMmHfDhvVcf/13+Pznv8hBBx3M9753DU1NzVx44afo7+/j5ptvYP/9D2DXrl2cfPKpnHDCiTz00Bqee66VhoYGXn75Jb71re9OLc+jqPgSESlRVirFqjWbaG2L0NVrUl8boKU5xMrlTXjcOnEhzjn55GX85je/5oMfPJf58xcAsGrVrwF461vfxo9/fAsXXvip4en3/fnJJx/nhRfWc+ONPwLgc5/7JMcf/2YikQhbt27lhhtuAeDhhx8a0+9xx53Ae9/7ASoqKgHwer1cdtlXCQaDvPOd7+aUU07Dsiw+/OFzxhRfS5YcwSGHNANQV1fPKaecNnw68pe/vI2FC/fn/PM/hmkOct55/8aqVffy17/+mVNPPZ13vOPdrF+/rmD5U/ElIlKiVq3ZxOq17cM/d/aawz+ft6LZqbCkxPRfcfWER6kmEgoF6YpECxzRxF55ZSODg4P88pe3ATB37lx27+7m1VdfYf/99x+ebsGCxozzv//9Z3PxxZ/hnHM+xK5dr7NgQSPJZJLXXnuVtrYwfn+A3bu7c46ptnbWcEyLFzcRjfbyuc99gV/96jbuumsVJ554EkuWHInLNeE7sydFxZeISAkyExatbZGMba1tHZy9bLFOQUrJcrvdpNNpNm3ayCGHNI/4uampmQ0b1nP++R8F4JlnnqaxcX/i8QTPPPPGBfbbt2/LuOw5c+ayePEhfO9713Dmme8E4PHHH+Hpp5/kllt+AsDdd6/KOG9VVRWxWD8Ar7++c/jzpqZm6utnc845HwLg/vv/Qm3tLF54YQNf+crlJJNJPve5T3LKKadhGG/KLzmA54orrsh7IcUWi8WvcDqGUlBdHSAWizsdRllTDgtDecxfthx29Q7yp8c2Z2wz40lOPmI+1ZW+YoVXFrQd5m+qOXzssUdYvfp+YrEYjY0LWb36fh5//FHmzp3H4sVNbNvWztq1T7Fp00ZOOOHEET+fddY5bN26hSeeeIwXX9zAv/71CsuWnc4BByzihRfW889/PkhbW5jt29vp6IhwzDHH4vP5R/RfXz+bP/zhHj73uUsAqK2tZc2avxMOv8Qrr2zk2WfXst9+dXR0dLB69f2Ypslxx53ArFmzWLXq/9i9ezc7duygre1l3vSmw3jLW05m9er7aWsLs3btU/h8Pg47bAn33fcnnnnmaV566QV8Pj/vfOd78HhG/tGTKYfV1YErJ8qfK50ee6FdqYlEoqUfpA1CoSARmw8PTzfKYWEoj/nLlkMzYXH5z56gs9cc0za7toKr/+OEGX/kS9th/pTD/GXKYSgUnPDcpK7YFBEpQQGfh5bmUMa2luaGGV94iZQzXfMlIlKiVi5vAoau8eqODlIXrKCluWH4cxEpTyq+RERKlMft5rwVzZy9bLGe8yUyjaj4EhEpcQGfhzl1VU6HISIFomu+RERERGyk4ktERETERiq+RERERGyk4ktERETERiq+RERERGyk4ktERETERiq+RERERGyk4ktERETERiq+RERERGyk4ktEcmImLHZ09GMmLKdDEREpS0V7vZBhGPOAq4GjwuHwcaPaPgz8CgiGw+G+YsUgIoVjpVKsWrOJ1rYIXVGT+mCAluYQK5c34XHr7zgRkckq5jfmycDvAde+HxqGcShwWBH7FZEiWLVmE6vXttPZa5JOQ2evyeq17axas8np0EREykrRiq9wOHwXEN33M8MwqoAvA1cWq18RKTwzYdHaFsnY1trWoVOQIiI5KNppx3FcA1wVDofjhmFMeqa6uiq8Xk/xoiojoVDQ6RDKnnKYux0d/XRFzYxt3dFBPH4foYZqm6Mqf9oW86cc5k85zF+uObSt+DIMY3+gDli5T+H1RcMw/hIOh9dONG93d6zY4ZWFUChIJBLNPqGMSzmcGithUR8M0Nk7tgCrC1ZgxRPKa460LeZPOcyfcpi/TDnMVozZVnyFw+GtwEf3/mwYxneA7+uCe5HSF/B5aGkOsXpt+5i2luYGAj4dmRYRmayiXfNlGMYy4HxgvmEYlxuGUbnn85BhGJfvmezLhmE0FisGESmclcubWHHsQmbXVuB2wezaClYcu5CVy5ucDk1EpKy40um00zFkFYlESz9IG+jwcP6Uw/yZCQuP34cVT+iIVx60LeZPOcyfcpi/cU47usaZHNBDVkUkRwGfh/kN1Sq8RESmSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIiIiI1UfImIiIjYSMWXiIjkxUxY7OqOYSYsp0MRKQveYi7cMIx5wNXAUeFw+Lg9n30FmAfsAI4F/jscDr9czDhERKTwrFSKVWs20doWoavXpL42QEtziJXLm/C49be9yHiKvXecDPwecO3zWQ3wxXA4/D3gbuC6IscgIiJFsGrNJlavbaez1yQNdPaarF7bzqo1m5wOTaSkFbX4CofDdwHRUZ99IxwOp/fpv6+YMYiISOGZCYvWtkjGtta2Dp2CFJlAUU87TsQwDD/wEeCz2aatq6vC6/UUP6gyEAoFnQ6h7CmHhaE85q+cc7ijo5+uqJmxrTs6iMfvI9RQXfQ4yjmHpUI5zF+uOXSk+NpTeP0Y+Ho4HH4l2/Td3bHiB1UGQqEgkUg0+4QyLuWwMJTH/JV7Dq2ERX0wQGfv2AKsLliBFU8UfXzlnsNSoBzmL1MOsxVjtl8RaRhGFfBT4PvhcPgZwzDOtjsGERHJT8DnoaU5lLGtpbmBgE9nK0TGU+y7HZcB5wPzDcO4HLgB+BWwBDjIMAyAaoYuvBcRkTKycnkTMHSNV3d0kLpgBS3NDcOfi0hmRS2+wuHwQ8BDoz4+q5h9ioiIPTxuN+etaObsZYvp6TOZVRPQES+RSXDsgnsREZkeAj4Pc+qqnA5DpGzoKXgiIiIiNlLxJSIiImIjFV8iIiIiNlLxJSIiImIjFV8iIiIiNlLxJSIiImIjFV8iIiIiNlLxJSIiImIjFV8iIiIiNlLxJSIiImIjFV8iIiXGTFjs6o5hJiynQxGRItC7HUVESoSVSrFqzSZa2yJ09ZrU1wZoaQ6xcnkTHrf+VhaZLlR8iYiUiFVrNrF6bfvwz5295vDP561odiosESkw/SklIlICzIRFa1skY1trW4dOQYpMIyq+RERKQE+fSVevmbGtOzpIT1/mNhEpPyq+RERKwKyaAPW1gYxtdcEKZtVkbhOR8qPiS0SkBAR8HlqaQxnbWpobCPg8NkckIsWiC+5FRErEyuVNwNA1Xt3RQeqCFbQ0Nwx/LiLTg4ovEZES4XG7OW9FM2cvW0xPn8msmoCOeIlMQyq+RERKTMDnYU5dldNhiEiR6JovERERERup+BIRERGxkYovERERERup+BIRERGxkYovERERERup+BIRERGxkYovERERERup+BIRERGxkYovERERERup+JIZw0xY7Ojox0xYTodSssyExa7umHI0w2k7ECmuor5eyDCMecDVwFHhcPi4PZ9VANcD24BDgO+Gw+G2YsYhM5uVSrFqzSZa2yJ0RU3qgwFamkOsXN6Ex62/P2BUjnpN6muVo5lI24GIPYq9N50M/B5w7fPZJcCWcDj8HeBG4OdFjkFmuFVrNrF6bTudvSbpNHT2mqxe286qNZucDq1kjMgRytFMpe1AxB5FLb7C4fBdQHTUx+8CHt/Tvh44yjCM2mLGITOXmbBobYtkbGtt69BpFZQjGaLtQMQ+RT3tOI45jCzIevd81jveDHV1VXi9nmLHVRZCoaDTIZSVHR39dEXNjG3d0UE8fh+hhmqboyotU82RtsX8lVIOy3VfKaUclivlMH+55tCJ4msXsG+UtXs+G1d3d6yoAZWLUChIJDL6QKJMxEpY1AcDdPaO/aVSF6zAiidmfE6nkiNti/krtRyW475SajksR8ph/jLlMFsx5sQVlH8GTgQwDOMIYF04HB73qJdIPgI+Dy3NoYxtLc0NBHw6oqocCWg7ELFTse92XAacD8w3DONy4AbgZuD6PT83ARcWMwaRlcubgKHrVrqjg9QFK2hpbhj+XJQjGaLtQMQernQ67XQMWUUi0dIP0gY6PJwfM2Hh8fuw4gn9FT8OM2HR02cyqyYwYY60LeavlHM42e3AaaWcw3KhHOZvnNOOrnEmB/SQVZlBAj4P8xuqS/qXidMCPg9z6qqUoxlO24FIcan4EhEREbGRii8RERERG6n4EhEREbGRii8RERERG6n4EhEREbGRii8RERERG6n4EhEREbGRii8RERERG6n4EhEREbGRii8RERERG6n4EpGcmAmLHR39mAnL6VBkDzNhsas7NqPWibZDKWdepwMQkfJgpVKsWrOJ1rYIXVGT+mCAluYQK5c34XHr7zgnjFgnvSb1tdN/nWg7lOlAxZeITMqqNZtYvbZ9+OfOXnP45/NWNDsV1ow2E9fJTByzTD/6M0FEsjITFq1tkYxtrW0dOvXjgJm4TmbimGV6UvElIln19Jl09ZoZ27qjg/T0ZW6T4pmJ62QmjlkKx/3qv6i+5krc27c5HYqKLxHJblZNgPraQMa2umAFs2oyt0nxzMR1MhPHLAWQSFB58w3UL3szVTffgPe5VqcjUvElItkFfB5amkMZ21qaGwj4PDZHJDNxnczEMUt+vE8/Sd2KU6i55krSwVp6f/oL4u94l9Nh6YJ7EZmclcubgKFra7qjg9QFK2hpbhj+XOw3E9fJTByz5M7Vs5vqq6+k4vZf4EqnGTj/Y/R/4wrS+9U5HRoArnQ67XQMWUUi0dIP0gahUJBIJOp0GGVNOcyfmbDw+H1Y8YSONOShkNuimbDo6TOZVROYMetE22FhTLvvxHQa/x/vpea/voxn1+skjTcRve5mkm8+sWhdZsphKBR0TTSPTjuKSE4CPg/zG6r1C6+EBHwe5tRVzah1ou1QRnNv3ULtv3+QWZ/4CO6e3fR/7Rt0P/BIUQuvqdJpRxERESlfySSV/+/HVH/vGlyxGPFTltF33Y1YB5fuqWgVXyIiIlKWvM89S82lF+Nbv45UfT3Ra7+P+cFzwTXhWT/HqfgSERGRsuLqi1L13aupvPWnuFIpBleeR98V15CePdvp0CZFxZeIiIiUDf99f6bma5fh2b6N5MGL6bv+ZhInn+p0WDlR8SUiIiIlz71jOzVf+xKBv/yRtM9H/xe/TOySy6CiwunQcqbiS0REREqXZVFx261UX3MV7r4oiRNOJHr9zVjGm5yObMpUfImIiEhJ8mxYT/Cyz+N79hlSs/YjesMtDH74AnCX95OyVHyJiIhIaenvp/r671L5kx/isiwGz/o3+q76Luk5c5yOrCBUfImIiEjJ8K35O8EvfxHPls1YBxxI9Hs3kFj+NqfDKigVXyIiIuI41+uvU/ONr1Bx7+9IezzELvoC/Zd+BaqqnA6t4FR8iZQQu9/RNxPeCViOYxwdczmOodjMhMWOjn6shKWclLtUiopf/S/V3/om7p7dJI45luj1t2AdvsTpyIrGkeLLMIwvAQcCHcAhwIXhcHjAiVhESoGVSrFqzSZa2yJ09ZrU1wZoaQ6xcnkTniJcWGp3f04oxzFmirmqwkf/QJzuaLwsxlBsI3IUNakPKiflzBN+meCln8f31BOkaoJEv3M9gx+9EDzTu6C2vfgyDGMe8DWgIRwOpwzD+D1wFvBru2MRKRWr1mxi9dr24Z87e83hn89b0Vz2/TmhHMeYKebOXnPEz6U+hmIrx/UqGQwOUnXTdVT94CZciQTmu99H3zXXkpq/wOnIbOHEnwkxIA7U7vm5BnjBgThESoKZsGhti2Rsa23rwExYZd2fE8pxjBPFPFqpjqHYynG9yli+hx+ibtmbqf7+daTmzKXnl6vo/cUvZ0zhBQ4c+QqHw717TjuuMgxjB9AObJponrq6Krze6X0IcrJCoaDTIZS9Usvhjo5+uqJmxrbu6CAev49QQ3XJ9VdqedyX3Tmdqn1zOFHMo5XSGOxULuu13Ni2L3d0wKWXwu23Dz2n65JL8Fx1FbOCpftdMlm55tCJ045HA18ClobD4aRhGDcA/w18ebx5urtjdoVX0kKhIJFI1Okwylop5tBKWNQHAyNOL+1VF6zAiicKGnMh+ivFPO7L7pxOxegcThTzaKUyBruVw3otN7bsy+k0gVV3UHPF13F3dZE44ij6vn8LyaNaYBAYLO91limH2YoxJ047NgJd4XA4uefnHUD5vZhJpEACPg8tzaGMbS3NDQW/k8vu/pxQjmOcKObRSnUMxVaO63Wm87yykVlnv4faz38G16BJ31XfZvf9Dw4VXjOYE3c7/hV4554jXruBJcAlDsQhUjJWLm8Chq5b6Y4OUhesoKW5Yfjzcu/PCeU4xrEx773bMcHuPrMsxlBs5bheZyTTpOqHN1F10/W4TBPzjLfT953rSe1/gNORlQRXOp12OoasIpFo6Qdpg1I/1VMOSj2H5fKcr1LP475K9RlZE+VQz/nKzkxYePw+rHhCOclDMfZl3xOPUXPZxXjbwlhz59H37e8Rf/f7wOUqaD+lYpzTjhMOdsqnHQ3DuGiq84pIZgGfhzl1Vbb9MrG7PyeU4xhHx1yOYyi2gM/D/IZq5aSEuHZ3U/PFi9jvvW/Hs7GNgY9eSPejTxN/z/unbeE1VZM67WgYxqHADQw9ENUDuIA64AfFC01ERERKXjpN4J67qLn8q7g7IiQPPYzo9TeTPO4EpyMrWZO95utyhu5I/AxwFbAI+FCxghIREZHS5978GsEvfwH/gw+Qrqig7/IrGPjMReDzOR1aSZvsacfN4XB4LRANh8Obw+HwPwG9DkhERGQmSiSovOVG6k89Af+DDxBfdjpdDz3BwOe/qMJrEiZ75OtAwzCqgWrDMD4AdAKnFC8sERERKUXeZ54meOnFeF/cQKqhgej3f4B51jm6risHky2+fs9QsXUDcC9Qjx4PISIiMmO4enuo/vZVVPzPrbjSaQY+fAH9/30V6bp6p0MrOxMWX4ZhLAGagd+Hw+FBwzBOB8LAFuA+G+ITERERJ6XT+P/0B2q+/mU8O3eQPKSZvutvJnHiSU5HVrayXfN1LXAsgGEYdcA9wHqgF93pKCIiMq2527dSe8GHmHXh+bi7Oun/8n/RveZRFV55ynbasTMcDv8XgGEY5wNPhcPhy/f8/D/FDk5EREQcYFlU3voTqr9zNa5YP/GTTqHvupuwmg5xOrJpIVvxte8jW98O3DVOm4iIiEwD3uefo+bSi/GtayVVV0f0O9dhfujDuqC+gLIVXyHDMALAQuAMhp7ztVdj0aISERERe/X1UX3tNVT+7Me4UikGz/kQfVd+m3RDg9ORTTvZiq//BTYDQeC74XB4l2EYR+/5/PFiBycihVWK7wgsxZgKbSaMUcqb/2/3UfPVy/C0b8U68CCi191EYtnpTodVcKWyL05YfIXD4T8bhnEAUBEOh3v3fPYccJQdwYlIYVipFKvWbKK1LUJXr0l9bYCW5hArlzfhcU/5Fa/TLqZCm2iMIiVh+3ZqP/1ZAn+8l7TXS/8llxH7wpegstLpyAqq1L5vsj7nKxwOx4G4DbGISJGsWrOJ1Wvbh3/u7DWHfz5vRbNiKpKJxnjxucc4FZYIpFJU3PZz+PaVBHp7SRx3AtHrb8Y69DCnIyuKUvu+mR5/XorIuMyERWtbJGNba1sHZsKyOaLSjKnQso1xMJ60OSKRIZ4XX2C/d72N4FcvBZeL6HU3sfuP90/bwqsUv29UfIlMcz19Jl29Zsa27uggPX2Z24qpFGMqtGxj7B6nTaRoYjGqv/VN6lacgu+Zpxl831nw0ksMfuTjME1O9WdSit830zfbIgLArJoA9bWBjG11wQpm1WRuK6ZSjKnQso2xbpw2kYJLpQjcezf1p76Zqh/cSGr+Anru+C3Rn90G8+c7HV3RleL3jYovkWku4PPQ0hzK2NbS3ODIHT+lGFOhZRtjhX+yr9YVmTrfww+x39tPp/aTH8O9vZ3YZy6i659PEl9xptOh2aYUv2+094vMAHvvrmtt66A7OkhdsIKW5gZH77orxZgKbSaMUUqTZ8N6ar713/gffACAwfefRf9Xv0Hq4MUOR+aMUtsXXel02pGOcxGJREs/SBuEQkEiEb1YIB8zPYeFesZNIfNYKs/dKaZMY5zp22IhKIdjubdspvq7VxO4+05c6TTxU5bR/40rSR69NOP0My2Hxfi+yZTDUCg44esAdORLZAYJ+DzMqatyOowRSjGmQpsJYxRnuV5/naof3kTl//wMVzxOYsmR9H/jShKnLddrgfZRKvuiii8REZEy5d7WTtUPb6Li17fjGhzEOmAR/V+9HPOsc6b1HYzlTsWXiIhImXH/6xWqfnAjFXf+H65EAmv/A4hd9AUGz/13COhO2lKn4ktERKRMeF5+iaqbbyBwz124UimSi5uIXXwp5tkfBJ/P6fBkklR8iYiIlDjv889RdeP1BP78BwCShy0h9oXLMN/9PvBMzxtVpjMVXyIiIiXK++QTVN10HYEH/g5AYukxxL7wZeJnvF0X0pcxFV8iIiKlJJ3G989/UHXT9fgffRiA+FtOJvaFL5E49TQVXdOAii8REZFSkE7j/9tfqbrpOnzPrAUgvnwF/Zd8ieSbT3Q4OCkkFV8iIiJOsiz8f/4D1Tdej/eF9QCY73wPsS9cRvKoFoeDk2JQ8SUiIuKERILA735L1c034N20kbTbzeBZ5xC7+FKsQw9zOjopIhVfMm048ZqafPscPX+hlzcTRWNx2nf1sXBODcEqf87zj85htuVNJue5rhetx2nONKn4za+p+sGNeLZsJu31MvDhC4hd9IUZ++7Fmcax4sswDAM4FxgAlgFXhMPhp5yKR8qXlUqxas0mWtsidPWa1NcGaGkOsXJ5E54iPeE53z5Hz18X9FNd6Sc2mCjI8uzIQamJJ5Ncc/uzbIv0kUqD2wWNoRq+fsFS/N7sX3WZcmgmLPoHkqQZu7zJ5HyiaSYbw0xbj9Nafz+Vv7qNyh/dgmfnDtKBAAMXfpLYZy8mtXB/p6MTGznyYm3DMDzAH4D3hMPhlGEY84FkOByOZJpeL9YeMtNegDpZd6xuY/Xa9jGfrzh2IeetaB7xWaFymEufucxf6OVNdv5cleK2+M1fPMXWXX1jPt9/Tg1Xfvz4rPNPdp3sXd5kcj7RNBefe8yYHNq9HstdKW6HmbiivVT84mdU/eSHuDs7SVdVM/CxTxD79OdIz53raGzlksNSNpUXazv1p9RxgAu4yDCMrwHvATocikXKmJmwaG3LWLPT2taBmbBKrs+J5i/08oqVg1ITjcXZFhlbeAFsi/QRjcUnnD+XdbIt0kdnz0DWnGdbL4Px5KRjmCnrcbpxdXVS9d2rqV+6hJprroSkRf+lX6Hz2Q30f/Nbjhde4hynTjsuAk4Ezg2Hwz2GYfwKiAO3ZZq4rq4Kr1fXPcBQhS1v2NHRT1fUzNjWHR3E4/cRaqge8Xm+OZxKn5Odv9DLm8z8U1VK2+L2jRFS4xwfT6UhGk9x8KLx481lnaTSsH23mTXnwITTdPeazN8nh06tx3JXStvhsJ074YYb4Mc/hv5+CIXgO9/B/Z//SXVtLaW2Fksyh2Um1xw6VXz1Ai+Hw+GePT8/ApzGOMVXd3fMnqhKnA4Pj2UlLOqDATp7x/7SqgtWYMUTI3JWiBzm2mcu8xd6eZOZfypKbVsM+t24XWQswNyuofZCrRO3CxbsF8iac2DCaepqAyNicmI9lrtS2w7d7Vup+uFNVPz6dlymiTVvPgNf+wYD//5RqKoCEyiheKH0cliOxjntOOE8Tp12fBKYvefaLxg6EtbmUCxSxgI+Dy3NoYxtLc0NRblTLN8+J5q/0MsrVg5KTbDKT2OoJmNbYyj7XY+5rJPGUA2zZ7VqfYMAACAASURBVFVmzXm29VLhH/m3r9Zj+fL8axM1l3yW+uOPovIXPyM1dz7R62+m6+nnGfjkfw4VXiL7cOTIVzgc7jIM4yvATYZhRIAQcJUTsUj523vnWGtbB93RQeqCFbQ0N4x7R1kp9Dl6/v1qAlRX+ogNJuiOmnkvz44clJqvX7B03LsdJ2NsDse/2zHz9GNznut60XosL56XXqTq5usJ3Ps7XKkUyUOaiV18KeZZ58Ak7rCVmcuRux1zpbsdh+jw8MQm82ykQudwpj7nq5S3xXJ5ztdEOdRzvibHqe3Q+9yzVN14PYH7/gRA8vAj6P/CZcTf9V7wlNf6KuV9uVxM5W5HleYybQR8HubU2Xt4P98+R89f6OXNRMEqP4ceWD/l+UfnMNvyJpPzXNeL1mNp8j7xONU3XYd/zWoAEsccR+yLXyK+4ky97FpyouJLRERkPOk0vocepOqm6/E/9ggA8ZNPJfaFL5E4+VQVXTIlKr5ERERGS6fx/+2vVN34PXzPPgOAueIMYpd8ieTxJzgcnJQ7FV8iIiJ7WRaBP95L1U034H1xAwDmu99H7JJLSR55tMPByXSh4ktERCSRIHD3nVTdfAPeVzaRdrsZ/LeVxC6+FMt4k9PRyTSj4ktERGYsV1cnFb+6ncr/+Rmebe2kfT4Gzv8osc9dQuqgg50OT6YpFV8iIjLjeF7YQOXPf0rFXatwDQ6Srqom9olPMfDZi0k1LnQ6PJnmVHyJiMjMkEziv/8+Km/9Cf5HHwbAWnQgA5/4FIPn/jvp2lkOBygzhYovERGZ1lydnVTc8Usqb7sVz9YtAMSXnc7Af3ya+FvPKLsHo0r5U/ElIiLTTzqNd+1TVP7PrQT+eC8u0yRdVcXARy9k4MJP6SJ6cZSKLxERmT76+qj43W+p/J9b8b6wHoDk4iYGP3ohgyvPI71fncMBiqj4EslLtnfw5fqOvkK/6zHf+MebZ0dHP1bCciTGYr//Mt91Zofp+O7HqeR93+3QE36ZyttuJXDnb3BHe0l7PJjvfh8DH72QxCnL9CR6KSkqvkSmwEqlWLVmE61tEbp6TeprA7Q0h1i5vAmP2521Pdvy6oJ+qiv9xAYTk5q/0PFnnSdqUh8MUFXho38gTnc0XvQY881JsddZocc/lTGUo3zy3ru7j7e2P8v7XrifhS89O9Q+bz79n/4sg//+EVLzF9g9HJFJ8VxxxRVOx5BVLBa/wukYSkF1dYBYLO50GGWtUDn8zQMbWb22nQHTAmDAtPjX9l4GzCRHHDw7a3vW5cUtevvjk56/0PFPdp7e/jgDcZtizDMnRV9nOcYzlW0x3z5L0VTy/vwDz/D2x+7mi/fdxKnrVlPbsYMthx8P115L3/duJHHKMtLBoN1DKUv6vZK/TDmsrg5cOdE85fmnkoiDzIRFa1skY1trWwfRWHzCdjNhTXp5k5k/V9niz7T8UopxKv3Zuc4KMf5MnOiz2HIaUzyO+3d389bL/4Of/fzTrHzyt/isBL9f+h4+/dEf8d/nXEX0zHeBz2dT9CJTp9OOIjnq6TPp6jUztnVHB2nf1Tdhe0+fyZy6qkktbzLz5ypb/JmWX0oxTqU/O9dZIcafiRN9FttkxjS/o52KX/0vFXfegbuzk9nAiwsO5f4j3sajzSdh+gIAuMs0BzIzqfgSydGsmgD1tQE6M/zSqAtWsHBOzYTts2oCk17eZObPVbb4My2/lGKcSn92rrNCjD8TJ/ostvHGFEiYnLn1KZrO/y6Bp54AIDV7NtFPfY5ve49gQ2DumGWVaw5kZtJpR5EcBXweWppDGdtamhsIVvknbB99J9dEy5vM/LnKFn+m5ZdSjFPpz851VojxZ+JEn8U2ekwH7/oXn3ngJ/zvTz/Gf9x7A4GnniC+7HR6bv1fOp97mcFvfZt5Jx2TcVnlmgOZmXTBfRnRhZH5K1QODzuwjgEzSU9fHDOepL62gpOOmMfK5U24Xa6s7dmWVxcM0LBfJT6PCzNuZZ2/0PFPbp4ADbMq8bpdmInix5hvToq9znKNZyrbYr59lqLDa9M0/f13nPf7Wzj3sf/jkNc3Ea/dj/hn/pPoD37C4H98Bss4FLxDJ2qmYw6cpN8r+ZvKBfeudDpd1KAKIRKJln6QNgiFgkQiUafDKGuFzuFMfc6Xx+/Diif0nK8pTL9XPtti2T/nKx7Hv/pvVKy6A//q+3ElEqQ9HqKnvY34+R8hfcaZw8XWeEZvhzI1+r2Sv0w5DIWCE/4loOKrjGgnyZ9yWBjKY/5mXA7TabzPP0dg1R1U3HMX7s5OAJKHLWFw5XkMnv1B0nPm5LTIGZfDIlAO8zeV4ksX3IuISNG4d+4gcNedVNx5B96XXwIg1RAi9qnPMrjyPKwlRzgcoYj9VHyJiEhhDQwQuO9PVKy6A99DD+JKpUj7/ZjveT+DK88lfvoKPY9LZjQVXyIikr90Gu+TT1Bx5x0Efn8P7mgvAIljjmXwg+dhvv8s0nX1DgcpUhpUfImIyJR5NrYR+N1vqbj7TjyvvQqA1biQ/gs/ifnBc7GaDnE4QpHSo+JLRERy4t7WTuDe3xH43W/xrV8HQLqqisFzPsTgyvNInHwqlOmLvkXsoOJLRESycnV2EvjjvQTuuQv/448CkPZ6Mc94O+YH/g3zzHdCTY3DUYqUBxVfIiKSWV8fgb/+eajgevABXMkkaZeL+FtOxjzrHMx3v5d0/WynoxQpOyq+RETkDYOD+NesJnDvXQTuvw/XwAAAiSOPHiq43n8WqQWNDgcpUt5UfImIzHR7C64/3IP//vtw9/cBkDx48VDBddY5unBepIBUfImIzEQDA/gffIDAH36H//6/Dhdc1gEHEvvYJzDf9wGSRx4Nel+iSMGp+JrBiv2OvEL3V46isTjtu/pYOKeGYJU/6/R2r5OpLM9MWOzo6MdKWAXpI9ccOa0QOS50DidtnIIrecAius79CNbZ/4Z76VLMZGpojMlUQfbtqb5DdKZ9X+RLOSsfjhVfhmFUAk8CfwuHw5c5FcdMZKVSrFqzida2CF29JvW1AVqaQ6xc3oRnEreH5zp/vv2Vo3gyyTW3P8u2SB+pNLhd0Biq4esXLMWf4YXBdq+TqSzv6EMaSAPrNnbQFTWpD+bXR645clohcjxiGQXI4aT09w8VXH+6d8wRrr6PfYI/zjuG+xOz6YrGqXtigOp1TxMbTBRk355Kzmbi90W+lLPy47niiisc6fiHP/zhtUAvELvooov+NtG0sVj8CluCKnHV1QFisXjey/nNAxtZvbadAdMCYMC0+Nf2XgbMJEccnP3OpVznz7e/QipUDrO56ra1bN3Vx943wqeB3v446zZ1cnrL2IuV7V4nU1neqzuivLojWrA+cs2R0wqRY7v2BVd3F4E/3EP1dd8h+KVLqLj7TrwvvUhq/gIGL/gYfddcS/83ruT25EJ+/4rJQHxPPHGL3v54wfbtqYw3l3ns2p9LXT7blXKYv0w5rK4OXDnRPI6UxIZhnA88CrzqRP8zmZmwaG2LZGxrbevATFgFnT/f/spRNBZnW6QvY9u2SB/RUTup3eskm4mWV6g+cs2R0wqR42LvC+4d26n4+f9j1tnvZfZhi6m96NME7vsT1v4H0H/JZXT//SG6nl5H/zeuJHlUC2YyNen1PJV9OxqL5zzemfh9kS/lrDzZfmzfMIzDgEPD4fB/GYZx5GTmqaurwuvV+WuAUCiY1/w7OvrpipoZ27qjg3j8PkIN1QWbP9/+iiHfHGazfWOEVDpzWyoN0XiKgxe9EYPd6ySbiZZXqD5yzZHTCpHjouwLbW1wzz1D/z355BufH388fOAD8IEP4DUMvMDoJeeynqeyb0fjqZzHO5UcFXt/LnWF2K5meg4LIdccOnFhxQeAQcMwvgqcDPgNw7gkHA7fNN4M3d0x24IrZaFQkEgkmtcyrIRFfTBAZ+/YnbUuWIEVT0zYR67z59tfoRUih9kE/W7cLjIWF27XUHshc1ToHE+0vEym0keuOXJaIXJckPWUSuFtfQb/3+4j8Jc/4Q2/DEDa4yFxyjLMd76b+DvePfI5XOMsM5f1PJV9O+h35zzeXHNkx/5c6vLdrpTD/GXKYbZizPbTjuFw+JpwOHxVOBz+LvAI8NREhZcUVsDnoaU5lLGtpbkh6x0yuc6fb3/lKFjlpzGU+TUrjaGxd/TZvU6ymWh5heoj1xw5rRA5nvIy+vrw//mP1Fz8n8xecgh173gr1Tdej2fza5hvfxe9t/yYzhc20XP3Hxm88FOTfgBqLut5Kvt2sMqf83hn4vdFvpSz8uTYBfeGYZwNnAMs+uEPfzhw0UUXbRhvWl1wP6RQF0YedmAdA2aSnr44ZjxJfW0FJx0xj5XLm3BP4pk+uc6fb3+FZNfFpScdMZd1mzrpi8VJM3Q0Z+GcoTv5Mt19ZPc6mcry3rJkLgctqKW3P1GQPnLNkdMKkePJLsO9dQuBu1ZRfe3VBL/yRSp+91t8G54nXV3D4PvPInbpV4ledxPmynOxlhwJlVUFGVNdMEDDfpX4PC7MuJX3vj2VnOUyjy4WH5LPtqkc5m8qF9y70ulxLrwoIZFItPSDtEGhDw/PxOd82X2Ifbo+58vj92HFE3rOVx7P+RqRQ8vC++xa/H+/n8D99+F96YXhaZOHH4F55tuJn/EOkkcvhSIUpqPHVOh9u1jP+dIps5GmkmflMH/jnHacsOpV8VVGtJPkTzksDOUxf6EK6Pnt7wn87T78D/wNd0cHAOlAgPjJpxI/4x3Ez3g7qcaFDkdaurQd5k85zN9Uiq/Se5KhiMh0lE7j2diGf83f8a9ZDY89wqz40KmKVGgOAx++YKjgOvU0qLb3DmARsZeKLxGRInH19uD750P4H3wA/4Or8bRvfaPx6KPpX34G8TPfQfKolqKcThSR0qTiS0SkUFIpvOvX4X/wAXxrVuN7+klc1tBDLlP77cfg+84ivnwFidOWM/tIg5hO94jMSCq+RETy4OrowP+PB/CvWY3/Hw+8ce2Wy0Vy6THET19B/PS3kmw5BkrwnZUiYj99E4iI5CKRwPvMWvwP/h3/mgfwPv8crj03Lllz5jL4oQ8TP/2txJedTrre3neXikh5UPElIjKRVArPSy/if/gf+B5+CN9jj+LuH3ovZdrnI/GWk4eObi1fgXX4ErD52XUiUn5UfImIjOLe/Br+hx/C9/A/8D/80PCpRIDk4ibMU5YRX/42EiefQrpG78UTkdyo+BKRGc/V0YH/kYfwPfwQ/n/+A8/m14bbrLnzGDznQ8RPWUbi1NMm/foeEZHxqPgSkZmnrw//E48OPQbi4YfwvrB+uClVOwvzHe8mfuoyEqechnVIs04likhBqfgSkWnP1RfF+9ST+B9/FN9jj+BtfQZXMgnseaL8KaftKbaWkTzyaN2VKCJFpW8YkRLi9PsvJ9N/NBZn+8YIQb+bYJV/zDxOv6sxGouz45VtLH5tA7XPPInviUfxrntu+HlbaY+HgcOPJHnqaaROO53EcSdAZeXw/GbCoqc7VrT3HBZCru9iLMUx2K0Y758UmSoVXyIlwEqlWLVmE61tEbp6TeprA7Q0h1i5vAmPDU8+n0z/8WSSa25/lm2RPlJpcAHVlV4CPs/wPGbCon8gSRpwu6AxVMPXL1iKv8hHklxdnbgefYTnfnEPB258jhN3vYqbocc/pL1eki3HYJ54EqurDuKv7ka2m56hMcZDrAwE8IyTg6oKH/0Dcbqj8azrxI51OLqPuqCf6ko/scFExj5zjcnp7bAYso1pOo5ZSp9erF1G9ALU/JVqDu9Y3cbqte1jPl9x7ELOW9FcEv1/8xdPsXVXX87L3n9ODVd+/Pi8Y9yXKxLB98Sj+B97BN9jj+J96YXhtoTHS3heMxsWHs6GhUvoO2opX//0sqxjHK99vOlHy3UdTmVbzDXGXGNyejvM1WRyONX1XqpjLrRS/U4sJ3qxtkgZMhMWrW2RjG2tbR2cvWxxUU+DTKb/eMJiWyT3wgtgW6SPaCw+9VOQ6TSeTRvxPf0k3qefxPfUE3g3tr3RXFHBwIkncw8LWb/wcNrmNRP3BYbb3T0WnT0DE47xPW85cNz2TNOPXid2rMOJ+sjU50RjcmoMdss2plxzJFIoKr5EHNbTZ9LVa2Zs644O0tNnMqeuytH+O3sGSU3x+HMqDe27+jj0wPrJzTAwgO+5Z4cLLd/ap3B3db2xvOoa4stOH3q46Yknk2xZyks7+vm/3zw3bv/hLbsnHGP7rr5x2zNNP3qd2LEOJ+ojU58TjcmpMdgt25hyzZFIoaj4EnHYrJoA9bUBOjP8EqgLVjCrJpBhLnv7rwx4cbuYUgHmdsHCOTXjt7++E+9Tewqtp5/A+/y64TsRAawDFjF42ltJHP9mEsedgHXY4eAZeTRi4RzXuPG5XWAcsN+EY1w4p2bc9kzTj14ndqzDifrI1OdEY3JqDHbLNqZccyRSKLqaUMRhAZ+HluZQxraW5oain/aYTP/BKj+NofELqIk0hva569Gy8GxYT8UvfkbwM5+g/tgjmX1EM7MuPJ+qn/4I7/PrSB55FLFPfZaen99O5/NhutauJ/qTnzP48f/AOuLIMYUXMGF8jaEaZs+qnHCMwSr/uO2Zph+9TuxYhxP1kanPicbk1Bjslm1MueZIpFB05EukBKxc3gQMXWfSHR2kLlhBS3PD8Oel0P/XL1g67t2O3VGTuuDYux2bauBri3qo/N63h67ZemYt7r43LkxN1dVhnvF2Ese/meRxJ5A4eumIxz7kYnR8+95tOZkxjm3fe7djgt19ZtZ1Ysc6HN3HfjUBqit9xAYTe9ZBtjE5Pwa75b7ey3/MUvp0t2MZ0V0p+Sv1HDr9rKHJPucrGk+Nfc6X30X1xpewnnwK6/EnqHv5efz/2jRi3uQhzSSOO4Hk3lOITYcU/Onx2Z4zluvznor1jKx8tkU952tILjnUc74yK/XvxHIwlbsdVXyVEe0k+VMOCyPUUEPn2vX4nl2Lt/UZfM+sxbt+HS7zjWtnUsFakkcvJbn0GBLHHU/i2ONJ1892MOrSom0xf8ph/pTD/OlREyJSFO7Xd+Jd14r3uVa861qh9Rlmd3QMt6c9HpKHLSG59FgSS48hufTYoXci6iGVIiJjqPgSkRFckQi+dc8OFVrPP4f3uVY8O3eMnGjRIgZPOpVkyzEklh5L8sijoEq35IuITIaKL5EZzNXZiXddK769R7Wefw7PtpFP+7bmzcd8+ztJHtVC8qijSRzZQsPhi4nqVIWIyJSo+BKZIVyRCN716/BueB7fntOHnq1bRkyTCs3BfNuZQ4XW0S0kj2ohNXeeQxGLiExPKr5EpptUCvdrr+Ld8DzeDeuH/l3/PJ7Xd46crKEB861v23NEa6jYSs2bX/C7D0VEZCQVXyLlzDTxhl8aKrLWr8O7YT2eFzaMeJYWgLWgEfPMd5A8/AiSS44cKrQaF6rQEhFxgIovkTLh6uzE+9ILeF9Yj3f90FEtT9vLI17Fk3a7sQ5pJn74ESSPOIrkkqFiKz1bj3gQESkVKr5ESs3gIN6NYTwvbMD70ot4X3oBz4sv4Nn1+ojJ0lVVQ6cLlxxJ8ogjhwqtQw+f8hPiRUTEHiq+RJySSuHeumWowHpxA569hdYrm3BZ1ohJrYX7Y77tTKxDDyd52OEkjzgK6+DFGd9zKCIipU3Fl0xbU3ldSFFeMZJO49q1C2/by3jaXsb74p4i66UXcff3jZg0FQwSO+oY0ocfTnrPkSzr0ENJz9rPvnjzlO31PrkqxTFK6dN2M5ZyUjocKb4Mw1gMXA08CywEOsPh8FVOxCLTj5VKsWrNJlrbInT1mtTXBmhpDk34otyJ5vFM9int6TTunTvwhF8eKrTCL+MNDxVc7t27R07q9WI1HUL8sMNJHno4CeNQ7tldwz873XRF40P97x9i5fGZ+y9IvAUWTybHfbG135v7V00pjlFKn7absabynSjF5dSRr3rgN+Fw+PcAhmG8aBjGn8Ph8DMOxSPTyKo1m1i99o0HhXb2msM/X3zuMTnPc96K5pETp1K4t7XvKbDCQ0ezwi/jaQvjjvaOmDTt8WAddDCJt5xC0jCwmt9E0jh06NU7gcDwdHesbmP1a5PsP9d4bXLN7c+yddcbR/JSadi6q49rbn+WKz9+fM7LK8UxSunTdjPWVL4TpbgcKb7C4fDToz5yA/1OxCLTi5mwaG2LZGxrbetgMJ4c8/l487hTFu1PrMcV30jFKxvfOG0YDuOKjdxc014v1uImEs3LSTYbWG86lGTzm4auy9qnyJpKzGcvWzziFEGu09shGouzLdKXsW1bpI9oLJ7TKchSHKOUPm03Y03lO1GKz/FrvgzD+ABwfzgcfnm8aerqqvB6Z9YOM55QKOh0CCVtR0c/XVEzY1t3dJDuXpP5o3K4c+NWQuF1HN21jcau7TR2b6Oxezvzd+/AbyVGLsTvB8OAww4b+u/ww+Gww3A1NeH1+fACE5daucfs8fsINVRPefpi2Xdb3L4xQiqdebpUGqLxFAcvmvy2WypjLDbtz/nbN4czZbvJxVS+EyV3ue7LjhZfhmGcDpwOXDLRdN3dMXsCKnGhUJCI3qc3ISthUR8M0Nk78svGayV4U7Kb2Q/+lb7nX8Dzyka8mzbieWUj87q6uG7UcmL+SjY3LCIyZ38Of8db4LDDsQwDa9FBkOn6pd2DwGBBYwaoC1ZgxRMj1nuu0xfD6G0x6HfjdpGxAHO7htpziakUxlhs2p/zNzqHM2G7yVW2nNTVBmZcTgot076crRhzrPgyDONdwCnAxcB8wzAWhcPhx52KR6aBdJrKrgjvTm7m9efX0di1begoVtc25vbuwpNOwY9g78mvtMeDdcAiEsccxwu+Bh6L17KtrpFtdQvorq4Dl4sVxy7kkCJfJxLweWhpDo24JmOvluaGMadJcp3eDsEqP42hmhHXfO3VGMr9rsdSHKOUPm03Y2XLSYXfi0ov+zl1t+MxwCpgLfAgUA38CFDxJROzLNzb2vG8+i88r7069N8+/++K9fNvo2bprZrFzuYjqVu6hKqjj6Bn3gFYTYdgLTpw6DQisCCVwlqziR1tHfREB5kdrKClucG2u4H29tPa1kF3dJC6LP3nOr0dvn7B0nHvdpyKUhyjlD5tN2MpJ6XHlU6Pc6FGCYlEoqUfpA1mzGkK08SzZTOe1/6F59V/4d63wNqyGVciMWaWVHUNqQMPwjroYJJNh2AdvJjBAxfTOe8Ago1zh//izZZDp5+Dk2v/TsU7UR71nK/JmTH7cxFNlMPput3kI1NOtB3mb5zTjhO+ONfxC+5lZnL1RXG/us+Rq81v/L97WzuuDH8UpGbPJnnkUVgHHox14EFD/x20GOvAg0iHQmNeEu0CGnKMK+DzMKeuauoDy1Ou/TsdbybBKj+HHlhfsOWV4hil9Gm7GUs5KR0qvqQ44nHc7VvxtG/Fs3UL7q2b8WzZgmfza0MFVkfmW5+t+QtInHjSnsLq4OGjWdaBB5GunWXzIERERApPxZdMjWni2bYV95YteNq3vlFcbd2Ce+sW3Dt3ZDx6lfZ4SO1/APEjjhwusKwDDx7694BFeim0iIhMeyq+JLPBwVHF1Zah67D2Flev7xy/uFrQSOLEk0jtfwDW/gdgHbBo+P9T8xeAz+fAgEREREqDiq+ZKhbDs30b7i2bx5wadG/dguf1nRlnS3s8pBoXknjLyW8UV/sfQOqARW8UV1N4j5+IiMhMod+S05Fp4t6+bai42ta+599tuLe349m2DfeObbi7uzPOmvZ6SS1YSPzkU4eKqdHF1bz5Kq5ERETyoN+i5SaRwL1zB+5t2/Bsbx/x71DB1Y67o2Pc2VPVNaQaG0keeTRW48LMxZVHt2WLiIgUi4qvUmFZuDo68Ly+A/frO3Hv3LnPv0Of8fpOGnZmvtYKIF1RgbWgkeSbDiO1oBGrsZHUgoWkGhux9vybrp015pEMIiIiYh8VX8WWqajauQP3668PF1XunTtxR3bhsqxxF5MOBGD+fBInnEhqQSOpxoVYe/7dW1yl6+tVWImIiJQ4FV9TtW9RtbeY2reo2nvkahJFVWrufJLHHEdq7jysefNIzd3nv3nzSc2dS3q/OkJzaunRk4hFRETKmoqv0SwLd0dkzxGpPIqqigpSc+eRPOY4rD0FVGrunn/nzd9TWM0jPWs/Ha0SERGZQVR8AaTTBD/9cXyPP4Z71+u4UqnxJ91bVB17PNaeAio1Z8+/+x6pKoOiqtDvPpuJ71Ir9HsMS62/Qsi2XTi93TgRn5mw2NHRj5WwCtJnoWPs7BkgvGU3xgH7MXuWHnwsUmgqvgDSady7dkEgQPK4EyY+UjUNLli3UilWrdlEa1uErl6T+toALc0hVi5vwuN2O768chBPJrnm9mfZFukjlQa3CxpDNXz9gqX4i/AoDrv7K4Rs24XT240T8Y1YZtSkPphfn4WOcSCe4Cs/fpy+geTwZzWVXq79zIlU+vVwZJFCcaXHuXOulEQi0dIP0gaFevv8HavbWL22fcznK45dyHkrmh1fXjEVKoff/MVTbN3VN+bz/efUcOXHj897+U73l81k8phtu3B6u3EivkL3WegYP3/zP0cUXnvVVHq55eJTc15esRVqf57JlMP8ZcphKBSc8CjN9DwsIeMyExatbZlfat3a1oGZGP86NjuWVw6isTjbImMLIYBtkT6isXhZ91cI2baLaCzu6HbjRHyF7rPQ+15nz0DGwgugbyBJZ89ATssTkfGp+JphevpMunrNjG3d0UF6+jK32bW8ctC+a+jUXyap9FB7OfdXCNm2i/ZdfY5uN07EV+g+C73vhbfszqtdRCZPxdcMM6smQH1tIGNbXbCCWTWZ2+xaXjlYOKcG9zgHlN2uofZy7q8Qsm0XC+fUOLrdOBFfofss9L5nHLBfXu0iMnkqaX4HrQAADiFJREFUvmaYgM9DS3MoY1tLc0POd0oVennlIFjlpzGUueBpDBX+LkS7+yuEbNtFsMrv6HbjRHyF7rPQ+97sWZXUVGa+eaOm0qu7HkUKyHPFFVc4HUNWsVj8CqdjKAXV1QFiBbi+57AD6xgwk/T0xTHjSeprKzjpiHmsXN6Eewp3chZ6ecVUqByedMRc1m3qpC8WJ80bR6C+fsHSotypZ3d/2Uwmj9m2C6e3GyfiK3SfhY7xtKUL+Odz24kn33jczt67HX0l+M7XQu3PM5lymL9MOayuDlw50Ty627GMFPqulJn4nK9C53CmPucrlzzqOV+Z+/T4fVjxhJ7zlQfdqZc/5TB/U7nbUcVXGdFOkj/lsDCUx/wph/lTDvOnHOZPj5oQERERKXEqvkRERERspOJLRERExEYqvkRERERspOJLRERExEYqvkRERERspOJLRERExEYqvkRERERspOJLRERExEYqvkRERERspOJrDzNhsas7hpmwnA4lIzNhsaOjv6Dx5TvmaCzOS691Ed3zQtHRy+vsGeCx9Tvo7BmYUjzZ4hu9/NHxZFr+vjnMtPzRn2VbZrYxZsvRZGLOZR1lWwej23d09PHXJzazo6MPgM07e7jrwU1s3tkzqf4yGb3MbHLNQa7TF3o7z2Yy/eW6P2dbZr7tImIvrxOdGoaxAjgL2AWkw+HwhG//LiYrlWLVmk20tkXo6jWprw3Q0hxi5fImPG7na9MR8UVN6oP5x5fvmOPJJNfc/izbIn2k0uACqiu9BHweunpN6oJ+emMJktYbr+SsqfRy7WdOpNLvyxpPXdBPdaWf2GAiY3wD8QRf+fHj9A0kh5fhAvb25nZBY6iGr1+wFL/XOyaHdTVjl3/0IQ2kgXUbO4ZjiCdTxAaTpNJjl5kphn3HmC1H9bUBzIRF/0CSdLaYJ7GORk8/q9pL36A1Yh14PS5qKjz09CeprfLQG7OGc3bnP14Zsby/PLkFF/D9z5/ErKpA1m0CoG8wzhd/8Ohwn3f+4xW8Hhffv+gkairGvgR8dI4mk4NccpZtO8omW3zZ1kGm/nLdn7MtM992EXGG7S/WNgyjCngeODwcDpuGYdwN/H/hcPiB8eYp5ou171jdxuq17WM+X3HsQs5b0VysbietGPHlu8xv/uIptu6a3FGNfdVUernl4lMnHc948X3+5n+OKHrGs/+cGq78+PGTXv5k7F3meDHsHeNUc5Qt5vHWUSHHuC8X8POvLh/zeaYXyX7yugdHFHt7eT0u/t+XTh/z+Xg5ynW95Tp9vtv53v5Gm8w6K9R63Tt9vu3lTi+Fzp9ymL9yebH2icDmcDhs7vn5UeBdDsSBmbBobYtkbGtt63D8EH0x4st3mdFYnG2R3IsKgL6B5JjTcxPFkym+HR19kyq8ALZF+ujsGZj08ie7zM07e8aNoW8gyeadPVPOUbaYM62jXHKYqzRM6hTkjo6+jIUXQNJKjzkFOdF2lOt6y3X6fLfzbZG+MacgJ7Nf5brvZZs+Govn1e7095vITObEacc5wL4lYu+ez8ZVV1eF1+speCA7OvrpipoZ27qjg3j8PkIN1QXvd7KKEV++y9y+MUIqj+OQ23ebvKnpjdU9UTyZ4tu0s3/SfaXSQ/1NdvmTXeaGzb0TTrNhc++Uc5Qt5kzrKJccTsWGzb0ce8TCMZ+HQsHh/394w+sTLmPTzn6OPHT+8M8TbUe5rrdcp893O0+lIRpPcfCiN8Y/mf0KKNh67Y4OEo2n8mp3+vutUPbdDmVqlMP85ZpDJ4qvXcC+Udbu+Wxc3d2xogRiJSzqgwE6e8d+QdUFK7DiCUcPxxYjvnyXGfS7cbuYcnGxYL/AiOVPFE+m+JrmTf6Xhds11N9klz/ZZS5ZVMvdE0yzZFEt90wxR9lizrSOcsnhVCxZVJvpkPqIz7Ktl6Z51SOmn2g7ynW95Tp9vtu52zXUPtl1sLc/oGDrte7/b+9eY+SqyziOf2dnZ7d772W33WWr20rhCRUkFPBCQ0JNkArhhZpIolETTQRMihIUG9oXlVTFlAiJa0J4oyY2EogSlFKa6osSScVSFWPAR0EDCC1I1W232+zszRfnLExnz9w6u+ecYX6fV7v/c+Y/zzznzPyf+Z/L9Cyjp62lruVJf74tBh0yq59yWL8Shx3LPiaJw46HgREzmz+LdzOwL4E4aM9luezCgchll13YT3tu8WfbarEU8dXbZ09nG8MD3TU/LwTnQ63q66g6nqj4hvq76e6o7jvD8EA3q/o6qu6/2j5HBvtKxtDd0crIYN8556hSzFHbqJYc1ioDjAz2VVxvqL+b1mz0KQ6t2QxD/Wfno9x+VOt2q3X9evfz4YFuejrPvoCgmvdVre+9Suv3dLbVtTzpzzeRZpbdtWtXrE+4bdu2qdHR0ReAO0ZHRz8MHHP3H5V7zMREftdSxbNx3QrOTE4zNp5nMj/Nyt5lbL5kkJs+uoGWTNnz5WKxFPHV2+fmS9bw3IsnGJ/IM0cwQHd3tNLTmWMyP8OKnjamZmbPmjWYvxIwl134gV8cz4qedvqXd5DLZpjMzyyI75pN5/HUn14nPz37dh+FUbdkYO3q4Kq0bEtLVf1fdfEa1p/Xy8nTU+E6bWSzGaZnZt++sq6wz6gYCl9jpRyt7G0n25JhKnx8pZgrbaPi9fu6ckzPzp21DVqzGXo7c+SnZ+ntbCU/Nbugn0LzVzsuyy0sNLu62pkoOu9py+XDHDzy6oLn/P62zbRFnDZQnKPKOagtZ5X2o0oqxVesmm1W73YtXr/e5Y0uaj+U2iiH9YvKYVdXe9m7OMR+teO5WMqrHedNTs0wNj5JX3d7Kr8RTk7NkG3LMZOfWrT46n3Npyby/OvNcdauDmYCivs7MXYGf+V/2HuXL5jxqiaeSvEV918cT1T/hTmM6r+4rVKflV5jpRxVE3Mt26jSNihefuytcZ578QSXbljFUH83Lx8f48gL/+bKiwbKzniVO1RR3Gclteag1vUXez+vpJrnq/X9XKnPepc3Kh0yq59yWL9zudpRxVcD0Zukfsrh4lAe66cc1k85rJ9yWL9GudWEiIiISNNS8SUiIiISIxVfIiIiIjFS8SUiIiISIxVfIiIiIjFS8SUiIiISIxVfIiIiIjFS8SUiIiISIxVfIiIiIjFqiDvci4iIiLxbaOZLREREJEYqvkRERERipOJLREREJEYqvkRERERipOJLREREJEYqvkRERERi1Jp0ABLNzM4HdgN/ANYCJ9z9bjO7D5gAxoFLga+5+/HkIk2vUjksWL4DuN3d+xMKMfXK7IdtwB0E++H7w/YdyUWaXmVyeDmwHXgW+BCwx90PJxdpeplZC/Ar4BmgDTgf+CLQAdwD/AO4ALjL3d9IKs40K5PD76AxpSqlcujuZ8LlVY8pKr7SayXwkLs/BmBmz5vZPuC0u+8M274J7AC2JRdmqkXm0N2Pmtk14XIpr9R+eD1wyN2fCts/kGCMaVcqh7uBB939UTP7BHA3cG2CcabdYXffDWBmjwGfBK4Gfu3uD5vZjcC9wOcSjDHtonKoMaU2UTncW+uYouIrpdz9SFFTCwVvkoK28fiiaiylcmhma4CbgO8BX4g9sAZSKofAZ4BXzGwTsAr4QdyxNYoyOXwDGAjbBoCjccbVSNx9lqBYxcxaCWYQnWDW69vhak8DP0kkwAZQKofuvrdgNY0pZZTK4bmMKTrnqwGE34oPuPtfC9qWAx8D9iQWWAOZzyHwN4Jp9ruSjajxFO2H64A5d78fOAQ8nGRsjaIohzuBz5vZvQQf2I8kGlwDMLPrgMeBx939WWA1cCpcfBJYEQ6KUkJEDufbNaZUqTCHBKcT1DymqPhKOTPbAmwBbi9o6wN+SHCs+T9JxdYoinK4CZgCbgZuBTrMbLuZXZBgiKkXsR+eJDjvAeC3wNVmlk0itkYRkcNfAt9w968DtwBPmFkmqfgagbsfcPetwHoz+wrwJtATLu4F/uvu04kF2AAicqgxpUaFOST4ElXzmKJvCClmZjcQnNPwVWDIzEaAvwP3A3e6+2tm9il3/3mScaZZcQ6BnLvfEi5bB3zJ3e9JLsL0K7Ef/gZ4H8GhnxHgJXefSS7KdCuRw/cAx8JVjgHtCYWXema2EVjv7vvCpn8S7H/7gI8ArwKbw/8lQqkcmlk/GlOqUiKHvecypuiHtVMqvBLqEMGVUABdBN9MbiMomue/nZxy9xvjjzD9SuXQ3X9sZhsIZhtuBb4L3Ofup5OJNL3K7IcHgW8BLwEXAaPu/vtEgky5MjkcAz4N/BnYCPzC3R9NJMiUC68Y3UNwiCdHsM/dBuQJzrN5meDKs+262jFamRw+gcaUqpTKobsfr3VMUfElIiIiEiOd8yUiIiISIxVfIiIiIjFS8SUiIiISIxVfIiIiIjFS8SUiIiISI93nS0Sagpk9QHAjxAfd/eaw7SHgL+6+O7wz+l5gGPiyuz+fXLQi8m6mmS8RaQrhjRBfB34KEN6R/wpga7h8muAmnXeq8BKRpaTiS0SayZPAx8O/Pwg8AFwR/q7dfNszUQ8UEVksKr5EpJnsJ5zpAq4DfgYcBa41sxYgo59JEpGlpuJLRJrJQeBiMxsEht39Nd6ZDbuSd34CSERkyaj4EpGm4e5jwO+AzwJvhc1PEsyCbQ3/FhFZUiq+RKTZ7Ad2AgfC/48CbcBV7n4ssahEpGmo+BKRZrMfyAJPA7j7LEEh9sckgxKR5pGZm5tLOgYRERGRpqGZLxEREZEYqfgSERERiZGKLxEREZEYqfgSERERiZGKLxEREZEYqfgSERERiZGKLxEREZEYqfgSERERidH/AS9geVGF1eavAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "aux = df.copy(deep=True).sort_values('pred').reset_index(drop=True)\n",
    "f, ax = plt.subplots(1,1, figsize=(10, 6))\n",
    "ax.scatter(df.W, df.Sa)\n",
    "ax.plot(aux.W, aux.pred, color='red', label='fitted values')\n",
    "ax.set_xlabel('W')\n",
    "ax.set_ylabel('Sa')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the context of Logistic regression (same applies to Poisson): https://online.stat.psu.edu/stat504/node/162/\n",
    "\n",
    "\"*When a logistic model fitted to n binomial proportions is satisfactory, the residual deviance has an approximate χ2 distribution with (n – p) degrees of freedom, where p is the number of unknown parameters in the fitted model. Since the expected value of a χ2 distribution is equal to its degree of freedom, it follows that the residual deviance for a well-fitting model should be approximately equal to its degrees of freedom. Equivalently, we may say that the mean deviance (deviance/df) should be close to one. Similarly, if the variance of the data is greater than that under under binomial sampling, the residual mean deviance is likely to be greater than 1.*\"\n",
    "\n",
    "The deviance is the sum of squared **deviance** residuals. The latter, in turn, are defined as (see https://bookdown.org/roback/bookdown-bysh/ch-poissonreg.html )\n",
    "\n",
    "$\\text{sign}(Y_i - \\hat{\\lambda}_i)\\sqrt{2\\left[Y_i\\log\\left(\\frac{Y_i}{\\hat{\\lambda}_i}\\right) - (Y_i - \\hat{\\lambda}_i)\\right]}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "567.8785724519456\n",
      "567.8785724519457\n",
      "567.8785724519457\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/pandas/core/series.py:679: RuntimeWarning: divide by zero encountered in log\n",
      "  result = getattr(ufunc, method)(*inputs, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "logdiff_y_yhat = np.log(endog_var / df.pred).values\n",
    "logdiff_y_yhat[logdiff_y_yhat==-np.inf] = -99999\n",
    "\n",
    "deviance_resid = np.sign(endog_var - df.pred) * np.sqrt(2 * (endog_var * logdiff_y_yhat - (endog_var - df.pred) ))\n",
    "\n",
    "print(m.deviance)\n",
    "print(np.sum((m.resid_deviance)**2))\n",
    "print(np.sum((deviance_resid)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.3209273242803836"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.deviance / m.df_resid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In contrast, the Pearson residual is defined as\n",
    "\n",
    "$\\frac{Y_i - \\hat{\\lambda}_i}{\\sqrt{\\hat{\\lambda}_i}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "544.1570112030288\n",
      "544.1570112030288\n",
      "544.1570112030288\n"
     ]
    }
   ],
   "source": [
    "pears_resid = (endog_var - df.pred) / np.sqrt(df.pred)\n",
    "print(m.pearson_chi2)\n",
    "print(np.sum((m.resid_pearson)**2))\n",
    "print(np.sum(pears_resid**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accounting for overdispersion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen above that $\\frac{\\text{deviance}}{\\text{residual degrees of freedom}}$ = 3.321. I.e. there is overdispersion. To correct for this, the quasi-Poisson model includes a dispersion parameter that scales the variance: $\\mathbb{V}(\\mu) = \\phi\\mu$. The estimator of $\\phi$ is precisely the Pearson $\\chi^2$ scaled by the residual degrees of freedom. ( https://stats.stackexchange.com/questions/314948/quasi-likelihood-quasi-poisson ). To account for over-dispersion, we have to multiply each standard error by $\\sqrt{\\hat{\\phi}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \"old\" standard errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.54224156, 0.01996535])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(np.diag(m.cov_params()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### over-dispersion adjusted standard errors\n",
    "\n",
    "Note that these are exactly the same as those by the `R` `glm()` with `family=quasipoisson(link=log)` (see `_2_overdispersion.R` in this repository)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.96729036 0.03561565]\n"
     ]
    }
   ],
   "source": [
    "# print(np.sqrt(m.deviance / m.df_resid) * np.sqrt(np.diag(m.cov_params())))\n",
    "print(np.sqrt(m.pearson_chi2 / m.df_resid) * np.sqrt(np.diag(m.cov_params())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Robust variance-covariance estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A better way than simply assuming over-dispersion is to estimate a robust covariance matrix. The discussion here is worth reading: https://github.com/statsmodels/statsmodels/issues/4942. For the different robust covariance matrices, see here: https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.RegressionResults.html#statsmodels.regression.linear_model.RegressionResults (scroll to bottom). Note that the robust standard errors that we get here are slightly **smaller** than the ones we get from the quasi-poisson model above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = sm.GLM(endog = endog_var,\n",
    "           exog = exog_vars,\n",
    "           family = sm.families.Poisson()).fit(cov_type='HC0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>y</td>        <th>  No. Observations:  </th>  <td>   173</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>   171</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>         <td>Poisson</td>     <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>          <td>log</td>       <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -461.59</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Sun, 12 Apr 2020</td> <th>  Deviance:          </th> <td>  567.88</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>12:11:59</td>     <th>  Pearson chi2:      </th>  <td>  544.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>5</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>        <td>HC0</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>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   -3.3048</td> <td>    0.840</td> <td>   -3.933</td> <td> 0.000</td> <td>   -4.952</td> <td>   -1.658</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.1640</td> <td>    0.030</td> <td>    5.381</td> <td> 0.000</td> <td>    0.104</td> <td>    0.224</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   No. Observations:                  173\n",
       "Model:                            GLM   Df Residuals:                      171\n",
       "Model Family:                 Poisson   Df Model:                            1\n",
       "Link Function:                    log   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -461.59\n",
       "Date:                Sun, 12 Apr 2020   Deviance:                       567.88\n",
       "Time:                        12:11:59   Pearson chi2:                     544.\n",
       "No. Iterations:                     5                                         \n",
       "Covariance Type:                  HC0                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         -3.3048      0.840     -3.933      0.000      -4.952      -1.658\n",
       "x1             0.1640      0.030      5.381      0.000       0.104       0.224\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}