chapter12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING *** OLE2 inconsistency: SSCS size is 0 but SSAT size is non-zero\n"
     ]
    }
   ],
   "source": [
    "nhefs_all = pd.read_excel('NHEFS.xls')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just a look at a couple basic details of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1629, 64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['seqn', 'qsmk', 'death', 'yrdth', 'modth', 'dadth', 'sbp', 'dbp', 'sex',\n",
       "       'age', 'race', 'income', 'marital', 'school', 'education', 'ht', 'wt71',\n",
       "       'wt82', 'wt82_71', 'birthplace', 'smokeintensity', 'smkintensity82_71',\n",
       "       'smokeyrs', 'asthma', 'bronch', 'tb', 'hf', 'hbp', 'pepticulcer',\n",
       "       'colitis', 'hepatitis', 'chroniccough', 'hayfever', 'diabetes', 'polio',\n",
       "       'tumor', 'nervousbreak', 'alcoholpy', 'alcoholfreq', 'alcoholtype',\n",
       "       'alcoholhowmuch', 'pica', 'headache', 'otherpain', 'weakheart',\n",
       "       'allergies', 'nerves', 'lackpep', 'hbpmed', 'boweltrouble', 'wtloss',\n",
       "       'infection', 'active', 'exercise', 'birthcontrol', 'pregnancies',\n",
       "       'cholesterol', 'hightax82', 'price71', 'price82', 'tax71', 'tax82',\n",
       "       'price71_82', 'tax71_82'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"We restricted the analysis to NHEFS individuals with known sex, age, race, ...\" (pg 149, margin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "restriction_cols = [\n",
    "    'sex', 'age', 'race', 'wt82', 'ht', 'school', 'alcoholpy', 'smokeintensity'\n",
    "]\n",
    "missing = nhefs_all[restriction_cols].isnull().any(axis=1)\n",
    "nhefs = nhefs_all.loc[~missing]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 64)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to add some columns to help calculate Table 12.1, and a `constant` column, which will be useful for modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs['constant'] = 1\n",
    "nhefs['university'] = (nhefs.education == 5).astype('int')\n",
    "nhefs['inactive'] = (nhefs.active == 2).astype('int')\n",
    "nhefs['no_exercise'] = (nhefs.exercise == 2).astype('int')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average weight gains in quitters and non-quitters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average weight gain\n",
      "      quitters: 4.5 kg\n",
      "  non-quitters: 2.0 kg\n"
     ]
    }
   ],
   "source": [
    "ave_gain_quit = nhefs[nhefs.qsmk == 1].wt82_71.mean()\n",
    "ave_gain_noquit = nhefs[nhefs.qsmk == 0].wt82_71.mean()\n",
    "\n",
    "print(\"Average weight gain\")\n",
    "print(\"      quitters: {:>0.1f} kg\".format(ave_gain_quit))\n",
    "print(\"  non-quitters: {:>0.1f} kg\".format(ave_gain_noquit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a simple linear model to get a confidence interval on weight difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>    1.9845</td> <td>    0.229</td> <td>    8.672</td> <td> 0.000</td> <td>    1.536</td> <td>    2.433</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    2.5406</td> <td>    0.451</td> <td>    5.632</td> <td> 0.000</td> <td>    1.656</td> <td>    3.425</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = sm.OLS(nhefs.wt82_71, nhefs[['constant', 'qsmk']])\n",
    "res = ols.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            estimate   95% C.I.\n",
      "difference     2.54   (1.7, 3.4)\n"
     ]
    }
   ],
   "source": [
    "est = res.params.qsmk\n",
    "conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
    "lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
    "\n",
    "print('            estimate   95% C.I.')\n",
    "print('difference   {:>6.2f}   ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create Table 12.1 in the margin of pg 149."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style  type=\"text/css\" >\n",
       "</style><table id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aead\" ><thead>    <tr>        <th class=\"index_name level0\" >qsmk</th>        <th class=\"col_heading level0 col0\" >1</th>        <th class=\"col_heading level0 col1\" >0</th>    </tr></thead><tbody>\n",
       "                <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row0\" class=\"row_heading level0 row0\" >Age, years</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow0_col0\" class=\"data row0 col0\" >46.2</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow0_col1\" class=\"data row0 col1\" >42.8</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row1\" class=\"row_heading level0 row1\" >Men, %</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow1_col0\" class=\"data row1 col0\" >54.6</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow1_col1\" class=\"data row1 col1\" >46.6</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row2\" class=\"row_heading level0 row2\" >White, %</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow2_col0\" class=\"data row2 col0\" >91.1</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow2_col1\" class=\"data row2 col1\" >85.4</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row3\" class=\"row_heading level0 row3\" >University education, %</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow3_col0\" class=\"data row3 col0\" >15.4</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow3_col1\" class=\"data row3 col1\" >9.9</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row4\" class=\"row_heading level0 row4\" >Weight, kg</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow4_col0\" class=\"data row4 col0\" >72.4</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow4_col1\" class=\"data row4 col1\" >70.3</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row5\" class=\"row_heading level0 row5\" >Cigarettes/day</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow5_col0\" class=\"data row5 col0\" >18.6</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow5_col1\" class=\"data row5 col1\" >21.2</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row6\" class=\"row_heading level0 row6\" >Years smoking</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow6_col0\" class=\"data row6 col0\" >26.0</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow6_col1\" class=\"data row6 col1\" >24.1</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row7\" class=\"row_heading level0 row7\" >Little or no exercise, %</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow7_col0\" class=\"data row7 col0\" >40.7</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow7_col1\" class=\"data row7 col1\" >37.9</td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                        <th id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadlevel0_row8\" class=\"row_heading level0 row8\" >Inactive daily life, %</th>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow8_col0\" class=\"data row8 col0\" >11.2</td>\n",
       "                        <td id=\"T_9643120e_2c11_11ea_8ea2_dbf823a5aeadrow8_col1\" class=\"data row8 col1\" >8.9</td>\n",
       "            </tr>\n",
       "    </tbody></table>"
      ],
      "text/plain": [
       "<pandas.io.formats.style.Styler at 0x7f51e850fd10>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summaries = OrderedDict((\n",
    "    ('age', 'mean'),\n",
    "    ('sex', lambda x: (100 * (x == 0)).mean()),\n",
    "    ('race', lambda x: (100 * (x == 0)).mean()),\n",
    "    ('university', lambda x: 100 * x.mean()),\n",
    "    ('wt71', 'mean'),\n",
    "    ('smokeintensity', 'mean'),\n",
    "    ('smokeyrs', 'mean'),\n",
    "    ('no_exercise', lambda x: 100 * x.mean()),\n",
    "    ('inactive', lambda x: 100 * x.mean())\n",
    "))\n",
    "\n",
    "table = nhefs.groupby('qsmk').agg(summaries)\n",
    "table.sort_index(ascending=False, inplace=True)\n",
    "table = table.T\n",
    "\n",
    "table.index = [\n",
    "    'Age, years',\n",
    "    'Men, %',\n",
    "    'White, %',\n",
    "    'University education, %',\n",
    "    'Weight, kg',\n",
    "    'Cigarettes/day',\n",
    "    'Years smoking',\n",
    "    'Little or no exercise, %',\n",
    "    'Inactive daily life, %'\n",
    "]\n",
    "\n",
    "table.style.format(\"{:>0.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to be modeling with squared terms and some categorical features. Here we'll explicitly add squared features and dummy features to the data. In later chapters we'll use Statsmodels' formula syntax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Squared features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
    "    nhefs['{}^2'.format(col)] = nhefs[col] * nhefs[col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dummy features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "edu_dummies = pd.get_dummies(nhefs.education, prefix='edu')\n",
    "exercise_dummies = pd.get_dummies(nhefs.exercise, prefix='exercise')\n",
    "active_dummies = pd.get_dummies(nhefs.active, prefix='active')\n",
    "\n",
    "nhefs = pd.concat(\n",
    "    [nhefs, edu_dummies, exercise_dummies, active_dummies],\n",
    "    axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to be creating a lot of IP weights from logistic regressions so a function will help reduce the work. The following function creates the denominators of the IP weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logit_ip_f(y, X):\n",
    "    \"\"\"\n",
    "    Create the f(y|X) part of IP weights\n",
    "    from logistic regression\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : Pandas Series\n",
    "    X : Pandas DataFrame\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Numpy array of IP weights\n",
    "    \n",
    "    \"\"\"\n",
    "    model = sm.Logit(y, X)\n",
    "    res = model.fit()\n",
    "    weights = np.zeros(X.shape[0])\n",
    "    weights[y == 1] = res.predict(X.loc[y == 1])\n",
    "    weights[y == 0] = (1 - res.predict(X.loc[y == 0]))\n",
    "    return weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.535408\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "X_ip = nhefs[[\n",
    "    'constant',\n",
    "    'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2'\n",
    "]]\n",
    "\n",
    "denoms = logit_ip_f(nhefs.qsmk, X_ip)\n",
    "weights = 1 / denoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IP weights\n",
      "   min:  1.05   expected:  1.05\n",
      "   max: 16.70   expected: 16.70\n",
      "  mean:  2.00   expected:  2.00\n"
     ]
    }
   ],
   "source": [
    "print('IP weights')\n",
    "print('   min: {:>5.2f}   expected:  1.05'.format(weights.min()))\n",
    "print('   max: {:>5.2f}   expected: 16.70'.format(weights.max()))\n",
    "print('  mean: {:>5.2f}   expected:  2.00'.format(weights.mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFlCAYAAAApldtwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAASuUlEQVR4nO3dfaykZ3nf8d9VNiSBNLHBCyW7Vhcah4SipFhb4gQ1anGSYkDYfwSJKA0raslSRQkJaYNppCK1UmXaKLyoFZWFHRsFQZDj1lZCXixDElUqbtZAeHNSrwi1Fzt4I4OTBqWJlat/nMftYfesd/fM8Zm9zn4+0mrmueeemfvx2fX3PM/MmVPdHQDg/PY31r0AAODMBBsABhBsABhAsAFgAMEGgAEEGwAG2LfuBTyZSy65pA8dOrTuZQDArrn33nv/pLv3nzx+Xgf70KFDOXr06LqXAQC7pqr+11bjTokDwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwADn9W/reiocuv7XdvTxvnjDq3b08QBgK46wAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAHOGOyqurmqHqmqz24ae1ZV3VVV9y+XFy/jVVXvqapjVfXpqrp8032OLPPvr6ojT83uAMDedDZH2LckecVJY9cnubu7L0ty97KdJFcluWz5c12S9yYbgU/y9iTfl+SlSd7+ROQBgDM7Y7C7+3eTPHrS8NVJbl2u35rkmk3j7+8NH09yUVU9L8k/TnJXdz/a3V9JcldO/SYAADiN7b6G/dzufjhJlsvnLOMHkjy4ad7xZex046eoquuq6mhVHT1x4sQ2lwcAe8tOv+msthjrJxk/dbD7xu4+3N2H9+/fv6OLA4CpthvsLy+nurNcPrKMH09y6aZ5B5M89CTjAMBZ2G6w70zyxDu9jyS5Y9P465d3i1+R5LHllPlvJvmRqrp4ebPZjyxjAMBZ2HemCVX1wST/MMklVXU8G+/2viHJh6vq2iQPJHntMv0jSV6Z5FiSryV5Q5J096NV9W+T/N4y799098lvZAMATuOMwe7uHzvNTVduMbeTvPE0j3NzkpvPaXUAQBKfdAYAIwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAOsFOyq+umq+lxVfbaqPlhV31RVz6+qe6rq/qr65ap6+jL3G5ftY8vth3ZiBwDgQrDtYFfVgSQ/meRwd784ydOSvC7JO5K8s7svS/KVJNcud7k2yVe6+zuSvHOZBwCchVVPie9L8s1VtS/JM5I8nOTlSW5bbr81yTXL9auX7Sy3X1lVteLzA8AFYdvB7u4vJfn5JA9kI9SPJbk3yVe7+/Fl2vEkB5brB5I8uNz38WX+s7f7/ABwIVnllPjF2Thqfn6Sb0/yzCRXbTG1n7jLk9y2+XGvq6qjVXX0xIkT210eAOwpq5wS/6Ekf9TdJ7r7r5LcnuQHkly0nCJPkoNJHlquH09yaZIst39bkkdPftDuvrG7D3f34f3796+wPADYO1YJ9gNJrqiqZyyvRV+Z5PNJPpbkR5c5R5LcsVy/c9nOcvtHu/uUI2wA4FSrvIZ9TzbePPaJJJ9ZHuvGJG9N8paqOpaN16hvWu5yU5JnL+NvSXL9CusGgAvKvjNPOb3ufnuSt580/IUkL91i7l8kee0qzwcAFyqfdAYAAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMsFKwq+qiqrqtqv6gqu6rqu+vqmdV1V1Vdf9yefEyt6rqPVV1rKo+XVWX78wuAMDet+oR9ruT/EZ3f1eS701yX5Lrk9zd3ZcluXvZTpKrkly2/LkuyXtXfG4AuGBsO9hV9a1JfjDJTUnS3X/Z3V9NcnWSW5dptya5Zrl+dZL394aPJ7moqp637ZUDwAVklSPsFyQ5keQXq+qTVfW+qnpmkud298NJslw+Z5l/IMmDm+5/fBn7OlV1XVUdraqjJ06cWGF5ALB3rBLsfUkuT/Le7n5Jkj/P/z/9vZXaYqxPGei+sbsPd/fh/fv3r7A8ANg7Vgn28STHu/ueZfu2bAT8y0+c6l4uH9k0/9JN9z+Y5KEVnh8ALhjbDnZ3/3GSB6vqhcvQlUk+n+TOJEeWsSNJ7liu35nk9cu7xa9I8tgTp84BgCe3b8X7vynJB6rq6Um+kOQN2fgm4MNVdW2SB5K8dpn7kSSvTHIsydeWuQDAWVgp2N39qSSHt7jpyi3mdpI3rvJ8AHCh8klnADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAwg2AAwwMrBrqqnVdUnq+pXl+3nV9U9VXV/Vf1yVT19Gf/GZfvYcvuhVZ8bAC4UO3GE/eYk923afkeSd3b3ZUm+kuTaZfzaJF/p7u9I8s5lHgBwFlYKdlUdTPKqJO9btivJy5Pctky5Nck1y/Wrl+0st1+5zAcAzmDVI+x3JfnZJH+9bD87yVe7+/Fl+3iSA8v1A0keTJLl9seW+QDAGWw72FX16iSPdPe9m4e3mNpncdvmx72uqo5W1dETJ05sd3kAsKescoT9siSvqaovJvlQNk6FvyvJRVW1b5lzMMlDy/XjSS5NkuX2b0vy6MkP2t03dvfh7j68f//+FZYHAHvHtoPd3W/r7oPdfSjJ65J8tLt/PMnHkvzoMu1IkjuW63cu21lu/2h3n3KEDQCc6qn4Oey3JnlLVR3LxmvUNy3jNyV59jL+liTXPwXPDQB70r4zTzmz7v7tJL+9XP9CkpduMecvkrx2J54PAC40PukMAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGGDbwa6qS6vqY1V1X1V9rqrevIw/q6ruqqr7l8uLl/GqqvdU1bGq+nRVXb5TOwEAe90qR9iPJ/mZ7v7uJFckeWNVvSjJ9Unu7u7Lkty9bCfJVUkuW/5cl+S9Kzw3AFxQth3s7n64uz+xXP+zJPclOZDk6iS3LtNuTXLNcv3qJO/vDR9PclFVPW/bKweAC8iOvIZdVYeSvCTJPUme290PJxtRT/KcZdqBJA9uutvxZezkx7quqo5W1dETJ07sxPIAYLyVg11V35LkV5L8VHf/6ZNN3WKsTxnovrG7D3f34f3796+6PADYE1YKdlV9QzZi/YHuvn0Z/vITp7qXy0eW8eNJLt1094NJHlrl+QHgQrHKu8QryU1J7uvuX9h0051JjizXjyS5Y9P465d3i1+R5LEnTp0DAE9u3wr3fVmSn0jymar61DL2r5LckOTDVXVtkgeSvHa57SNJXpnkWJKvJXnDCs8NABeUbQe7u/9btn5dOkmu3GJ+J3njdp8PAC5kPukMAAYQbAAYQLABYADBBoABBBsABhBsABhAsAFgAMEGgAEEGwAGEGwAGECwAWAAwQaAAQQbAAYQbAAYQLABYIBt/z5sNhy6/td29PG+eMOrdvTxANgbHGEDwACCDQADCDYADCDYADCAYAPAAIINAAMINgAMINgAMIBgA8AAPunsPOOT0wDYiiNsABhAsAFgAKfE97idPsWeOM0OsA6OsAFgAMEGgAEEGwAG8Bo258yPngHsPkfYADCAYAPAAIINAAMINgAMINgAMIBgA8AAgg0AAwg2AAzgg1NYOx/EAnBmjrABYABH2Ow5fqUosBc5wgaAAXb9CLuqXpHk3UmeluR93X3Dbq8BzpXX2YF129Uj7Kp6WpL/lOSqJC9K8mNV9aLdXAMATLTbR9gvTXKsu7+QJFX1oSRXJ/n8Lq8D1soRO3CudjvYB5I8uGn7eJLv2+U1wJ7zVLzRbqddiN9U+Mbs/DP5a7Lbwa4txvrrJlRdl+S6ZfN/V9UfPuWr2lmXJPmTdS9iB9mf89uY/al3nNW0MftzlnZ0f87yv+FTxddmC0/R1+RvbzW428E+nuTSTdsHkzy0eUJ335jkxt1c1E6qqqPdfXjd69gp9uf8Zn/Ob3tpf/bSviQz92e3f6zr95JcVlXPr6qnJ3ldkjt3eQ0AMM6uHmF39+NV9c+T/GY2fqzr5u7+3G6uAQAm2vWfw+7ujyT5yG4/7y4aezr/NOzP+c3+nN/20v7spX1JBu5PdfeZZwEAa+WjSQFgAMHeIVV1aVV9rKruq6rPVdWb172mVVXV06rqk1X1q+tey6qq6qKquq2q/mD5Gn3/ute0iqr66eXv2Wer6oNV9U3rXtO5qKqbq+qRqvrsprFnVdVdVXX/cnnxOtd4Lk6zP/9h+fv26ar6L1V10TrXeC622p9Nt/2LquqqumQda9uO0+1PVb2pqv5w+bf079e1vrMl2Dvn8SQ/093fneSKJG/cAx+7+uYk9617ETvk3Ul+o7u/K8n3ZvB+VdWBJD+Z5HB3vzgbb+B83XpXdc5uSfKKk8auT3J3d1+W5O5le4pbcur+3JXkxd39PUn+Z5K37faiVnBLTt2fVNWlSX44yQO7vaAV3ZKT9qeq/lE2Pmnze7r77yb5+TWs65wI9g7p7oe7+xPL9T/LRhAOrHdV21dVB5O8Ksn71r2WVVXVtyb5wSQ3JUl3/2V3f3W9q1rZviTfXFX7kjwjJ32ewfmuu383yaMnDV+d5Nbl+q1JrtnVRa1gq/3p7t/q7seXzY9n43MnRjjN1ydJ3pnkZ3PSB16d706zP/8syQ3d/X+WOY/s+sLOkWA/BarqUJKXJLlnvStZybuy8Q/zr9e9kB3wgiQnkvzicor/fVX1zHUvaru6+0vZOBp4IMnDSR7r7t9a76p2xHO7++Fk4xvgJM9Z83p20j9N8uvrXsQqquo1Sb7U3b+/7rXskO9M8g+q6p6q+p2q+vvrXtCZCPYOq6pvSfIrSX6qu/903evZjqp6dZJHuvveda9lh+xLcnmS93b3S5L8eWadbv06y2u7Vyd5fpJvT/LMqvon610Vp1NVP5eNl8w+sO61bFdVPSPJzyX51+teyw7al+TibLyE+S+TfLiqtvr47POGYO+gqvqGbMT6A919+7rXs4KXJXlNVX0xyYeSvLyqfmm9S1rJ8STHu/uJMx63ZSPgU/1Qkj/q7hPd/VdJbk/yA2te0074clU9L0mWy/P+FOWZVNWRJK9O8uM9+2do/042vkH8/eX/CweTfKKq/tZaV7Wa40lu7w3/IxtnE8/rN9IJ9g5ZvjO7Kcl93f0L617PKrr7bd19sLsPZePNTB/t7rFHcN39x0kerKoXLkNXZvavdH0gyRVV9Yzl792VGfwmuk3uTHJkuX4kyR1rXMvKquoVSd6a5DXd/bV1r2cV3f2Z7n5Odx9a/r9wPMnly7+tqf5rkpcnSVV9Z5Kn5zz/5SaCvXNeluQnsnE0+qnlzyvXvSj+nzcl+UBVfTrJ30vy79a8nm1bzhTcluQTST6TjX/Hoz61qao+mOS/J3lhVR2vqmuT3JDkh6vq/my8E/mGda7xXJxmf/5jkr+Z5K7l/wf/ea2LPAen2Z+xTrM/Nyd5wfKjXh9KcuR8Pwvik84AYABH2AAwgGADwACCDQADCDYADCDYADCAYAPAAIINAAMINgAM8H8BQcwKzJLjJ5EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.hist(weights, bins=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the main model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = nhefs.wt82_71\n",
    "X = nhefs[['constant', 'qsmk']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Weighted least squares gives the right coefficients, but the standard error is off."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>    1.7800</td> <td>    0.288</td> <td>    6.175</td> <td> 0.000</td> <td>    1.215</td> <td>    2.345</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    3.4405</td> <td>    0.408</td> <td>    8.434</td> <td> 0.000</td> <td>    2.640</td> <td>    4.241</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wls = sm.WLS(y, X, weights=weights) \n",
    "res = wls.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GEE gives the right coefficients and better standard errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>    1.7800</td> <td>    0.225</td> <td>    7.920</td> <td> 0.000</td> <td>    1.340</td> <td>    2.220</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    3.4405</td> <td>    0.525</td> <td>    6.547</td> <td> 0.000</td> <td>    2.411</td> <td>    4.470</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gee = sm.GEE(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    groups=nhefs.seqn,\n",
    "    weights=weights\n",
    ")\n",
    "res = gee.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate   95% C.I.\n",
      "theta_1       3.44   (2.4, 4.5)\n"
     ]
    }
   ],
   "source": [
    "est = res.params.qsmk\n",
    "conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
    "lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
    "\n",
    "print('           estimate   95% C.I.')\n",
    "print('theta_1     {:>6.2f}   ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a simple check that there is no association between `sex` and `qsmk`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>763.607760</td>\n",
       "      <td>763.623497</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>801.748892</td>\n",
       "      <td>797.200691</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk           0           1\n",
       "sex                         \n",
       "0     763.607760  763.623497\n",
       "1     801.748892  797.200691"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(nhefs.sex, nhefs.qsmk, weights, aggfunc='sum')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(This matches the R output, but the Stata output is different.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_indices = (nhefs.race == 0) & (nhefs.sex == 1)\n",
    "subset = nhefs.loc[subset_indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now a check for positivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "crosstab = pd.crosstab(subset.age, subset.qsmk).sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAF5CAYAAABk5qjLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZhc5XXm36/2rWvtXb1IAgktSEKNjBcwZrcxqw0EnIGBGG8zcWJskwlOACs2iR2PncRMZozxJjt2vBBjgxHCW2wIBoNQa0egvfe9q6prX7/549atrq6u5d6uW/v5PU8/3XWr6tbXpVa995zvnPcwzjkIgiAIgqhPVNVeAEEQBEEQK4eEnCAIgiDqGBJygiAIgqhjSMgJgiAIoo4hIScIgiCIOoaEnCAIgiDqmIoJOWOslzH2O8bYMcbYUcbYJ1LHdzLGxhhjB1Jf763UmgiCIAii3mGV6iNnjHUB6OKcDzLGWgDsA3AzgD8B4Oecf7kiCyEIgiCIBkJTqRfinE8AmEj97GOMHQOwaiXnam1t5atXr1ZwdQRBEARRu+zbt2+Wc96W676KCXkmjLHVALYDeAXAxQA+zhj77wBeA/Bpzrm70PNXr16N1157rdzLJAiCIIiagDE2lO++ihe7McYsAH4K4D7O+QKArwE4B8AFECL2r+R53kcYY68xxl6bmZmp2HoJgiAIopapqJAzxrQQRPwHnPMnAYBzPsU5T3DOkwC+AeCiXM/lnD/OOd/BOd/R1pYzu0AQBEEQTUclq9YZgG8BOMY5/6eM410ZD3sfgCOVWhNBEARB1DuV3CO/GMBdAA4zxg6kjv0NgA8wxi4AwAGcBfDRCq6JIAiCIOqaSlatvwiA5bjr2UqtgSAIopFIJpOYnZ2Fx+NBIpGo9nKIEjEYDOjp6YFWq5X1vKpUrRMEQRClMzo6CsYYVq9eDa1WC2EHk6hHOOeYm5vD6Ogo1qxZI+u5ZNFKEARRpwQCAaxatQo6nY5EvM5hjMHlciEcDst+Lgk5QRBEHaNS0cd4o7DSizH6CyAIgiCIOoaEnCAIgmharr32Wnz3u9+t9jJKgoScIAiCaFr27NmDu+++GwCwa9cuXHLJJUvuv+eee/Dggw9WY2mSaXohPzsbwO/emK72MgiCIIgGJB6Pl/01ml7Ifzo4inu/uxfJZGXGuRIEQTQLq1evxpe//GVs3boVNpsNt99+e7oq+xvf+AbOPfdcOJ1O3HjjjRgfH08/jzGGxx57DOvWrYPD4cCf//mfo9DI7V//+tfYsGEDbDYbPv7xj+Nd73oXvvnNbwIAdu7ciTvvvDP92LNnz4IxlhbYyy67DN/85jdx7NgxfOxjH8PLL78Mi8UCu92Oxx9/HD/4wQ/wpS99CRaLBTfccAMAYHx8HLfccgva2tqwZs0aPProo+nz79y5E7feeivuvPNOWK1W7Nq1C6+++ip27NgBq9WKjo4OfOpTn1LuTQb1kcNm1CLJAX80DqtBXhM+QRBELfF3vziK18cXyvoam7qt+OwNmyU//ic/+Qmee+45GAwGXHzxxdi1axfWr1+Pz3zmM/jVr36FzZs34/7778cdd9yBF154If28Z555Bnv37sXCwgIuvPBC3HDDDXjPe96z7Pyzs7O45ZZb8O1vfxs33XQT/vVf/xWPPfYY7rrrLlm/18aNG/HYY4/hm9/8Jl588cX08Zdeegk9PT145JFHAAgmPDfccANuuukm/PCHP8To6CiuuuoqnHfeeXj3u98NAHjqqafwxBNP4Hvf+x4ikQiuuOIKfOITn8Bdd90Fv9+PI0eUdSJv+ojcahTE2xuMVXklBEEQjcdf/uVforu7G06nEzfccAMOHDiAH/zgB/jgBz+IgYEB6PV6fOELX8DLL7+Ms2fPpp/3wAMPwG63o6+vD5dffjkOHDiQ8/zPPvssNm3ahFtvvRVarRb33XcfOjs7y/b77N27FzMzM3j44Yeh0+mwdu1afPjDH8aPfvSj9GPe/va34+abb4ZKpYLRaIRWq8XJkycxOzsLi8WCt73tbYquiSJyUchDMfRWeS0EQRClICdSrhSZomoymTA+Po65uTkMDAykj1ssFrhcLoyNjWH16tU5n+f3+wEAmzdvxtCQMJp7z549GB8fR2/v4qc3Y2zJbaUZGhrC+Pg47HZ7+lgikcA73/nO9O3s1//Wt76Fhx9+GBs2bMCaNWvw2c9+Ftdff71iayIhzxBygiAIovx0d3enxRgQHOrm5uawatWqos89evToktunT5/GyMhI+jbnfMlts9mMYDCYvj05OZn33LkMWbKP9fb2Ys2aNThx4oTk86xbtw4//OEPkUwm8eSTT+LWW2/F3NwczGZz3nPIoelT6yTkBEEQleVP//RP8Z3vfAcHDhxAJBLB3/zN3+Ctb31rOhqXw3XXXYejR4/iySefRDwex6OPPrpErC+44AK88MILGB4ehtfrxRe+8IW85+ro6MDo6Cii0eiSY6dPn07fvuiii2C1WvGP//iPCIVCSCQSOHLkCPbu3Zv3vN///vcxMzMDlUqVjuTVarXs3zUfTS/kdhMJOUEQRCW58sor8fnPfx633HILurq6cOrUqSV7zHJobW3FE088gQceeAAulwsnTpzAxRdfnL7/6quvxu23346tW7fiwgsvLJjSvuKKK7B582Z0dnaitbUVAHDvvffi9ddfh91ux8033wy1Wo1f/OIXOHDgANasWYPW1lZ86EMfgtfrzXve5557Dps3b4bFYsEnPvEJ/OhHP4LBYFjR75sLVqikv1bZsWMHf+211xQ5VzAax6aHf4kHrt2Aj73rHEXOSRAEUQmOHTuGjRs3VnsZNcdll12GO++8Ex/60IeqvRTZ5Ps3ZYzt45zvyPWcpo/IjVo1tGoGD1WtEwRBEHVI0ws5Yww2o5ZS6wRBEERd0vRV64DQS75AQk4QBNEQ/P73v6/2EipK00fkACgiJwiCIOoWEnIAdhJygiAIok4hIQdF5ARBEET9QkIOQcg9wWjxBxIEQRBEjUFCDkHIfZE4jTIlCIIg6g4ScghV65wDvnD5B8ATBEEQhfnYxz6Gz3/+89VeRt1AQg7yWycIgqglHnvsMTz00EMAhFaynp6eJffv3LkTd955ZzWWVpOQkAOwm3QASMgJgiCagXi8sbKvJOSgiJwgCKIc7N+/HwMDA2hpacHtt9+OO+64Aw8++CB27dqFSy65ZMljGWM4efIkAOCee+7Bgw8+iEAggGuvvRbj4+OwWCywWCz493//d/zDP/wDfvzjH8NisWDbtm0AAK/Xi3vvvRddXV1YtWoVHnzwQSQSCQDArl27cPHFF+OTn/wknE4ndu7ciZMnT+Jd73oXbDYbWltbcfvtt1f2zVEQcnbDopB7QlS5ThBEHbPnAWDycHlfo3MLcO0Xiz4sGo3i5ptvxn333YePf/zjeOqpp/CBD3wAf/3Xfy35pcxmM/bs2YM777wTo6Oj6ePHjx/HyZMn8f3vfz997O6770ZHRwdOnjyJQCCA66+/Hr29vfjoRz8KAHjllVdwxx13YHp6GrFYDB/84AdxzTXX4He/+x2i0SiUGsRVDSgiB0XkBEEQSvPHP/4RsVgM9913H7RaLW699Va85S1vKctrTU1NYc+ePfiXf/kXmM1mtLe345Of/OSS0ajd3d34i7/4C2g0GhiNRmi1WgwNDWF8fBwGg2FZhqCeoIgcJOQEQTQIEiLlSjE+Po5Vq1aBMZY+1t/fX5bXGhoaQiwWQ1dXV/pYMplEb29v+nbmzwDwpS99CQ899BAuuugiOBwOfPrTn8YHP/jBsqyv3JCQAzBoVdCpVSTkBEEQCtHV1YWxsTFwztNiPjw8jHPOOQdmsxnBYDD92MnJybznybwQyHest7cXer0es7Oz0Ghyy1r2czo7O/GNb3wDAPDiiy/iqquuwqWXXopzzz1X2i9YQ1BqHalRpiaagFaMn+0fxZnZQLWXQRBEHfD2t78dGo0Gjz76KOLxOJ588km8+uqrAIBt27bh6NGjOHDgAMLhMHbu3Jn3PB0dHZibm4PX611y7OzZs0gmkwCEi4ZrrrkGn/70p7GwsIBkMolTp07h+eefz3veJ554Ir3v7nA4wBiDWq1W4DevPCTkKchvvTDJJMf9TxzCd186W+2lEARRB+h0Ojz55JPYtWsXHA4HfvzjH+P9738/AGD9+vV4+OGHcdVVV2HdunUF96c3bNiAD3zgA1i7di3sdjvGx8dx2223AQBcLhcGBgYAAN/73vcQjUaxadMmOBwO3HrrrZiYmMh73r179+Ktb30rLBYLbrzxRnz1q1/FmjVrFHwHKgfjvP5sSXfs2MGVrjC85WsvQa9R4d8//DZFz9soeIMxbPvcr/DuzR34+l07qr0cgiAAHDt2DBs3bqz2MiRzzz33oKenB4888ki1l1Kz5Ps3ZYzt45zn/PCliDwFReSFEd+bCW+4yishCIIgMiEhT0FCXhgScoIgiNqEqtZTkJAXRnxvZv0RRONJ6DR0DUgQhDx27dpV7SU0JPRpnMJq1MIXjiNBo0xzIgo558DUAkXlBEEQtQIJeQp7yhTGF6aoPBeZ2QpKrxNE7SC2YBH1z0qLz0nIU6T91oMk5LlYKuShKq6EIAgRs9mMsbExRKPRFYsAURtwzjE3NweDwSD7ubRHnoJsWgvjCUWhYkCSU0ROELVCT08PZmdnMTQ01HCjOZsRg8GwbPa6FEjIU9hMJOSFWAjF4DTrEYknMOGhiJwgagGVSoX29na0t7dXeylEFSEhT0EReWG8oRjsJi3UTIdxisgJgiBqBhLyFCTkhfGGYrAZtbDoNZgkIScIgqgZqNgtBQl5YUQh77YbqNiNIAiihiAhT2HQqqHX0CjTfIhC3mk1YtYfRSSeqPaSCIIgCJCQL8Fm1MJL7Wc58QQFIe+yC60RU95IlVdEEARBACTkSyCb1twkkhy+cBxWoxbdNiMAYJzS6wRBEDUBCXkGJOS5Ed3u7EYtOm1CRE4FbwRBELUBCXkGJOS5Ed8Tm1GLrpSQU0ROEARRG5CQZ2AzkZDnIlPIzXoNrAZqQSMIgqgVSMgzoIg8N2khT7nfdduNGPeQkBMEQdQCJOQZ2Ixa+CNxxBM0TSgTcZCM2GvfZaNecoIgiFqhYkLOGOtljP2OMXaMMXaUMfaJ1HEnY+zXjLETqe+OSq0pG1GoFsI0fCCTzNQ6AHTajJRaJwiCqBEqGZHHAXyac74RwNsA/DljbBOABwD8lnO+DsBvU7erArm75SZbyLttBswFogjHyBSGIAii2lRMyDnnE5zzwdTPPgDHAKwCcBOA76Ye9l0AN1dqTdmQkOdmIRSDXqOCQasGgHQL2tQCReUEQRDVpip75Iyx1QC2A3gFQAfnfAIQxB5A1ebx2cs0yvQ7fziDR555XdFzVhLRnlWk254yhaGCN4IgiKpTcSFnjFkA/BTAfZzzBRnP+whj7DXG2GszMzNlWZsoVp5gVLFzJpIc//d3p/Afg6OKnbPSZAu52EtOBW8EQRDVp6JCzhjTQhDxH3DOn0wdnmKMdaXu7wIwneu5nPPHOec7OOc72trayrI+q1jspmBE/uqZecz6I/AEY3Wbshd91kW6UjatE1TwRhAEUXUqWbXOAHwLwDHO+T9l3PU0gLtTP98N4KlKrSmbcuyR7z48nv55ZD6o2HkrSXZEbtSpYTdpKSInCIKoASoZkV8M4C4AVzDGDqS+3gvgiwCuZoydAHB16nZV0GvUMGiVG2WaSHI8d2QSa9vMAIDhBhFyAOi0GqgFjSAIogbQVOqFOOcvAmB57r6yUusohpLubq+cmcOsP4r/9e4N+F8/PVS3Qr4QiqVd3UTI3Y0gCKI2IGe3LOxGnWJCvvvQBIxaNa7f1gWnWYehufoT8kSSwxeJL4vIyd2NIAiiNiAhz8Jm1KYtSUshnkjiuSOTuGJjO0w6DXqdprrcI1/IMoMR6bIZ4A7GyBSGIAiiypCQZ2FVKLX+ypl5zAWiuH5LFwCgz2mqy9R6tqubCFWuEwRB1AYk5FnYjFpF2s92H56ASafGZecJ/jb9ThPGPCHE6mwgi6dARA5QLzlBEES1ISHPQoliNzGtfuXGDhh1gq1pn9OERJJjos4KxPJG5Cl3t3r7fQiCIBoNEvIsbEYtAtFESZHzH0/PYz4QxXVbOtPHep0mAPXXgiYKud1EETlBEEQtUrH2s3rBZhTekoVQDC6LfkXn2H14fElaHQD6XYKQD80HcAlaS19ohRCF3JoVkRu0ajhM2prbIx+ZD+LktB+Xb6iaZX9FODntx9MHxsBLOIeKMdxxUW+63oEgiPqEhDwLu0kHQNgbXomQi2n1qzZ2pKeFAUCH1QCdWlV3EXm+qnVAKHirNSH/1otn8JPXRvD6595T7aWUlW+8cBo/fm0EqnzODBJIcoAx4L6r1iu3MIIgKg4JeRal2rS+fHoO7mAM123tWnJcrWLocRjrrgXNG4rBoFVBr1Evu6/bbsBYje2RL4RjCEYTCMcSSy6kGo0pXxhbe2x4+uOXrPgc2z/3K8z5lRsQRBBEdaA98iysJQr57kMTMOvUeNf65YNd+lymujOF8QSjOaNxQJhLXmt75IFIHADgVnCCXS0yvRBB2wq3fkRcFj3mAhGFVkQQRLUgIc/CVsIEtFgiieeOTuKqTR05o8E+pwnDc0FwXsrOZmXJ5bMu0mUzwhOMIRStHVOYQERYiztQn5PmpDLjj6DdWpqQO806isgJogEgIc+ilNT6S6fm4AnGcN2Wrpz39zlN8EXidTXO1BuKwW7U5byvFivXfamIXMmZ8rVGIskx5y89Im+16DAXaNz3iSCaBRLyLNJCvgKb1mcPTcCi1+DSHGl1QBByAHWVXveG4ssq1kVq0d1tMbVePxdLcpkLRJDkQJvVUNJ5nGYd5knICaLuISHPQqdRwaRTpx3NpCKm1a/Ok1YHhD1yoL56yRcKpNa77WJEXotC3rgCNb0g7GuXGpE7zXq4g1EkkvWz1UMQxHJIyHOwEne3P5ychTeUP60OAL2O+hPyQnvkHamIcMJTO6l1fxOk1md8gpCXukfeatGB88a+6CGIZoCEPAcrEfLdhybQotfgnevzm72Y9Rq0WvQYrpPUeiyRhD/HCFMRg1YNl1mH8RqJyDnn6Yh8voGL3UQhLz0iF2ofKL1OEPUNCXkO5E5Ai8aT+GUqrZ6r3zqTPqexbiLyRTOY/HYDnTYDJmuk2C0US0DMEjdyRD7tEy6c2lpKbD8zC8+f9VMLGkHUMyTkOZA7Ae0Pp2axEI4vM4HJRT2NM00PTDHljsiB2nJ3E9PqQGOni2d8EVgNmpINb1wWisgJohEgIc+B3NT67kMTaDFocMm64h7qfS4zJrwhROO1P840PTAlT/sZIBS81YqQiz3kQGNXrU/7ImgvsWIdAFyp1PpKe8mTycWtDIIgqgcJeQ7sRi08EoWAc47fHJvC1RuLp9UBISJPcmCshgrE8pFvYEomnTYDvKEYgtHqf6CLomLUqhs6tT7jK72HHBDmCjCGFfeS/8fgKN7xxf9EJF47hkAE0YyQkOfAZtQiFEtIiponF8LwBGPY3meXdO6+Ohpnmm8WeSbdqV7y8RrwXBdT6z0OYxNE5KULuVrF4DTpMLfCPfI3J33whmKUmieIKkNCngNxT1hKev34lB8AsK6jRdK560nIC00+E+lMubtN1kB63R8WhLzXacJCONaQ/dGcc8UicqA0Uxixer7R7XAJotYhIc+BHJvWE1M+AMB6iULe3qKHXqPC8Fxg5QusEOL2gqSIvAYq1wPRxYic85UPvqll/JE4QrGEIhE5UJrfulg938jbGARRD5CQ50DOBLTjUz60WnTpntxiqFQMvXVSue4NxWDUqqHT5P8z6bAJglITEXlGah1ozMr1dA95ia1nIq0lTEAT1zLfgO8zQdQTJOQ5kDMB7fiUH+vapUXjIkILWvUj2GJ4QzHYC7SeAYBeo0arRVcTg1MCaSEXti8aMVKcFl3dWkqvWgdKS62La2nkegSCqAdIyHNgTwm5J1T4A45zjpPTfqzvsMg6vzDONFDz40wL2bNm0mUz1kixWwKMAd32VETegHu3SkfkLosO7mAM8YS8dshwLAFfqibBQ8VuBFFVSMhzIHUC2rg3DH8kLrnQTaTPaUIgmqj5al9vKFaw9UxEcHerASEPx2HWaeA0CdscjZhaX4zIFRJys/heybvoES8oVvJcgiCUhYQ8B4t75IV7o4/LLHQTqZfKdakRebfNUBvFbpE4zHo17OZURqUBBWbGF4FOrZL07yIFV6r6Xe4++XSGkDfiFgZB1BMk5DnQqlUw69RFi91Oiq1n7fJS6/11Ms5UcmrdboQvHF9ikVoN/NE4zHoNWvQaaFSsQSPyMNpa9GCMKXK+9OAUmZXrM6mKda26Md9ngqgnSMjzIMWmVahY18MhsWJdRCzGqvUpaNL3yMVe8upG5YFIHC16DRhjsJu0DSkwM74IWhVKqwPCKFMAmJW5zSOm1s9ps1BqnSCqDAl5HqRMQDu+gkI3ADDq1Ghv0dd0RB5LJBGMJtKFf4XoqhF3NyG1Lkxqs5t0DVvsptT+OAA4UxPQ5mW6u037IlAxQcgptU4Q1YWEPA9CRJ7/A4pzjpNTPtn74yL9LhOGaljIpUw+E+mqEXc3X3hRyB0NHJErVbEOCB0aqhX4rc/4InCa9emq92oxMh/Ex/5tX014/RNEtSAhz4PdVDgiH/OEEIgmsG4FETkg2IiO1IOQS4jIO6wGMFZ9d7dANA5LRkTeaMVusUQSc4GoohG5SsUEdzeZQj6dygzYTbqq2uG+cGIGzx2dxBuTvqq8PkHUAiTkeSi2R34iVei20oi8z2nC5EIY4VhtTo6SMvlMRKdRodWir3pEHogkYNYLE+gaMSIXrVSVMoMREWxa5aXWxcyA06Stqh3uRGo7R26xHkE0EiTkeSgm5OnWM5mubiL9LhM4B0bd1W/byoVXgs96Jl02A8arLOT+SBwWvbBeRyoir3XTHTmI3uZKptYBwGXWy/Y0mPaF0d6yWOhZrYumidTf3EptZgmiESAhz4PNqEU4lsw7a/n4lB/tLXpJe8i5EHvJazW9Lie1DghCPlHFGevReBLReBIWMSI36xBNFew1CjMKm8GIOC3yUuvJJMesP4q2VGodqF4vuWgNvNKZ6gTRCJCQ56HYBLQT0ysvdAOEPXKgdnvJ5Qu5saqpddFnPbPYDWgsd7dphe1ZRVplTkCbD0aRSHIhIhff5yp1CIh/cyud4EYQjQAJeR6sBWxak0mOE1P+FRe6AUCbRQ+jVo2hGu0lX0lE7ovE4QtX5wPdnyXki5Fi4xS8iRF5q0KzyEWcZj28oRhiEv3W05kBqwGOKtrhcs7TBZa1bndMEOWEhDwPohDkisjHPCGEYomSInLGWGoKWu0KuVmnhlYt7U+kKzWopFpRuTiL3JKOyBvPb33aF4bDpC04VnYluFKmMG6JYpiZGbBXMfPhCcYQjgkXH7Myi/UIopEgIc9DodT6iWmh0E2uNWs2tdyCJtXVTUTsJa9WwZuYWrcsS603VkSudMU6sDg4Reo+8/SC8G/c3qKHJW2HW/n3WSx0Y4wicqK5ISHPQyEhPy56rJcQkQNIR+S1WFntCUqbfCYiCnm1Ct78EaGobXlqvXE+4KcVNoMRSQ9OkbjPPONfjMgFO1xdVd5nsdDtnDYL7ZETTQ0JeR4KC7kPHVZ9yROo+l0mhGKJ9AdjLbEgMyIXTWEmqhSR+8NLI3J7lYuwyoHS9qwiznRELu3vcHohAoteA5Muw0WvCu+z+Ld2frcV84FoTV4QE0QlICHPg9UgfEjlTK1P+UvaHxep5RY0ual1rVqFNos+HSVVmsWqdXV6PS16TcPskXPOyxeRi0IuIyLPXIfDpKvK+zzhDUGtYjiv04poIglflafvEUS1ICHPg0atgkWvWVb1nExynJz2Y90KjWAyEVvQarFy3RuKpaNaqXTZjdWLyLP2yAHAbtaWJeX75qQPR8a8ip+3EAvhOKLxZFmE3GbUQq1ikveZZxaWCrndpK1Kd8CEJ4yOFn06S0HubkSzQkJeAJtRi4WsiHzULVasl1boBgA9DiMYq81ecrkROQD0Oow4PRMo04oKk91HDgiR4nwZBOYzTx7Ch7/3GpIV9BefKZOrGyD4rTtMOsmp9dqJyMPoshvhtMgr1iOIRoOEvAC5bFpFa9ZSC90AwKBVo9NqqDkhj8aTCMUSsoX8gl47xjyhdFVzJfFH49BrVEva5cpRhBWOJXB4zIsJbxj7R9yKnrsQ0wuiq5vyVeuAMJdcamp9eiG8ZK9eyHxU3g53whtCp82AVrNYrFd7tSYEUQlIyAuQU8jF1jMFInKgNlvQ5JrBiAz0OwAAg8OVEzgRfzi+JK0OlGdwytFxL2IJQbCeOTSh6LkLkVkpXg6kTkALROIIRBPLIvJK2+FyzjHhDaPbZkhH5NSCRjQrJOQFyCXkJ6b86LIZYDWUVrEu0u801dweuTiHXU77GQBs7rZCp1ZhcNhTjmUVJBCJL0mrA6nBKQpXUw8OCb/bQJ8dzx6eqFh6PR2RW8sj5C6LtMEpi37vi5mBatjhuoMxROJJdNmMsvvgCaLRICEvQE4hn/YpklYX6XOaMO2LIFRDwz1WGpHrNWqcv8qKwaEqROSRxDIht5u08EXikq1HpbBvyI0+pwl3v2M1phYiFcs+zPgj0GuESvxy4JI4yjRXZqAadrjjKb+CLpsBBq0aZp2aesmJpoWEvAA2kxaeDCEXK9bXl+jolkmfS6hcH3XXTlQuCrn4AS2HgT4HDo15EY0rJ55SCETiy0RO7I9WSmA45xgcdmOgz44rN3ZAp1FVLL0+vRBGu1UwYCkHLrMuXRlfeB3LJ7A5qzDKVLQCFq2BhQlutEdONCck5AWwGbWIxpMIx4RoecQdRDiWVGx/HFjsJa+l9PpKI3JA2CePxpM4Ol7Z9qxANJ7uIRdR2t1tzBPCtC+CgX4HLHoNLj+vrWLp9Rl/BG0KD0vJRNxnLibGuarnq2GHK/oViI6CK5mpThCNAgl5AbLd3ZSyZs2krwbHmYoT31Yi5Aj03uEAACAASURBVBemC94qu0/uD+faI1dWYMTfaaBP+B3fu6UL074IXqvAVsL0Qnl81kXEfeZiw0emfRGoVQzOjGxNNexwJ7xhaFQsPQnOZdZhllLrRJNCQl6A5UKuzLCUTJxmHcw6dW0JeUjoyRbd7eTQYTVgld1Y8cp1fyRX1bqyKd/BITeMWjU2dAoXcldu7IBeo8Kzh8ufXs/u3VYa0W+9WFQ744ug1aKDSrWY4rcbK2+HO+ENo8NqgDq1DpdFh3lKrRNNSsWEnDH2bcbYNGPsSMaxnYyxMcbYgdTXeyu1HilkC/mJKR+6bQa0KFSxDqTGmbrMNSXknlBUmGolcYRpNtv77Nhf4YK3XFXrojOdUpHi4LAb23pt6fdFSK+349nDE0iUMb0eiSfgCcbK4rMu4pRo0zqdYwKbRq1Ci6GydrjjnlA6rQ4IM9XJb51oVioZke8C8J4cx/+Zc35B6uvZCq6nKGkhDy6m1pVMq4v0OY01JeQrcXXLZKDPgXFvuGK+68kkRyCaKBCRlx4phmMJvD6+kE6ri1y3NZVePztf8mvkQ0wZlzMiT5uqSIjIc62j0u5ukwvhdKEbIBjaxBIcC2HyWyeaj4oJOef8BQDl+7QrA6KYeUIxJJIcp2b8ilizZtOXMoWppOVnIRZC8kaYZpM2hhmqzD55MFWMmC3kJp0aOrVKEYE5NOpFPMmXCfkVG9ph0Kqwu4zp9fT87zL1kAOA1SjMFS/WgjadZwKbYL5TmdS6aAazNCInUxiieamFPfKPM8YOpVLvjuIPrxxiatYbimF4PohIPFmWiLzfZUYknsThCg/iyIc3FEvve66ETV1W6DWqiu2TiyNMs1PrwqxsrSKmMOLvsr3PvuS4OZ1enyxbel00YWmzlK/YjTEGp1lXUAgTSY75QG4hr+RM8vlAFNF4comQL85Up31yovmotpB/DcA5AC4AMAHgK/keyBj7CGPsNcbYazMzMxVZnLgX7g3F0oVuSowvzeb6rV1oa9Hjr396qOL917koNbWu06iwtcdWOSHPGmGaiVIp331Dbqx2mdKCkcl1W7sw649gb5nS69O+8rq6iTiLVH7P+SNI8twp/nLY4eZDnLDXZVtMrZO7G9HMVFXIOedTnPME5zwJ4BsALirw2Mc55zs45zva2toqsj61iqHFoMFCKIaT06nWMwUr1kXsJh3+4X1b8MakD//6nycUP79cShVyQNgnPzq2kO7BLyeBHCNMRewKCAznHPuH3ektg2zS6fUymcPM+CJgbFGsykWrRV+w8lu8oGjL0QZnL4Mdbj4yXd1EpBbrEUQjUlUhZ4x1Zdx8H4Aj+R5bLUSb1uNTPqyyG5elb5Xi6k0deP/AKvzf35/C4dHqptg9wRhsMmeRZ7O9z4FoojLGMLlGmIoIEXlpAjMyH8KsP7psf1zEpNPgyg0d2HOkPNXr074IXGbdirsIpFIstZ5O8ecpdlPaDjcfkwuiq1uuPXJKrRPNRyXbz34I4GUA5zHGRhlj9wL4EmPsMGPsEIDLAXyyUuuRyqKQ+xV1dMvFZ6/fjFaLDvc/cRCReHW818OxBCLxZOkReb+wl1yJgjdfgYjcYdaWvHcrbhHkE3JATK9H8cqZuZJeKxdC73Z50+pAagJagYh2cWBKDiE3i61+5Y/Kxz1haNUsXWkPCCOBLXoNmcIQTUklq9Y/wDnv4pxrOec9nPNvcc7v4pxv4Zxv5ZzfyDmv3FxIidiMWswFoqmKdeX3x5e8lkmLL7x/C96c8uHR31Ynxb6Q6pkvpWodEKZj9TorYwxTOLWuK3lW9uCwG2adGud15v/3v/y8dhi16rKk12d8YbRby1foJtJqEaLqfBeR0znsWUUq6e426Q2hw2pYYkoDiKYwJORE81HtYreax2bU4o2JBUTjybLsj2dzxYYO3HphDx57/jQOjlR+HGh6YEqJQg4IEezgsLvsJh2FU+taxJM8HbWvhH1DblzQZ0+7iOXCqFPjio3t+OXRScQVTi/P+Mrrsy7iNBd2d5vxRWA1aGDQ5ioqrJzf+rg3jO6MQjcRYaY6pdaJ5oOEvAh2kxaRVCV5uSNykYeu34Q2ix73P3GwIsVimZQyMCWbgT4HphYiGPOU1xjGH8ndRw4smsKstBArGI3jjUlfwbS6yPVbhPT6q2eUq17nnGPGHyl7xTogRLRA/oKx6TxmMIDydriFmPCG0GlbnqFwmfVU7EY0JSTkRchMMZ9bgYgcEET0i7dswYlpP75a4RS70kIOlH+ASiASh4oBBu3yP+dSBebgiBeJHEYwubjsvHaYdGo8o6A5jCcYQyzBKxKRF2vhmslhzyqitB1uPpJJjilvZEmhm4irSLEeQTQqJORFEAWtx1G+ivVcXHZeO27f0YuvP38K+ys4gMRTwuSzbDZ0tcCgVWGwzL7r4sCUXLO6xSKslQp5PiOYXBh1aly5sQPPHVEuvV6pHnIgc3BK7vR0oYh8cSZ5eVPrc4EoookkunLUDIh75OS3TjQbJORFEAWtUmn1TP72+o3otBoqmmJXMiLXqlXY1mMv+4VIrslnIotFWCsTmMEhN9a2mdPnKcZ1WzoxH4jij6eVSa8vurpVpmodyJ1a55ynIvLc6zBq1dBplLHDLcSkaAZjz71HHk9yLITIb51oLkjIiyAKWrlbz3JhNWjxxVu24tRMAP/86+MVeU2vQlXrIgP9DhwdL68xTK7JZyKlpNY559g/4sGFEtLqImJ6XSnvdbFSvBJV61aDBlo1y5la90fiCMUSeSNyxhgcCtnhFmI8NYgnV7GbuMc/SwVvRJNBQl6EdETeXvmIHAAuXd+GD1zUh8f/6zQOja58r/mPp+fwhWePFX2cNxRDi0FTsEJbDgN9DsSTHIckmNz4wjH81RMHcXY2IOs1/AWE3GbUgrGVpXzPzgUxH4jmdXTLhUGrxlUbO/DckQlF0uuFTFiURvRbz+VXPiMhxV+JCWgTqcLJfMVuAA1OIZoPEvIibOu146YLuvGu8ypjC5uLv3nvBqgZw54jkys+x7+9PISvv3A67RmfjwUF7FkzEfeWpfST//3uY3hi3yheOCHPSz9QILWuVjFYDSszhRH39qUUumVyybpWuIMxRar1p30RmHTqvL+f0ohzvXOtAyg8uMVu0pbdEGZiIQydWpXTrnZxa4AicqK5ICEvgtWgxVfv2F4RZ618tBi02NxtLaloTBTSZ4oYlijhs55Jq0WPfpep6Np//+Y0frR3BIB8v+xCe+TAykdsDg670aLXyPYP6HeaAECRGfP55n+Xi1aLLmdqXWpEPl/2iDyMDpt+mRkMgPT/URqcQjQbJOR1wvY+Bw6OelbkZT3uCWHCGwZjwO5D4wWrepUWcgC4sM+BwWFP3tf1hmJ44KeHsa7dAqtBI9vUIxBJFOwoWOmITdEIJpdoFKLPJQj50FzpQj7tC+ctMCsH+WxapyUU3VVilOmkN7xk6lkmYocC9ZITzQYJeZ0w0O9AOJbEGxOFU+O5EKPxP7mwF6dmAjg+5c/7WE8ZhHx7vwOz/ghG5nOnmh955nXM+CP4yp9sQ1tL7tRuIYSIfLnbmMhKRmz6I3Ecn5JmBJNNR4sBOo0KI3UYkbvypNZnfBHo1Kp0v3guHKnUejnbv8a9IXTn2B8HAL1GjRaDhvbIiaaDhLxOuLBfNFeRn14fHPLAoFXhk1evhyoVleejHBH5QIF98v98YwpP7BvF/3jXOdjaY5ftzsU5L1i1DqSKsGRWUx8c8SDJIavQTUSlYuh1GBVJrU8XMGEpBy6LDv5IfFmXwbQvjLYWfc5efRGHSVeyHW4hkkmOqYUwOvNE5IBgCkOpdaLZICGvE7ptBnRY9SsT8mE3tvbY0Wkz4G1rXXjm8ETBNLfSQn5eRwtMOvWytXuDMXzmycPY0NmCv7jyXACCkMj5II7Ek4gnedHUutyIXNzTv6C3uBFMLvqcppJT6+FYAr5wvMIRuTgOdOn7NeOLoLXIOtLubmVqQZsNRBBLcHTncHUTcVn0VOxGNB0k5HUCYyw9hEQO4VgCR8e96RTxdVu7cHomgDcml6fow7EEovFkybPIs9GkjGGy1/53zxzFrD+KL9+2DXqNkBovNhM7G38q+msxFC52C0YTskbD7ht2Y127ZcUXNX1OE0bmgyWlmSvZeiaSzxSmkBmMSLn91ic8Qk99Z4Geerl/PwTRCEgWcsZYH8uRV2MCfcoui8jFQJ8DI/OhtEmIFI6MeRFL8HR6+z2bO1Pp9eXV60q6umVzYb8DxyZ8CEYF4f3N61N4cnAMf375uTh/lS39OJdFD3cwikRSmgCmJ5/pCkTkZnnubskkx/5hz4r2x0X6XGb4IvGS2rGmqyDkrnTl99KoVspefal2uMWYSLm6dedwdRNptehoJjnRdMiJyM8AyNVM7UzdR5SZgf7UXvOQdGMYMQoW93pdFj3efo4Lu3Ok18sp5AP9diSSHAdHvPAEo/jMzw5jY5cVH7/83CWPc5l14Fy6GPgLjDAVWRyxKe2cp2cD8IZi6bqEldCXakEbKmGffEZ0datCaj0zIo8lkpgLRIuuo1Q73GJMpFzduvIUuwFCRO4ORpGUeCFIEI2AHCFnAHL977AAkB4iEitmc7cNOrVKlnf54JAHfU7Tkj7467Z048xsAMeyKuCVHJiSzfbexWK9nU8fhTsQxZdv2wqdZumfYLFRmtkECowwFXGKKV+Je7eLFz8r2x8HFoW8lIK3akTkTsvyPXLx36JoRF7u1Lo3DJ1GlU7/58Jp1iOR5OmLUoJoBoraRTHGHk39yAF8gTGW+cmkBnARgANlWBuRhUGrxuZVVsn75Jxz7Bt245JzW5ccf/fmDjz01BHsPjyOTd3W9PFyRuQOsw5rW834t5eHMLkQxievWo/N3bZlj0vv0QYiAIrb4qZT6wXazxYjRWkCs3/YDatBg7WtK/fXF4W8lBa0GV8EKrZoPVoJWvQa6NSqJX7lab/3ItXzpdjhSmHCG0aXzVCwcr7VsjiK1VFA8AmikZASkW9JfTEAGzNubwFwLoBBAPeUaX1EFgN9Dhwc9SIaL24MM+oOYcYXSe+Pi7gserzjHBd2H1qaXi+nkAOCqc3kQhibu634n5efk/MxaXcuiRG5T0qxW3rvVprA7BtyY3ufQ7YRTCZGnRptLXoMzcnzjc9keiECl0WvmO+9FES/9fmM919q0V0pdrhSmPCECha6AWTTSjQnRYWcc3455/xyAN8FcK14O/X1bs75RznnJ8q/VAIQhDwaT+L1iYWij12cpb18r/e6LV04OxfE0fHF84hCbjeWJ5K5dH0rjFo1vnzbNmjVuf/0nHnan/IRkLRHLj3luxCO4cS0v6RCN5E+p6mk1PqMv3ileDkQ53qLpGeiS1iLsEddvoi8UKEbQINTiOZE8h455/zPOOfF1YMoK4sFb8XT6/uHPTDp1NjQuTxFfc3mTqhVbMm4TW8oBsYKR7elcOO2bux/+Gps7LLmfYzDpANj0iMqKUJu0Kph0KokRYoHhj3gHCUVuon0O0153eykIJqwVBqnWYfZwPKIXMq8AWFwivIimkiZwRQqdAMyR5mSkBPNg5z2MwNj7K8ZY79ijB1gjB3K/CrnIolFumxGdNsMkvbJBSMYGzQ5ol+nWbcsvb4QiqFFrykppVwIxhgM2vx72YCQnnWYpJvC+CW0nwHiiM3ikeLgsBuMAdt6l+/fy6XXacK4NySrfz0TKb3b5aDVosd81h65w6RdVpiYi3KNMp31RxBP8qJCLmZf5qkFjWgi5FSt/z8ADwA4C+DnAH6a9UVUiO39DuwfLtyCFoom8Pr4QsHI8vqtXRieD+LImJBo8QSjipvBrAQ5ph6BSBxGrbroPrLUgR6Dwx6c19GCFkPp70Of0wTOgTG3/Kg8keSY9UerFpFn1ihML0j3e7ebtLLtcKUg9pDnG5giotOoYDVollyIEESjIyeHejOA2zjnvynXYghpDPQ5sPvQBCa9YXTmiVAOjXoQT/KCe73XbOrE3/7sCHYfnsCWHltZ7FlXgivPBK5c+CNxWCRsBUgZZSoYwbhx/dZuSa9djH7XYgva2jZ5FfCiKU4lfdZFXBYdgtEEQtEEjDp1aq9e2jrKFZFPpGa75/t7z6TVoqfUOtFUyInIgwBGyrUQQjqFhpCIDKYi9lyFbiIOsw4Xn9uK3YeF0aY1I+QWneRRpv5IomAPuYgUgTk544cvHFdkfxworZd8ekF6gZnSuJa0AMqLyFdihysFKa5uItlV9wTR6MgR8i8B+BRjjPzZq8zmbht0GlXBgrfBYTfWtJoLmmcAQvX6yHwIh8e88IZiZatYl4PLrJe8Ry5MPiu87w6IKd/C5xTfz+x2vZXS1qKHQavC8AqGp8z4K28GI5JZ+c05l1U9Xy53twlvCHqNKu3SVwinWfqFIEE0AnJE+WoAtwM4wxjbwxh7OvOrTOsjcqDTqLB1lS1vRM65kCLeLkGQrtncAY2KYfehCXhDcVhrICJ3mnXwBGOIJ4r3yvsj8aKFboAQkXtDsYLWnYPDbjhMWqxpNctabz4YYytuQZtekGbCUg6cGe56C6E4ovGkjIi8PO5u4xLMYERcFvkz7QminpEj5LMAfgbgPwFMApjL+iIqyEC/A0fGFnKmMEfmQ5j1RyWliO0mHS5Z14pnDk1goUZS66I717wEMQhE4pJS63aTFkku9InnY3DYg+19DkliIZWVCnl1I/JFd7QZf1jWOtK+9goXvE16w0UL3URcqWJJ8lsnmgXJxW6c8z8r50IIeQz02fH4C0kcGVtemb5veD71GGl7vddt6cJfvSl0ENaCkDszUrvFIlLpxW5ipBhLp38z8QSjODntx80XKFPoJtLrNOGlU3PgnMu6QJheiKBFr4FRV3zbQGnSE9D8kYy9emmZAbl2uFKZ8ITwtrUuSY91WXRIcsATihXdWiKIRoD2u+sUUaRzDVAZHPLAotdgfUdxr3JAqF7XqgWRqQUhlzM4RdgjlyDkRUZs7h8RigOVcHTLpM9pQjCakLznLzLjl15gpjRmnRo6jQrzgajszIBcO1wpJJIcU74IuuzSLibIppVoNiRH5Iyxw8g9/QwAwDnfqsiKCEm0Ww1YZTfm3CcfHHZjW69Nske3zaTFJee24ndvztSGkGekdovhl5hadxSJFPcPuaFiwLZeZQrdRMQWtKG5oCRnNJEZGZXiSsMYQ6tZMOVJR+TW6u2Rz/giSCS55NR62q8/EMU6xVZBELWLnIj8P7DUAOZpAMMAelM/ExXmwn7HstnkwWgcb0z6ZEeW16V6p2shFZmZ2i1EPJFEOJaUXOwG5N+7HRz2YEOnVVJ0L4eVTkGbXKiOPauI06LDnD+CGX8Eeo0KLRLfFzl2uFIZlzCHPBNnjpnqBNHIyNkj/7tcxxljfwWgX7EVEZIZ6LPj6YPjGPeE0v21B0e8SCQ5BmT2Qt98QTc0KoaL1jjLsVRZ2I1aqFjxwReBqFDoJ6X9rFCkmEhyHBjx4Obtyu6PA0CPQ34vuScYxfB8ELe/pVfx9UjFZRYqv6cXwmi36mXt70u1w5XKpERXNxFXeqY6pdaJ5kCJPfInAfw3Bc5DyEQU630Z/eRiqn2gV56Qa9Qq3Lx9VUVHZuZDlfJbny0SUfkljDAVaTFooGK5+5uPT/ngj8QV3x8HhAi102rAkIxe8nLt18vBZRbe/xl/BG0ytgQA6Xa4Uhn3yIvIxYs2uXUJBFGvKCHkl0JwfSMqzMYuKwxa1ZJ98sEhN85pM9eEZ3opCKM0C0dUUiafiahUDPY87m7pi58yCWef0yQrtT445IZaxRQZ3LJSxFGm0wvS7VlFpNjhymHSG4ZBq4Jd4t+0Vq2Czail1DrRNMgpdsveB2cAugBsB5Az7U6UF61aha2r7Gk7Vs459o94cOWG9iqvrHSkDE7xyxByQByxuVxgBoc8cJl16cI0pel1mvCHk7OSHz847MaGzhaYJOz9lwunWY9QLIERd1By25eIw6zDsQnlJh5PeMPothllpfezZ6oTRCMjJyLPNoCZBvAbANdyzj9XhrUREhjod+D1cS/CsQTOzgUxH4jK3h+vRVwWfdGISozIpVStA/n91gUXPGWNYDLpd5kwuRBGOFbcfzyR5Dgw7KlqWh1Y7BwIx5Ky/d4deS6YVsq4NyRpWEomwtYA7ZETzQEZwtQ5A312PJbgODLmTe/DKjX0o5q4zMVnkgckziIXcZi0GPOElxxzB6I4PRvArTt6VrZQCYiV66PuIM5tL9zbf3zKh0A0gYF+Zdvg5CIWjAHy3eUcqT3yZJIrMtt+0hvGO85plfUcl1mPUzP+kl+bIOoB2XvkjLG1jLHrGWPXMcbWlmNRhHQyC972DbvRYtDgXJkjM2sRl1kPbyiGWAG/dV9YerEbkLsIa/9IeffHASG1DkirXBcLFy/sq273QGYbotQechG7SXBWE/99SiGeSGJqISy50E3EaSl+IUgQjYJkIWeMWRljTwA4CeDnAJ4CcIIx9hPGmDQLMUJxWi169DlNGBx2Y3DIjQt67YpEQdVGHNxRaGKZnGI3QCzCWnq+fanCsq095SssyzSFKcbgsButFh16ndJarcpFpnlNm0V+sRsgzSu/GDP+CJIckl3dRFrNuvRMd4JodORE5F8FsBXA5QCMqa8rU8f+RfmlEVIZ6LPj1TPzOD4l3wimVmlNRYSFWtDk9JEDQqQYjiURii7uVQ8OebCpy1rWwjKXWQeTTi0pIt9fhsEtK6GUiFxJd7fx1FZIt8QechGnWQfOlZ/CRhC1iBwhvxHAhzjnz3POY6mv3wP4CICby7I6QhID/Q64gzEkORqi0A1YFJJClcf+SBxaNYNeI03IswUmnkji4KhHsfnj+RDHmRZrQZsPRHFmNlATF2MmneDQxthi4ZtUxDYxJXrJJ1KubrKL3SyLg3cIotGRI+RG5B5XOg+g8kOTiTTiBz9jwAUKe4VXi7RNa4FecqkDU0TSIzZTAvPmlA/BaKIiFz99TlPR1Pr+dD979f8NGWNwmfVwmXXQqOWV0hSzw5WD6OomNyJ3kU0r0UTIySf+AcDnGWN3cc6DAMAYM0PoIX+pHIsjpCH0HKvR4zDWxNATJZDyQewPSxuYIrI4YlMQGLH/vhIRcJ/ThOePzxQcZ7pvyA2NimFrT/WFHBAq16Px/MWG+ZCTWv/a70/h4Ign7/3Hp30watWwGuVtfUi5ECSIRkHO/45PAXgOwBhj7BCESWjbILi6XVOGtRES0ahV+NA716JD5l5mLWMzaqFWsYIfxFInn4lkjzIdHHKj1aJHj6P8hWV9LhMi8SRmfBG0W3MnsAaH3djUba3KDPJcvG/7qhUVixWyw81k2hfGl375BjpaDHkvQLUqFW7b0SO7ZkDK1gxBNApy+sgPM8bOBXAngA0QnN2+D+AHnPNQmdZHSORTV6+v9hIURfRbL/RBHIjKS60705GiGJG7cWG/vSKFZWIv+dB8MKeQxxNJHBzxVnVQSjZ/dvGaFT2vkB1uJr88MgnOge/dexHWdyjb+CJuoxTz6yeIRkCORevfAxjhnD+WdfxjjLFVnPOHFF8d0dSIgzvy4Y8kZG0lpFPrgShm/REMzQXxpxf1lbxOKYhCPjwXxFtWL+8Rf2PSh1Asge01sD+uBPnscDN55tAE1rVbFBdxQMhSOUxamoBGNAVyqljuArA/x/FBAP9dmeUQxCLF/LIDkTgsElvPAECnUcGsU8MdjGG/uD9eoSr/VQ4jGMtvCrO/zINbKk0+O1yR6YUwXj07j+u2dpVtDU6zjordiKZAjpC3A5jJcXwWQIcyyyGIRYoNTpFb7AYsursNDguFZVtWVWbCmF6jRrfNmFfI9w250d5Smf36SlBsAtqeVFr9ui3lE3KXRU/ubkRTIEfIhwG8M8fxSwGMKrMcglik1aIvOPhCbvsZIBS8uYNR7BtyY3O3FQZt5QrLep35hXwwNSil2kYwSlFsJvnuQxM4r6MF68qQVhdxmXWYo8EpRBMgR8i/DuCfGWMfZoydk/r6CICvAHi8PMsjmhmnWQdfOJ6zBYpzjkBUfkTuMAn77odGPRU3z+lzmnIK+aw/guH5YNUHpShJLjtckamFMPYOzeO9ZYzGARplSjQPcqrWv8IYawXwKADR6ikK4Kuc8y+VY3FEcyNO4JoPRJc5e4ViCSS5dJ91EbtJh5dOzSGR5BXfj+53mTHjG0UwGl9iCTs41Fj748CiHW44lliW9dhzeEJIq2/tLOsanGY93MEY4omkbFMbgqgnZP11c84/A6AVwNsAvB1AG+f8gXIsjCDSpjA5Ko/9MgemiDhM2nRvdKUjcnEK2sj80m7NwWEPtGqG8yu0X18JxD7uXFH57sMT2NDZUnSka6m0Wpa2GxJEoyL7MpVzHuCc7+Wcv8o5lzzwlzH2bcbYNGPsSMYxJ2Ps14yxE6nvjROSECWTdufKUXnsF0eYriAiB4AOqx7dMv27S6UvzzjTwSE3NnfbKrpfX27SdrhZNq2T3jD2nnWXtchNhExhiGahkvmmXQDek3XsAQC/5ZyvA/Db1G2CAFD4gzgQESefyY/IAVSlsKxfNIWZC6SPxRJJHBrzNFRaHci0w136b/fs4QkAwHvL2HYm4kzb/FLBG9HYVEzIOecvQBiwkslNAL6b+vm7oClqRAatZiEiz1W5vphalxfFij7gF1ZhSpzdpEWLXrNkCtqxiQWEY8mGKnQDFt/n7JnkYlr9nDZL2dfQmvZbp4icaGyqXQHSwTmfAIDU9/Z8D2SMfYQx9hpj7LWZmVzt7ESjYTVqoFGxPBG5IORyq9ZXt5qhYsA7zmlVZI1yYIyhN6tyvREL3YDMSXOLqfVxTwj7hty4vgLROEAROdE8VFvIJcM5f5xzvoNzvqOtra3ayyEqAGMMjjzuXIHoyordLui1Y9+DV2NTt1WRNcql32XCmWARfgAAIABJREFUUKaQD3vQaTWg294YRjAimXa4InuOTAJA2dvORBwmHRijPXKi8am2kE8xxroAIPV9usrrIWoMl1mXMzXqW2GxGwA4zLriDyoTfU4TRudDSKYq5/cNuauS5i83mXa4IrsPjWNTlxVrK5BWBwB1avDOLAk50eBUW8ifBnB36ue7ATxVxbUQNYjLosvZfhZYYftZtel1mhBNJDHlC2N6IYwxT6hhBqVkk+nuNuYJYXDYU1Zv9Vy4zDrMk9860eBU7FOQMfZDAJcBaGWMjQL4LIAvAvgJY+xeCBawt1VqPUR94DLrMer2LDseiMTBGGCqkdndUul3iZXrwbTIVbqfvVKIdriAYAIDVC6tLuI0574QJIhGomJCzjn/QJ67rqzUGoj6w5knovJHEjDrNHXnTZ7ZS35y2g+dWoXNVdqvLzfCBDQhtf7MoQls7rZiTau5omtotehxbHKhoq9JEJWm2ql1gihIq0UHXySOSDyx5LgwMKW+onEA6LYboVYxjMwHMTjkxvmrrNBr6u/3kIKYWh91B3FgpPJpdaD4BD2CaARIyImaxpnqJc/+MPZH5A9MqQW0ahW67QacnPbj0Ji3IQvdRMRRpnsOC9XqlXBzy8Zl0cETjCGWWD54hyAaBRJyoqYRB6dkt6DVq5ADQnr9+eMziMaTDdc/nondpMNCOIanD47j/FVW9Lsqm1YHFv36801iI4hGgIScqGkWB6cs/SBeySzyWqHPaUIwKmwVNGqhGyBE5JwDh8e8uG5Ld1XWIGZ0cnkREESjQEJO1DSLg1OWVh7761rIhch0ld2IDmtlB7dUEtGmFahOWh1YOgqXIBoVEnKipsk3OCUQre/UOoCG7R8XsadsWrf22NCXarurNGJGJ5dfP0E0CiTkRE1jNWigVTPMZu+Rh+tXyMVe8kbeHwcEDwCgzL3jE4eA/7MDCMzlXkMqo3NmNpDzfin8+yvDuOH/vAg3RfVEjUJCTtQ0jLFUC9HSiCoQSdRtan1ztxWP3Hw+/uQtvdVeSlnZ3G3F52/ajDvf1l++Fzn7IjB3Apg4kPNup1mHS9e34evPn8bwXDDnYwpxctqPnb84isNjXuz8xdFSV0sQZYGEnKh5nGb9kmKlaDyJaCIJSx32kQPCxcmdb+uv24yCVFQqhrvevrq8v6f77NLvOfji+7dAo2L4q/84mPa4l0IiyXH/Ewdh0qlxzztW46kD43juyERp6yWIMkBCTtQ8rZalg1Pq1WedKAMShLzbbsRD12/CK2fm8b2X8z8um2/812kcGPHgczedj7+9biPOX2XFgz8/QoVzRM1BQk7UPNl+2X4SckLEfWbp9zzctqMHl53Xhn987k2clbBffmLKh3/61XG8Z3MnbtjaBa1ahS/ftg3eUAwPP3VEiZUThGKQkBM1j8usX+K3Lgr5SkaYEg1EMgm4h4SfC0TkgLCd8YX3b4FGXTzFHk8kcf8TB2ExaPDI+85P+/lv6LTivqvW45lDE9h9iFLsRO1AQk7UPC6LDoFoAuGYYKJCqXUCAOCfBBIRQGsSBJ0X3v/ushnx8PWbsPesG9956Wzex339hdM4OOrF5286H62pqneRj166Flt7bHjoqSPU0kbUDCTkRM2T7e5GqXUCADCfSqf3XwxEFoDgfNGn3HphD67Y0I7//cs3cHrGv+z+Nyd9+OpvTuC6LV05h7xo1Cp85bZt8IfjeOjnR8CLXDwQRCUgISdqnrQpTCq9HogIkXmjV30TRRDT6edcvvR2AcQUu06twl/9xyEkMlLssVRKvcWgwedu2pz3HOs6WvDJq9djz5FJPEMpdqIGICEnah7R1GM2VfC2mFqvz/YzQiHcZwGmAlZfkrpduOBNpMNqwM4bN2PfkBvf+cPic77+/CkcHvPikZvPT//N5ePD71yDbb12PPTUEUz7wiv9DQhCEUjIiZrHlRWR+9LFbtqqrYmoAdxnAFsP4Dp38bZE3rd9Fa7a2IH//cs3cWrGj2MTC/jqb0/ghm3duFaCE52QYt+KYDSBB39GKXaiupCQEzVPepQpReREJu6zgGMNoDMDlg5JqXURxhj+4f3nw6hT4/4nDuL+Jw7CZtThczfmT6lnc257C+6/Zj1+9foUnj44Ln/9BKEQJOREzWPRa6BTq9LFboFIHHqNCho1/fk2Ne6zgGO18LNj9WIrmkTaWwz4uxs3Y/+wB0fHF/D37zsfDrOu+BMzuPeStRjos+Php45ixkdV7ER1oE9CouYR/dZFm1Z/pH4HphAKEfEDgZmlQj4vPbUucuO2btzzjtX46KVr8e7NnbKfr1YxPHLzFnhDMfzm2JTs5xOEEpCQE3WBy6JLW2MG6ngWOaEQYho9U8gXxoC4vKiYMYadN27GZ967ccVL2djVAqdZh8Eh94rPQRClQEJO1AVCRC58SFNETqSF3LlG+O5YA4ADnpGKL4UxhoE+OwaHSciJ6kBCTtQFrRb9EkMYEvImJ1dEnnm8wmzvc+DUTACeIA1UISoPCTlRF2TukQuzyKlivalxnwEMNsDoEG6nhVz+PrkSDPQJ69g/7KnK6xPNDQk5URe4LDqEYgmEognaIycWW89EWjoBjaFqEfm2XhvUKkbpdaIqkJATdcGi33qEUuvE0tYzAGAs1YJ2tirLMek02NjVQkJOVAUScqIucJkFy8w5f5SEvNlJJoSe8UwhB1bcgqYUA30OHBj2LPFvJ4hKQEJO1AXOlLvbrD+CYDRBqfVmZmEcSMYWK9ZFHGuEiLxKdqkDfQ4Eogm8OemryusTzQsJOVEXtKYi8pH5IACafNbUZFesizhWA7EAEJit8IIExII3Sq8TlYaEnKgLxIh8KCXkFJE3MWJlei4hz7y/wvQ6jWi16EnIiYpDQt4sxEJVMctQCrNODZ1GheE5Ucip/axpcZ8FVBrA2rP0uJhql1vwFpgDgvMlL0s0hqnXFrRwLIFxT6jayyBWAAl5s/CbvwMeu0QoFKpDGGNoNevSEXmLgSLypsV9FrD1AuqsvwF73+L9cvjxncDP/4cSK8NAvwNnZgNpO+F64vEXTuOaf34B0Xiy2kshZEJC3gwkE8CRnwJhT9Xac5TAadGl98jNOhLypmX+zPK0OgBojUBLl7zK9VgIGH0VmHpdkaWl98nr0Hf90KgH/kgcZ2YD1V4KIRMS8mZg6CUgMC38PH2sumspAZdZj0gqWqA98ibGfXZ5xbqIWLkulfH9QDIOLIwC8dKj6K09Nmjq1BjmxLQ/9Z2q7usNEvJm4OjPAI1R+HnmjequpQRcGbOiqWq9SQl7gdB87ogckG8KM/Kq8J0nAW/pNSQGrRqbu611J+ShaALDqWzX8Sl/lVdDyIWEvNFJxIFjTwPnvQew9dW3kFsWhZwi8iYlX+uZiGM14BsXUuZSGN0LgKXOrUy1+/Y+Bw6OeBFP1M9e86kZf7r9/sQUReT1Bgl5ozP0ByAwA2y6GWg7D5iuXyF3pnrJASp2a1rSQp4ntS6m3D3Dxc/FOTDyCrDm0qXnLpGBfgdCsQTeqCNjmOMp8V7bZk7/TNQPJOSNzus/B7QmYN01QPsGYPZ43VauixG5WsWg19CfblOSFvL+3PfLGWfqPpu6yL0RUOuVE/I+O4D6MoY5PuWHVs1w9aYOnJ0LIhKvz8+IZoU+DRuZRBx4/Wlg/XsAnQlo2wgkIlX1oy4FcY/crFODMVbl1RBVYf4MYHQKI0xzIQq5lL/x0b3C9963ChcGCgn5KrsR7S36uqpcPzHlw9pWCzZ1WZFIcqpcrzNIyBuZoReB4Cyw+X3C7fYNwveZ+qxcd1mE1DoVujUxhSrWAcDcBmjN0kR55FVAZwHaNwmp+nkJz5GAYAzjwGAdGcMcn/ZhXYcF69pbhNtU8FZXkJA3Mkd/JnyorbtauN16nvC9TvfJ0xE5CXnzkj2+NBs540xHXwVWDQAq9eJzFBq4cmG/A8PzQcz4Ioqcr5wEo3GMzIewvqMFa9vMUDEqeKs3SMgblUQcOPYLoVpdm2o901sE96u6jcgFIbdQoVtzkogLLWKFhBxIiXKR1Ho0AEweAXouEm471wBRnyJWrQAw0F8/++QnU/3j6zssMGjVWO2igrd6g4S8UTn7X0BwbjGtLtK2sW4jcpNOA4NWRan1ZmVhVDBvyVexLuKUMM50bBDgCWF/HFB84Mrmbhu06vowhhHT6Os6WlLfLWlzGKI+ICFvVI7+TNj/O/eqpcfbNwBzJ4Topg5xmfUk5M1KsR5yEcdqIB4G/FP5HzOaMoLp2bH0nAoVvAnGMDbsH6rMPvmsP4L3/78/4KVT8ke4npjyQadWod9pAgCs72jBUI1WrvvCMdz6tZewb0iZzEmjQELeiCRiqbT6tYtpdZG2jUAiCsyfrs7aSuThGzbhw5eurfYyiGogR8iBwpXrI3sB1zrA5BRu21PtbAqOQB3oc+DQmAexMhvDcM7x0M+PYHDYg92HJmQ//8S0H2vbzNCoBTk4t92CRJLj9EztVa7/+vUpvDbkxn++MV3tpdQUJOSNyJkXBBvL7LQ6UPeV6+/e3JkeTEE0GfNnALUOsHYXfpyjyDhTzoWIvPeixWM6E2DpVKxyHRAK3sKxJI5NLCh2zlw8c2gCe45MQq9RrahS/viUL51WB4SIXDxea4gXKlRVvxQS8kbk6M8AXQtwzpXL72tdL3yfebOyayKIUnGfFYo1VUVm0dt7AbD8Qj5/Wqgf6XnL0uNyfdqLIBa87StjP/mML4KHnzqCbb12fPida/Hm5AL8EenbZoFIHKPuENa3W9LH1raZoVYxnKgxsfSGYnjhxAwAqqrPhoS80UjEgDeeSaXVDcvv15mFNGIdT0EjmpRirWciGj1gXZVflMVBKWKhm4hT5uS0InTZjOiyGcrWT845x4M/P4xANIEv37oVO1Y7kOTAoRHprydWrGdG5HqNGv0uU81F5L95fQqxBMel69swNB9EOFZ7e/jVgoS80TjzPBBy506ri7RvrOvhKUST4j5TvGJdxLkm/3736KuA3gq0bVh63LEaWBgD4sr1fg/0Ocrm8Pb0wXH88ugUPn31eqzraMH23tQsdBmV8qJYr++wLDm+vr2l5irXdx+ewCq7Ebdd2APOFy9CCBLyxuPoz4QPqf/f3pnHR1Wf+//9zR6yQIAkBAgkQEIgsggI7oIbS1CJW9t7+2tvrVdbly63VmtrrbWt1dpFa62tba29Xe1V2RFQVNyQXYGw7wFCAkkICZD9/P545pBJmEnOLMlkyPN+vXglM2eZ7xxyzvN9Pt9nGX61931S8+D4LvHeFSUcOFMpLUydeOTQfsnV4rUwaCJEtHn8pWQBlrOGKw65cEgfDp84Q+nJ2qCdE6DsZC2Pzi/iwiF9uPMKCf7s3SuaEWmJPkn5u8pqiImKYGi/hFbv56YncqD8VLfxeqtON/D+rmMUjM1g5ABRD7RvegtqyM8nmhpg2yIYOcuzrG6TNgqaG8I2cl3pgTiNWLdJyZL0s/rTrd+vq4ayotaBbmeP6SBIzg8mDnV5yUH0yi3L4rtzt1Db0MTPbxtHZERL34GJQ1LYWHwCy2GFup2l1QxPTWx1DhCpvdmS9qbdgeVbj9LQZDFrTAZZ/RKIijAa8OaGGvLzib0rofYE5M9pfz9bUtR1ciVcsFPJHBtyL0b58HqwmlsqurU6Jqv1ZwWB/IG9iYmKCGphmHmfHOatbaU8cP1Ihqe2lsQnDO3DidMN7HXY9GRXac05sjq0RK53l4A3W1YfN1iuZ3b/BA14c6NbGHJjzH5jzGZjzCfGmHWhHk/Y4kRWB1fkutF1ciV88Nkj92LIi10dz+xCMO4kpknL3yB65DFREYwZ1DtoAW+lJ2v5wfwiJg5N4Y7Lz40XsFMznSgANXWNHD5x5qzRdie7v3i93UG+rjrdwAe7jjN7bMbZroe56d1vDT+UdAtD7mKaZVnjLcvycIcpHdJYD9sXQl6BRO22R0wveSCqR66EC5X7pbNZ7Lneo0e8VWo7tEYUqfg+5x7jS8MVH5gwpA+bD1dR3xhYYRjLsnj49c3UNzXz9K1jz5HDAYanJpIcF+Vo4mB7tDlp517TmKgIsvondAv5etnWozQ2WxSMzTj7Xk56IgcrTnOmvnus4Yea7mTIlUDY+64EA7UXre6ORq4r4UTlPufeOEjFtpik1pHrliU9yNvmj7vjpOGKj0wYkkJ9YzNFR6oCOs9rGw7z9vYyvj09j2Gpnic0ERGGCx1Gyu8qPTf1zJ2ctMRuIV8v3lRCZt94xgxq6UGfm56E1Y3W8ENNdylabQHLjTEW8HvLsl4M9YDO0lgHa/8Ik+44t9ypUw6tg+oSGHWDf8c3N8GaF6HqkPd9DnwEsb1h2DRn50zNg13LxZOPinF2zPqXJdq9PYZeBnmznJ1P8c62hVL8JGNcqEfSPajcD5kXO9/fGOib1dq7Lt8t0e+eAt1sUrJkUmxZco4gMMEV8Part3Yx0sN6tBMsC15ZV8zkrL586dIszzusewnyCpgwJIVnVuzkZG0DyXHRXs+5q6ya2KgIhrhqrLclJz2JpUVHqW1oIi66gyI8nUTlqXo+3H2cL1+RfVZWh5Z0uZ2l1VzgZuB7Kt3FkF9mWdYRY0wa8KYxZrtlWe+572CMuQu4C2DIkCFdN7Itr8Gy74qsN/Z2/87x5g/g4Cq4801Je/GVVc/Dm9+X9TvaebhMucu5UU7Nk05SFXvEO++IEwdh4dchMhYivPzZNNXDxr/CA7udj0M5lwOr4JX/Bwn94Z7VkNAv1CMKLU0NMokdm+XbcSlZrSsY2oVgPAW6nT0mGxpOw6ljsmYeBNKT47hkWD/W7a9g3X7/m32kJcXys1vHEuFBUufYDlj8P3B8JxNGPIBlwafFJ7giJ9Xr+XaW1niMWLfJTU88m68dKmO53CWrzx7Tuizv0H4JREdq5LpNtzDklmUdcf0sM8bMBSYD77XZ50XgRYBJkyY5y60IBkXz5GfxGv8MeVMDHHG1TJx3D9y1sv3UsLYc2wFv/xjyZsNn/hY0L+FszfWybc4M+db58vPej6Gvl6YlO96Af35WitLkXBeccfY06k/B/HsgKUOMyRvfhltfCvWoQsuJgxJp7ou0DrL/zuXQ3Cw548WrIa53S5lib8eAePJBMuQA/7zLBzXBH+xublvnM37q4xgDGw60b8h3lVYzObuv1+1nI9fLQuf1Lt58lCF9e3HBoORW70dHauS6OyFfIzfGJBhjkuzfgeuBLaEdlYszlbDnbfm9eLV/5yjdIjP88Z+XNemVTzo/tqkR5n1VyqrO/lXwjDjIw8xEOF8nL5orMq83Iw4SLR+bLPsq/rHiccnvv/lFmPqQKEL2JKqnYsvjfR1WdbNJyYamOlnWgpb18baFYFodkyU/g5iC1iXYz6fqEpLKNpCblsT6dlLeqmsbOFJV63V9HAh5vrYtqxe4Rau7k5OexM5uEFXfHQi5IQfSgQ+MMZ8Ca4DFlmUtDfGYhO1LpHDK8GugtEi8JV+x012mfgcmfAE+fFbWzJ2w6jnJe531dFC9A0DW+1OynBnyygMyjo4C6aJiJWp+2yJZe1d8Y/8HsPp3MPluyL4CLvsmZIyHRf8Dp3zvM33e4GvqmY27d11bJepTe7I6SFxCew1XuivFayU+JSoOiuYyYWgKGw9W0tzsWby0U7c8pZ7ZhDpfe1nRUZqaLQrGZHjcnpuWRHHFGU7XO28Sc74SckNuWdZey7LGuf7lW5b1k1CP6SxFc+XGnnK3SOOHN/h+jkNrRCbtPRiu/wkkDRQvu6GDko1l2+CdJ2DUjXDBLf6NvyNSR0GZA0O+1bW8MLqDQjMgxr6uCva+E9jYehr1p2D+veJFXvsDeS8yCua8AHUnYfG3Qju+UFK5X2IzEgf4dtxZQ75PJqJYkNlOxDrIslfywPAy5Gcq4fgOCXTNuQ62zmdCZhLVtY1eo7p3eamx3pbc9KSQeeSLN5cwtF8v8gcme9xuj11rrncDQ95tOV0hxii/sCVdxV6H8oXiNXK8MRCXDDc9B8d3wjvtzFdsST02CQp+GVxJ3Z20PAl268h7LpoHAy90Jm0OmybR8yqv+8Zbj4nyMee3spRikz5a1Jyt82DL6yEbXkixU8/ak8Q90WeILB9V7ncpYwYGOShT0QkpaJ3KofXyM3OyPK9qjnJpzG7AewOVnaU1xEVHkJniOWLdJic9keLKrs/XrjhVz0d7yikY41lWl7F1r+pzoUQNuTd2LJGo7tFzJCe1f25L1KtTqkvhxIHW7RKHXw0TvwQfPef9fB8+A0c2QsEvINF7sErApI6S71i+2/s+lfslWM9pfnpUDIyaLcsSQewidV6z7z1JL5zyFRh66bnbL/06DJwgXnlNWdePL9Q4bV/alshoUcIq98sactoomUx3REpw25l2OofWyIRl0ETImQ5RcQw8tJQ+vaLZcMBzYZidpdWMSEv0HAHvRqjytc/K6mM9y+oAWf16ERMZoevkqCH3TtFc6ds98EJ5PXiyBMs4bEYAtHjwbfNWr/8R9M50SexnWm8rLYJ3nxTD6dR4+osduX6snQpvdtT+6Jucn3f0HJHX96i83iF1NSKp9x0G1zzqeR9bYq+vkRQjX/4Gwx3LEqXCH0MOclzFXolLaS9/vO0x1SXn3pvdleLVkJYvVe9iEyHnesy2BUzMTPYa8LartIbcNO/r4zbu+dpdyeJNJWT3T2B0hveJV1RkBMNSE9QjRw25Z05XSFGI/MIWWTvzIjhd7lvHsOI1EBlzblGP2CSR2Mt3S2qZTVODGPf4PjDrFwF/jQ7plyMz+fbWyYvmijfoy4N02FRJ81F5vWPefBROFIuhjmlH5kzLg2nfk0IxW17ruvGFmtMVEiPga8S6TUq2qFt1VR0Hup09Jkt+Vh7w7zO7kuYmkdbd1/7zC6GmlILe+9ldVkPV6dbtiqvONHD0ZPsR6zahyNcur6njoz3H25XVbXLSk7p8ktEdUUPuie2LRHJ294jth4Av8vqhtWLEPdU+HzYVJn1Zir0c/Fje++BXUPKprIt3RRGQ6DjxBL155BX7oOQT35WBqBjIu0GWJzoK6uvJ7H0X1v0JLrkXhjjIM770fom3WPKALNv0BPyNWLdJyZIcdHDukduThnCQ149th/rq1pOU3OkQFc8lte8DsLG4tVe+u8x7jfW2hCJfe2nRUZot2pXVbXLTEjlUeYZTdT07cl0NuSeK5spM3t2TTs2THGmnAW+N9eIJtOcFXPc49HFJ7MVrYOVTcMGtMPrGwMbvC6l53j1yO1q9o7aonsgvFE/KzsNXWlN7EubfB/1GwNWPODsmIhJu+q302F70zZ4hsdtBZ4EYcoD4FLnWvhwTDoa82MPyXUwC5E4n/fByokzzOQ1UbCm6vdQzd7o6X3vxphKG9U8gb0DH48vRyHWgm1R261acrpC+3pd9rXW0eESEtD506pEf3QyNte17AbGJ8mD+y2x4uQDi+0rOeFeSmicV2RrrzlUOiuZKlG8fP0riDrsK4vrIZKAra6831omX0p6Ri+/jv2EIFm9+H04ehjuW+1bDPzUXrvk+LH8EVv/emSfvjchoSBvtf1aEZUnqUy/v1cECxjbkfYb6d7ztXduZI07o1Q9iEn2LXC/fA3UBGrvemb4rcYfWynjbFmrKn0PE1nnc3O8gGw+2rkGxs7SG+OhIBqc4+7vLTUti8aYSTtc30iumfZNRcaqexNgoYqL88xGP19Tx8d5y7p02okNZHVoi13eWVjMu00NHux6CGvK2bFsoOeOe5OTBk+G9n8kNG9vBbNFboFtbsq+AyXdJ1PINz3TuQ9ETaaPk+x7fBQMuaHm/fI/I/Nf7mdYfGS1NYormibzuS1naQHjzUSmq0hG3vdz5wYTe2PO2NKC59Gsd5zV74uJ75O906UOBj+WiOyU7wh+WPgzr/wx3rmj9txMsLEsm1cmD2o8faI+UbIiIlmIpTjHGt8j1Q+vhj1f7NbxWpGTD1zb6NrEqXiPPpbbH5FwP0b24JXYNXz44gqZm62xN9V1lziLWbdzztccO9m4sD1WeZuaz7/OZSZk8Mnu08+/gxtItzmV1gKF9JXK9p/cmV0PelqK5MrsdMPbcbZkXyXrb4fWyxt0exWsgebAUl+iIGU/KerkdRd6VpNqR69tbP4y3+hGt3pb8OdJEZc8KqfjW2TQ1SiBY9lWSyuWN934m1dKGXhb8inkdUVsF8++XdMZp3/PvHBGR8PnXYf/7gcnrOxZLZ7+8AkmL9IV978HqF+T3eV+B/35HJm/BZP3L8h0Lfun/OeL7wFfeb7+0sCdShraflunO5v+TgjW3/gmMn13CDq2FD34pz5bBDnLdQdTD8l0w/nPnbnPJ6+N2vceZulvYVVZN3gCJAN9ZWs1lI/o7HlqL1+vdkFuWxUOvbaK6tpH5nx7h4VmjvDZjaY/Fm0oYnprASIeyvx253tMD3tSQu3PquDygLv+G51mxXUyieG3HhvzQWufeVkRkaIw4QP8cefiUtQl4K5oncmSfTP/PnX2VrE0Wze0aQ37gQ2k0MumO9uX8vsPg91dIKtftf+28gjueWPY9qD4CX34rMJUiNhFGzgxsLMOnwcHVMrG4Z5WzHGtwS5kbDlMfhtfvhPd/IYVrgkXlAVk+GDZV/j8DwUlToLakZMHut1oarnijuVkmvTnX+d+mGKR+wKrfyL3i1JAfcpV/9haHk19IXNFcJkdsZ8OB8eQNSKbqTAOlJ+scr49DS752ewFvf199kA93lzN1ZCrv7jjG2v0VXDzMt2WCsupaVu8r576rcxzJ6ja56Umsd9B//XxGg93caU9WB5ndp47qOODt5BGoKnae7hJKomJdketuAW/le+DopsClZ1te3/FG1+TkFs2VVq8517e/X6hSuXa9JQrFZV+HwX60sw020fGS9lZ9RIymU86mzP0Wxt4GY26H955kcX8pAAAft0lEQVSGkk3BGVdzMyy4DzBw43NdO9Gy6ZstMS41HWQHFK+WnPNA75X4PtLToWiefH8nFK+RSfigCZ63j7gOKzqBW2LXnK3w5rQ0qztn87W9yNfFFad5Ysk2Lh/Rn+f/YwJx0REs3lTi+Pw2y2xZ3UttdW/kpidy+ETPjlxXQ+7O1nniZaS3s96XeZHcQO3dbGcjSad436c7kZbX2pDb+d+ByOo2+YVSyGT3isDP1R5NjbBtAeTOcLae2tWpXGdOwIL7ZSI49eHO/zynZF4ka/Ub/iIeaEfseefclLmZT0nA1byvBqdZzvqXRBmb/mP/Ai2DgdPI9aK50qgkd3rgn5lfCCcPwWGHTZUOrZHlMPeSvu7E9MLkTud6s5ZPDkjTHTsfPMdBMRh3vOVrNzdbPPjqJiKM4clbxpAQG8XVeWm8sUUqs/nC4s0ljEhL9GmSYY8N6NHr5GrIbWxZ3b0IjCcGT4baE+2vnx1aKzf3gDHBH2dnkDpKCt3YOd9F82QS0ntw4OfOulKi8Tu7OMyBD6Rgj9NUua5O5Vr2XfHu5vzWc12BUDL1YYmVWPA1WcP3Ru1JmYz0y2mdMterL8x+Rlr2vhdg1kXFPlj+qKzZT/hiYOcKhBQ7l7ydyPXmZmkxO+LajoNfnTByhhSQcnKvNDVKE6eOVL/8QpKbT5BWsY7KU/XsLK2mV0wkg/r4kCmB93ztv68+wKq95TxSMIrBrrrtBWMGcrymjjX7KhyfX2T1CkdFYM4Zm1vkek9FDbnNtgUSyNaRRGZHobcnrxevkfaTUTHBG19nkpYn3718FxzfDaWbgxfRHRklefGdLa8XzYXoBBhxnfNj7FSuHYslYKmz2LkMPvk7XP5N7zJoKImOkwlG9VGZcHjDTpmb88K5KXN5s2DsZ2Wt/Mgn/o2juVly6yMiQyep2/TObGm44o3ij6HmaPDulbjeMinYOr9jeb1sqyhdHWXF5FxHU1QvZkd8zMbiSnaVVZPjQ8T62dO4jKV7vvbB8tM8sWQ7V+am8pmLWmJppuWlEh8dyeLNRxyff+mWo1g+RKu7M6RvL2Kj2l/DP99RQ25TNFc8jfT89vfrlyM3nLd88sY6qYbmT1pRqEh1BQOVbYetQZTVbfILoeEU7HozeOd0p6kRti4Qj8bXNKWL7xH1Ycm3xZAFmzOV4ummjYarHgz++YPFoImydr/xb7Bz+bnbd69wpczd7/1ve+aTkJDqktj9aJiz9o+irEx/IjhqUCBExUjWSXuG/KysPiN4n5tfKJMlO5DNG7YjMbijtqzxWLkzmRG5ho37j7OztMZRada2tK253txs8cCrnxIVYXjy5jGtvOheMVFcPSqNpVuO0tjkbL1/0aYSctMTfQrCs4mMMAxPTQxZu9XugBpygJpjsP+DjmV1cBWGmez9RivZBE314bM+DtBvuATNHNvmktUvdpY255Shl8saamfJ6/vfgzMV/nlGtsTeWAsLvxF8iX3pwxJJP+eF7iept2Xqd2RSt/BrMgGxqa0SSb3/SJjajscenwI3/lq8xZU/8+2zy/fAWz8QReXCz/s3/mCTMtS7IW9uEs8553rJIAgWuTMkla2je6V4LSSkOSpsFDXmZvqaGiq2rOBYdZ2j0qxtGdImX/t/V+1nzb4Kvn/DaAZ6kOkLxmRwvKbekbxedrKWtfsrmOVjkJs7OemJ6pH3eM7K6g7XVzMnS7qWp/XE4tXyMxwi1m2iYsWYb18s65zBLpQSGQWjbhSJuf50cM8N8tCLSRRZ0h/6j5DOYzvfgE//FbxxbV8Cn/4TrvgWDBwfvPN2FlGxUPiCtEpd6mawl31XIrPnvNBxylzudBj/n9I34PB6Z5/b3CzpbBHRMhEIpaTuTt9sWbP3xMFVEvPgT/ni9ohLllS2rR1Erx9aI88hJ9dqxDXURcSTf0LKJfvj9brna+8/foonl25n2shUbpvoWTmZNjLNJa93HL3+hi2rB2DIc9OTOFJVS3VtQ8c7n4eoIQcxBP1zRf50wuCLAEtaI7bl0BqJtE1KD+oQO51UO3LddE6td1te3x1keb2pAbYtkpxqX0qdtmXKV2DIJfDGQ5I+GCinK2DRNyQD4spvB36+rmLghXDF/8Cn/4AdS2U5ZOPffEuZm/4EJKbDvHucNc1Z/TsxjDOfDK4SFCgpWXCqDOpPnbutaB5ExUv/72CTXygTJ9spaMup4xKc2pGsbhMdT9nAq5kRuZYoGs/WJ/eV3PQkdh6t5tuvfkp0ZAQ/vXms18C0+JhIrnEory/eVMLI9CS/JH8bW2XoqZHrasirS6WQiBNZ3WbQRMB4lteL14aXN25jF80YcknnPEyHXibrp8GW1/e5ZPXRAXpGEZFw0/OyLLLw64FL7G88JFH0c14In6BHmysflAnIwq/L+r6vKXPxfSRY7dh2WPlk+/se3w0rHhdJeZyHCmWhxFs7U1tWzw2yrG6TO719ed1+7jjt5gbEjbuVvqaGaTHbfY5YPzus9ESOVNWydn8lj92Qz4De7aszs8dmUH6qntXtyOtHq2pZe6DCryC31mNzBeP10HVyrezmNFrdnbhk8d7bBrxVHZLiGj7cYN0Gu1RrZ9Uft+X1T/8pHo633FdfCVRWd6ffcLj2Malh/usL/S85ajVLeuLUhyHDQ6nf7k5UjESx/+FqmdB89u++r+/nXAsX/j/48FlZsvHGqePyebOf6T6Suo17Clq6m1p34CPx1DvrXolNcsnr86V8c9vKcsWrISJK1BOH9B8/k5rF8dzea53P6V02tsd8TV4aN08Y1OH+U0em0SsmkkWbSryWhH1jSwmWRUDr4wCZrsj1npqCpoY8oT9ccIvvZRwzJ8OW11uXcPTUUjBcyLkOLrkPxn2m8z4jv1CKiexaHpyHYFODVGcbOSt4TVkm3yUe/vGdgZ0nr0DWxsOVjHFw84vy9+1vytz0J2S549SxdnYyMOlLkBzYg7xT8FYUxmkFwUDIL4TtiyTFbeilrbcVr5VeED4sJZnoeCoHX8vUsvflvvFjknr5iP7ccVk2X5k6zNFkIC46kmtGpbOs6Cg/uimfqMhzBeDFm0rIG5DECD8C8NyJjDCMSEtkZw+V1tWQ5xf6Z1QyJ0vnp+M7WiYBxWtk3ay9ynDdldgkmO5npzOnDL1UIm2L5gXHkO9dKcV5gukZRUTAtHYis3sSF9wS2PFxyV3fljeYxKdAbO/WhtyuIJhzffBUJU/kzpDUtqK5rQ15UyMc2QATvuDzKTMv/w/410LYt9IvBSshNopHb/Ctq1nBmAwWfnqEj/dWcHlOa6/8aFUt6w5U8q3rcn0eiydy05P4eG95UM4Vbugaub/Y6+Du8vqhNeK9BLsL1PlCRKQE0u1c5jmAyFe2zoXYZN87dymKE4yBvlmtI9ftxjyd3QI3NlEmC1vny5q8TekWaDjtPNDNneFXy/3S2VUW3Zg6MpWEGM/FYZa4ItpnBbg+bpOTnkhJVS0ne2Dkuhpyf+k3XEqP2oUZGmolh9yfG6wnkV8IjWfEmAdCY70rWj2IsrqitCUlq7VHvnVe58vqNvmFkuJ2cFXLe34Eup0lOk7ul22LglMT3wFx0ZFcOzqdpVuO0tAmen3x5hJGZSQzPDU4AYO5rvrxu3pgwJsacn8xRox2sevGKvkEmhvCc328KxlyiaQmBeoV7LNl9SDn8SqKOylZcOKAxArYFQSdNuYJlNzpslTnfq8Ur4bEAVJC1h/yC+W+2bcyOGN0QMGYDCpPN7BqT4vsfeTEGdYfqGR2kLxxaIlc74mFYdSQB0LmZFkjP10RnoVgQkFEpJR/3bVc+lr7S5HK6koXkJItKYnVR1yNeY53vqxuE5MgKW5bF7TI68U+FILxxPBpsu7fhfL6lbmpJMZGtWptelZWDzBa3Z3BKfHERUf0yFKtasgDwfa+D6+XGywlGxJTQzumcGD0HCmJustPeb2xXiJ68wq6f9lTJbxxj1y3G/Pk+NCYJ1DyCyXV7cBHUnHvxIHAVL+oWGlws72L5fVRaSzb2iKvL9lcwuiMZLL7By9gMMIVub6rTD1yxRcGTpAOScVrZO1KZXVnDLlY5EF/vYK970p53K7yjJSei23Iy3e7Uh1nBFZB0Fdyrpc1+aK5LYG1gap++YVy/+x9N+DhOaVg7EBOnG7goz3lHD5xhg0HTwRcBMYTuWlJukau+EhsonRLK3pdglI00M0ZZ+X1N6HOj9lz0VyRB4dNC/7YFMWd3pnSUGjDX1397rt48hiTIGvl2xZI0FtEtOT4B8KwrpfXr8jpT1JsFIs3HeENl6weSG11b+SkJ3H0ZC1VZ3pW5Loa8kDJnCKzdVCP3BfyC0Ve9zV6vbFOKoWNmh1+pU+V8CMyCvpkwuF1wasg6Cv5hZLytuGv0nwn0CyNqBi5f7Yv9q/drB/ERUdy3eh0lhWVMu+Tw1wwKJmsIMrqNna71d09TF5XQx4otswVnQBpHfQyV1rInAJJGb57BXvegbqqwGurK4pTbHk9t4tldZsR14m8XlcVvGDa0XPkfHveCc75HDBrTAZVZxrYcvhkUIPc3LEj13tawJsa8kDJdMnpgybI7F1xRkREi7xee9L5cUVzIa43DJvaWSNTlNbYNddDFZMR00smEdDyvAmUYVPlPvJlIt1YB+v/4rcXf0WuyOsQgKze3AxbXpNMIQ8M6hNPfHRkQDXX68+cYvVrz1JbXen3OboaNeSBkpItQW+jOqH15/lOfiE01TmX149uhi2vQv7NKqsrXUfW5dLmeMQ1oRvDxP8SBWvo5cE5X1QM5N0AO5Y4N8xv/wgWfg0+/ZdfHxkbFclnLspk6shUhvbzU1Zf9yd49Q54/S6PHQojIgyjByazdr/3jmsdceCVB5iy+VHOPHcpHN7g93m6EjXkgWIM3PUOTLkr1CMJPwZPhqSBzryCxnqY91Wpf3319zt/bIpiM+ZWuG9taGR1m2FXwbe2Bze9Nb8Q6k7Cnrc73vfgavjoN/L71nl+f+Qjs0fz8pf8XB6o2Adv/kAKSu1+Ez75u8fdZuQPYMvhkxwo96MM9P4PyNn/D5Y1TeJ0XR3Nf7oePn4h8LbGnYwaciV0RERIZbbdDuT1938hHvnsX0FCv64Zn6Kczwy7CuL6dDyRrj8tk+jemdIdcO9KONXFzUmam2H+fZLxcudbMPQyWPqwtI5uw8wxAwApAesTdTU0zb2HA1Y6q8c/yV29fskqMx6Wfgde+Tyc6b5SuxpyJbSMniOVs3a84X2fkk/h/Z/DmNth1A1dNzZFOZ+JjHZFry+RXhHeePvHULEHbvqN9Ji3mmD7wq4bJ8DaP0hlvelPQJ8hcNPz0NwIC752jrc8OKUXFw7p06qSnCPeeoyIqoM8UH83c6bk8MitV/Cfp77BskH3w86l8Lsr4dC6IH6p4KGGXAktgy+C5EHevYLGeph3D/TqBzOf6tqxKcr5Tn4h1FfDnhWetx9YBR//Fi66Uzz4AWOg77AuzUGnfA+89ZhE71/4eXmvbzZc97iMe8P/nnNIwZgMio6cZN9xh/L63pWw9g8sTyrkaMqFjBnUm0uG9+O/Ls3m7j2XsGX6v2W/l6bDR891O6ldDbkSWiIixCvfs0KqTbXlvaeldePsZ6BX364fn6Kcz2RfJXEnRR7WvetPwfx7xAO+9ofynjFi/Pe9D6eOd/74mpth/r1SCOfGX7euMT/py5B1BSz7HpwobnXYTFdU/BIn8npdNcy/j6aUYXyr/EYKxgzEuD7nwRkjGdqvF/esjOD0He9I9sDyR+Cfn/UaOR8K1JAroSe/0LO8fmSjrI2P+5zUh1YUJbhERsty1Y4l0HCm9bYVj0PFXpGxY91ajeYXiry+rQvk9dW/k4p2M5+E5IGtt0VEiNxvNcOC+1p5yYP6xDPBqbz+5qNQVcy7eY9R0xzTqiNbr5gonr51HMWVp3nq3aPwmb/BzJ9JgODvrpAgwG6AGnIl9AyeJIE07nJdYx3M/SokpsGMn4ZubIpyvpNfCPU1sNtNXt//gRjRyXdB9hWt90+/APqN6Hx5/fhumUzkzpDJvCdSsuD6H0nd+PV/brWpYOxAtpacZO+xdorD7HkH1r0El9zLy4cGMLRfL/IHJrfaZXJ2X750aTZ/WXWAj/aWw5S74cvLpW7In2fCB78S5SCEqCFXQo8xUhxm9wo4c0LeW/kUHNsGNzwr0p+iKJ1D1pUQ37fFMNfViJydkgXXPnbu/ra8vv99qDnWOWNqbhJZPypGltXaa9s66Q5ZIlj+fag8cPbtWa7oda/yeu1JWHA/9MuhYsqDfLSnnFljMs7K6u58e/pIsvsn8OCrmzhV1wgDL4S73xM1463H4B+3d81SgxfUkCvdg/xCaG4Qie/wepnljv+8NIxQFKXziIxyyetviLz+1mPStvWm30rTFk/kF4qkvW1B54zp499C8WqY+TQkd1AFzhiR2EEkdpd3nNE7nklDU1jkTV5f/gicPAxzXmD5ziqami2vFefiYyJ5+taxHD5xhp++sU3ejOsNt70MBb+Afe/B7y6XdrMhQA250j0YNFHk9U2vSJR64gCY/pNQj0pRegb5hdBwSgqurP0DTPkqZF3mff+00dAvp3Pk9WM7JeVtZAGMvd3ZMX2GyPNi33uw/qWzb88ak8H2o9XsLmsjr+9eARv+ApfeD5kXsXhzCVkeZHV3JmX15c7Ls/nbxwf5YJfL+zZGIvrvfEvq4b9cIAG6XSy1qyFXugfGSHGYve/Cse1w43MQ3yfUo1KUnkHWFZLiueb3kl52zaPt72/L6wc+hJqy4I3DltSj46X4U3uSelsmfBGGXw3LH5V878r93DCknsGmjA/Wymsq98tEYcH90H8kTP0u5TV1fLSnnIKxnmV1d751/UiGpSbw0GubqK51a5WaMRbuXinlo9/+Mfzt5uBelw5QQ650H/Jvlp8TvgA5IWgXqSg9lcgoiVPBuCT1Xh0f0xny+qrfwKG1MOvnkJTu27HGiAMQEQl/vAaeHUfqS5P5IPYb/Nfam+DZcfLv+YugugTmvADRcSwrKnXJ6gM7/Ii46Eh+fts4SqrO8MSS7a03xibBLX+UuJ6Dq+DPs6Cp0bfv4CfarkvpPgyaAF9YIC1OFUXpWq7+vkSHZzqshZ42SrzaonkiLwdK2XZ4+yeQNxsuuMW/c/QeDHeukP7xLt7fdYx5G4/w4IyRpCe7ern3HwmDJwKwePMRsvsnMCojydFHTBiSwn9fOYzfr9zLzAsGcGWuW/17Y6TBzeCLpHxsF3XEVEOudC+GXRXqEShKz6RXX+jlQ0MTW15f+RRUl/ruQbvT1Cj13GMSfJfU25KaK/9c5A6r5fUNKxhSn8vXx+e02rW8po5Ve8q5Z+qIDmV1d755bS4rtpXx0GubWPbNK0mOi269Q3q+/OsiVFpXFEVR/CN/DmAFLq9/9Gs4skEiwBPTgjI0m/TkOC7K6svizUfO2ba06CjNFhSM9a0/ui2xl56s5SeLtgVrqH6jhlxRFEXxj7RRkJoXWPR66VZ496eyRn/BzcEbmxsFYzLYWVrDztLqVu8v3lTCsNQE8gY4k9XdGZ/Zh7uvGs4r64p5Z0fXBbZ5Qg25oiiK4j/5hZI/XX3U92ObGkRSj02Ggl8Gf2wuZl4wAGNoVbL1WHUdH+8tZ7aXIjBO+Ma1OeSmJ/Kd1zZRdaah4wM6CTXkiqIoiv+MdsnrW/2Q1z98Bko+gdm/hIT+QR+aTVpyHJOz+rbqUW7L6rN8lNXdiY0Sif14TT0/WrQ1GEP1CzXkiqIoiv+k5UmBGF/l9aNb4N2nJO109E2dMzY3Zo/NYHdZi7y+ZFMJw1MTGJnuu6zuztjBffjqVcN5df0h3t5eGoyh+owackVRFCUw8gsld/rkuQFlHrEl9fg+kjPeBUy/YAARBhZtKqGsupbV+8opGDvQb1ndnfuvGUHegCS+89pmqk53vcTeLQy5MWaGMWaHMWa3MeY7oR6PoiiK4gO+yuvv/xKObpJUs4R+nTo0m7SkOKZk92PxpiMs3SKy+uwAZHV3bIm9/FQ9P1xYFJRz+kLIDbkxJhJ4HpgJjAY+Z4wZHdpRKYqiKI5JzYW0fGfyeskmeO9nMOY2adbShcwam8GeY6f4/cq9jEhLJDdAWd2dCwb15t5pI3h942He3Nq1EnvIDTkwGdhtWdZey7LqgX8Bnb9goiiKogSP/EIo/hiqDnvfp7FeJPVe/WDmz7pubC5m5Iu8fvjEGa+dzgLhvmkjGJWRzHfnbqbyVH3Qz++N7lDZbRBQ7Pb6ENBujc7y8nJefvnlzhyToiiK4gPJDc3cDFQ9N5W6SM+12qOba0lpOMqKtDsp/ncntUDtgKHxvdl3OoaGvat5+eXgtx29MjaSP5Sk8IVnF3LLwOqODwgC3cGQe4o0sM7ZyZi7gLsABg0a1NljUhRFUXzgZHQ6W5KnkVLvpf830BAZx57EyRT3GtOFI2vNVf1PM7CmkbTYpk45f0ZcE1f1P01NYwTNFkQEHkvXIcayzrGZXYox5hLgMcuyprtePwxgWdZPvR0zadIka926dd42K4qiKErIsCwrKNHw7hhj1luWNcnTtu6wRr4WyDHGZBtjYoDPAqHRXBRFURQlQIJtxDsi5NK6ZVmNxpj7gGVAJPCSZVldH7+vKIqiKGFIyA05gGVZS4AloR6HoiiKooQb3UFaVxRFURTFT9SQK4qiKEoYo4ZcURRFUcIYNeSKoiiKEsaoIVcURVGUMEYNuaIoiqKEMWrIFUVRFCWMUUOuKIqiKGGMGnJFURRFCWPUkCuKoihKGBPy7mf+YIw5BhwI9ThCQH/geKgHcR6g1zF46LUMHnotg8f5eC2HWpaV6mlDWBrynooxZp23NnaKc/Q6Bg+9lsFDr2Xw6GnXUqV1RVEURQlj1JAriqIoShijhjy8eDHUAzhP0OsYPPRaBg+9lsGjR11LXSNXFEVRlDBGPXJFURRFCWPUkHdDjDFxxpg1xphPjTFFxpgfut7PNsasNsbsMsa8YoyJCfVYwwVjTKQxZqMxZpHrtV5LPzDG7DfGbDbGfGKMWed6r68x5k3XtXzTGJMS6nGGA8aYPsaYV40x240x24wxl+i19A1jzEjX36L976Qx5hs97TqqIe+e1AFXW5Y1DhgPzDDGXAw8BfzKsqwcoBL4cgjHGG58Hdjm9lqvpf9MsyxrvFt6z3eAFa5rucL1WumYZ4GllmXlAeOQv0+9lj5gWdYO19/ieGAicBqYSw+7jmrIuyGWUON6Ge36ZwFXA6+63v8LMCcEwws7jDGDgQLgj67XBr2WweQm5BqCXktHGGOSgSuBPwFYllVvWdYJ9FoGwjXAHsuyDtDDrqMa8m6KSwr+BCgD3gT2ACcsy2p07XIIGBSq8YUZzwAPAs2u1/3Qa+kvFrDcGLPeGHOX6710y7JKAFw/00I2uvBhGHAM+LNryeePxpgE9FoGwmeBf7p+71HXUQ15N8WyrCaXXDQYmAyM8rRb144q/DDGzAbKLMta7/62h131WjrjMsuyJgAzgXuNMVeGekBhShQwAXjBsqwLgVOc5/JvZ+KKcbkR+L9QjyUUqCHv5rjktneBi4E+xpgo16bBwJFQjSuMuAy40RizH/gXIqk/g15Lv7As64jrZxmyFjkZKDXGZAC4fpaFboRhwyHgkGVZq12vX0UMu15L/5gJbLAsq9T1ukddRzXk3RBjTKoxpo/r93jgWiQQ5h3gVtduXwTmh2aE4YNlWQ9bljXYsqwsRHp727Ks/0Svpc8YYxKMMUn278D1wBZgAXINQa+lIyzLOgoUG2NGut66BtiKXkt/+Rwtsjr0sOuoBWG6IcaYsUiARiQy2fq3ZVmPG2OGIV5lX2Aj8HnLsupCN9LwwhgzFXjAsqzZei19x3XN5rpeRgH/sCzrJ8aYfsC/gSHAQeA2y7IqQjTMsMEYMx4JwIwB9gJfwnW/o9fSMcaYXkAxMMyyrCrXez3qb1INuaIoiqKEMSqtK4qiKEoYo4ZcURRFUcIYNeSKoiiKEsaoIVcURVGUMEYNuaIoiqKEMWrIFUVRFCWMUUOuKIqiKGGMGnJFURRFCWPUkCtKD8cYM8MY874xptIYU2GMWWaMGeW2fYoxZoMxptbVqWuWMcZyVcqz9xltjFlsjKk2xpQZY/5pjBkQki+kKD0MNeSKoiQgjWQmA1OBKmChMSbGGJMILAK2AxORdrBPux/sakrxHlJ3fTLSGyARWGCM0WeMonQyWqJVUZRWuBqinASuAvKBnwKDLMs649r+H8DfgWmWZb1rjHkcaW96jds5UoAKYIplWWu6+jsoSk9CZ8uK0sMxxgw3xvzDGLPHGHMSKEWeDUOAPGCLbcRdrG5zionAlcaYGvsf0sQCYHhnj19RejpRHe+iKMp5zkLgMHC362cj0lIzBjBAR7JdBLAYeMDDtlIP7ymKEkTUkCtKD8bV7nEUcK9lWe+43ptAy7NhG/AFY0y8m1c+uc1pNgC3Awcsy2rogmEriuKGSuuK0rOpBI4D/22MGWGMuQr4HeKVg6yFNwF/cEWmXwt817XN9tSfB3oDr7gi3IcZY641xrxojEnquq+iKD0TNeSK0oOxLKsZ+AwwFok6fx74PlDn2l4D3IAEvW1EItYfcx1e69rnCHAZ0AwsBYpc56mzz6MoSuehUeuKoviEMeYmYC6QZlnW8VCPR1F6OrpGrihKuxhjvgjsRSLRL0ByzheqEVeU7oEackVROiId+CGQARxFItQfCumIFEU5i0rriqIoihLGaLCboiiKooQxasgVRVEUJYxRQ64oiqIoYYwackVRFEUJY9SQK4qiKEoYo4ZcURRFUcKY/w+unIa56gyN9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.axhline(0, c='gray')\n",
    "ax.plot(crosstab.index, crosstab[0], label='non-quitters')\n",
    "ax.plot(crosstab.index, crosstab[1], label='quitters')\n",
    "ax.set_xlabel('age', fontsize=14)\n",
    "ax.set_ylabel('count', fontsize=14)\n",
    "ax.legend(fontsize=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that there are actually a few ages with zero counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>63</th>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64</th>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65</th>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>6</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>71</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>74</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk  0  1\n",
       "age       \n",
       "63    3  3\n",
       "64    7  1\n",
       "65    3  2\n",
       "66    4  0\n",
       "67    2  0\n",
       "69    6  2\n",
       "70    2  1\n",
       "71    0  1\n",
       "72    2  2\n",
       "74    0  1"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "crosstab.iloc[-10:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a discussion on ages with zero counts, see Fine Point 12.2, pg 155."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"The effect estimate obtained in the pseudo-population created by weights $0.5 \\, / \\, f(A|L)$\n",
    "is equal to that obtained in the pseudo-population created by weights $1 \\, / \\, f(A|L)$.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>    1.7800</td> <td>    0.225</td> <td>    7.920</td> <td> 0.000</td> <td>    1.340</td> <td>    2.220</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    3.4405</td> <td>    0.525</td> <td>    6.547</td> <td> 0.000</td> <td>    2.411</td> <td>    4.470</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gee = sm.GEE(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    groups=nhefs.seqn,\n",
    "    weights=(0.5 * weights)\n",
    ")\n",
    "res = gee.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Second, we need to estimate Pr[A=1] for the numerator of the weights. We can obtain a nonparametric estimate by the ratio 403/1566 or, equivalently, by fitting a saturated logistic model for Pr[A=1] with an intercept and no covariates.\" pg 154"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "qsmk = (nhefs.qsmk == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.25734355044699875"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# option 1\n",
    "qsmk_mean = qsmk.mean()\n",
    "qsmk_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.570260\n",
      "         Iterations 5\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>   -1.0598</td> <td>    0.058</td> <td>  -18.335</td> <td> 0.000</td> <td>   -1.173</td> <td>   -0.947</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# option 2\n",
    "lgt = sm.Logit(qsmk, nhefs.constant)\n",
    "res = lgt.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "lgt_pred = res.predict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check for equivalence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "equivalent: True\n"
     ]
    }
   ],
   "source": [
    "equivalent = np.all(np.isclose(lgt_pred, qsmk_mean))\n",
    "print('equivalent: {}'.format(equivalent))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create stabilized IP weights. Shortcut: modify the IP weights already calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_weights = np.zeros(nhefs.shape[0])\n",
    "s_weights[qsmk] = qsmk.mean() * weights[qsmk]    # `qsmk` was defined a few cells ago\n",
    "s_weights[~qsmk] = (1 - qsmk).mean() * weights[~qsmk]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stabilized weights\n",
      " min   mean    max\n",
      "------------------\n",
      "0.33   1.00   4.30\n"
     ]
    }
   ],
   "source": [
    "print('Stabilized weights')\n",
    "print(' min   mean    max')\n",
    "print('------------------')\n",
    "print('{:>04.2f}   {:>04.2f}   {:>04.2f}'.format(\n",
    "    s_weights.min(),\n",
    "    s_weights.mean(),\n",
    "    s_weights.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refit the model from the last section, using the new weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>    1.7800</td> <td>    0.225</td> <td>    7.920</td> <td> 0.000</td> <td>    1.340</td> <td>    2.220</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    3.4405</td> <td>    0.525</td> <td>    6.547</td> <td> 0.000</td> <td>    2.411</td> <td>    4.470</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gee = sm.GEE(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    groups=nhefs.seqn,\n",
    "    weights=s_weights\n",
    ")\n",
    "res = gee.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate   95% C.I.\n",
      "theta_1       3.44   (2.4, 4.5)\n"
     ]
    }
   ],
   "source": [
    "est = res.params.qsmk\n",
    "conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
    "lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
    "\n",
    "print('           estimate   95% C.I.')\n",
    "print('theta_1     {:>6.2f}   ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate is the same as in the previous section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check again for no association between sex and qsmk in the the pseudo-population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "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>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>567.098228</td>\n",
       "      <td>196.513582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>595.423986</td>\n",
       "      <td>205.154456</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk           0           1\n",
       "sex                         \n",
       "0     567.098228  196.513582\n",
       "1     595.423986  205.154456"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(nhefs.sex, nhefs.qsmk, s_weights, aggfunc='sum')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subset the data to subjects that smoked 25 or fewer cigarettes per day at baseline. In this case, we can either obtain the subset from the original dataset, or we can obtain it from the reduced dataset that we've been using. I'll get it from the reduced subset, since it already contains dummy features we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1162, 64)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from original dataset\n",
    "\n",
    "intensity25 = nhefs_all.loc[\n",
    "    (nhefs_all.smokeintensity <= 25) & ~nhefs_all.wt82.isnull()\n",
    "]\n",
    "intensity25.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1162, 83)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from reduced dataset\n",
    "\n",
    "intensity25 = nhefs.loc[nhefs.smokeintensity <= 25]\n",
    "intensity25.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the stabilized IP weights $SW^A = f(A) \\, / \\, f(A|L)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"we assumed that the density f(A|L) was normal (Gaussian) with mean $\\mu = E[A|L]$ and variance $\\sigma^2$.  We then used a linear regression model to estimate the mean $E[A|L]$ and variance of residuals $\\sigma^2$ for all combinations of values of L.\" pg 156\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = intensity25.smkintensity82_71\n",
    "X = intensity25[[\n",
    "    'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = sm.OLS(A, X)\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_pred = res.predict(X)   # i.e., E[A|L]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The denominator is the distribution, $N(\\mu, \\sigma)$, evaluated at each point of $y = A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "fAL = scipy.stats.norm.pdf(\n",
    "    A,                        # A\n",
    "    A_pred,                   # mu = E[A|L]\n",
    "    np.sqrt(res.mse_resid)    # sigma\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"We also assumed that the density f(A) in the numerator was normal.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAFlCAYAAAA+t0u5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXLklEQVR4nO3df5Bd9Xnf8fcTIFjDpggKbBWhVnSiZIzZmIQdwoz7x67xFBk8Ec6EDB7qSDGt0hnisafqJML+w05TpsqkmDat4xkleKw0jteMfwQN4DZEYYd4JpggmyCw6lqxVSzQiHEM2GtTOouf/nGPyvVqtXu1e+/eR+e+XzM7e8/3nHvv8+zZ1Uffc889NzITSZJUw48NuwBJkvQ6g1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpELOHXYBAJdccklu3rx52GX0zfe//30uuOCCYZex5ka1bxjd3u179Ixq74Po++DBg9/OzEsXjpcI5s2bN/PEE08Mu4y+mZ2dZWpqathlrLlR7RtGt3f7Hj2j2vsg+o6I/73YuIeyJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqpMSnS0nVbd794JLrd03Ms2ORbY7uuWlQJUlqqWVnzBHxhoh4PCL+NiKeiYjfbsaviIgvRcTXI+LTEfHjzfj5zfKRZv3mwbYgSVJ79HIo+1XgrZn5ZuBqYGtEXAf8LnBPZm4BXgRub7a/HXgxM38KuKfZTpIk9WDZYM6OuWbxvOYrgbcCn2nG9wE3N7e3Ncs066+PiOhbxZIktVhk5vIbRZwDHAR+Cvgo8HvAY82smIjYBHwhM6+KiKeBrZl5rFn3d8AvZOa3FzzmTmAnwPj4+DUzMzP962rI5ubmGBsbG3YZa67NfR967uUl14+vgxOvnDo+sfHCAVVUQ5v3+VJGtW8Y3d4H0ff09PTBzJxcON7TyV+Z+RpwdUSsBz4PvHGxzZrvi82OT0n/zNwL7AWYnJzMqampXko5K8zOztKmfnrV5r4XO7Gr266Jee4+dOqf09HbpgZUUQ1t3udLGdW+YXR7X8u+z+jtUpn5EjALXAesj4iT/xJdDjzf3D4GbAJo1l8IfKcfxUqS1Ha9nJV9aTNTJiLWAW8DDgOPAL/cbLYduL+5vb9Zpln/l9nL8XJJktTToewNwL7mdeYfA+7LzAci4qvATET8e+ArwL3N9vcC/y0ijtCZKd86gLolSWqlZYM5M58Cfm6R8W8A1y4y/n+AW/pSnSRJI8ZLckqSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIcsGc0RsiohHIuJwRDwTEe9rxj8cEc9FxJPN141d97kzIo5ExNci4oZBNiBJUpuc28M288CuzPxyRPwEcDAiHm7W3ZOZ/7F744i4ErgVeBPwk8BfRMRPZ+Zr/SxckqQ2WnbGnJnHM/PLze3vAYeBjUvcZRswk5mvZuY3gSPAtf0oVpKktovM7H3jiM3Ao8BVwL8BdgDfBZ6gM6t+MSL+K/BYZv5Jc597gS9k5mcWPNZOYCfA+Pj4NTMzM6vtpYy5uTnGxsaGXcaaa3Pfh557ecn14+vgxCunjk9svHBAFdXQ5n2+lFHtG0a390H0PT09fTAzJxeO93IoG4CIGAM+C7w/M78bER8DfgfI5vvdwHuAWOTup6R/Zu4F9gJMTk7m1NRUr6WUNzs7S5v66VWb+96x+8El1++amOfuQ6f+OR29bWpAFdXQ5n2+lFHtG0a397Xsu6ezsiPiPDqh/MnM/BxAZp7IzNcy84fAH/L64epjwKauu18OPN+/kiVJaq9ezsoO4F7gcGZ+pGt8Q9dm7wSebm7vB26NiPMj4gpgC/B4/0qWJKm9ejmU/Rbg3cChiHiyGfsA8K6IuJrOYeqjwK8DZOYzEXEf8FU6Z3Tf4RnZkiT1Ztlgzswvsvjrxg8tcZ+7gLtWUZckSSPJK39JklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiHLBnNEbIqIRyLicEQ8ExHva8YvjoiHI+LrzfeLmvGIiN+PiCMR8VRE/Pygm5AkqS16mTHPA7sy843AdcAdEXElsBs4kJlbgAPNMsDbgS3N107gY32vWpKkllo2mDPzeGZ+ubn9PeAwsBHYBuxrNtsH3Nzc3gb8cXY8BqyPiA19r1ySpBaKzOx944jNwKPAVcCzmbm+a92LmXlRRDwA7MnMLzbjB4DfyswnFjzWTjozasbHx6+ZmZlZZSt1zM3NMTY2Nuwy1lyb+z703MtLrh9fBydeOXV8YuOFA6qohjbv86WMat8wur0Pou/p6emDmTm5cPzcXh8gIsaAzwLvz8zvRsRpN11k7JT0z8y9wF6AycnJnJqa6rWU8mZnZ2lTP71qc987dj+45PpdE/PcfejUP6ejt00NqKIa2rzPlzKqfcPo9r6Wffd0VnZEnEcnlD+ZmZ9rhk+cPETdfH+hGT8GbOq6++XA8/0pV5KkduvlrOwA7gUOZ+ZHulbtB7Y3t7cD93eN/2pzdvZ1wMuZebyPNUuS1Fq9HMp+C/Bu4FBEPNmMfQDYA9wXEbcDzwK3NOseAm4EjgA/AH6trxVLktRiywZzcxLX6V5Qvn6R7RO4Y5V1SZI0krzylyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVsmwwR8THI+KFiHi6a+zDEfFcRDzZfN3Yte7OiDgSEV+LiBsGVbgkSW3Uy4z5E8DWRcbvycyrm6+HACLiSuBW4E3Nff4gIs7pV7GSJLXdssGcmY8C3+nx8bYBM5n5amZ+EzgCXLuK+iRJGimreY35NyLiqeZQ90XN2EbgW13bHGvGJElSDyIzl98oYjPwQGZe1SyPA98GEvgdYENmviciPgr8dWb+SbPdvcBDmfnZRR5zJ7ATYHx8/JqZmZm+NFTB3NwcY2Njwy5jzbW570PPvbzk+vF1cOKVU8cnNl44oIpqaPM+X8qo9g2j2/sg+p6enj6YmZMLx89dyYNl5omTtyPiD4EHmsVjwKauTS8Hnj/NY+wF9gJMTk7m1NTUSkopaXZ2ljb106s2971j94NLrt81Mc/dh079czp629SAKqqhzft8KaPaN4xu72vZ94oOZUfEhq7FdwInz9jeD9waEedHxBXAFuDx1ZUoSdLoWHbGHBGfAqaASyLiGPAhYCoirqZzKPso8OsAmflMRNwHfBWYB+7IzNcGU7okSe2zbDBn5rsWGb53ie3vAu5aTVGSJI0qr/wlSVIhBrMkSYUYzJIkFWIwS5JUiMEsSVIhBrMkSYUYzJIkFWIwS5JUiMEsSVIhBrMkSYUYzJIkFWIwS5JUiMEsSVIhBrMkSYUYzJIkFWIwS5JUiMEsSVIhBrMkSYUYzJIkFXLusAtQXZt3P3hG2++amGdqMKVI0shwxixJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhfh2qT5Y+LaiXRPz7FjmrUZH99w0yJIkSWcpZ8ySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIcsGc0R8PCJeiIinu8YujoiHI+LrzfeLmvGIiN+PiCMR8VRE/Pwgi5ckqW16mTF/Ati6YGw3cCAztwAHmmWAtwNbmq+dwMf6U6YkSaNh2WDOzEeB7ywY3gbsa27vA27uGv/j7HgMWB8RG/pVrCRJbXfuCu83npnHATLzeERc1oxvBL7Vtd2xZuz4ykuUFrd594Mrut/RPTf1uRJJ6p/IzOU3itgMPJCZVzXLL2Xm+q71L2bmRRHxIPAfMvOLzfgB4Dcz8+Aij7mTzuFuxsfHr5mZmelDO8Nx6LmXf2R5fB2ceGXp+0xsvHCAFfXHwr6WM74OLrt47fo60/pOWsnPfrnnOt0+Pxv282rMzc0xNjY27DLW3Kj2DaPb+yD6np6ePpiZkwvHVzpjPhERG5rZ8gbghWb8GLCpa7vLgecXe4DM3AvsBZicnMypqakVljJ8OxbM3HZNzHP3oaV/tEdvmxpgRf2xsK/l7JqY51fWcD+eaX0nreRnv9xznW6fnw37eTVmZ2c5m/92V2pU+4bR7X0t+17p26X2A9ub29uB+7vGf7U5O/s64OWTh7wlSdLylp0xR8SngCngkog4BnwI2APcFxG3A88CtzSbPwTcCBwBfgD82gBqliSptZYN5sx812lWXb/ItgncsdqiJNV28sS7XRPzPb+k4El3Um+88pckSYUYzJIkFWIwS5JUiMEsSVIhBrMkSYUYzJIkFWIwS5JUyEovySlpQPxwDmm0OWOWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpELOHXYB0lrbvPvBYZcgSafljFmSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEIMZkmSCvGSnNIAeflPSWfKGbMkSYUYzJIkFWIwS5JUiMEsSVIhqzr5KyKOAt8DXgPmM3MyIi4GPg1sBo4Cv5KZL66uTEmSRkM/ZszTmXl1Zk42y7uBA5m5BTjQLEuSpB4M4lD2NmBfc3sfcPMAnkOSpFZabTAn8OcRcTAidjZj45l5HKD5ftkqn0OSpJERmbnyO0f8ZGY+HxGXAQ8D7wX2Z+b6rm1ezMyLFrnvTmAnwPj4+DUzMzMrrmPYDj338o8sj6+DE68sfZ+JjRcOsKL+WNjXcsbXwWUXr11fZ1rfIPWyzwdtLX+nTv7sz6Tvs+F3vldzc3OMjY0Nu4yhGNXeB9H39PT0wa6Xgf+/VQXzjzxQxIeBOeBfAVOZeTwiNgCzmfkzS913cnIyn3jiib7UMQwLr+60a2Keuw8tfV7d0T03DbKkvjjTq1btmpjnvbdtG1A1p6p0Va1e9vmgreXv1Mmf/Zn0fTb8zvdqdnaWqampYZcxFKPa+yD6johFg3nFh7Ij4oKI+ImTt4F/DjwN7Ae2N5ttB+5f6XNIkjRqVvNf/HHg8xFx8nH+NDP/e0T8DXBfRNwOPAvcsvoyJUkaDSsO5sz8BvDmRcb/Hrh+NUVJkjSqvPKXJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBViMEuSVIjBLElSIQazJEmFGMySJBUy3I/D0RlbyScqtelTfaRB6vXva9fEPDu6tvVvTP3kjFmSpEIMZkmSCjGYJUkqxGCWJKkQg1mSpEI8K1slrORsc0lqI2fMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIH/soqbSVfiTo0T039bkSaW04Y5YkqRBnzOqrlc5uJEkdzpglSSrEGbM0wjzCIdVjMI8A//EdDe5nqR1aGcyexSnV438cpN74GrMkSYUYzJIkFdLKQ9lnAw/rSZIWYzBLkha12ARi18Q8O5aYWJwN5+qsZGL0ia0XDKCSxQ0smCNiK/CfgXOAP8rMPYN6LklaqK1HpVbS19kQlnrdQF5jjohzgI8CbweuBN4VEVcO4rkkSWqTQZ38dS1wJDO/kZn/F5gBtg3ouSRJao1BHcreCHyra/kY8AsDei5J0ghq68sVkZn9f9CIW4AbMvNfNsvvBq7NzPd2bbMT2Nks/gzwtb4XMjyXAN8edhFDMKp9w+j2bt+jZ1R7H0Tf/yQzL104OKgZ8zFgU9fy5cDz3Rtk5l5g74Cef6gi4onMnBx2HWttVPuG0e3dvkfPqPa+ln0P6jXmvwG2RMQVEfHjwK3A/gE9lyRJrTGQGXNmzkfEbwD/g87bpT6emc8M4rkkSWqTgb2POTMfAh4a1OMX18pD9D0Y1b5hdHu379Ezqr2vWd8DOflLkiStjB9iIUlSIQZzn0TE70XE/4yIpyLi8xGxvmvdnRFxJCK+FhE3DLPOQYiIWyLimYj4YURMLljX9t63Nr0diYjdw65nkCLi4xHxQkQ83TV2cUQ8HBFfb75fNMwaByEiNkXEIxFxuPk9f18z3ureI+INEfF4RPxt0/dvN+NXRMSXmr4/3Zzg2zoRcU5EfCUiHmiW16xvg7l/HgauysyfBf4XcCdAcynSW4E3AVuBP2guWdomTwO/BDzaPdj23kfw0rOfoLMfu+0GDmTmFuBAs9w288CuzHwjcB1wR7Of2977q8BbM/PNwNXA1oi4Dvhd4J6m7xeB24dY4yC9DzjctbxmfRvMfZKZf56Z883iY3Teuw2dS5HOZOarmflN4AidS5a2RmYezszFLhDT9t5H6tKzmfko8J0Fw9uAfc3tfcDNa1rUGsjM45n55eb29+j8Y72RlveeHXPN4nnNVwJvBT7TjLeub4CIuBy4CfijZjlYw74N5sF4D/CF5vZilyfduOYVDUfbe297f70Yz8zj0Akw4LIh1zNQEbEZ+DngS4xA783h3CeBF+gcFfw74KWuSUhbf+f/E/CbwA+b5X/IGvbt5zGfgYj4C+AfLbLqg5l5f7PNB+kc+vrkybstsv1Zdyp8L70vdrdFxs663pfQ9v7UJSLGgM8C78/M73YmUe2Wma8BVzfnzHweeONim61tVYMVEe8AXsjMgxExdXJ4kU0H1rfBfAYy821LrY+I7cA7gOvz9fehLXt50rPBcr2fRit6X0Lb++vFiYjYkJnHI2IDnZlV60TEeXRC+ZOZ+blmeCR6B8jMlyJils5r7Osj4txm9tjG3/m3AL8YETcCbwD+AZ0Z9Jr17aHsPomIrcBvAb+YmT/oWrUfuDUizo+IK4AtwOPDqHEI2t67l57t9Lu9ub0dON3Rk7NW8/rivcDhzPxI16pW9x4Rl558d0lErAPeRuf19UeAX242a13fmXlnZl6emZvp/E3/ZWbexhr27QVG+iQijgDnA3/fDD2Wmf+6WfdBOq87z9M5DPaFxR/l7BQR7wT+C3Ap8BLwZGbe0Kxre+830vnf9MlLz9415JIGJiI+BUzR+ZSdE8CHgD8D7gP+MfAscEtmLjxB7KwWEf8M+CvgEK+/5vgBOq8zt7b3iPhZOic5nUNnEndfZv67iPindE50vBj4CvAvMvPV4VU6OM2h7H+bme9Yy74NZkmSCvFQtiRJhRjMkiQVYjBLklSIwSxJUiEGsyRJhRjMkiQVYjBLklSIwSxJUiH/D44FeXkmHv3+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "A.hist(bins=30, ax=ax);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.057659208261618, 10.467830908151853)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.mean(), A.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "fA = scipy.stats.norm.pdf(A, A.mean(), A.std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the stabilized IP weights are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw = fA / fAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stabilized weights\n",
      " min   mean    max\n",
      "------------------\n",
      "0.19   1.00   5.10\n"
     ]
    }
   ],
   "source": [
    "print('Stabilized weights')\n",
    "print(' min   mean    max')\n",
    "print('------------------')\n",
    "print('{:>04.2f}   {:>04.2f}   {:>04.2f}'.format(\n",
    "    sw.min(),\n",
    "    sw.mean(),\n",
    "    sw.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now fit the marginal structural model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = intensity25.wt82_71\n",
    "X = pd.DataFrame(OrderedDict((\n",
    "    ('constant', np.ones(y.shape[0])),\n",
    "    ('A', A),\n",
    "    ('A^2', A**2)\n",
    ")))\n",
    "\n",
    "model = sm.GEE(\n",
    "    y,\n",
    "    X,\n",
    "    groups=intensity25.seqn,\n",
    "    weights=sw\n",
    ")\n",
    "res = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>    2.0045</td> <td>    0.295</td> <td>    6.792</td> <td> 0.000</td> <td>    1.426</td> <td>    2.583</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A</th>        <td>   -0.1090</td> <td>    0.032</td> <td>   -3.456</td> <td> 0.001</td> <td>   -0.171</td> <td>   -0.047</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>A^2</th>      <td>    0.0027</td> <td>    0.002</td> <td>    1.115</td> <td> 0.265</td> <td>   -0.002</td> <td>    0.007</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the estimate and confidence interval for \"no change\", you can read off the values in the `constant` row above (because `A` and `A^2` will be zero).\n",
    "\n",
    "Getting Statmodels to calculate the estimate and confidence interval for when smoking increases by 20 cigarettes / day will take a couple extra steps. In Chapter 11, the regression result had a `get_prediction` method. The GEE result doesn't (yet?) have that _method_, so we'll use the hidden `get_prediction` _function_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.regression._prediction import get_prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "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>mean</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>2.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.9</td>\n",
       "      <td>-1.7</td>\n",
       "      <td>3.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mean  mean_ci_lower  mean_ci_upper\n",
       "0   2.0            1.4            2.6\n",
       "1   0.9           -1.7            3.5"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred_inputs = [\n",
    "    [1, 0, 0],       # no change in smoking intensity\n",
    "    [1, 20, 20**2],  # plus 20 cigarettes / day\n",
    "]\n",
    "pred = get_prediction(res, exog=pred_inputs)\n",
    "summary = pred.summary_frame().round(1)\n",
    "summary[[\"mean\", \"mean_ci_lower\", \"mean_ci_upper\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can relabel the rows and columns to make this table a little nicer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "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>estimate</th>\n",
       "      <th>CI lower</th>\n",
       "      <th>CI upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>no change</th>\n",
       "      <td>2.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>2.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+20 per day</th>\n",
       "      <td>0.9</td>\n",
       "      <td>-1.7</td>\n",
       "      <td>3.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             estimate  CI lower  CI upper\n",
       "no change         2.0       1.4       2.6\n",
       "+20 per day       0.9      -1.7       3.5"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = summary[[\"mean\", \"mean_ci_lower\", \"mean_ci_upper\"]]\n",
    "summary.index = [\"no change\", \"+20 per day\"]\n",
    "summary.columns = [\"estimate\", \"CI lower\", \"CI upper\"]\n",
    "summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: since the `get_predictions` function wasn't attached to the GEE regression result, it might not work correctly with other versions of the GEE model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"if interested in the causal effect of quitting smoking A (1: yes, 0: no) on the risk of death D (1: yes, 0: no) by 1982, one could consider a _marginal structural logistic model_\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>   -1.4905</td> <td>    0.079</td> <td>  -18.881</td> <td> 0.000</td> <td>   -1.645</td> <td>   -1.336</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    0.0301</td> <td>    0.157</td> <td>    0.191</td> <td> 0.848</td> <td>   -0.278</td> <td>    0.338</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sm.GEE(\n",
    "    nhefs.death,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    groups=nhefs.seqn,\n",
    "    weights=s_weights,\n",
    "    family=sm.families.Binomial()\n",
    ")\n",
    "res = model.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Odd ratio is $\\exp(\\hat{\\theta}_1)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate   95% C.I.\n",
      "odds ratio    1.03   (0.8, 1.4)\n"
     ]
    }
   ],
   "source": [
    "est = np.exp(res.params.qsmk)\n",
    "conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
    "lo = np.exp(conf_ints[0]['qsmk'])\n",
    "hi = np.exp(conf_ints[1]['qsmk'])\n",
    "\n",
    "print('           estimate   95% C.I.')\n",
    "print('odds ratio  {:>6.2f}   ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the numerator of the IP weights. Reuse the basic `weights` for the denominator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.567819\n",
      "         Iterations 5\n"
     ]
    }
   ],
   "source": [
    "numer = logit_ip_f(nhefs.qsmk, nhefs[['constant', 'sex']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_AV = numer * weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stabilized weights\n",
      " min   mean    max\n",
      "------------------\n",
      "0.29   1.00   3.80\n"
     ]
    }
   ],
   "source": [
    "print('Stabilized weights')\n",
    "print(' min   mean    max')\n",
    "print('------------------')\n",
    "print('{:>04.2f}   {:>04.2f}   {:>04.2f}'.format(\n",
    "    sw_AV.min(),\n",
    "    sw_AV.mean(),\n",
    "    sw_AV.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 83)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th>        <td>    1.7844</td> <td>    0.310</td> <td>    5.752</td> <td> 0.000</td> <td>    1.176</td> <td>    2.393</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>            <td>    3.5220</td> <td>    0.658</td> <td>    5.353</td> <td> 0.000</td> <td>    2.232</td> <td>    4.811</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sex</th>             <td>   -0.0087</td> <td>    0.449</td> <td>   -0.019</td> <td> 0.985</td> <td>   -0.890</td> <td>    0.872</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk_and_female</th> <td>   -0.1595</td> <td>    1.047</td> <td>   -0.152</td> <td> 0.879</td> <td>   -2.212</td> <td>    1.893</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs['qsmk_and_female'] = nhefs.qsmk * nhefs.sex\n",
    "\n",
    "model = sm.WLS(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk', 'sex', 'qsmk_and_female']],\n",
    "    weights=sw_AV\n",
    ")\n",
    "res = model.fit(cov_type='cluster', cov_kwds={'groups': nhefs.seqn})\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going back to the original dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1629, 64)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll add features that were added to the reduced dataset that we've been using"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add constant feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['constant'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add dummy features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "edu_dummies = pd.get_dummies(nhefs_all.education, prefix='edu')\n",
    "exercise_dummies = pd.get_dummies(nhefs_all.exercise, prefix='exercise')\n",
    "active_dummies = pd.get_dummies(nhefs_all.active, prefix='active')\n",
    "\n",
    "nhefs_all = pd.concat(\n",
    "    [nhefs_all, edu_dummies, exercise_dummies, active_dummies],\n",
    "    axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add squared features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
    "    nhefs_all['{}^2'.format(col)] = nhefs_all[col] * nhefs_all[col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also add a feature to track censored individuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['censored'] = nhefs_all.wt82.isnull().astype('int')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the IP weights for treatment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ip = nhefs_all[[\n",
    "    'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.542264\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "ip_denom = logit_ip_f(nhefs_all.qsmk, X_ip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.575901\n",
      "         Iterations 5\n"
     ]
    }
   ],
   "source": [
    "ip_numer = logit_ip_f(nhefs_all.qsmk, nhefs_all.constant)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_A = ip_numer / ip_denom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stabilized weights\n",
      " min   mean    max\n",
      "------------------\n",
      "0.33   1.00   4.21\n"
     ]
    }
   ],
   "source": [
    "print('Stabilized weights')\n",
    "print(' min   mean    max')\n",
    "print('------------------')\n",
    "print('{:>04.2f}   {:>04.2f}   {:>04.2f}'.format(\n",
    "    sw_A.min(),\n",
    "    sw_A.mean(),\n",
    "    sw_A.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the IP weights for censoring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# same as previous, but with 'qsmk' added\n",
    "\n",
    "X_ip = nhefs_all[[\n",
    "    'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
    "    'qsmk'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.142836\n",
      "         Iterations 8\n"
     ]
    }
   ],
   "source": [
    "ip_denom = logit_ip_f(nhefs_all.censored, X_ip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.161989\n",
      "         Iterations 7\n"
     ]
    }
   ],
   "source": [
    "ip_numer = logit_ip_f(\n",
    "    nhefs_all.censored,\n",
    "    nhefs_all[['constant', 'qsmk']]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_C = ip_numer / ip_denom\n",
    "sw_C[nhefs_all.censored == 1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stabilized weights\n",
      " min   mean    max\n",
      "------------------\n",
      "0.94   1.00   1.72\n"
     ]
    }
   ],
   "source": [
    "print('Stabilized weights')\n",
    "print(' min   mean    max')\n",
    "print('------------------')\n",
    "print('{:>04.2f}   {:>04.2f}   {:>04.2f}'.format(\n",
    "    sw_C.min(),\n",
    "    sw_C.mean(),\n",
    "    sw_C.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create the combined IP weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_AC = sw_A * sw_C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stabilized weights\n",
      " min   mean    max\n",
      "------------------\n",
      "0.35   1.00   4.09\n"
     ]
    }
   ],
   "source": [
    "print('Stabilized weights')\n",
    "print(' min   mean    max')\n",
    "print('------------------')\n",
    "print('{:>04.2f}   {:>04.2f}   {:>04.2f}'.format(\n",
    "    sw_AC.min(),\n",
    "    sw_AC.mean(),\n",
    "    sw_AC.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now model weight gain using the combined IP weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "wls = sm.WLS(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    weights=sw_AC[nhefs_all.censored == 0]\n",
    ") \n",
    "res = wls.fit(cov_type='cluster', cov_kwds={'groups': nhefs.seqn})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<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>constant</th> <td>    1.6620</td> <td>    0.233</td> <td>    7.136</td> <td> 0.000</td> <td>    1.206</td> <td>    2.118</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    3.4965</td> <td>    0.526</td> <td>    6.648</td> <td> 0.000</td> <td>    2.466</td> <td>    4.527</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}