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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "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, ...\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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 11."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style  type=\"text/css\" >\n",
       "</style>  \n",
       "<table id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774\" > \n",
       "<thead>    <tr> \n",
       "        <th class=\"index_name level0\" >qsmk</th> \n",
       "        <th class=\"col_heading level0 col0\" >1</th> \n",
       "        <th class=\"col_heading level0 col1\" >0</th> \n",
       "    </tr></thead> \n",
       "<tbody>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row0\" class=\"row_heading level0 row0\" >Age, years</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row0_col0\" class=\"data row0 col0\" >46.2</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row0_col1\" class=\"data row0 col1\" >42.8</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row1\" class=\"row_heading level0 row1\" >Men, %</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row1_col0\" class=\"data row1 col0\" >54.6</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row1_col1\" class=\"data row1 col1\" >46.6</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row2\" class=\"row_heading level0 row2\" >White, %</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row2_col0\" class=\"data row2 col0\" >91.1</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row2_col1\" class=\"data row2 col1\" >85.4</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row3\" class=\"row_heading level0 row3\" >University education, %</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row3_col0\" class=\"data row3 col0\" >15.4</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row3_col1\" class=\"data row3 col1\" >9.9</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row4\" class=\"row_heading level0 row4\" >Weight, kg</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row4_col0\" class=\"data row4 col0\" >72.4</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row4_col1\" class=\"data row4 col1\" >70.3</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row5\" class=\"row_heading level0 row5\" >Cigarettes/day</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row5_col0\" class=\"data row5 col0\" >18.6</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row5_col1\" class=\"data row5 col1\" >21.2</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row6\" class=\"row_heading level0 row6\" >Years smoking</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row6_col0\" class=\"data row6 col0\" >26.0</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row6_col1\" class=\"data row6 col1\" >24.1</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row7\" class=\"row_heading level0 row7\" >Little or no exercise, %</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row7_col0\" class=\"data row7 col0\" >40.7</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row7_col1\" class=\"data row7 col1\" >37.9</td> \n",
       "    </tr>    <tr> \n",
       "        <th id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774level0_row8\" class=\"row_heading level0 row8\" >Inactive daily life, %</th> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row8_col0\" class=\"data row8 col0\" >11.2</td> \n",
       "        <td id=\"T_c02a4804_09c1_11e8_855c_74c63b4fe774row8_col1\" class=\"data row8 col1\" >8.9</td> \n",
       "    </tr></tbody> \n",
       "</table> "
      ],
      "text/plain": [
       "<pandas.io.formats.style.Styler at 0x7f192a45edd8>"
      ]
     },
     "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAD9BJREFUeJzt3X+snmV9x/H3Z1RUcFp+HBi2zYqz\n/prRQTpEycxm1QkYyh+SsDhtXJMmC1MUN60zGcmWLLgZUbOFhVCkZgQlyEajTG0KziwZzILKD6tr\ng6w9Uukx/NBJnBK/++O5KofT057yPOV5DlzvV3Ly3Pd1X/dzf89zfnzOdd33/ZxUFZKk/vzapAuQ\nJE2GASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq1JJJF3AoJ554Yq1cuXLSZUjS\nM8odd9zxo6qaWqjfog6AlStXsn379kmXIUnPKEn+53D6OQUkSZ0yACSpUwaAJHXKAJCkThkAktQp\nA0CSOmUASFKnDABJ6pQBIEmdWtR3Ao9q5cYvDb3v/ZedewQrkaTFxxGAJHXKAJCkThkAktQpA0CS\nOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlT\nBoAkdcoAkKROGQCS1CkDQJI6tWAAJLk6yb4k98xqOz7J1iQ72+NxrT1JPp1kV5K7kpw+a591rf/O\nJOuenk9HknS4DmcEcA3wtjltG4FtVbUK2NbWAc4GVrWPDcAVMAgM4FLgdcAZwKX7Q0OSNBkLBkBV\nfR14aE7zWmBzW94MnD+r/bM1cBuwNMkpwB8CW6vqoap6GNjKgaEiSRqjYc8BnFxVewHa40mtfRmw\nZ1a/6dZ2sHZJ0oQc6ZPAmaetDtF+4BMkG5JsT7J9ZmbmiBYnSXrCsAHwYJvaoT3ua+3TwIpZ/ZYD\nDxyi/QBVdWVVra6q1VNTU0OWJ0layLABsAXYfyXPOuCmWe3vblcDnQk82qaIvgK8Nclx7eTvW1ub\nJGlClizUIcl1wO8DJyaZZnA1z2XA9UnWA7uBC1r3m4FzgF3AY8B7AKrqoSR/A3yj9fvrqpp7YlmS\nNEYLBkBV/dFBNq2Zp28BFx3kea4Grn5K1UmSnjbeCSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6\nZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMG\ngCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROjRQAST6Q5N4k9yS5\nLsnzkpya5PYkO5N8PsnRre9z2/qutn3lkfgEJEnDGToAkiwD3gesrqpXA0cBFwIfAy6vqlXAw8D6\ntst64OGqeilweesnSZqQUaeAlgDPT7IEOAbYC7wJuKFt3wyc35bXtnXa9jVJMuLxJUlDGjoAquoH\nwMeB3Qx+8T8K3AE8UlWPt27TwLK2vAzY0/Z9vPU/Ye7zJtmQZHuS7TMzM8OWJ0lawChTQMcx+Kv+\nVODFwLHA2fN0rf27HGLbEw1VV1bV6qpaPTU1NWx5kqQFjDIF9Gbg+1U1U1W/AG4E3gAsbVNCAMuB\nB9ryNLACoG1/EfDQCMeXJI1glADYDZyZ5Jg2l78G+A5wK/CO1mcdcFNb3tLWadtvqaoDRgCSpPEY\n5RzA7QxO5t4J3N2e60rgw8AlSXYxmOPf1HbZBJzQ2i8BNo5QtyRpREsW7nJwVXUpcOmc5vuAM+bp\n+zPgglGOJ0k6crwTWJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKn\nDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoA\nkKROGQCS1CkDQJI6ZQBIUqcMAEnq1EgBkGRpkhuSfDfJjiSvT3J8kq1JdrbH41rfJPl0kl1J7kpy\n+pH5FCRJwxh1BPAp4MtV9QrgtcAOYCOwrapWAdvaOsDZwKr2sQG4YsRjS5JGMHQAJHkh8EZgE0BV\n/byqHgHWAptbt83A+W15LfDZGrgNWJrklKErlySNZJQRwEuAGeAzSb6Z5KokxwInV9VegPZ4Uuu/\nDNgza//p1iZJmoBRAmAJcDpwRVWdBvyUJ6Z75pN52uqATsmGJNuTbJ+ZmRmhPEnSoYwSANPAdFXd\n3tZvYBAID+6f2mmP+2b1XzFr/+XAA3OftKqurKrVVbV6ampqhPIkSYcydABU1Q+BPUle3prWAN8B\ntgDrWts64Ka2vAV4d7sa6Ezg0f1TRZKk8Vsy4v7vBa5NcjRwH/AeBqFyfZL1wG7ggtb3ZuAcYBfw\nWOsrSZqQkQKgqr4FrJ5n05p5+hZw0SjHkyQdOd4JLEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjpl\nAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaA\nJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1auQASHJUkm8m+WJbPzXJ7Ul2\nJvl8kqNb+3Pb+q62feWox5YkDe9IjAAuBnbMWv8YcHlVrQIeBta39vXAw1X1UuDy1k+SNCEjBUCS\n5cC5wFVtPcCbgBtal83A+W15bVunbV/T+kuSJmDUEcAngQ8Bv2zrJwCPVNXjbX0aWNaWlwF7ANr2\nR1t/SdIEDB0ASd4O7KuqO2Y3z9O1DmPb7OfdkGR7ku0zMzPDlidJWsAoI4CzgPOS3A98jsHUzyeB\npUmWtD7LgQfa8jSwAqBtfxHw0Nwnraorq2p1Va2empoaoTxJ0qEMHQBV9ZGqWl5VK4ELgVuq6p3A\nrcA7Wrd1wE1teUtbp22/paoOGAFIksbj6bgP4MPAJUl2MZjj39TaNwEntPZLgI1Pw7ElSYdpycJd\nFlZVXwO+1pbvA86Yp8/PgAuOxPEkSaPzTmBJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjpl\nAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaA\nJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1augASLIiya1JdiS5N8nFrf34\nJFuT7GyPx7X2JPl0kl1J7kpy+pH6JCRJT90oI4DHgQ9W1SuBM4GLkrwK2Ahsq6pVwLa2DnA2sKp9\nbACuGOHYkqQRDR0AVbW3qu5syz8BdgDLgLXA5tZtM3B+W14LfLYGbgOWJjll6MolSSM5IucAkqwE\nTgNuB06uqr0wCAngpNZtGbBn1m7TrW3uc21Isj3J9pmZmSNRniRpHiMHQJIXAF8A3l9VPz5U13na\n6oCGqiuranVVrZ6amhq1PEnSQYwUAEmew+CX/7VVdWNrfnD/1E573Nfap4EVs3ZfDjwwyvElScMb\n5SqgAJuAHVX1iVmbtgDr2vI64KZZ7e9uVwOdCTy6f6pIkjR+S0bY9yzgXcDdSb7V2v4SuAy4Psl6\nYDdwQdt2M3AOsAt4DHjPCMeWJI1o6ACoqv9g/nl9gDXz9C/gomGPJ0k6srwTWJI6ZQBIUqcMAEnq\nlAEgSZ0yACSpUwaAJHXKAJCkThkAktSpUe4EflZbufFLQ+97/2XnHsFKJOnp4QhAkjplAEhSpwwA\nSeqUASBJnTIAJKlTBoAkdcrLQJ8GXkIq6ZnAEYAkdcoRwCLj6EHSuDgCkKROGQCS1CmngJ5FRpk+\nAqeQpN44ApCkThkAktQpA0CSOmUASFKnDABJ6pRXAelXvAlN6svYAyDJ24BPAUcBV1XVZeOuQUee\n4SE984w1AJIcBfwj8BZgGvhGki1V9Z1x1qHFxfCQJmPc5wDOAHZV1X1V9XPgc8DaMdcgSWL8U0DL\ngD2z1qeB1425Bj2LjHr386Q8U0cujtaemsX+eo07ADJPWz2pQ7IB2NBW/zfJ9572qhZ2IvCjSRdx\nENY2nInWlo8dcvNifd1GqmuBz3lUi/U1gyFrG/H1+s3D6TTuAJgGVsxaXw48MLtDVV0JXDnOohaS\nZHtVrZ50HfOxtuFY21O3WOsCaxvWuM8BfANYleTUJEcDFwJbxlyDJIkxjwCq6vEkfwZ8hcFloFdX\n1b3jrEGSNDD2+wCq6mbg5nEfd0SLakpqDmsbjrU9dYu1LrC2oaSqFu4lSXrW8b2AJKlTBsAhJFmR\n5NYkO5Lcm+TiSdc0W5KjknwzyRcnXctsSZYmuSHJd9tr9/pJ17Rfkg+0r+U9Sa5L8rwJ1nJ1kn1J\n7pnVdnySrUl2tsfjFlFtf9++pncl+ZckSxdLbbO2/XmSSnLiYqotyXuTfK997/3dJGqbjwFwaI8D\nH6yqVwJnAhcledWEa5rtYmDHpIuYx6eAL1fVK4DXskhqTLIMeB+wuqpezeBChAsnWNI1wNvmtG0E\ntlXVKmBbW5+Eaziwtq3Aq6vqNcB/Ax8Zd1HNNRxYG0lWMHibmd3jLmiWa5hTW5I/YPCOB6+pqt8G\nPj6BuuZlABxCVe2tqjvb8k8Y/CJbNtmqBpIsB84Frpp0LbMleSHwRmATQFX9vKoemWxVT7IEeH6S\nJcAxzLkPZZyq6uvAQ3Oa1wKb2/Jm4PyxFtXMV1tVfbWqHm+rtzG4j2fsDvK6AVwOfIg5N5eO00Fq\n+1Pgsqr6v9Zn39gLOwgD4DAlWQmcBtw+2Up+5ZMMvtl/OelC5ngJMAN8pk1PXZXk2EkXBVBVP2Dw\n19duYC/waFV9dbJVHeDkqtoLgz9AgJMmXM/B/Anwb5MuYr8k5wE/qKpvT7qWebwM+L0ktyf59yS/\nO+mC9jMADkOSFwBfAN5fVT9eBPW8HdhXVXdMupZ5LAFOB66oqtOAnzK5aYwnafPpa4FTgRcDxyb5\n48lW9cyT5KMMpkevnXQtAEmOAT4K/NWkazmIJcBxDKaR/wK4Psl8b4szdgbAApI8h8Ev/2ur6sZJ\n19OcBZyX5H4G76j6piT/PNmSfmUamK6q/SOlGxgEwmLwZuD7VTVTVb8AbgTeMOGa5nowySkA7XHR\nTBcAJFkHvB14Zy2ea8h/i0Gof7v9TCwH7kzyGxOt6gnTwI018F8MRu0TOUk9lwFwCC2lNwE7quoT\nk65nv6r6SFUtr6qVDE5i3lJVi+Iv2ar6IbAnyctb0xpgsfy/h93AmUmOaV/bNSySE9SzbAHWteV1\nwE0TrOVJ2j9z+jBwXlU9Nul69ququ6vqpKpa2X4mpoHT2/fiYvCvwJsAkrwMOJpF8sZ1BsChnQW8\ni8Ff2N9qH+dMuqhngPcC1ya5C/gd4G8nXA8AbVRyA3AncDeD7/+J3aWZ5DrgP4GXJ5lOsh64DHhL\nkp0MrmiZyH/MO0ht/wD8OrC1/Sz80yKqbVE4SG1XAy9pl4Z+Dli3WEZP3gksSZ1yBCBJnTIAJKlT\nBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq1P8DgbsPmQtEQZEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f192a52fba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(weights, bins=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the main model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "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>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "crosstab = pd.crosstab(subset.age, subset.qsmk).sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11833b518>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8lNXZ939n9mRmssxkQhISSMIi\nJCyBAKK4VgXBHUHrVrV9a21rl6eV6tPtsevja23tW61abau2glXAXdS6FAQlYAJhCwgkBLJCMlkm\nyWT28/5x5p5MJrPcsw8z5/v5zCfJzD1zn0wy57qv7XcRSik4HA6Hk7lIkr0ADofD4SQXbgg4HA4n\nw+GGgMPhcDIcbgg4HA4nw+GGgMPhcDIcbgg4HA4nw+GGgMPhcDKcuBkCQkgZIeQ/hJAmQsghQsj3\n3Pc/RAjpIIQ0um+r4rUGDofD4YSGxKuhjBBSDKCYUrqHEKIF0ADgegA3ARimlD4alxNzOBwOJyxk\n8XphSmkXgC7390OEkMMAJkfyWgUFBbS8vDyGq+NwOJz0p6GhoZdSagh1XNwMgTeEkHIACwDsArAM\nwH2EkK8AqAfwQ0ppf7Dnl5eXo76+Pt7L5HA4nLSCEHJSzHFxTxYTQjQANgP4PqXUBOApANMA1IB5\nDL8P8Lx7CCH1hJD6np6eeC+Tw+FwMpa4GgJCiBzMCKynlL4KAJTS05RSJ6XUBeBZAEv8PZdS+gyl\ndBGldJHBENKz4XA4HE6ExLNqiAD4G4DDlNI/eN1f7HXYDQAOxmsNHA6HwwlNPHMEywDcAeAAIaTR\nfd+PAdxCCKkBQAG0AvhGHNfA4XASgN1uR3t7OywWS7KXkpGoVCqUlpZCLpdH9Px4Vg3tAED8PLQl\nXufkcDjJob29HVqtFuXl5WDBAE6ioJTCaDSivb0dFRUVEb0G7yzmcDhRY7FYoNfruRFIAoQQ6PX6\nqLwxbgg4HE5M4EYgeUT73qe1Ifjo8Gk8ufV4spfB4XA4KU1aG4Ltx3rx1NbmZC+Dw+FkEKtWrcLA\nwAAGBgbw5JNPeu5vbW3Fhg0bkriywKS1IcjPVmDI4oDd6Ur2UjgcToawZcsW5OXlxcQQOByOWC/P\nL+ltCNSslGrAbE/ySjgcTjxpbW3F7Nmz8fWvfx3V1dVYvnw5RkdH0djYiKVLl2LevHm44YYb0N/P\n1GwuueQSPPDAA1iyZAlmzpyJ7du3+33dhoYGzJ8/H/Pnz8e6deswZ84cAMDzzz+P++67z3Pc1Vdf\nja1btwJgkji9vb148MEH0dzcjJqaGqxbtw4PPvggtm/fjpqaGjz22GNwOp1Yt24dFi9ejHnz5uEv\nf/kLAGDr1q248MILce2116KqqgojIyO46qqrMH/+fMyZMwcvv/xyzN+/hGgNJYu8bAUAYMBsg0Gr\nTPJqOJzM4BdvHUJTpymmr1lVkoP/uaY66DHHjh3DSy+9hGeffRY33XQTNm/ejEceeQSPP/44Lr74\nYvz85z/HL37xC/zxj38EwK62d+/ejS1btuAXv/gFPvzwwwmveffdd+OJJ57ARRddhHXr1oW15ocf\nfhgHDx5EYyNro9q6dSseffRRvP322wCAZ555Brm5ufj8889htVqxbNkyLF++HACwZ88eHDx4EBUV\nFdi8eTNKSkrwzjvvAAAGBwfDWocY0toj0LkNQd+ILckr4XA48aaiogI1NTUAgNraWjQ3N2NgYAAX\nX3wxAODOO+/EJ5984jl+9erVnmNbW1snvJ4Q57/ooosAAHfccUdM1/vvf/8b//jHP1BTU4Nzzz0X\nRqMRx44dAwAsWbLE0xMwd+5cfPDBB3jggQewfft25ObmxnQdQNp7BCw01M9DQxxOwgh15R4vlMox\nr18qlWJgYEDU8VKp1BOLv/vuu7F3716UlJQEjefLZDK4XGO5x0hq+CmlePzxx7FixYpx92/duhVq\ntdrz88yZM7Fnzx5s2bIFP/3pT3HZZZfh5z//edjnC0ZaewT56rHQEIfDySxyc3ORn5/vif//85//\n9HgHgXjuuefQ2NjoSfjm5eVhx44dAID169d7jisvL0djYyNcLhfa2tqwe/fuCa+l1WoxNDQU8OcV\nK1bgqaeegt3OLlSPHj2KkZGRCa/T2dmJ7Oxs3H777Vi3bh327NkTxrsgjrT2CDyhIW4IOJyM5IUX\nXsC9994Ls9mMyspKPPfcc2E9/7nnnsNXv/pVEEI88XsAWLZsGSoqKlBVVYXZs2dj4cKFE56r1+ux\nbNkyzJkzBytXrsRvf/tbSKVSzJ8/H3fddRe+973vobW1FQsXLgSlFAaDAa+//vqE1zlw4ADWrVsH\niUQCuVyOp556Kvw3IgRxG1UZSxYtWkQjHUxzzk/fxZ3nl+PHq2bHeFUcDkfg8OHDmD07vT9jra2t\nuPrqq3HwYGoKJvv7GxBCGiili0I9N61DQwDrJejnyWIOh8MJSPobArWCJ4s5HE7UlJeXp6w3EC3p\nbwiy5ejnOQIOh8MJSAYYAgU3BBwOhxOE9DcEajmXmOBwOJwgpL8hyFZgwGyDy5X61VEcDoeTDNLe\nEORlK+CigMnCvQIOh8Oor6/Hd7/7XQCsk/ezzz7zPPb666+jqakpWUtLCmlvCPLTWGbi/753BN98\nsSHZy+BwzjoWLVqEP/3pTwBiYwgSJRcdL9LfEKjTV3iu8dQAdrYYk70MDicl+M1vfoOZM2figgsu\nwC233IJHH30Ul1xyCYRm1N7eXpSXlwNgm//VV1+N1tZWPP3003jsscdQU1ODbdu24c0338S6detQ\nU1OD5uZmNDc348orr0RtbS0uvPBCHDlyBABw11134d5778W5556LH/3oR9i2bRtqampQU1ODBQsW\njJOTSHXSWmICYDkCID31hvrNNgyY7Riy2KFVyZO9HA6H8e6DQPeB2L5m0Vxg5cMBH25oaMC//vUv\nNDY2wuFwYOHChaitrQ35suXl5bj33nuh0Whw//33AwCuvfZaXH311VizZg0A4LLLLsPTTz+NGTNm\nYNeuXfjWt76Fjz/+GADQ3t6Ozz77DFKpFNdccw3+/Oc/Y9myZRgeHoZKpYrBL54YMsAQpG9oSCiL\n7RgYxawibgg4mcv27dtxww03IDs7GwDbzGPB8PAwPvvsM6xdu9Zzn9Vq9Xy/du1aSKVSAEx/6Ac/\n+AFuu+02rF69GqWlpTFZQyJIf0PgDg2lm8wEpdRj3Nr7RjGrKCfJK+Jw3AS5ck803nLRkUhFu1wu\n5OXleYbL+OItF/3ggw/iqquuwpYtW7Bs2TK8//77mDVrVmQLTzBpnyPQKmWQSUjaNZWN2p2wOdg/\neFu/Ocmr4XCSy0UXXYTXX38do6OjGBoawltvvQWAhX4aGlhBxaZNm/w+N5hcdE5ODioqKrBx40YA\n7AJs3759fl+nubkZc+fOxQMPPIDFixd7cglnA2lvCAghyMuWp11oyDv53d4/msSVcDjJZ+HChbj5\n5psxf/58rFy5EosXLwYA3H///XjqqaewYMEC9Pb2+n3uNddcg9deew01NTXYvn07vvzlL+N3v/sd\nFixYgObmZqxfvx5/+9vfMH/+fFRXV+ONN97w+zp//OMfMWfOHMybNw9yuRwrV66M2+8ba9JehhoA\nrvjDNkwv1OCp20Mnj84WDnYM4urH2cCMFdWT8Jc7QirNcjhxI9VkqB966KFxCeBMgMtQhyA/W5F2\n5aNCqEurknGPgMPhREXaJ4sBNrv4pDG94uiCYZtTkoumLlOSV8PhpBYPPfRQspdwVpERHoFOnX4K\npIKQ3tzSXAyO2rmEBifpnA1h5nQl2vc+IwxBnluKOp3+UQXDVl3Cykbb+3h4iJM8VCoVjEZjWn3G\nzhYopTAajVE1sGVEaCg/Ww67k2LE5oRGmR6/8oDZjhyVDOV6Vsfc3m9GVQnvJeAkh9LSUrS3t6On\npyfZS8lIVCpVVA1s6bErhsC7qSxdDEHfiA35agVK87MA8BJSTnKRy+WoqKhI9jI4ERK30BAhpIwQ\n8h9CSBMh5BAh5Hvu+3WEkA8IIcfcX/PjtQYBQW8onfIE/WYb8rIV0KkVyFZIuSHgcDgRE88cgQPA\nDymlVQCWAvg2IaQKwIMAPqKUzgDwkfvnuJKOekMDZjt02XIQQlCan4V23l3M4XAiJG6GgFLaRSnd\n4/5+CMBhAJMBXAfgBfdhLwC4Pl5rEMiLgQLpa3vbcfzMcKyWFDV9IzaPp1Oan4027hFwOJwISUjV\nECGkHMACALsATKKUdrkf6gYwKcBz7iGE1BNC6qNNQOminEnQ1mfGf728D3/d3hLVOmLJgDs0BIB7\nBBwOJyribggIIRoAmwF8n1I6rvOJslozv/VmlNJnKKWLKKWLDAZDVGvIzZKDkMhDQ5v3tAMAWnpG\nolpHrLA6nBixOaFTs5BXaX4WhiwODI6mT+iLw+EkjrgaAkKIHMwIrKeUvuq++zQhpNj9eDGAM/Fc\nAwBIJQQ5KnlEoSGXi44Zgt7UMARCM5ngEZTlMw127hVwOJxIiGfVEAHwNwCHKaV/8HroTQB3ur+/\nE4B/Kb8Yw7qLw79i3nWiD219o6gqzkHvsDUlOniF6ifvHAHAS0g5HE5kxNMjWAbgDgBfIoQ0um+r\nADwM4ApCyDEAl7t/jjt52fKIhtNsamiHVinDNy6uBJAa4aH+EWaM8r1CQwDLZXA4HE64xK27ilK6\nAwAJ8PBl8TpvIPKzFThtCm9C0bDVgS0HunD9gskeKYcTvcOoKcuLxxJF4+sR5GXLoea9BBwOJ0Iy\nQmsIYJvmQJihoS0HujBqd2JNbSmm6NSQkBTxCHwMAeslyOaGgMPhRER66C2IID9bHnb56Kb6dlQa\n1Fg4JQ+EEJTpslMiYTyWLB4bWJ/oEtL+ERuGLA5M0Wcn7JzRUNdiDDs0WF2Se9b8fhxONGSOIVAr\nMGp3wmJ3QiWXhjy+tXcEu1v78KMrzwHLewOVBerU8AhGbMhWSMf9HmW6bOw+0QdKqWe98eTRf3+B\nz5qN+M/9l8T9XNFyymjGl5+pC/t5tVPzsfmb58dhRRxOapE5hsDTXWxHUW5oQ7B5TzskBFi9YEzR\nr9Kgwc4WI1wuCokk/pttIPrMY13FAqX5WRiyOmAadSDXy1OIF6dNVnQMjCbM8ETDKXcS/Xdr5mFu\naa6o5zz2wVHsbx+M57I4nJQhgwwB2xz7Rmwoyg2u2+1yUWxuaMeFMwzjjq0oUMNid6HbZEFJXlZc\n1xuMAbN9XFgI8Koc6jcjN1vcZhcNw1Y7bA4Xhq0OaFXxNzzR0O0uElhcrkN5gVrUcyoNGnx85EzS\njT6HkwgyJlkcjt7QZ81GdA5asHbReH3vSgPbRJIdHuo32zyyGQKlCW4qG7I4AEQu25FIugdZEj3U\nBYA3hVol7E7Ku7U5GUHGGAKh5l5MU9mmhjbkqGS4fPZ4GaRpBg0AoKU3ueJz/SNjOkMCiZ5LIBiC\n3uGzwBCYLMjLlovKDQkYtEoAwJkha7yWxeGkDBljCHTujbMvhEdgstjx7sFuXFczecLGUahVQq2Q\npoBHYPeEugRys+TQKmUJNATMoJ4dHoEFRTnhjfEzaJgh6OGGgJMBZIwh8ISGQmxcb+/rgtXhwpra\niWPfCCGoMKiTWkLqcLpgstgnJIsJIZicoBJSSqlXaCj1N8quQUtYYSFgzCPoGQ6vCZHDORvJGEOg\nkEmgVkhDhoY2NbRh5iQN5gWoLqko0OBEEkNDg6N2UIoJHgHgnkuQgCH2FrsLDhcTjT0bQkOnTRYU\nR2oIwvAI+OB2ztlKxhgCgPUSBEsWt/WZsefUAG5cWBqwJLKyQI32/lFY7M54LTMogiHL90kWA2NN\nZfHekIasY8Y01UNDVocTvcM2TAozNKRRyqCSS0QbgoMdg5j1s/e4AiznrCSzDEG2ImiOQKgbXza9\nIOAxlQY1KAVOGpPzgRcMmW+yGGCGYMTmDFtKI1yEsBAAGIdTOzR0xsTWF65HQAiBQasMyxBYHS40\np0DDIYcTLhllCPKy5UFDQ01dg5BJCGZM0gQ8prKAPZas8JBwBa7zawgSI0c9zhCkuEcg9BAU5Ybf\n92HQKNEj0tB1DbLzpLph5HD8kVGGQBciNNTUacL0Qg2UssBlhhXuXoJkXfn50xkSKNMJJaTx9VaE\niiGtSgZjiucIut0bdLhVQwDC8ggEZdtUD5VxOP7IKEOQn60I+kE91GlClVtuOhAapQyTcpQ4kaTK\nIY/yqN8cQWI9gooCdcpvfB5DEGZoCAAKtSrRfQQejyDF3w8Oxx8ZZQjysuUYsjjgcLomPNYzZMWZ\nISuqioMbAoBtgC09yQkN9ZvtUEhZBZQvuVlyaFUytCXIIyjXq2EcsaZ0tUzXoAVZcilyVOGrqRi0\nSgyY7bA6QhcGCB4BDw1xzkYyyhAIsgwDfmQDDneZACCkRwAwHZpk9RKwrmJ5wKqmRMwlEDyCcn02\n7E6KIasjxDOSh1A6GokwnlBCKib8JXgEqe4hcTj+yChDIFTa+NOlb3Ibguri0IJtlQVqDJjtEY2+\njJZ+P8qj3iRiLoFgCMp0LBSVynmCrsHRsEtHBcR2F4/anB5NorOhr4LD8SWjDIHQhOWvcuhQpwmT\n87JESTh7xOeSUDk0YLZ7dJP8Ueb2COIZrhmyOKBWSFHo3mBTubv4tMkadumogNimMqEySSYhcfUI\n6lv7cO8/G+B0pW4ojnN2kmGGwO0R+KkcauocFBUWAsZKSJNROeRvFoE3pflZMNucosT1ImXIYodW\nJYfeHWpLVY/A6aI4bQpfXkJgTGYiuCHocqubTi/UxDVH8OHhM3jvUDfODHHZC05sySxDoPYfGjLb\nHGjpHRGVKAbYZiuXkqRUDg2YJyqPeuOZS9AXv/DQkMUBrUrmybmkaqWMcdgKh4tGbAj0Gvb7hfII\nhERxdUkuRmzOuHWdCyE/LoTHiTWZZQgChIaOdA+BUqBapEcgk0owRZed8MohSqlf5VFvElFCOmS1\njzMEqZog9TSTRZgjUMqkyMuWh9x4hUSx8P8TL8Mo/E25IeDEmowyBFlyKRQyyYSmsqZO8RVDApUG\nTcI9ApPFAaeLThhK483k/Pg3lTGPgOn7a5Qy9KZoyWRXFD0EAoUimsq6By3IUck83lhfnEJl3BBw\n4kVGGQJCCHTZigk5gqYuE3Kz5JgcxvjJygI1Wo3mhCbugukMCeRmyZGjiu9cAiE0BLDwSap6BKdN\n0RsCg1YZMibf7Za51rurjHrjkDwftTk9BpcbAk6syShDALCmsr6R8aGhpk4Tqopzwqo1rzSoYXO4\n0JGgQTCAl/JoiMqmMl12AjwCZgh0akXKJou7Bi2QSQgK1MqIX0OM3lC3yYKi3CxP8jweHkHHwNjf\nU6z+EYcjlowzBPnZ4/WGnC6KI92hpSV8qUzC2Eohye1PXsIb1ksQT4/A7hlYr1crUzZZfHrQgkk5\nqqiGzwt6Q8HKcbsHLSjOUXmSy/HwkNq8/p6CoiqHEysyzhDo1ONDQyd6h2Gxu0RXDAlUFCR+kL1H\nZyhIaAhwD6iJ01wCm8MFq8MFrdIdGlIroi6Z/PuOE/iieygWyxtHJJPJfDFolbDYXRgO0D1td7rQ\nM2zFpFwVNEoZFFJJXEJDgmGfZlBzj4ATczLOEPhKUR+KIFEMsA0wRyVLaMJYbGhomkEDi90Vl7V5\nK48CYzmCSI1O1+Aofvl2Ex5570jM1ijQbQp/VrEvoZrKmLcAj4yFXqOIS2iovd8MhUyCqpJcniPg\nxJyMMwRCaMjlTvI2dZmgkEowvTDwDAJ/sPnFmoSGhgbMNkgIkKMKbgjOrdQBAOpa+mK+BkFeQggN\n6dQKOFwUptHI9IZ2ude49WhPTBulKKWeJG40GDTs+YE23y4fmWudWhGXUFl73yhK87JQlBM6VMXh\nhEvmGQK1Ai46tqE1dZows0gDuTT8t2JagTqhoaG+EdZMFirmXVmgRqFWiZ0txpivYcwQMI+gwF0p\nY4wwHLKz2QilTAKni+K1PR2xWSQA06gDo3Zn7DyCAOEYX5nruBmCfjMm52fBoFVi1O7EiC05o1I5\n6UnmGQJ3WKXPzMIZQsVQJFQa1OgatMBsS4z65oDZ7ncgjS+EECyt1KOuxRjzK0dhXrHGq2oIiDxB\nWnfCiItmGrBwSh42NbTHbL3dMSgdBVgfARDYIxDOI+gZFWiUcZGZaO8fRWl+tmj9Iw4nHDLQEIzp\nDZ0ZssI4YovYEFR4xlYmxisIpTzqzXnT9OgZssZcLlvwCHK8QkNAZKqbnQOjOGk0Y2mlHmsXleHY\nmWHsc8+NjhZB/ydSwTmB3Cw55FIScEBN9+AolDIJcrPG3o9YVw2ZbQ4YR2wozc8KGaricCIhboaA\nEPJ3QsgZQshBr/seIoR0EEIa3bdV8Tp/IITSywGzDYc62aZTPTm09LQ/PCqkCQoP9Y2INwRLK/UA\nWOgllgQKDUWy+dW5Q1dLK3W4al4xVHIJNjW0xWSdQsgmUglqAYmEoEATuLu4261uKvSg6DUKmG1O\njMYwdCP0qpS6Q0MANwSc2BJPj+B5AFf6uf8xSmmN+7Yljuf3iyc0NGL3SEvMKtJG9FpCCWmiPIKB\nEDpD3pTrszEpR+nZbGPFWNUQW4cgiR1JOKSuxYjcLDlmF+UgRyXHldVFeLOxMyaibULIJlpDAASf\nXdztM+/Ao8gawxLSdo8hGAsNcQVSTiyJmyGglH4CIPZlK1EiyDMMmG1o6jJhqj7bs6mFi0ouxeS8\nrISIzzHBOVvIZjIBQgjOq9SjrqUvpnkCX49AKZOyIfYReQR9OLdC50l+r6ktg8niwAdNp6NeZ/eg\nBQUaBRSy6P/FDUE8gq5By7jwk14duYcUCGH0aJkuC3lZcsgkhHsEnJiSjBzBfYSQ/e7QUX6ggwgh\n9xBC6gkh9T09PTE7eY5KBqmEoN9siypRLFBVkoOPj5yJ+1SwUbsTVodLdGgIYOGh3mErmmNoqIYs\ndqjkknFVVvoIKmU6BkZxqs/sCWEBwPnT9Jicl4WNDe1Rr7M7ijkEvhi0/mUmXC6KMyYrinLHNKp0\nmtjPaGjvZ3kIg0YZMlTF4URCog3BUwCmAagB0AXg94EOpJQ+QyldRCldZDAYYrYAQgjys+Vo6xtF\nq9EsWno6ED9ZNRuUAt/esBc2hytGq5yI2GYybzx5ghj2EwjKo97oNcqwp5TVuXMX500bMwQSCcGN\nCydjx7EeT4w/UroHLSjKES8iGAyDllUC+QoM9pltsDldKMoZ0zLSx2FGg1A6KuQhCnNC6x9xOOGQ\nUENAKT1NKXVSSl0AngWwJJHnF8jLVnhq7MPtKPalvECNR9bMw762ATz8buy7YwUEnaFgyqO+TNVn\nozhXFdM8wZDV4ZGXEIhEeK6uxYi8bDnOmTQ+P3NjbSlcFNi8JzqvgHkEkYvNeWPQKuGiE8M9Yz0E\nYwZHUCCNZQmpUDrqWQ/3CDgxJqGGgBBS7PXjDQAOBjo2nuR7DRupEjGsPhQr5xbjrvPL8fdPT+C9\ng91Rv54/BJ2hYLMIfBH6CXbFsJ/AW3lUoEATfmhoZ4txXH5AYKpejSUVOmyOoqdg1ObEgNmO4tzY\neASBegl8m8kAQK1gMy9imSNghmDsdwmWvOZwIiGe5aMvAdgJ4BxCSDsh5GsAHiGEHCCE7AdwKYD/\nitf5gyHE2XVqBSblxOaq8b9XzcL80lys27QPp4yxzxdEEhoCWGlm77ANx8/EJk/grTwqINTOu0TO\nZmjrM6O9fxTneeUHvFlTW4qW3hHsOdUf0RpjWTEEIGClTpdPMxnAjG+BWhFRX4U/RqwO9Ll7CLzX\nYxyx8SH2nJgRz6qhWyilxZRSOaW0lFL6N0rpHZTSuZTSeZTSaymlXfE6fzAEQ1BdEt4MgmAoZVI8\ncetCEADf3rAHVkdsJQDEDKXxx3mVBQAQs/CQP49Ap1bC6aIwWewBnjWeXSdYzmLpNP+G4Kq5xchW\nSLEpwqSxcKUebTOZQKAmrtODFkjdyVtvdBpF2DmTQAilo2XeoSEte79TdSAQ5+wj4zqLASDPXfse\nbcWQL2W6bPxu7Xwc6BjEb985HNPX7ncP0xEjMTF+TVkoyVXFTHeIeQQTQ0OA+ATpzmYj8rPlmFno\nv39DrZRh5ZxivLWvK6LGrG4T2zxjVTVUoHUPsfeJ+3cNWlCoVULqE96K5YwGoRptnEeg4U1lnNiS\nkYZA576qjjZR7I8V1UX42gUVeGHnSbyzPzyHZ8hiDxgX7zfboFXJwhbHG9MdCt5PMGx1iIrJ+6sa\nEvIWYhPGdS1GLK3UBxXPW7uoFMNWB94/FH7OpXuQbZDRCs4JZCtk0ChlEz2CACWq+hhObfNuJhMI\nJYTH4YRLRhoC4cM7N0JpiVA8cOUsLJiShx+/dgB2p7iS0jNDFiz+zYcBwyHh6Az5snSaHn0jNhwL\nkCc42DGIxb/+EK/UB5d3cDhdMNuc0PhUDY01UYXemNr6zOgYGB3XP+CPJeU6lOmy8Na+zpCv6Uv3\n4Ci0KhnUPuuMBn8J2q7BUb/GJpZznNv7zVDKJB6vCwAKtVxviBNbMtIQrJxTjM3fPN8zbjLWKGQS\n3HNhJQZH7dgvUkRtZ7MRFrsLL+465ffx/jDkJXw5L4ju0JDFjvs27MGo3YnmEJpJI1YWpvENDQkj\nGsUkSHd69IWCGwKJhKCmLB9fnA5/clnXYPQDaXzxV7IZaN6BTs2komOhSitUDHnnsjyhKm4IODEi\nIw2BQiZB7dSATc0xYUmFMBxGXGxeOG5f2wCO+dn8+kfEy0v4Upqfhcl5WRPWQinFg5sPoK1/FNkK\naciNRUgG+w7GETwVMVfBdS1G6NQKzBAxCKiyQI2OgdGwtYcChWyiwbe7eMhix4jN/7wDfZihsmC0\n9ZvHhYWAwKEqDidSMtIQJAK9RolzJmlFG4KdzUYsnJIHmYT4DQ9FExry9BOc6BtX4vnPupN450AX\n7l9+Ds4p0obcWHx1hgQUMglyVLKQTVSUUtQ1G7G0cmL/gD8qDWpQCpwMsxw3Lh6BT2jIXw+BgD7M\n5Hkw2vtHUaab2A8RSPaCw4kEbgjiyHnT9Khv7Q8pPdE1yOQuVs0txqWzCvHq3g44fHILYofSBGJp\npQ59IzYcPcO8jQPtg/j124eawX0qAAAgAElEQVRx6TkGfOOiSlHdqr7Ko94UaEJXyrT1jaJz0BIy\nLCRQ6Zn3IL4HQhgmH6vSUQGDVokhi8PjnYwNpJm4SY8N64luox6y2DFgtk/wCAAhVMUVSDmxQZQh\nIIR8JOY+zniWVuowanfiQMdA0OOEub1LK/VYW1uKniErth0dE9qzOVwYtjo81U6RrYVtvnXNRgyO\n2vGtDQ0o0Cjwh5tqIJEQUVeYgTwCQJzMRJ3I/IBAhXveQ6jchTfCMPmiGHUVC/jOAfCdVeyN0FcQ\nbVNZx8DYHAJ/6wk0LIfDCZeghoAQoiKE6AAUEELyCSE6960cwORELPBsZkmFuOEwO5uNyFHJMLs4\nB5fOKoRerRgXHvI0k0WYIwBYj0NpfhZ2thjxo0370DVgweO3LvTkHQq1KvSN2IJWOQljKv0ZAjGV\nMjtbjNCLzA8AgEYpQ6FWGdbgn7ERlbHpGBcY6y5mm68QGir005ke7fhOgfa+iaWj3uvhOQJOrAjl\nEXwDQAOAWe6vwu0NAE/Ed2lnPzq1ArOKtKgLof5Zd8KIcyv1kEoI5FIJrl8wGR8ePu3ZSCKVl/Bl\naaUe/246jfcPncYDV84alzAXNrreIF7BsNsj0Pj1CJRBh7FQSj39A+F0c1ca1GGFhjyx+xgpjwr4\nNnF1myzQqRVQyaUTjs1WSKGSS6IWnvPXTOZZj0+oisOJhqCGgFL6/yilFQDup5RWUkor3Lf5lFJu\nCESwtFKP+pN9AfME3nN7BdbUlsLupHizsQPA2JVlpMli77VQClw+exL+z4UV4x4TMwLR5DOv2Bu9\nWoF+sz2g3tCpPjO6Bi1YWqkLa82VBk1Yc5e7giRxo6HQp4mrO0hCmhASk+7itv5RZMmlniokb/jI\nSk4sEZUjoJQ+Tgg5nxByKyHkK8It3otLB5ZW6mGxu7Cv3X+ewHtur8Ds4hzMmZzjGdAihIaiNQQr\n5xTh+5fPwO/Xzp9wVS5mYxmyOCCXEij9TP3SaxRwuigGR/3rDQnhMbH5AYHKAjUGzHaPDHcoTpss\nUMgkUXtPvujUChDi5RH4TCbzRa+Jvru4vd88oYdAgHcXc2KJ2GTxPwE8CuACAIvdt0VxXFfasLRS\nB0LGBrH44j2315s1C0txqNOEpk7TWGhIHd3mplbK8P3LZyLXzyYpzhAw5VF/G5MuxKzeuhYjCjRK\nTBeZHxCodCeMW0SGh4TRkbESExSQSSXQqxXjQkOTghgCQZE1Gnzlp73hekOcWCK2fHQRgGWU0m9R\nSr/jvn03ngtLF/KyFZhVlIO6E4EMQZ9fXf7raiZDIZVgU0O7ZxZBtB5BMAQJg1Aegb9EMXu+MJBl\n4ubH8gN9bqMY3gZd4S4hFVs5dHrQEjP5aV+EEZEWuxN9IzYUBzlP7AzBxEQxEHhGAocTCWINwUEA\nRfFcSDpzXiXrJ/CVphbm9p7nR445X63A5VWFeL2xAz1DVmTJpX4Tk7FCKZMiN0seNNTgT3lUQBdk\nRGOr0Yxuk/j+AW/K8rMglxKcEJkn6DKNxryHQKAwR4WeYSvOmNh7FMwjKNAo0TtsjXi4jslix+Co\nPaBHIISqeAkpJxaINQQFAJoIIe8TQt4UbvFcWDqxtFIHq8OFfW3jdYfqQsTN19SWom/Ehrf3d8U8\n5u2PUCWJQxbHBME5gWDdtOH2D3gjk0owRZeNlp7QoSFKKU4PWmPeVSxg0CjRY7Kga5CVdQYzODq1\nAlYHE+mLhGClo8DEUBWHEw1i5Rkfiuci0p1zK/QghCVMBQ0igNXV5/uZ2ytw0QyDZ3OO9ewEfxSG\nMATDVgfKdP43JiFs5a9kcmezEQatEtPc8f5wqSjQiPII+kbcw+Tj5BEITXfdfiaT+eKtNxSJCmqw\n0lGBAj67mBMjxFYNbfN3i/fi0oXcbDmqinMm6A7VtRhxbkVgXX6ZVILVC1jfXjiziiMlVLdqsByB\nXCpBXrZ8Qlw80v4Bb6YZ1Gg1mkOOZgzW7RsLDFol7E6KI91MpiNYLmLMQ4pso/ZMJgtgeIX18Koh\nTiwQWzU0RAgxuW8WQoiTEGKK9+LSiaWVeuw51e9pABLm9oaqq19TWwog/MlkkSDoDQWKa5ssdr89\nBAL+ZCZO9I7gzJA14HxiMVQUqGFzuNDh3hwDcdoUnx4CAaGy6mDHIDRKmV/NJQFhRkOkJaTtbkXY\nYCHBQq0Kvdwj4MQAsR6BllKaQynNAZAF4EYAT8Z1ZWnG0ko9rA4XGttYP4Enbh5gbq/AjEla3Hne\nVCyvjn+u3qBlOvojfuLaLhfFsDWwRwC4J3P5XAHv9NMnES7C3IhQJaRNnezapKIgshBUKISSzf3t\ng5jkR1rCm2hlJoL1EHjWow1uuDkcsYStPkoZrwNYEYf1pC1LKtz9BO6Nsa6lDzq1IuDcXm9+cd0c\nXDu/JN5LDNpLMGJzgFL/OkMCerVywhVwXUsfCrXKqDZnTy9BiBLSuhNGzC7OQV6cymyF92dw1O5X\nddQbz7CeKEJDgRLF3uuxOV0wjUY/AIeT2YjKYhFCVnv9KAHrK+AauGGQmyVHdUmOlyEw+u0fSCbe\nhsB34x5THg0SGtIo8HnrmCEQ8gPnT4s8PwAwT0OrkgVNGFsdTtS39uPWc6dEfJ5QCO8PEDr8lK2Q\nIUsuRV+EoaG2fjMWlQcfnjTWXWzx2yTI4YhFrEdwjddtBYAhANfFa1HpytIKPfacGsDxM0Oi5vYm\nmmAegWAIApWPAkCBWoF+s82T1G3uGUHPkDXq35MQ4tYcChwa2tc2CKvDFVUuIhQ5KplHXkNMQjrS\nprLBUTuGLI6gFUPAWKhK6GvgcCJFlEdAKb073gvJBM6bpsdfd5zAk1ubPT+nEmND0Sc6e8NBJKgF\ndGoFXJRpI+k1yqj6B3yZVqD25Bv8UddiBCEYV54bawhhcxva+0dFJaQLNAr0RmAIhNLRMhGhIYDr\nDXGiR2zVUCkh5DVCyBn3bTMhpDTei0s3FlfoICHA63s7wtLlTxR5WXLIJMTvxmISERrSu69Qhavg\nuhYjinJUKNcH39DEUFGgRtegJeBA+J3NRswuil9+QEDYfMV7BOFv0kLpqJgcAcBlJjjRIzY09ByA\nNwGUuG9vue/jhEGOSo45k3Phooiqrj5eSCQEBRql31DDkEeCOnjVEMAmc0WjL+QPoXLIX57AYndi\nz6n+hHhYQjhGjEeg10xMnothzBAEDw3lqGRQyCTcEHCiRmzLo4FS6r3xP08I+X48FpTuLK3UY3/7\nYFTllPEkUJNSsHnFAt4eQXPPMHqHrTHbnIXk9YneEVSX5I57bF/bAKwOV0JyLsJVuBg9I1ZOy4yi\nP2NY12L0lLx689Hh01ArpCF7RwghIbvBORwxiDUERkLI7QBecv98C4Dg8xc5flleNQkv1p3ExTML\nk70Uvxi0Sk9jljfB5hULeEtR72xhV8Kx2pwFQ+CvhHSnkB8oj79xnVeaix3Hs0V1eus1Cs+8aV8D\nOmJ14KvPfx5Qi0hspRXvLubEArGG4KsAHgfwGAAK4DMAd8VpTWnNonIdmn55ZbKXERCDRomDHYMT\n7h+y2CEhbAxjIPKz5SCEddMePzOM4lwVpgSRSAiHLIUUk/Oy/IrP1bUYUV2Sk5ASypsXT8HNi8WV\nqOrUYx6SryHYcqALZpsT//jqEswvzZvwXH/jQP1h0Chxqs8s6lgOJxBicwS/BHAnpdRAKS0EMwy/\niN+yOMnCoGUjFn11fQTl0WBXqTKpBHlZcvQOW1HXYsR5Mc6DVBSoJ+QIWH5gAEsrUqsCC/BqKvOT\nJ9jU0I5yfTYunFGA3Gz5hJtUZH9JKH0oDkcMYg3BPEppv/ADpbQPwIL4LImTTApzlHC6qGcYjsCw\nZWJ4wx86tQK7T/TBOGKLecy+0qBGS8/IOEmFvacGYEtQfiBc9AFkJk4Zzdh1og9rakujNpQGrRJ9\nIzbYnf5nYnM4YhBrCCSEEE+bIyFEB/FhJc5ZRKARiKYgyqPe6DVKHDvDwjexruKpKFBjyOoYd4Vd\n12KEhLDS3FRjTG9o/Hu5aU87CAFWL4y+AltIXkc7H5mT2Yg1BL8HsJMQ8itCyK/AcgSPBHsCIeTv\n7p6Dg1736QghHxBCjrm/Bu+h5yScQLXpQyGURwWEq+DJeVkhyx/DxSM+55UnYPmBXORmpZ7EgqBA\n6m24XC6KzQ3tuGB6AUryon9/+OxiTiwQqz76DwCrAZx231ZTSv8Z4mnPA/DNij4I4CNK6QwAH7l/\n5qQQgiHwjTsHm0XgjRAXPzdG/QPeVAqVQ+48gcXuxN5TAylbipulkCJbIR0XGqprMaJjYNQjLx4t\n3npDHE6kiA7vUEqbADSFcfwnhJByn7uvA3CJ+/sXAGwF8IDY1+TEn4IAV5hDVjtmqkJ3QguVMvHQ\n/CnJy4JCJvEkjPec6ofN6Uo5qQ5v9BrFuKltGxvaoVXJsCJGsuKFOYIsCPcIOJETtgx1lEyilHa5\nv+8GMCnB5+eEQK2UQa2Q+gkNOUSVNJbmZ0EqIXHZnKUSggq92hMaqmvpg4SwktxURadWeuY4D1ns\nePdgF66ZXwKVPHAZbjgUuD0wbgg40ZC0hC+llBJCAk7UIITcA+AeAJgyJX7SwpyJ+DYpUUpFVw1d\nXzMZC8ryQurkREpFgRpHz7BRkXXNRsydnCsqd5EsCtQKzwjNd/Z3wWJ3YW2MwkIAoJRJkZsl54aA\nExWJ9ghOE0KKAcD99UygAymlz1BKF1FKFxkMhoQtkCNMvhqLOVvsLjhcVFSOQCGTYMak0MN2IqXS\noMYpoxlDFjsa2wZSsmzUG28p6k0N7ZhmUKOmbGIDWTTwXgJOtCTaELwJ4E7393cCeCPB5+eIoFCr\nGneFKUZnKFFUGjRwuCjeaOyEzZma/QPe6DRsfGdLzzDqT/Zj7aKymCfRhVnTHE6kxM0QEEJeArAT\nwDmEkHZCyNcAPAzgCkLIMQCXu3/mpBgGHyEzkwjl0UQhaA5t2HUKUgkJOcUr2RSolbA7KZ7/rBUS\nAtywYHLMz8H1hjjRErdPNqX0lgAPXRavc3Jig0GrhMnigMXuhEou9fIIkm8IprnnFzd1mTC/LC8l\nvJRgCE1lr9S34eKZBkwSMccgXHwNN4cTLokODXHOAnyblMTMK04UedkK5LvF5VK1f8Aboa/CYndh\n7aKyuJyjUKuE2ebEiJUPsedEBjcEnAn4jkAUM684kQgdxvGcTxwrhO7ivGw5LpsdH+lxPqmMEy3c\nEHAm4LuxiJlXnEimGdTu/EDqewQFWuYRXDe/BEpZlL0Dxz4AnrkEcIzXFRL+Xke6Jw658Yfd6cLX\nnv8cP37twDgBP07mkhqfbE5K4WsIUik0BADfuHgaLj2nMGU8lGAU52bhtzfMxfLqGPROHvsA6NwL\n9DUDhbM9dy+aqkNlgRo/e+MQaqfqPH+/QPzu/S/w0RFWuT27SIs7ziuPfm2csxruEXAmoFcrQMiY\nITClWGhomkGDlXOLk70M0dx67hSPdEdUGI+P/+omSyHFn29bCNOoHd9/ee+EWRLefNh0Gs980oJb\nz52CS88x4FdvH8aB9omDiDiZBTcEnAnIpBLo1QqvHIEdGqVM9LAUTpwwHmNfe49NeGh2cQ5+eV01\nPj1uxOMfT3wcANr7zfjhxn2oLsnBz6+uwu9vqoFeo8C3N+yByV0ZxslMuCHg+KXAq0lJrPIoJ47Y\nLcBAG/vexyMQuGlRGVYvmIz/99ExfHq8d9xjNocL396wFy4XxZO3LYRKLoVOrcATty5A58AofrRx\nP88XZDDcEHD84l2bLngEnCTS1wI2LhwBDQEhBL++YQ6mGTT43r/24oxpTCbk4XePYF/bAB5ZMw9T\n9WrP/bVTdfjRlefgvUPdeP6z1jj+ApxUhhsCjl/GGwLuESQdISxUssBvaEggWyHDk7ctxIjVie+8\ntBcOpwvvHezG3z89gbvOL/ebW/n6hZW4fHYhfrvlMBrbBuL1G3BSGG4IOH4RDAGlFMNWccqjnDgi\neAEzrwRG+wBzX8BDZ07S4lfXz8GuE334yWsHsW7TPswvzcV/r5rl93hCCB5dOx+FWhW+vX4PBs08\nX5BpcEPA8YtBo4TN6YJp1ME9glSg9zigLQaKa9jPAcJDAmtqS7G2thQv17eBAHji1oVB+xjyslm+\n4MyQBT9+7UAMF845G+Cfbo5fPJOvhi0Ysti5R5BsjMcB/XR2A1h4qGxJ0Kf88ro5oACuqylBmS70\nfIgFU/Jx97IK/HV7C/+bZxjcI+D4RdAbOjNkhcniSAnl0YzGeIwZgfypgEQW0iMAWH/Bo2vn48IZ\n4ud5XDLTABcF6lv7o1kt5yyDGwKOX4Tu1I7+UdgcLh4aSibmPmC0HyiYAUjlQH7FWPI4xiycmg+F\nVIKdLca4vD4nNeGGgOMXwRC0uAfF8/LRJCJUCQlhIf10ljOIAyq5FDVT8lDHDUFGwQ0Bxy85KhkU\nMglO9DBDwOPFScToYwgKprO+ApczLqdbWqnHwY5B3m2cQXBDwPELIQQGjRItvcMAUkd5NCMxHgck\nciBvKvtZPwNwWoHBtricbmmlDi4KfH4icIkqJ73ghoATEINWidZeMwDuESSV3mOArgKQuo2x4BmI\nSBhHwsIp+VDIJDw8lEFwQ5AqUArsfwUYTZ3OzkIt6yUAuEeQVIzNzAsQKHB/LyZPcPgtYPhMWKdT\nyaVYUJaXMglji92JzQ3tXAspjnBDkCqc/BR49evAnheSvRIP3rr2OdwjSA4uJ8sH6KeN3ac2AMrc\n0JVDg+3Ay7cDu58J+7TnTdPjUKcJg6PJzxO8vrcDP9y4D3tOpc5FUrrBDUGqsHc9+9qdOl2d3oaA\newRJYuAUywcUeHkEhDDDECo01LqDfe09GvZpl1bqQVMkT3Cok01ea+oSN4GNEz7cEKQC1mGg6Q32\nfYoaAg03BMnB2My+CnkBgYIZoUNDrdvZ1whKTWvK8qCQpUY/gWAAmjq5IYgX3BCkAk2vA/YRoOJi\ndvVmMyd7RQDGuotVcgnkUv6vkhQ8paMzxt+vnw6Y2gHbSODnnnAbgr5mwOUK67QquRS1U/KTnjB2\nuigOewwBn6QWL/inOxVo3MA+2Iu/BlAXcOZwslcEYMwj4BVDScR4HFDlAuqC8fcLHkJfi//nDZwC\nBk4ChtmAw8KMRpgsrdSjqcuUVDXSk8YRmG1OFGgUONI9BIczPIPGEQc3BMmmr4UlimtuBYrns/u6\n9yd3TW7GDAEPCyWNXrfGEPEZE+qpHAqQMG79lH1ddHfw44KwtFIHSoFdJ5LnFQhhoetrJsPqcOFE\nbxAPiBMx3BAkm8aXACIB5n2ZNQwpc1ImTyAMXOceQRIxHp8YFgIAXaX78Wb/z2vdDmTpgNnXBj8u\nCDVT8qCUSVDXkryEcVOnCTIJwXU1kwGMJY45sYUbgmTicgH7XgIqLwVyJ7OrvqK5KeMRqORS5GbJ\nufJosrCNAKaOiYliAFCogZzSwCWkrduB8mWAtghQaCMSqVPKpKidmh/zhLFx2Iq+EZuoYw91mjC9\nUINZxVooZJKkVA5lghfCDUEyaf2EyQTU3Dp2X9E84PShuOnIhEt5gRrFuapkLyMzEa7iC/wYAuF+\nfyGf/pMsR1B+ofhS0wAsrdTjSLcJA2ZxG3coBsw2XPvEp/jW+gZRxzd1mVBdkgu5VIJZRdqEVw69\nf6gblz66FfvSfIQnNwTJZO961hg06+qx+4rmAnZz4CRggnnursX42dVVyV5GZiJs3v5CQwDzFIzN\nrCvdG6F/oPxC9lVMqWkAzpumd+cJog8PuVwUP3xlHzoGRtFwsh9mmyPo8WeGLOgZsqKqJAcAUFWc\ng6YuU0I7jP+1+xQAYO+p9J7PwA1BsrAMsvb/uTcCcq8r7qK57GuKhId0agXPESQLwRAI+QBf9DMA\n6yAw0jP+/tYdQLYeMMwaO26wDbCPhr2EeaW5UMkl2NkcfXjo2e0t+OjIGSyvmgS7k6LhZPDN9XDX\nEABmAACgqiQHfSM2dJssUa9FDKdNFmw7yt7bdG9m44YgWRx6DXCMAjW3j7/fMIspTXalhiHgJJHe\nY0BuGaAIMGbSe2ylN607gKnLAIn7462fBoBG5GUqZVIsmqqLup+gvrUPj7z/BVbNLcIfbq6BVEJC\nGpdD7r4Bb48ASFxj2Wt7O+CiQEWBmhsCTpxo3AAUnANMXjj+fpkCKJyVMpVDnCRiPD5eY8iXAj8q\npP0ngcFTQMVFXseFKDUNwdJKHY50D6FfZILXl74RG+7bsBel+Vl4+MZ50ChlmFeaG9K4NHWaUJqf\nhdws5pHOKs4BIYkxBJRSbKxvQ+3UfCyvnoSj3cOwp3EPQ1IMASGklRBygBDSSAipT8YakkrvMaBt\nF7Dgton14QBLGHNDkNlQGrh0VCC3DJAqx1cECbIS5ReM3adzG5MIx1surdQDiKyfwOWi+K+XG9E3\nYsOfb13oES9cWqnH/vZBjFgD5wmaukweLwBgU/LK9Ym5Om9sG0BzzwjW1paiqjgHNqcLx88Mx/28\nySKZHsGllNIaSumiJK4hOTRuAIgUmHez/8eL5gIjZ4Ch04ldFyd1GOkBrKbxYnO+SKQsf+DdI+Cb\nHwAApQbQlkTUSwAA80rzkCWXRtRP8NS2Zmw72oOfXVOFOZNzPfefV6mHw0VRHyBPYLY5cKJ3xBMW\nEqgqzklIL8HGhnao5BJcNa8Y1SWJDUklAx4aEksklQpOx8Sbwwrs+xcw/XJW4+2Ponnsq5iEsb9z\nOB1ha8ukDemiWe+ZUxwkNASMLyGllBmC8gv8dCIHKDUVgUImwaLyfOxsNsLhdIm+7Ww24vf//gJX\nzyvG7edOGfeatVPzIZOQgOGhw11DoBSoLskdd39VSQ5O9ZnjOkbTYnfirX2dWDmnGFqVHBUFGqjk\nyelhSBTJ6hSiAP5NCKEA/kIpDV8wPRqOfwi8cR9w746JGi7+oBR4+gJg+mXAFb8Ud47//C+w7f+C\n/ap+uPJ/Az+3aA772r0fmHFF4OPe+SHw+V/9P6YtAe77nF0NZgr7NwIf/g9wy0tjch1nK6FKRwX0\nM4Av3gWcdjZ/YLANWPY9P8dNBw5uZv/L/sKRIVhaqcfv3v8C03/ybljPqyhQ439XzwXxPufr34Ka\nujC/7CsBDYGw6U7wCNw/H+kawpIKXVhrEcv7h7oxZHFgbW0pAEAqIZhVlONJXqcjyTIEF1BKOwgh\nhQA+IIQcoZR+4n0AIeQeAPcAwJQpU/y9RuTUPQ0MdQEtW4G5a0If33MEOH2Q3SbXAlXXBT/+i3eB\nbQ8DM1ey431Rasf3DviiymVyE8HyBJZBYO+LLClYftH4x0b7gLongcNvjm9WS2fOHAHe+i7rwXjl\nTuAb29j7eLZiPMbi/7llwY/TTwdcDtZAdvIzdp/QPzDuuBnsf8ZsFHfx48Nt506BhJCwRN8kEoLr\nakrGlx87bO6KOQsuW3g7fr/ThGGrAxrl+K2oqdOE3Cw5SnyaGavdOYNDnYNxMwSbGtoxOS/LkxsB\nmAF6e18nKKXjjVqakBRDQCntcH89Qwh5DcASAJ/4HPMMgGcAYNGiRbHz902dQPNH7PvWHeIMgdCg\no5/BPImiuYFruwdOAa/dy8I7a58f3yMQDkVzgxuCg68yVcnLH5pobCgFjr7PGtYywRDYRoCNdzLZ\nheufAjZ9FXjzO8DaFyK6+k0Jet0VQ5IQ0VvviqDWHWx6meGc4MdFYAjyshX45iUhwlRi6NzLjDWA\nlc5teMRVg/rWPlxyTuG4w1hHcc6ETdegVaJAo4hbvL5zYBQ7jvfiO1+aAYlk7NxVxTnYsOsUOgZG\nUZofoJz3LCbhOQJCiJoQohW+B7AcwMGELWD/y0zqubBqbIMPRet2dmV2+2a2sWy8C7D7aWpx2Nhj\n1AXc9ELkRgBgoQ1jMxta44/GDUxiuGThxMcIYQbg5A6g70TkazgboJSFyHq+AG78K1B9PXDZz9mg\nn93PJnt1kWM87l9jyBfPIPtjbn0hP/kBYCzXEKeB96IRqpoMszG17XXIpZiQhHY4XTjiUzEkQAjB\nbHeHcTx4bW8HKAXWLCwdd3+6J4yTkSyeBGAHIWQfgN0A3qGUvpeQM1PKrpKnnA/Mv4V9eExdwZ/j\nco0l4PKnAtc/DXTtA/79k4nHfvg/QEcDcN0TgT0GsRTNBUCZ7pAvPUeB9t1ssw90xTv/ywAIE7VL\nZ/a+yH7HSx4EKi9h953/XWDmlcD7PwY69iRzdZHhtAP9J8QZgmwdUxk99gETqPMuG/UmbyprVIyw\nhDRmtG4HCquBpd+ExHgUayadniBqd6J3BFaHa0J+QKC6JBfHTg/D5ohtQYTQO3BuhQ5T9OOv+mcV\n5UBC0lf9NOGGgFLaQimd775VU0p/k7CTt9ezD0LNrUCFO4568tPgz+k5wuKqwgds1irgvPtYkvbg\n5rHjDr/N4vJLvhE6hyCGYFITjeuDl58CQG4pMO1SJnOdrhVEpw8BW+5nBuCidWP3SyQsRKQtYh7a\n6FkmGDZwisX9g5WOelMwAzixjX3vLz8AjJWaRqg5FBMcNuDULvbZq74BkGXhZvl2HOwYxJBXFVCg\nRLFAVUl86vrrT/aj1WjGmtrSCY9lKaRp3WGcWeWjjS8C8mwWPiiax7T/T3wS/DmeBh2vD9jlDwGl\nS4A3v8fCN/2twOvfYmGa5b+KzVpzJrMrPd88gcvJwlszrgC0k4K/Rs1trMtU+B3SCesQSwqrcoHV\nz7KNzptsHbDmOXaV/Ma3z66yUk/pqEhDIHgO6kKgYGbg4wpmJDc01LmHyaqUXwCocoCqazGn7wPI\nXNZx/QRNnSYoZBJMM/ivePNITcR4U95U345shRSr5hb7fby6JJeHhs56bGaWYK26jlXtSKTA1PND\n5wlatwO5U1hYSEAqB9Y+B0hlbDPaeBdAwO6TKQO9UngEmk3Q/DGreKq5LfRrzLqKqZs2bojNmlIF\nSoG3/4vN4l3zd0BT6AKSKvsAABqQSURBVP+4ssWs3PfI20DdU4ldYzQYRfYQCAiGIFB+wPu4vhbW\nZ5IMTmwHQJgOEgDU3AqZfQgrZQ2o89IdOtRpwjmTtAHnZFcUqJEll8Z0UzbbHHh7fydWzS2GWum/\nhqaqJAcdA6NJHd0ZLzJn4siRd1inpvcGWn4hcPQ9VkmUUzLxOS4XG/k388qJj+WWAjc8A2xYy36+\neT2QXx7bNRfNZQlPp4MZHYCFhbJ0/tfkizwLmLOaNbCt+h27CosV+18Bdgdo/9AWA9c+DmTlxe58\n3jQ8DxzYCHzpp4Fj4gJLv8XKKj/4GXDo1fDOM/0K4JIHxB175gjrG7nil0BeiJLPUHQfZN3B2SLL\nI4UQUqj3Qj8dcNmZlxgoh3Vwc/hGU6oErv1TaMPVuh2YNGfs9yq/CMgtw12jn+J/Wtj/M6UUTV0m\nXDE7sLcrlRDMKtaiqct/Xb/TRfGbdw7jghl6fGlWCK/ZzXsHuzFic3p6B/wheCKHugZx/rTwK69S\nmczxCBpfBPKmjF2NAGMfnEBewZkmVpNfESDuOnM5cPVjwMrfAbOD9AVEStE8wGkdu0I09zGDNu8m\nJk4nhgW3M3e86fXYrcvlAj7+FTDYwbwr75tCA3yxJX7hmK79wLsPANMuAy74YejjCWHJ+zlrJq41\n2M3lBLb+Vpw3ZR0GXrmDGZpNX2XJ3khp2cpCf8H6THypuAhYcAeLuwfDU0IaIDxEKfDxr1ljWjjv\nVftuYNdfgp/bYQXado83VhIJMP8WzLftRU9HC4Ysdpw2sellgfIDAlXFOWjq9D+b4E8fHcPfPz2B\n37xzWPTsgo317Ziiyw7am1CVxpVDmeERDLQBLdtYZYl3XXbRXBZjbt3ONldfBAPhbTx8WfTV2K7V\nm2JBauIAUDibXa05beH1BkyuZXHjveuBhV+JzbpO7mAJzdXP+n/fPnsc+PdPgV1PA0u/GZtzAoDF\nxPoFsvXA6mdC19gLZOUDq0NsVL64nMA/rgPe/gFQXANMCjCcRwhTGY+zIoKdTwAfPgSsiKAGYqgb\n2Px/WB9AsM5zX1S5zNiFQsg5GI+BVW37cKqOhY6ue5IJIopl413AgVdYfixQaLSjgV2Q+F5U1dwK\n8skjuF6yHZ+3rvDcXR3KEJTkYP2uU2jvH0WZbqzCZ8exXvzp42Mo02WhuWcEe9sGsHBKftDXausz\nY2eLET+4YmbQZrECjRKTcpRpmTDODI9g378AUHdJpRcSKdvkA3kErduZF+GdH0gk+hnM7e7ax35u\nXA9MmhuefAIhLBzWVhe7ipHGDSzRHuiq9bz7gHNWAf/+GdAubiRhSChlTWL9J1leIIKmqLCQSFlf\nglLLjE+gfo49L7BN8JL/Zpv/4q8zY3BkS3jnczqYEbCNsEY4hTr638GXbB2gygucMG5cD8jV4Ve9\n1dwOjPazMGsgWncAIMCU88bfr6uAa8oy3CT7BHXNRhzqYJvsLD89BN74SxifNlnw/Zf3YrpBg83f\nPB8quQSbGtpDLn/znnYQAtwYJCzkfd509AjS3xBQyv7Byy/0H8Mvv4BdBQ12jL/f5WKlpb7yDYlE\nKmNXot0HgNNNrCszkk7heTcDRALsi0HS2DrEmrWqbwg8MIUQ4PonWa5g410spBUtn/+Vhbcu+xkw\n9bzQx8cCbRFw47OsiuedH0wMdXUfALb8CKi8FLjQHaZa8RvmQbx+LzNaYtn2MLvwuOoPbB5FPCDE\nPbbSTy+BbYRJP1RfH74+1bRL2d967/rAx7RuZxpafvIekgW3oZx0Y/DoDjR1mVCuz54gOeGLb12/\nw+nCd17aixGrE0/ethCFWhVWzinGW/s6YbEHnv/tclFs3tOO86fpMTkvK+SvWlWSg+NnhmF1pMZM\n8ViR/obg1E7WnBOoykYoC/X1Cs4cYlc5oRJw8UaYTdC4HpDI/IdiQpFTzNROG19iIY9oOPQ6kwhY\ncHvw47LymcTGUFf0+YLOvaw5bMZy4Hw/gmrxpPISFlLc/zKw5x9j91tMrGIsWze+fFWmZL83BbDp\nblY7H4rjHwKfPMre05pbYv87eKMPUEJ6+C3ANiyuGs0XiZR528c/YOEtXzz5gQAXVVXXwSbJwsK+\nd1B/sn+C4qg/shRSVBo0nqvzxz48it0n+vCbG+ZgxiQtAGBtbSmGLA68f8jPmtzsOtGHtr5RrK0V\nl+CvKs6Fw0Vx7HR6zSZIf0PQuJ4lMKuu9f/4pDnMXfattfcMAE+2IZjLEtYNL7BKoUhDIjW3AkOd\nLBkZDY3r2WZSujj0saW1LG78xRYWLomE0QG24aoLgRv+Ij4vEEsuWscMwrs/YhU9lAJvfY9dYNz4\nN0BjGH+8rgK4/s8sLv7Bz4O/9mAH8Oo9LAe08nfx+g3G0E9jxtk31LX3ReYxTz0/stetuY1Jq+x/\neeJj7fVMFyvQZ0mpwUDFKqyS7MLw0GDIRLFAVXEODneZsPWLM/jzf5px86IyrPaShlhaya7yg4WH\nNja0QauUYUV1AEl4H9JVaiK9DYFthF3BVl8fOOYqkbjzBD6G4MR29sGIthQwWoTZBLahyK7WBM5Z\nxQxeYxD3PRTGZuZhBZO28OXce4HZ17AEatvu8M5HKfMmTB2sR0NsOWWskUjZVb8qF3jlK8Bnf2IV\nQl/6KVAeoJBg9jXAud8Edj0FNL3p/xinA9j8NaZbtfaFwKG2WCJUDvmOt2zdzv6/IhXpK5jBmiz3\nrp/o/Qn5gSAhvdzz74aWjOJKyed+NYb8IdT1f//lRswq0uIX11WPe1wiIbixthQ7jveic2B0wvOH\nrQ68e6AbV88vRpZCOuFxf0zRZUOtkKZdwji9q4aa3nS7uyHCGBUXAl+8w6qL8srG8gPxKAkNl0nV\nAAjzBILNJgiFTAnMXcvCG6P9LHQTLo0bWK7BN+keDEKAa58A/nIRyxfMWS3+uaYu1gy2/NdA2ZKw\nlxtTNIUsSf3CNewqf/oVwLL/Cv6cK37JSivf+DYbTeq7yQqGdfVfAUOQjuBY4hGpOw6U1LDv970E\ngDD9rWhYcBvzlDr2MG9QoHU7q4AL8j+nrLwA3dJirJVuw7SSh0SdTjAYdocLf75tIVTyiZv5moWl\n+NNHx/Dqnnbc96Xxndpb9ndh1O70KykRCImEid6l22yC9DYE3fvZvNYpS4MfJ7isJz8F8r7M5g5Y\nBpKbKBZQaliHcNkS1tEcDQtuAz5/lnVYL/5aeM91OdmGMe1L/pvvgpGVx9RY/3U78PnfwntuzW2s\nCikVKL8AWPFbFv4QE6aSKZjMxYurgfq/+z9m2feAeWtjv9ZA6CoBkDGPwOViBr7ioui93+obgHcf\nZD07giGwW5gnuOTrwZ9LCHqn3Yjzjz4B6uwGUB7ydPPL8jDNoMYPrjgnoBzFFH02zq3QYVNDO759\n6fRx5aGbGtpRWaAOWV7qS1VJDl7d0wGXi46Tqj6bSW9DcOX/Apf+JLS7W1jNrlZObGdXux59oSD9\nA4nky1GEc7wprmHy240bwjcEJ7axEM3yX0d27pIFwA/8KKmebSz9Zni9EflTge/EqIQ2Fsiz2IYv\nVA6d/BQYOMnCXNGiymUhsQObgRX/y2TYO+pZU2QgMTwv5qz6BnD0zyD7XhbV0Z2bJcdHP7wk5HFr\nakuxbtN+1J/sx+JyFl5s7R3B7tY+rFtxTtiDZqpLcvCPnSfR1m/GVH0cynyTQHrnCABxpXC+eYLW\nHUB+BZORSCeEnoKOeqbfHw6NG9gH/ZxV8VkbJ3Hop495BKF6QsJlwW2AdZCF9AB2cUUkob1ygPXs\nVFzE8lgxVMxdNbcY2QopNtWPJY0372mHhAA3Lgz/M15VzKqa0kmSOv0NgVjKL2RXRn0nmL5QIFmJ\ns515NzEJ63CSxqMDrLxw7trohu1wUgOhhNQ6xHozgvWEhItbP8gjzdG6gxU8iNWdqrmNfQ5PfRab\n9QBQK2W4am4x3t7fCbPNAaeLYnNDOy6cYUBRbvj/zzMmaSCVkLSqHOKGQEDY+Hc9za5oRLiyZyWa\nQmDmCtZtLVaF0j1jNiPGXmYCBTNYEcWup8X1hISDWz8IzR+zZHj75+FdVM2+BlBogzenRcCa2lKM\n2Jx472A3djYb0TloCStJ7I1KLsV0gyatKoe4IRAwzGaqnvXPsZ+T3T8QT2puBYZPsw+rGBrXBx6L\nyTn7EFRCP31cfE9IONTcAsA9QlRkfsCDIhuYcwPrXg8k6xEBSyp0mKLLxsb6dmxsaEOOSoYrqsQp\nk/qjuiS9pCa4IRCQSNjm77SySqNwK2POJmasYMJtjS+GPrbnKLuqC6d3gJPaCOJz1sH4/F11lSzn\n1vIfd34gTEmQmtsB+0hMFXMJIVhTW4qdLUa8e6Ab19aU+C03FUtVSQ66TRYYh60xW2My4YbAG+HK\nJZ29AYCVNc67Gfji3dA6QGLGYnLOLnImA7Ks8HtCwkFofiyuCX8ORtkSltAWIwHe84XoxPKNtaUg\nBLA5XaIlJeByMYkXnya5SKakDfV140zrQdHHJxJuCLyZ9iW26c3wI9GbbtTcyiStD2wKfEzfCRYq\nm7ki9FhMztmDRMIUbGeujJ/nW3UdK8mOpAmSEPb/efJTJggZiP2vAH9ewnIdIpicl4VLZhpQVZyD\neaWh9YwAsJkUT1/Ausm9qCphone7WkQKKtpGYHpqOfTPXYCBN/+b9VekEETs4IZksmjRIlpfX5+Y\nkw2fAdSGzAiDPH0huyr8xraJjzmswN+WMz2db2xPnhQ3Jz5YTEzEMJ6yFuY+JuMdSSOkqRN4rBq4\n8H7gSz+Z+HjPUeCZS1gIqbAK+OZnoj6zQtWQViViTcc/Al68kZVNW4eAu7eMK4P9yt93o/nMMLb/\n6NLgjWWUYvSVe6Bs2ogPXLVYIa2HSz8DkuufYuNU4wghpIFSuijUcdwj8EVTmBlGAGDue1cjcNpP\no9f7P2GPXf8UNwLpiCon/tpG2brIu+FzSpi8976XJoZ+bGam+STPAi5+gE0S7NwrbkkKmTgjYOoE\nXv06EwP89m7W47DxbmBkbLbymtpSdAyMYmeLMcgLAdj7IrIOv4L/51iNriv/ijtsD2LQZAL+vpwN\ncLJP1EFKNNwQZDJz1wIS+cRY7MFXmRTFefcxeQsOJxksuA0YbANaPxl//5b7gZ4jbFbE0m8BMpW4\nfIJYnA5gk5cYoHYSkxY39wKv3eMxTMurJkGrkgUfftN9EHTL/aiXzkdd2Vdx17L/396dB1dV3QEc\n//6yCYQlJrLEgBp2BCQsRZBN2QShKFsVtQMuZZgCwrTWQqcytFWrTDsUtGqtLP0DSk1kqyiCSK04\nVTAggoDIniCb7ItJCPz6x7mRKNl5Lw/u/X1mMnn3vXvzzm8ezO+dc885v1RadhtMt9PPsvvm4a6a\n36tdXXW4CLJEEGTxSdCsn9s7p6DO7tGdsPQJN6Ww99RIts4EXbMBblim8JqCDfPcBIbuv3L39Kom\nuFXRm9JDN+6++hm3oO3HMy5tBnhjmtuyZsd78NF0wK0nGNTmRt7ZfIBTOUXUqc49DekjyY+twZiz\nYxjawfWsn+zblBa3pDBg1xD2D1rg7tXN7gfLJ7veTgRYIgi6tIfh7BH4aqX7j5Q+0lVGGzbnyje5\nM+ZKxFaBVsPcqvack65K37Jfutl9d066dF7bh9wmkV+WszxoUbavgDXTof2oyzcD7PAYtBwC7z/j\ndh8AhndoQM75iyz7/MD3z1WFf0+EY7uYk/w05+ISGdA6GYCY6CheHOF2S330g3i+fXwN/Ohx+Phl\neLUL7A3dquqyskQQdI17u6Ivn82D5ZPcVLnBf4t8HQZjwN3Hyv/W9QTSR7qbz0NnXaoIB5Daw02J\nvdLhoZPZbuinbmvo9/zlr4u4XsL1qZDxKJw5Qpv6tWhcp/rlw0OZc2BzBue7T2bmznr0b5VMfKHy\nm/VqVWH6/WlsP3yaKe/sgQF/gpFvuV1+59wD7/za1VOpJJYIgi46BtrcD9uWuX+8XSa46aLGXA1S\n2kHt5q5U6dEdMPT1y6cyR0V721qscjd5K+LCeXcz+EK+2zI9tpj6xVVqutdzTsDCnyHfHufh22qw\nY28Wu7Oy3EyprLVuO+5GvVhWawRncvMZ3uHy7Sx6NK3N2Dsbk56Z7RJJajf4+f+g42g3JfaVOy4v\noRsmlgiMt/hHoUEn6Pl0pFtjzCUFawpQ6DEJGvYo+ry0B4svlVkW7011RYQGzby0BUdx6rWG/i+4\nldPTUhm15i42VhlN6qxWMC0VZvVxK/eHvEb6+v00SKxKx1uKrq43sXcTbk9N5LeLN7H90GlXSfGe\naTDqbUBg7oDiK9yFkK0jMM7uD6Feq4pVLjMmnC6cd7W2G/X8/pDQD83uB2e/gXHryjcFfNvbsGCE\nG6cf8OeyXaPqetEn3ZDQ/LX7OHgyh4m9mxAlAk37ki316DZtNRN7NWVC7ybF/qnDp3K4Z+aHJFSL\nY8nYLpeGkPLOuVrfd4wvvodSCltHYMontZslAXN1io51K5RLSgLgegVHv4LscnxpPL4XFo9xK63v\nfq7s14m4UradxkCnMSTcNZ6ZZ3vx38Sh7rnEhixcvx9VGNIupcQ/VadmFWY80JadR87w9OLNfPfl\nPK4a9HiqwkmgPCwRGGP8oeVgiK1Wts0UAfLzIOMR9+1++FxX17uCerWoQ0K12O9uGqsqGZnZ3NEo\niQaJpS/c69L4Bib0asLCDft549OsCrejoiwRGGP84boa0GKQWxBZlvn4K6fA/ky49yWvlvMVvHVM\nNPelpbBiyyFOnjvP2t3H2HfsXLlqHozv2YSujW9gypIv2FrJtQ4ikghEpJ+IfCkiO0RkUulXGGNM\nGbR9CHJPufH7kmxZCp+8ArePcRvkhcCw9vXJy7/I0s+/JiMzm+rXxdC/VXKZr4+OEqbfn0bNqrGM\nnbeeM7llLBwVApWeCEQkGvgr0B+4FRghIrdWdjuMMT50c1e3L1BJw0PHdsOScZDSHvr8IWRv3fLG\nmjSvV4N5H+9l2aYDDLwtmapx5at5ULvGdbw4oi17jp7lNws3UVmTeSLRI+gI7FDVXaqaBywAQpOS\njTHBFhUFbR6EXR/AiSLG2vNzIX0UCG71fExcyN66oPjNtoOnOZd3ocKlMDs1TOIXfZqydOPXzF+7\nL2TtK0lM6aeEXApQ+BPKBm4v6YKjR48yd+7ccLbJGOMT1c/HMQzl7ItdyIv6/oyb2Iu5VL9wnFV1\nHidryeqQv3dOvhBFEgmxF9i0eimbK7iRcTWFRvG1mLJ4E/s2fEhylfAOE0UiEZSJiIwGRgOkpJQ8\n/coYYwqciU0iM2EgSXlF7wq6sWpfsqq1Dst7V49R+tc9Q0LsxSvazT5KYEjyKd78uibREv7hoUpf\nUCYinYGpqnq3dzwZQFX/WNw1tqDMGGPK72peULYOaCIiqSISBzwAhH8NtTHGmCJV+tCQquaLyDjg\nXSAamK2qRZTIMsYYUxkico9AVd8GQrB5uDHGmCtlK4uNMSbgLBEYY0zAWSIwxpiAs0RgjDEBZ4nA\nGGMC7pqoUCYiR4C9kW5HBdwAfBPpRkSAxR08QY39ao/7ZlWtXdpJ10QiuFaJyKdlWdXnNxZ38AQ1\ndr/EbUNDxhgTcJYIjDEm4CwRhNdrkW5AhFjcwRPU2H0Rt90jMMaYgLMegTHGBJwlghARkSoislZE\nNorIFyLyO+/5VBH5RER2iMi/vK23fUdEokVkg4i85R37Pm4R2SMim0TkMxH51HsuUURWishX3u/r\nI93OUBORBBHJEJFtIrJVRDr7PW4RaeZ9zgU/p0Rkol/itkQQOrlAT1VtA6QB/USkE/ACMF1VGwPH\ngcci2MZwmgBsLXQclLjvUtW0QlMIJwGrVLUJsMo79psZwHJVbQ60wX3uvo5bVb/0Puc0oD1wDliE\nT+K2RBAi6pzxDmO9HwV6Ahne8/8A7otA88JKROoDA4DXvWMhAHEX415cvODDuEWkFtAdmAWgqnmq\negKfx/0DvYCdqroXn8RtiSCEvOGRz4DDwEpgJ3BCVQsqT2cDfizA/BfgKeCid5xEMOJWYIWIZHo1\ntgHqquoB7/FBoG5kmhY2qcARYI43FPi6iMTj/7gLewD4p/fYF3FbIgghVb3gdR3rAx2B5hFuUtiJ\nyEDgsKpmRrotEdBVVdsB/YGxItK98IvqpuT5bVpeDNAOeEVV2wJn+cFwiE/jBsC71zUISP/ha9dy\n3JYIwsDrKq8GOgMJIlJQCa4+sD9iDQuPLsAgEdkDLMANCc3A/3Gjqvu934dx48UdgUMikgzg/T4c\nuRaGRTaQraqfeMcZuMTg97gL9AfWq+oh79gXcVsiCBERqS0iCd7jqkAf3E201cAw77SRwJLItDA8\nVHWyqtZX1VtwXeb3VfUhfB63iMSLSI2Cx0BfYDOwFBcv+DBuVT0IZIlIM++pXsAWfB53ISO4NCwE\nPonbFpSFiIjchrtZFI1LsG+o6u9FpCHum3IisAF4WFVzI9fS8BGRO4EnVXWg3+P24lvkHcYA81X1\nWRFJAt4AbsLtmPsTVT0WoWaGhYik4SYGxAG7gEfw/s3j77jjgX1AQ1U96T3ni8/bEoExxgScDQ0Z\nY0zAWSIwxpiAs0RgjDEBZ4nAGGMCzhKBMcYEnCUCY4wJOEsExhgTcJYIjCmFiCz2Npb7omBzORF5\nTES2ezUo/i4iL3nP1xaRN0VknffTJbKtN6Z0tqDMmFKISKKqHvO2DlkH3A18hNtj5zTwPrBRVceJ\nyHzgZVVdIyI3Ae+qaouINd6YMogp/RRjAu8JERnsPW4A/BT4oGArARFJB5p6r/cGbnUlGQCoKSLV\nC9WqMOaqY4nAmBJ4+yf1Bjqr6jkR+Q+wDSjuW34U0ElVcyqnhcZcObtHYEzJagHHvSTQHOgExAM9\nROR6b6vtoYXOXwGMLzjwNmgz5qpmicCYki0HYkRkK/A88DGutsJzwFrcvYI9wEnv/CeADiLyuYhs\nAcZUeouNKSe7WWxMBRSM+3s9gkXAbFVdVNp1xlyNrEdgTMVM9epTbwZ2A4sj3B5jKsx6BMYYE3DW\nIzDGmICzRGCMMQFnicAYYwLOEoExxgScJQJjjAk4SwTGGBNw/wfPip7CwK91/wAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\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')\n",
    "ax.set_ylabel('count');\n",
    "ax.legend()"
   ]
  },
  {
   "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>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\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": [
    "## 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.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\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": [
    "# 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": [
    "# 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.\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ols = sm.OLS(A, X)\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAE9hJREFUeJzt3X+s3XV9x/HnW0AkXNfCgLtaml0W\nOwNSRblBEvfHveBGAWNxkQ1CXFG2+gcajE1c0WRqHBHjkM1MMVUIdTKvTCA0FKZYuSMkA6SIFOgY\nVTtpaUoYULnK2Irv/XG+1dP23HtO7/n96fORnNzz/Xw/55zX/fb01W+/93u+NzITSVK5XtPvAJKk\n7rLoJalwFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqnEUvSYU7vN8BAI477rgcGxvrd4yGfvnL\nX3L00Uf3O8a8DGv2Yc0NZu+XQzX7pk2bnsvM45vNG4iiHxsb46GHHup3jIamp6eZmJjod4x5Gdbs\nw5obzN4vh2r2iPivVuZ56EaSCmfRS1LhLHpJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJUuKZFHxGv\ni4gHI+LHEfF4RHymGj8pIh6IiKci4tsR8dpq/MhqeWu1fqy734IkaS6tfDL2FeCszJyJiCOA+yLi\nLuBjwLWZORURXwUuA66rvr6QmW+MiIuAzwN/3qX8GnJjazYcMLZ62R4u3W9829Xn9yqSVJyme/RZ\nM1MtHlHdEjgL+E41vg64oLq/olqmWn92RETHEkuSDkpkZvNJEYcBm4A3Al8GvgDcn5lvrNYvAe7K\nzFMj4jFgeWZur9b9BHhHZj6333OuAlYBjI6Onj41NdW576qDZmZmGBkZ6XeMeRmG7Jt37D5gbPQo\n2PXyvmPLFi/oUaL2DMM2n43Z+6Od7JOTk5syc7zZvJYuapaZrwKnRcRC4Dbg5EbTqq+N9t4P+Nck\nM9cCawHGx8dzUC9IdKheLKlX9j9EA7VDN9ds3vetue2SiR4las8wbPPZmL0/epH9oM66ycwXgWng\nTGBhROz923gi8Ex1fzuwBKBavwB4vhNhJUkHr5Wzbo6v9uSJiKOAdwFbgHuA91XTVgK3V/fXV8tU\n63+QrRwfkiR1RSuHbhYB66rj9K8Bbs7MOyLiCWAqIv4W+BFwfTX/euCfImIrtT35i7qQW5LUoqZF\nn5mPAm9rMP5T4IwG4/8DXNiRdJKktvnJWEkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4i16SCmfR\nS1LhLHpJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSqcRS9JhbPoJalwFr0k\nFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqXNOij4glEXFPRGyJiMcj4opq/NMRsSMiHqlu59U95sqI\n2BoRT0bEOd38BiRJczu8hTl7gNWZ+XBEvB7YFBF3V+uuzcy/q58cEacAFwFvBt4AfD8i/jAzX+1k\ncElSa5ru0Wfmzsx8uLr/ErAFWDzHQ1YAU5n5Smb+DNgKnNGJsJKkgxeZ2frkiDHgXuBU4GPApcAv\ngIeo7fW/EBH/CNyfmd+sHnM9cFdmfme/51oFrAIYHR09fWpqqt3vpStmZmYYGRnpd4x5GYbsm3fs\nPmBs9CjY9fK+Y8sWL+hRovYMwzafjdn7o53sk5OTmzJzvNm8Vg7dABARI8AtwEcz8xcRcR3wWSCr\nr9cAHwSiwcMP+NckM9cCawHGx8dzYmKi1Sg9NT09zaBma2YYsl+6ZsMBY6uX7eGazfu+NbddMtGj\nRO0Zhm0+G7P3Ry+yt3TWTUQcQa3kb8rMWwEyc1dmvpqZvwa+xm8Pz2wHltQ9/ETgmc5FliQdjFbO\nugngemBLZn6xbnxR3bT3Ao9V99cDF0XEkRFxErAUeLBzkSVJB6OVQzfvBN4PbI6IR6qxTwAXR8Rp\n1A7LbAM+BJCZj0fEzcAT1M7YudwzbiSpf5oWfWbeR+Pj7nfO8ZirgKvayCVJ6hA/GStJhbPoJalw\nFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4i16SCmfR\nS1LhLHpJKpxFL0mFs+glqXAWvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSqcRS9JhWta9BGxJCLu\niYgtEfF4RFxRjR8bEXdHxFPV12Oq8YiIL0XE1oh4NCLe3u1vQpI0u1b26PcAqzPzZOBM4PKIOAVY\nA2zMzKXAxmoZ4FxgaXVbBVzX8dSSpJY1LfrM3JmZD1f3XwK2AIuBFcC6ato64ILq/grgG1lzP7Aw\nIhZ1PLkkqSUHdYw+IsaAtwEPAKOZuRNq/xgAJ1TTFgNP1z1sezUmSeqDyMzWJkaMAP8GXJWZt0bE\ni5m5sG79C5l5TERsAD6XmfdV4xuBj2fmpv2ebxW1QzuMjo6ePjU11ZnvqMNmZmYYGRnpd4x5GYbs\nm3fsPmBs9CjY9fK+Y8sWL+hRovYMwzafjdn7o53sk5OTmzJzvNm8w1t5sog4ArgFuCkzb62Gd0XE\noszcWR2aebYa3w4sqXv4icAz+z9nZq4F1gKMj4/nxMREK1F6bnp6mkHN1swwZL90zYYDxlYv28M1\nm/d9a267ZKJHidozDNt8Nmbvj15kb+WsmwCuB7Zk5hfrVq0HVlb3VwK3143/RXX2zZnA7r2HeCRJ\nvdfKHv07gfcDmyPikWrsE8DVwM0RcRnwc+DCat2dwHnAVuBXwAc6mliSdFCaFn11rD1mWX12g/kJ\nXN5mLklSh/jJWEkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4i16SCmfRS1LhLHpJKpxFL0mFs+gl\nqXAWvSQVzqKXpMJZ9JJUOItekgpn0UtS4Sx6SSqcRS9JhbPoJalwFr0kFc6il6TCWfSSVDiLXpIK\nZ9FLUuEsekkqXNOij4gbIuLZiHisbuzTEbEjIh6pbufVrbsyIrZGxJMRcU63gkuSWtPKHv2NwPIG\n49dm5mnV7U6AiDgFuAh4c/WYr0TEYZ0KK0k6eE2LPjPvBZ5v8flWAFOZ+Upm/gzYCpzRRj5JUpva\nOUb/4Yh4tDq0c0w1thh4um7O9mpMktQnkZnNJ0WMAXdk5qnV8ijwHJDAZ4FFmfnBiPgy8O+Z+c1q\n3vXAnZl5S4PnXAWsAhgdHT19amqqI99Qp83MzDAyMtLvGPMyDNk379h9wNjoUbDr5X3Hli1e0KNE\n7RmGbT4bs/dHO9knJyc3ZeZ4s3mHz+fJM3PX3vsR8TXgjmpxO7CkbuqJwDOzPMdaYC3A+Ph4TkxM\nzCdK101PTzOo2ZoZhuyXrtlwwNjqZXu4ZvO+b81tl0z0KFF7hmGbz8bs/dGL7PM6dBMRi+oW3wvs\nPSNnPXBRRBwZEScBS4EH24soSWpH0z36iPgWMAEcFxHbgU8BExFxGrVDN9uADwFk5uMRcTPwBLAH\nuDwzX+1OdElSK5oWfWZe3GD4+jnmXwVc1U4oSVLn+MlYSSqcRS9JhbPoJalwFr0kFc6il6TCWfSS\nVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcPO6TLH6a6zBpX0buXH50V1OImkYuEcvSYWz6CWpcBa9\nJBXOopekwvnD2Abqf9i5etmehr/XFGDb1ef3KpIkzZt79JJUOItekgpn0UtS4Sx6SSqcRS9JhbPo\nJalwFr0kFc6il6TCNS36iLghIp6NiMfqxo6NiLsj4qnq6zHVeETElyJia0Q8GhFv72Z4SVJzrezR\n3wgs329sDbAxM5cCG6tlgHOBpdVtFXBdZ2JKkuaradFn5r3A8/sNrwDWVffXARfUjX8ja+4HFkbE\nok6FlSQdvPle62Y0M3cCZObOiDihGl8MPF03b3s1tnP+ETVIWv2lJ14HSBockZnNJ0WMAXdk5qnV\n8ouZubBu/QuZeUxEbAA+l5n3VeMbgY9n5qYGz7mK2uEdRkdHT5+amurAt9MZm3fs/s390aNg18uN\n5y1bvKBHifZVn28uJy04jJGRkb68dqvbptHzNdrm/drWB2tmZqbj27xXzN4f7WSfnJzclJnjzebN\nd49+V0QsqvbmFwHPVuPbgSV1804Enmn0BJm5FlgLMD4+nhMTE/OM0nmX7nf1yms2N95M2y6Z6FGi\nfc12Nc393bj8aDq9XVt97Va3TaPna7TN+7WtD9b09HTHt3mvmL0/epF9vkW/HlgJXF19vb1u/MMR\nMQW8A9i99xCPVLqxNRvmvKz1Xh7WUq81LfqI+BYwARwXEduBT1Er+Jsj4jLg58CF1fQ7gfOArcCv\ngA90IbMk6SA0LfrMvHiWVWc3mJvA5e2GkiR1jp+MlaTCWfSSVDiLXpIK5y8HV1H8QJd0IPfoJalw\nFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4L2qmrmj1\n4mKSus89ekkqnEUvSYWz6CWpcBa9JBXOopekwln0klQ4T6/UUPB0TWn+3KOXpMK1tUcfEduAl4BX\ngT2ZOR4RxwLfBsaAbcCfZeYL7cWUJM1XJ/boJzPztMwcr5bXABszcymwsVqWJPVJNw7drADWVffX\nARd04TUkSS1qt+gT+F5EbIqIVdXYaGbuBKi+ntDma0iS2hCZOf8HR7whM5+JiBOAu4GPAOszc2Hd\nnBcy85gGj10FrAIYHR09fWpqat45Om3zjt2/uT96FOx6ufG8ZYsX9CjRvurzzeWkBYcxMjLSl9du\nx1zbvFO68We3ecfulrL3633TzMzMTMffL71yqGafnJzcVHfYfFZtFf0+TxTxaWAG+CtgIjN3RsQi\nYDoz3zTXY8fHx/Ohhx7qSI5OqD+Vb/WyPVyzufHPrLddfX6vIu2j1VMNb1x+NBMTE3157XbMtc07\npRt/dmNrNrSUvV/vm2amp6c7/n7plUM1e0S0VPTzPnQTEUdHxOv33gf+BHgMWA+srKatBG6f72tI\nktrXzm7TKHBbROx9nn/OzH+NiB8CN0fEZcDPgQvbjylJmq95F31m/hR4a4Px/wbObieUJKlz/GSs\nJBXOopekwln0klQ4i16SCmfRS1LhvB59D7T6IaNB/SCNytDsfbh62R4uXbPB92GB3KOXpMJZ9JJU\nOItekgrnMXoB/k5WqWTu0UtS4Sx6SSqcRS9JhbPoJalwFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuH8\nZKw0oLzqqTrFoi/Y5h27udRLG0iHPIteasLrAGnYWfQDxELpHbe1DiVDX/QH8xfWY5kaBP4jo17z\nrBtJKlzX9ugjYjnwD8BhwNcz8+puvVa/uGcm9cf+f/f2/r7b/fXzf/Gt9sONy4/ucpIuFX1EHAZ8\nGfhjYDvww4hYn5lPdOP1pEPZMOxweKpof3Xr0M0ZwNbM/Glm/i8wBazo0mtJkubQrUM3i4Gn65a3\nA+/o0mtJUk8Mw/+eGonM7PyTRlwInJOZf1ktvx84IzM/UjdnFbCqWnwT8GTHg3TGccBz/Q4xT8Oa\nfVhzg9n75VDN/vuZeXyzSd3ao98OLKlbPhF4pn5CZq4F1nbp9TsmIh7KzPF+55iPYc0+rLnB7P1i\n9rl16xj9D4GlEXFSRLwWuAhY36XXkiTNoSt79Jm5JyI+DHyX2umVN2Tm4914LUnS3Lp2Hn1m3gnc\n2a3n76GBP7w0h2HNPqy5wez9YvY5dOWHsZKkweElECSpcBZ9AxHxhYj4j4h4NCJui4iFdeuujIit\nEfFkRJzTz5yNRMSFEfF4RPw6Isb3WzfQ2aF26Ywq39aIWNPvPHOJiBsi4tmIeKxu7NiIuDsinqq+\nHtPPjLOJiCURcU9EbKneL1dU4wOdPyJeFxEPRsSPq9yfqcZPiogHqtzfrk4CGUgRcVhE/Cgi7qiW\nu57dom/sbuDUzHwL8J/AlQARcQq1M4jeDCwHvlJd7mGQPAb8KXBv/eAwZK+7dMa5wCnAxVXuQXUj\ntW1Zbw2wMTOXAhur5UG0B1idmScDZwKXV9t60PO/ApyVmW8FTgOWR8SZwOeBa6vcLwCX9TFjM1cA\nW+qWu57dom8gM7+XmXuqxfupfQ4AapdxmMrMVzLzZ8BWapd7GBiZuSUzG334bOCzM2SXzsjMe4Hn\n9xteAayr7q8DLuhpqBZl5s7MfLi6/xK14lnMgOfPmplq8YjqlsBZwHeq8YHLvVdEnAicD3y9Wg56\nkN2ib+6DwF3V/UaXdljc80TzMwzZhyFjM6OZuRNqZQqc0Oc8TUXEGPA24AGGIH916OMR4Flq//v+\nCfBi3c7ZIL9v/h74OPDravl36UH2of/FI/MVEd8Hfq/Bqk9m5u3VnE9S+y/uTXsf1mB+z09baiV7\no4c1GBu0U66GIWNRImIEuAX4aGb+oraDOdgy81XgtOpnZ7cBJzea1ttUzUXEu4FnM3NTREzsHW4w\ntePZD9miz8x3zbU+IlYC7wbOzt+eg9r00g690Cz7LAYiexPDkLGZXRGxKDN3RsQianudAykijqBW\n8jdl5q3V8NDkz8wXI2Ka2s8YFkbE4dWe8aC+b94JvCcizgNeB/wOtT38rmf30E0D1S9N+WvgPZn5\nq7pV64GLIuLIiDgJWAo82I+M8zAM2Uu4dMZ6YGV1fyUw2/+w+qo6Nnw9sCUzv1i3aqDzR8Txe8+C\ni4ijgHdR+/nCPcD7qmkDlxsgM6/MzBMzc4zae/sHmXkJvciemd72u1H7QeXTwCPV7at16z5J7Zjg\nk8C5/c7aIPt7qe0ZvwLsAr47LNmrjOdRO9PpJ9QORfU90xxZvwXsBP6v2uaXUTvmuhF4qvp6bL9z\nzpL9j6gdIni07n1+3qDnB94C/KjK/RjwN9X4H1DbcdkK/AtwZL+zNvk+JoA7epXdT8ZKUuE8dCNJ\nhbPoJalwFr0kFc6il6TCWfSSVDiLXpIKZ9FLUuEsekkq3P8D2uEiYqe845oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f19203df908>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A.hist(bins=30);"
   ]
  },
  {
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "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": [
    "### 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 1992, one could consider a _marginal structural logistic model_\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "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": 53,
     "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": 54,
   "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": 55,
   "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": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sw_AV = numer * weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "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": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 83)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "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": 59,
     "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": [
    "This matches the Stata regression. The R regression is a little different."
   ]
  },
  {
   "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": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1629, 64)"
      ]
     },
     "execution_count": 60,
     "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": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nhefs_all['constant'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add dummy features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "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": 63,
   "metadata": {
    "collapsed": true
   },
   "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": 64,
   "metadata": {
    "collapsed": true
   },
   "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": 65,
   "metadata": {
    "collapsed": true
   },
   "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": 66,
   "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": 67,
   "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": 68,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sw_A = ip_numer / ip_denom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "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": 70,
   "metadata": {
    "collapsed": true
   },
   "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": 71,
   "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": 72,
   "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": 73,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sw_C = ip_numer / ip_denom\n",
    "sw_C[nhefs_all.censored == 1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "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": 75,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sw_AC = sw_A * sw_C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "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": 77,
   "metadata": {
    "collapsed": true
   },
   "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": 78,
   "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": 78,
     "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}