chapter12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING *** OLE2 inconsistency: SSCS size is 0 but SSAT size is non-zero\n"
     ]
    }
   ],
   "source": [
    "nhefs_all = pd.read_excel('NHEFS.xls')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just a look at a couple basic details of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1629, 64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['seqn', 'qsmk', 'death', 'yrdth', 'modth', 'dadth', 'sbp', 'dbp', 'sex',\n",
       "       'age', 'race', 'income', 'marital', 'school', 'education', 'ht', 'wt71',\n",
       "       'wt82', 'wt82_71', 'birthplace', 'smokeintensity', 'smkintensity82_71',\n",
       "       'smokeyrs', 'asthma', 'bronch', 'tb', 'hf', 'hbp', 'pepticulcer',\n",
       "       'colitis', 'hepatitis', 'chroniccough', 'hayfever', 'diabetes', 'polio',\n",
       "       'tumor', 'nervousbreak', 'alcoholpy', 'alcoholfreq', 'alcoholtype',\n",
       "       'alcoholhowmuch', 'pica', 'headache', 'otherpain', 'weakheart',\n",
       "       'allergies', 'nerves', 'lackpep', 'hbpmed', 'boweltrouble', 'wtloss',\n",
       "       'infection', 'active', 'exercise', 'birthcontrol', 'pregnancies',\n",
       "       'cholesterol', 'hightax82', 'price71', 'price82', 'tax71', 'tax82',\n",
       "       'price71_82', 'tax71_82'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"We restricted the analysis to NHEFS individuals with known sex, age, race, ...\" (pg 149, margin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "restriction_cols = [\n",
    "    'sex', 'age', 'race', 'wt82', 'ht', 'school', 'alcoholpy', 'smokeintensity'\n",
    "]\n",
    "missing = nhefs_all[restriction_cols].isnull().any(axis=1)\n",
    "nhefs = nhefs_all.loc[~missing]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 64)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to add some columns to help calculate Table 12.1, and a `constant` column, which will be useful for modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs['constant'] = 1\n",
    "nhefs['university'] = (nhefs.education == 5).astype('int')\n",
    "nhefs['inactive'] = (nhefs.active == 2).astype('int')\n",
    "nhefs['no_exercise'] = (nhefs.exercise == 2).astype('int')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Average weight gains in quitters and non-quitters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average weight gain\n",
      "      quitters: 4.5 kg\n",
      "  non-quitters: 2.0 kg\n"
     ]
    }
   ],
   "source": [
    "ave_gain_quit = nhefs[nhefs.qsmk == 1].wt82_71.mean()\n",
    "ave_gain_noquit = nhefs[nhefs.qsmk == 0].wt82_71.mean()\n",
    "\n",
    "print(\"Average weight gain\")\n",
    "print(\"      quitters: {:>0.1f} kg\".format(ave_gain_quit))\n",
    "print(\"  non-quitters: {:>0.1f} kg\".format(ave_gain_noquit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a simple linear model to get a confidence interval on weight difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>constant</th> <td>    1.9845</td> <td>    0.229</td> <td>    8.672</td> <td> 0.000</td> <td>    1.536</td> <td>    2.433</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    2.5406</td> <td>    0.451</td> <td>    5.632</td> <td> 0.000</td> <td>    1.656</td> <td>    3.425</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ols = sm.OLS(nhefs.wt82_71, nhefs[['constant', 'qsmk']])\n",
    "res = ols.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            estimate   95% C.I.\n",
      "difference     2.54   (1.7, 3.4)\n"
     ]
    }
   ],
   "source": [
    "est = res.params.qsmk\n",
    "conf_ints = res.conf_int(alpha=0.05, cols=None)\n",
    "lo, hi = conf_ints[0]['qsmk'], conf_ints[1]['qsmk']\n",
    "\n",
    "print('            estimate   95% C.I.')\n",
    "print('difference   {:>6.2f}   ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create Table 12.1 in the margin of pg 149."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style type=\"text/css\">\n",
       "</style>\n",
       "<table id=\"T_62c15_\">\n",
       "  <thead>\n",
       "    <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>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row0\" class=\"row_heading level0 row0\" >Age, years</th>\n",
       "      <td id=\"T_62c15_row0_col0\" class=\"data row0 col0\" >46.2</td>\n",
       "      <td id=\"T_62c15_row0_col1\" class=\"data row0 col1\" >42.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row1\" class=\"row_heading level0 row1\" >Men, %</th>\n",
       "      <td id=\"T_62c15_row1_col0\" class=\"data row1 col0\" >54.6</td>\n",
       "      <td id=\"T_62c15_row1_col1\" class=\"data row1 col1\" >46.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row2\" class=\"row_heading level0 row2\" >White, %</th>\n",
       "      <td id=\"T_62c15_row2_col0\" class=\"data row2 col0\" >91.1</td>\n",
       "      <td id=\"T_62c15_row2_col1\" class=\"data row2 col1\" >85.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row3\" class=\"row_heading level0 row3\" >University education, %</th>\n",
       "      <td id=\"T_62c15_row3_col0\" class=\"data row3 col0\" >15.4</td>\n",
       "      <td id=\"T_62c15_row3_col1\" class=\"data row3 col1\" >9.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row4\" class=\"row_heading level0 row4\" >Weight, kg</th>\n",
       "      <td id=\"T_62c15_row4_col0\" class=\"data row4 col0\" >72.4</td>\n",
       "      <td id=\"T_62c15_row4_col1\" class=\"data row4 col1\" >70.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row5\" class=\"row_heading level0 row5\" >Cigarettes/day</th>\n",
       "      <td id=\"T_62c15_row5_col0\" class=\"data row5 col0\" >18.6</td>\n",
       "      <td id=\"T_62c15_row5_col1\" class=\"data row5 col1\" >21.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row6\" class=\"row_heading level0 row6\" >Years smoking</th>\n",
       "      <td id=\"T_62c15_row6_col0\" class=\"data row6 col0\" >26.0</td>\n",
       "      <td id=\"T_62c15_row6_col1\" class=\"data row6 col1\" >24.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row7\" class=\"row_heading level0 row7\" >Little or no exercise, %</th>\n",
       "      <td id=\"T_62c15_row7_col0\" class=\"data row7 col0\" >40.7</td>\n",
       "      <td id=\"T_62c15_row7_col1\" class=\"data row7 col1\" >37.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_62c15_level0_row8\" class=\"row_heading level0 row8\" >Inactive daily life, %</th>\n",
       "      <td id=\"T_62c15_row8_col0\" class=\"data row8 col0\" >11.2</td>\n",
       "      <td id=\"T_62c15_row8_col1\" class=\"data row8 col1\" >8.9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n"
      ],
      "text/plain": [
       "<pandas.io.formats.style.Styler at 0x7fb7d030d9e8>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summaries = OrderedDict((\n",
    "    ('age', 'mean'),\n",
    "    ('sex', lambda x: (100 * (x == 0)).mean()),\n",
    "    ('race', lambda x: (100 * (x == 0)).mean()),\n",
    "    ('university', lambda x: 100 * x.mean()),\n",
    "    ('wt71', 'mean'),\n",
    "    ('smokeintensity', 'mean'),\n",
    "    ('smokeyrs', 'mean'),\n",
    "    ('no_exercise', lambda x: 100 * x.mean()),\n",
    "    ('inactive', lambda x: 100 * x.mean())\n",
    "))\n",
    "\n",
    "table = nhefs.groupby('qsmk').agg(summaries)\n",
    "table.sort_index(ascending=False, inplace=True)\n",
    "table = table.T\n",
    "\n",
    "table.index = [\n",
    "    'Age, years',\n",
    "    'Men, %',\n",
    "    'White, %',\n",
    "    'University education, %',\n",
    "    'Weight, kg',\n",
    "    'Cigarettes/day',\n",
    "    'Years smoking',\n",
    "    'Little or no exercise, %',\n",
    "    'Inactive daily life, %'\n",
    "]\n",
    "\n",
    "table.style.format(\"{:>0.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to be modeling with squared terms and some categorical features. Here we'll explicitly add squared features and dummy features to the data. In later chapters we'll use Statsmodels' formula syntax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Squared features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
    "    nhefs['{}^2'.format(col)] = nhefs[col] * nhefs[col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dummy features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "edu_dummies = pd.get_dummies(nhefs.education, prefix='edu')\n",
    "exercise_dummies = pd.get_dummies(nhefs.exercise, prefix='exercise')\n",
    "active_dummies = pd.get_dummies(nhefs.active, prefix='active')\n",
    "\n",
    "nhefs = pd.concat(\n",
    "    [nhefs, edu_dummies, exercise_dummies, active_dummies],\n",
    "    axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to be creating a lot of IP weights from logistic regressions so a function will help reduce the work. The following function creates the denominators of the IP weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logit_ip_f(y, X):\n",
    "    \"\"\"\n",
    "    Create the f(y|X) part of IP weights\n",
    "    from logistic regression\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : Pandas Series\n",
    "    X : Pandas DataFrame\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Numpy array of IP weights\n",
    "    \n",
    "    \"\"\"\n",
    "    model = sm.Logit(y, X)\n",
    "    res = model.fit()\n",
    "    weights = np.zeros(X.shape[0])\n",
    "    weights[y == 1] = res.predict(X.loc[y == 1])\n",
    "    weights[y == 0] = (1 - res.predict(X.loc[y == 0]))\n",
    "    return weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.535408\n",
      "         Iterations 6\n"
     ]
    }
   ],
   "source": [
    "X_ip = nhefs[[\n",
    "    'constant',\n",
    "    'sex', 'race', 'age', 'age^2', 'edu_2', 'edu_3', 'edu_4', 'edu_5',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2', 'wt71', 'wt71^2'\n",
    "]]\n",
    "\n",
    "denoms = logit_ip_f(nhefs.qsmk, X_ip)\n",
    "weights = 1 / denoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IP weights\n",
      "   min:  1.05   expected:  1.05\n",
      "   max: 16.70   expected: 16.70\n",
      "  mean:  2.00   expected:  2.00\n"
     ]
    }
   ],
   "source": [
    "print('IP weights')\n",
    "print('   min: {:>5.2f}   expected:  1.05'.format(weights.min()))\n",
    "print('   max: {:>5.2f}   expected: 16.70'.format(weights.max()))\n",
    "print('  mean: {:>5.2f}   expected:  2.00'.format(weights.mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFpCAYAAABjxXptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEslJREFUeJzt3X+s3fdd3/HXezEBWgZJG7crdjS3IxS6CtbIK4FqaGuANW3V5A8iFTFqdZEiTV0plI2mQ1qlTZrSDdEf2tQpSkJSrWqpQrZEkAFRWoYmrRlOWvorsFilS9yExihpYFQMIt77434zLo4d2/fc3OP39eMhWed8P+dzzvl8fW0/7/d7zj2u7g4AMMtfW/cCAIAzJ+AAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADDQnnUv4NlcdNFFfeDAgXUvAwB2zH333feH3b33VPPO6oAfOHAghw8fXvcyAGDHVNX/Pp15TqEDwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEBn9f9G9lw4cN2vbuvjffn6N2zr4wHA6XAEDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AA50y4FV1c1U9VlWf3zT2gqq6u6oeXC4vXMarqj5YVUeq6rNVdemm+xxa5j9YVYeem90BgHPD6RyB35LkdceNXZfknu6+JMk9y3aSXJHkkuXXtUk+lGwEP8l7knxfklcnec/T0QcAztwpA97dv5Xk8eOGr0xy63L91iRXbRr/cG/4VJILquolSf5hkru7+/HufiLJ3XnmNwUAwGna6mvgL+7uR5NkuXzRMr4vycOb5h1dxk42DgBswXa/ia1OMNbPMv7MB6i6tqoOV9XhY8eObeviAGC32GrAv7qcGs9y+dgyfjTJxZvm7U/yyLOMP0N339DdB7v74N69e7e4PADY3bYa8DuTPP1O8kNJ7tg0/pbl3eiXJXlyOcX+60l+pKouXN689iPLGACwBXtONaGqPprk7ye5qKqOZuPd5Ncn+XhVXZPkoSRXL9PvSvL6JEeSfD3JW5Okux+vqn+d5LeXef+qu49/YxwAcJpOGfDu/rGT3HT5CeZ2kred5HFuTnLzGa0OADghn8QGAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADLRSwKvqp6vqC1X1+ar6aFV9U1W9tKruraoHq+qXqur8Ze43LttHltsPbMcOAMC5aMsBr6p9SX4yycHufmWS85K8Ocl7k7yvuy9J8kSSa5a7XJPkie7+jiTvW+YBAFuw6in0PUm+uar2JHlekkeTvDbJbcvttya5arl+5bKd5fbLq6pWfH4AOCdtOeDd/ZUkP5/koWyE+8kk9yX5Wnc/tUw7mmTfcn1fkoeX+z61zH/hVp8fAM5lq5xCvzAbR9UvTfLtSZ6f5IoTTO2n7/Ist21+3Gur6nBVHT527NhWlwcAu9oqp9B/KMnvd/ex7v7zJLcn+YEkFyyn1JNkf5JHlutHk1ycJMvt35bk8eMftLtv6O6D3X1w7969KywPAHavVQL+UJLLqup5y2vZlyf5YpJPJvnRZc6hJHcs1+9ctrPc/onufsYROABwaqu8Bn5vNt6Mdn+Szy2PdUOSdyV5Z1UdycZr3Dctd7kpyQuX8XcmuW6FdQPAOW3PqaecXHe/J8l7jhv+UpJXn2Dunya5epXnAwA2+CQ2ABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYaKWAV9UFVXVbVf1uVT1QVd9fVS+oqrur6sHl8sJlblXVB6vqSFV9tqou3Z5dAIBzz6pH4B9I8mvd/V1JvjfJA0muS3JPd1+S5J5lO0muSHLJ8uvaJB9a8bkB4Jy15YBX1bcm+cEkNyVJd/9Zd38tyZVJbl2m3ZrkquX6lUk+3Bs+leSCqnrJllcOAOewVY7AX5bkWJJfrKpPV9WNVfX8JC/u7keTZLl80TJ/X5KHN93/6DIGAJyhVQK+J8mlST7U3a9K8if5y9PlJ1InGOtnTKq6tqoOV9XhY8eOrbA8ANi9Vgn40SRHu/veZfu2bAT9q0+fGl8uH9s0/+JN99+f5JHjH7S7b+jug919cO/evSssDwB2ry0HvLv/IMnDVfXyZejyJF9McmeSQ8vYoSR3LNfvTPKW5d3olyV58ulT7QDAmdmz4v3fnuQjVXV+ki8leWs2vin4eFVdk+ShJFcvc+9K8vokR5J8fZkLAGzBSgHv7s8kOXiCmy4/wdxO8rZVng8A2OCT2ABgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIAEHgIEEHAAGEnAAGEjAAWAgAQeAgQQcAAYScAAYSMABYCABB4CBVg54VZ1XVZ+uql9Ztl9aVfdW1YNV9UtVdf4y/o3L9pHl9gOrPjcAnKu24wj8HUke2LT93iTv6+5LkjyR5Jpl/JokT3T3dyR53zIPANiClQJeVfuTvCHJjct2JXltktuWKbcmuWq5fuWyneX2y5f5AMAZWvUI/P1JfjbJXyzbL0zyte5+atk+mmTfcn1fkoeTZLn9yWU+AHCGthzwqnpjkse6+77NwyeY2qdx2+bHvbaqDlfV4WPHjm11eQCwq61yBP6aJG+qqi8n+Vg2Tp2/P8kFVbVnmbM/ySPL9aNJLk6S5fZvS/L48Q/a3Td098HuPrh3794VlgcAu9eWA97d7+7u/d19IMmbk3yiu388ySeT/Ogy7VCSO5brdy7bWW7/RHc/4wgcADi15+LnwN+V5J1VdSQbr3HftIzflOSFy/g7k1z3HDw3AJwT9px6yql1928m+c3l+peSvPoEc/40ydXb8XwAcK7zSWwAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADDQlgNeVRdX1Ser6oGq+kJVvWMZf0FV3V1VDy6XFy7jVVUfrKojVfXZqrp0u3YCAM41qxyBP5XkZ7r7u5NcluRtVfWKJNcluae7L0lyz7KdJFckuWT5dW2SD63w3ABwTttywLv70e6+f7n+x0keSLIvyZVJbl2m3ZrkquX6lUk+3Bs+leSCqnrJllcOAOewbXkNvKoOJHlVknuTvLi7H002Ip/kRcu0fUke3nS3o8sYAHCGVg54VX1Lkl9O8lPd/UfPNvUEY32Cx7u2qg5X1eFjx46tujwA2JVWCnhVfUM24v2R7r59Gf7q06fGl8vHlvGjSS7edPf9SR45/jG7+4buPtjdB/fu3bvK8gBg11rlXeiV5KYkD3T3L2y66c4kh5brh5LcsWn8Lcu70S9L8uTTp9oBgDOzZ4X7vibJTyT5XFV9Zhn7F0muT/LxqromyUNJrl5uuyvJ65McSfL1JG9d4bkB4Jy25YB393/PiV/XTpLLTzC/k7xtq88HAPwln8QGAAMJOAAMJOAAMJCAA8BAAg4AAwk4AAwk4AAwkIADwEACDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BAe9a9gOkOXPer2/p4X77+Ddv6eADsTo7AAWAgAQeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIJ/EdpbxyW4AnA5H4AAwkIADwEBOoe9y231KPnFaHuBs4AgcAAYScAAYSMABYCCvgXPG/KgbwPo5AgeAgQQcAAYScAAYSMABYCABB4CBBBwABhJwABhIwAFgIB/kwtr5YBiAM+cIHAAGcgTOruO/UAXOBY7AAWCgHT8Cr6rXJflAkvOS3Njd1+/0GuBMeZ0eONvs6BF4VZ2X5D8kuSLJK5L8WFW9YifXAAC7wU4fgb86yZHu/lKSVNXHklyZ5Is7vA5YK0f0wKp2OuD7kjy8aftoku/b4TXArvNcvHFvu52L32T4Ru3ss5u+Jjsd8DrBWP+VCVXXJrl22fw/VfV7z/mqttdFSf5w3YvYRvbn7DZmf+q9pzVtzP6cpm3dn9P8PXyu+NqcwHP0NfmbpzNppwN+NMnFm7b3J3lk84TuviHJDTu5qO1UVYe7++C617Fd7M/Zzf6c3XbT/uymfUl2x/7s9I+R/XaSS6rqpVV1fpI3J7lzh9cAAOPt6BF4dz9VVf80ya9n48fIbu7uL+zkGgBgN9jxnwPv7ruS3LXTz7uDxp7+Pwn7c3azP2e33bQ/u2lfkl2wP9Xdp54FAJxVfJQqAAwk4Nukqi6uqk9W1QNV9YWqese617Sqqjqvqj5dVb+y7rWsqqouqKrbqup3l6/R9697Tauoqp9e/px9vqo+WlXftO41nYmqurmqHquqz28ae0FV3V1VDy6XF65zjWfiJPvz75Y/b5+tqv9cVResc41n4kT7s+m2f1ZVXVUXrWNtW3Gy/amqt1fV7y1/l/7tuta3VQK+fZ5K8jPd/d1JLkvytl3wMbHvSPLAuhexTT6Q5Ne6+7uSfG8G71dV7Uvyk0kOdvcrs/GG0Devd1Vn7JYkrztu7Lok93T3JUnuWbanuCXP3J+7k7yyu78nyf9K8u6dXtQKbskz9ydVdXGSH07y0E4vaEW35Lj9qap/kI1PAv2e7v7bSX5+DetaiYBvk+5+tLvvX67/cTYCsW+9q9q6qtqf5A1Jblz3WlZVVd+a5AeT3JQk3f1n3f219a5qZXuSfHNV7UnyvBz3eQpnu+7+rSSPHzd8ZZJbl+u3JrlqRxe1ghPtT3f/Rnc/tWx+KhufezHCSb4+SfK+JD+b4z6A62x3kv35J0mu7+7/u8x5bMcXtiIBfw5U1YEkr0py73pXspL3Z+Mv6l+seyHb4GVJjiX5xeUlgRur6vnrXtRWdfdXsnG08FCSR5M82d2/sd5VbYsXd/ejycY3xEletOb1bKd/nOS/rnsRq6iqNyX5Snf/zrrXsk2+M8nfq6p7q+q/VdXfXfeCzpSAb7Oq+pYkv5zkp7r7j9a9nq2oqjcmeay771v3WrbJniSXJvlQd78qyZ9k1unZv2J5bfjKJC9N8u1Jnl9V/2i9q+JkqurnsvES20fWvZatqqrnJfm5JP9y3WvZRnuSXJiNlzz/eZKPV9WJPu77rCXg26iqviEb8f5Id9++7vWs4DVJ3lRVX07ysSSvrar/tN4lreRokqPd/fQZkduyEfSpfijJ73f3se7+8yS3J/mBNa9pO3y1ql6SJMvluFOax6uqQ0nemOTHe/bP7P6tbHzD+DvLvwv7k9xfVX9jratazdEkt/eG/5mNs41j3piXCPi2Wb5zuynJA939C+tezyq6+93dvb+7D2TjzVGf6O6xR3jd/QdJHq6qly9Dl2f2f2H7UJLLqup5y5+7yzP4TXmb3Jnk0HL9UJI71riWlVXV65K8K8mbuvvr617PKrr7c939ou4+sPy7cDTJpcvfran+S5LXJklVfWeS8zPsP2sR8O3zmiQ/kY2j1c8sv16/7kXx/709yUeq6rNJ/k6Sf7Pm9WzZcibhtiT3J/lcNv4ej/pUqar6aJL/keTlVXW0qq5Jcn2SH66qB7PxTufr17nGM3GS/fn3Sf56kruXfw/+41oXeQZOsj9jnWR/bk7ysuVHyz6W5NC0syQ+iQ0ABnIEDgADCTgADCTgADCQgAPAQAIOAAMJOAAMJOAAMJCAA8BA/w+DxArUXQnKyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.hist(weights, bins=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the main model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = nhefs.wt82_71\n",
    "X = nhefs[['constant', 'qsmk']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Without clustering, 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": [
    "With clustering, GEE and WLS give the right coefficients and better standard errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "WLS with clustering"
   ]
  },
  {
   "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.915</td> <td> 0.000</td> <td>    1.339</td> <td>    2.221</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>     <td>    3.4405</td> <td>    0.526</td> <td>    6.543</td> <td> 0.000</td> <td>    2.410</td> <td>    4.471</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wls = sm.WLS(y, X, weights=weights) \n",
    "res = wls.fit(cov_type=\"cluster\", cov_kwds={\"groups\": nhefs.seqn})\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GEE with clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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": 23,
     "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": 24,
   "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": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>763.607760</td>\n",
       "      <td>763.623497</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>801.748892</td>\n",
       "      <td>797.200691</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk           0           1\n",
       "sex                         \n",
       "0     763.607760  763.623497\n",
       "1     801.748892  797.200691"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.crosstab(nhefs.sex, nhefs.qsmk, weights, aggfunc='sum')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_indices = (nhefs.race == 0) & (nhefs.sex == 1)\n",
    "subset = nhefs.loc[subset_indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now a check for positivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "crosstab = pd.crosstab(subset.age, subset.qsmk).sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAF7CAYAAAAgxam6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXmcZGV5739v7VvX2vv0MgszzMLMMM2IIgRFFkUWUQZxwQsJbvfGBDUmFxW4EyVivGoi994EUeNoYqIQUJABXCMoKsxMz87A7L3vXVVd+/reP06d6urqWs7pOrU/38+nP91dXXX67Z6e+tWz/R7GOQdBEARBEI2BqtoHIAiCIAhCOUjYCYIgCKKBIGEnCIIgiAaChJ0gCIIgGggSdoIgCIJoIEjYCYIgCKKBIGEnCIIgiAaChJ0gCIIgGoiKCTtjrJcx9l+MsROMseOMsXtSt+9mjI0xxg6l3t5ZqTMRBEEQRKPBKuU8xxjrAtDFOR9kjLUAOADgFgDvBeDnnH+1IgchCIIgiAZGU6lvxDmfADCR+tjHGDsBYNVKrtXa2spXr16t4OkIgiAIonY5cODALOe8Tcp9KybsmTDGVgPYAeBlAJcD+ARj7L8B2A/grzjn7kKPX716Nfbv31/uYxIEQRBETcAYG5J634o3zzHGLACeAPBJzvkCgH8GsA7AxRAi+q/ledxHGWP7GWP7Z2ZmKnZegiAIgqgnKirsjDEtBFH/Aef8SQDgnE9xzhOc8ySAbwG4NNdjOeePcs53cs53trVJykYQBEEQRNNRya54BuA7AE5wzr+ecXtXxt3eDeBYpc5EEARBEI1GJWvslwP4EICjjLFDqds+B+D9jLGLAXAA5wF8rIJnIgiCIIiGopJd8b8DwHJ86dlKnYEgCKKRSCaTmJ2dhcfjQSKRqPZxiBIxGAzo6emBVqst6TpV6YonCIIgSmd0dBSMMaxevRparRZCxZOoRzjnmJubw+joKNasWVPStchSliAIok4JBAJYtWoVdDodiXqdwxiDy+VCOBwu+Vok7ARBEHWMSkVP442CUi/O6C+CIAiCaFquv/56fO9736v2MRSFhJ0gCIJoWp577jnceeedAIA9e/bgiiuuWPL1u+66C/fdd181jrZiSNgJgiAIokzE4/GKf8+mF/bzswH812vT1T4GQRBEw7F69Wp89atfxbZt22Cz2XD77benm8O+9a1v4YILLoDT6cTNN9+M8fHx9OMYY3jkkUewfv16OBwO/Pmf/zkKbSL9xS9+gY0bN8Jms+ETn/gE3vKWt+Db3/42AGD37t2444470vc9f/48GGNpwX3rW9+Kb3/72zhx4gQ+/vGP4w9/+AMsFgvsdjseffRR/OAHP8BXvvIVWCwW3HTTTQCA8fFx3HrrrWhra8OaNWvw8MMPp6+/e/du7Nq1C3fccQesViv27NmDV155BTt37oTVakVHRwc+/elPK/dLzkHTC/sTg6O4+3v7kExWZn0tQRBEM/HYY4/h+eefx7lz53DkyBHs2bMHv/71r/HZz34Wjz32GCYmJtDf34/3ve99Sx73zDPPYN++fTh8+DAee+wx/OxnP8t5/dnZWdx666148MEHMTs7i3Xr1uGll16Sfc5NmzbhkUcewWWXXQa/3w+Px4OPfvSj+OAHP4i/+Zu/gd/vx09/+lMkk0ncdNNN2L59O8bGxvCrX/0K//iP/7jkfE899RR27doFj8eDD37wg7jnnntwzz33YGFhAWfOnMF73/te2eeTQ9PPsduMWiQ54I/GYTWUZgpAEARRTf72p8fx6vhCWb/H5m4r/tdNWyTf/y//8i/R3d0NALjppptw6NAh7Nu3D3/2Z3+GgYEBAMBDDz0Eh8OB8+fPQ1zJfe+998Jut8Nut+Oqq67CoUOH8I53vGPZ9Z999lls3rwZu3btAgB88pOfxNe+lnOXmCLs27cPMzMzeOCBBwAAa9euxUc+8hH88Ic/xNvf/nYAwGWXXYZbbrkFAGA0GqHVanH69GnMzs6itbUVb3rTm8p2PoAidliNgph7g7Eqn4QgCKLx6OzsTH9sMpng9/sxPj6O/v7+9O0WiwUulwtjY2MFHwcAW7ZsgcVigcViwW9/+1uMj4+jt7c3fV/G2JLPlWZoaAjj4+PpFx12ux1f+tKXMDU1lb5P9vf/zne+g5MnT2Ljxo14wxvegGeeeaZs5wMoYodNFPZQDOX7UyAIgig/ciLpatLd3Y2hocX14oFAAHNzc1i1alXRxx4/fnzJ52fPnsXIyEj6c875ks/NZjOCwWD688nJybzXzjVHnn1bb28v1qxZg1OnTkm+zvr16/Ef//EfSCaTePLJJ7Fr1y7Mzc3BbDbnvUYpNH3EninsBEEQRPn5wAc+gO9+97s4dOgQIpEIPve5z+GNb3xjOg0vhxtuuAHHjx/Hk08+iXg8jocffniJeF988cV48cUXMTw8DK/Xi4ceeijvtTo6OjA6OopoNLrktrNnz6Y/v/TSS2G1WvH3f//3CIVCSCQSOHbsGPbt25f3uv/2b/+GmZkZqFQq2O12AIBarZb9s0qFhJ2EnSAIoqJcffXV+OIXv4hbb70VXV1dOHPmDH74wx+u6Fqtra14/PHHce+998LlcuHUqVO4/PLL01+/9tprcfvtt2Pbtm245JJLcOONN+a91tve9jZs2bIFnZ2daG1tBQDcfffdePXVV2G323HLLbdArVbjpz/9KQ4dOoQ1a9agtbUVH/7wh+H1evNe9/nnn0+XEO655x788Ic/hMFgWNHPKwVWaISgVtm5cyffv3+/Itea8IZw2UO/xkPv2Yr3X9qnyDUJgiAqwYkTJ7Bp06ZqH6PmeOtb34o77rgDH/7wh6t9FNnk+zdljB3gnO+Ucg2K2CliJwiCIBqIphd2o1YNrZrBQ13xBEEQRAPQ9F3xjDHYjFqK2AmCIBqE3/zmN9U+QlVp+ogdEGbZF0jYCYIgiAaAhB2giJ0gCIJoGEjYAdhJ2AmCIIgGgYQdFLETBEEQjQMJOwRh9wSjxe9IEARBEDUOCTsEYfdF4rS6lSAIogb4+Mc/ji9+8YvVPkbdQsIOoSuec8AXjlf7KARBEE3PI488gvvvvx+AMLrW09Oz5Ou7d+/GHXfcUY2j1QUk7CD3OYIgiGYiHm/sII6EHYDdpANAwk4QBKEkBw8exMDAAFpaWnD77bfjfe97H+677z7s2bMHV1xxxZL7MsZw+vRpAMBdd92F++67D4FAANdffz3Gx8fTO9j//d//HV/60pfwox/9CBaLBdu3bwcAeL1e3H333ejq6sKqVatw3333IZFIAAD27NmDyy+/HJ/61KfgdDqxe/dunD59Gm95y1tgs9nQ2tqK22+/vbK/nDJCwg6K2AmCIJQmGo3illtuwYc+9CHMz8/jtttuwxNPPCHrGmazGc899xy6u7vh9/vh9/vxgQ98AJ/73Odw++23w+/34/DhwwCAO++8ExqNBqdPn8bBgwfx85//HN/+9rfT13r55Zexdu1aTE9P4/Of/zzuv/9+XHfddXC73RgdHcVf/MVfKPrzV5Omt5QFFoXdE6LOeIIg6pjn7gUmj5b3e3RuBa7/ctG7/fGPf0QsFsMnP/lJMMawa9cufP3rXy/LkaampvDcc8/B4/HAaDTCbDbjU5/6FB599FF87GMfAwB0d3enxVuj0UCr1WJoaAjj4+Po6elZlkGoZyhiB0XsBEEQSjM+Po5Vq1aBMZa+rb+/vyzfa2hoCLFYDF1dXbDb7bDb7fjYxz6G6enp9H16e3uXPOYrX/kKOOe49NJLsWXLFvzLv/xLWc5WDShiBwk7QRANgoRIulJ0dXVhbGwMnPO0uA8PD2PdunUwm80IBoPp+05OTua9TuYLg3y39fb2Qq/XY3Z2FhpNblnLfkxnZye+9a1vAQB+97vf4ZprrsGVV16JCy64QNoPWMNQxA7AoFVBp1aRsBMEQSjEZZddBo1Gg4cffhjxeBxPPvkkXnnlFQDA9u3bcfz4cRw6dAjhcBi7d+/Oe52Ojg7Mzc3B6/Uuue38+fNIJpMAhBcR1113Hf7qr/4KCwsLSCaTOHPmDF544YW813388ccxOjoKAHA4HGCMQa1WK/CTVx8SdqRWt5pow1sxfnxwFOdmA9U+BkEQdYBOp8OTTz6JPXv2wOFw4Ec/+hHe8573AAA2bNiABx54ANdccw3Wr19fsL69ceNGvP/978fatWtht9sxPj6O2267DQDgcrkwMDAAAPj+97+PaDSKzZs3w+FwYNeuXZiYmMh73X379uGNb3wjLBYLbr75ZnzjG9/AmjVrFPwNVA/Gef25re3cuZPv379f0Wte8/UXsKHDgn/64CWKXrdRSCY51t/3HD70pn7svnlLtY9DEASAEydOYNOmTdU+hmTuuusu9PT04MEHH6z2UWqWfP+mjLEDnPOdUq5BEXsKwS+eIvZ8+MJxJJIcE95QtY9CEARBFICEPQVteCuM+LuZ8IarfBKCIAiiENQVn8Jm1OLklK/ax6hZSNgJgiiVPXv2VPsITQFF7CkoYi+M+LuZ9UcQjSerfBqCIAgiHyTsKaxGbbqOTCxHFHbOgakFitoJgiBqFRL2FPaUSY0vTFF7LjKzGZSOJ4jaQZzlJuofpabUSNhTpP3iqTM+J0uFnTrjCaIWMJvNGBsbQzQaVUwUiOrAOcfc3BwMBkPJ16LmuRRkK1sYTygKFQOSnCJ2gqgVenp6MDs7i6GhoYbfMd4MGAwG9PT0lHwdEvYUNhMJeyEWQjE4zXpE4glMeChiJ4haQKVSob29He3t7dU+ClFDkLCnoIi9MN5QDHaTFmqmwzhF7ARBEDULCXsKEvbCeEMx2IxaWPQaTJKwEwRB1CzUPJeChL0worB32w3UPEcQBFHDkLCnMGjV0GtodWs+RGHvtBox648iEk9U+0gEQRBEDkjYM7AZtfDSuFtOPEFB2LvswijGlDdS5RMRBEEQuSBhz4BsZXOTSHL4wnFYjVp024wAgHFKxxMEQdQkJOwZkLDnRnTjsxu16LQJETs10BEEQdQmJOwZkLDnRvyd2IxadKWEnSJ2giCI2oSEPQObiYQ9F5nCbtZrYDXQyBtBEEStQsKeAUXsuUkLe8qdr9tuxLiHhJ0gCKIWIWHPwGbUwh+JI56gbUmZiItxxFn/LhvNshMEQdQqFRN2xlgvY+y/GGMnGGPHGWP3pG53MsZ+wRg7lXrvqNSZshGFayFMyxQyyUzFA0CnzUipeIIgiBqlkhF7HMBfcc43AXgTgD9njG0GcC+AX3HO1wP4VerzqkDuc7nJFvZumwFzgSjCMTKpIQiCqDUqJuyc8wnO+WDqYx+AEwBWAXgXgO+l7vY9ALdU6kzZkLDnZiEUg16jgkGrBoD0yNvUAkXtBEEQtUZVauyMsdUAdgB4GUAH53wCEMQfQNX2D9rLtLr1uy+dw4PPvKroNSuJaCcr0m1PmdRQAx1BEETNUXFhZ4xZADwB4JOc8wUZj/soY2w/Y2z/zMxMWc4mipcnGFXsmokkx//7rzP4z8FRxa5ZabKFXZxlpwY6giCI2qOiws4Y00IQ9R9wzp9M3TzFGOtKfb0LwHSux3LOH+Wc7+Sc72xrayvL+axi85yCEfsr5+Yx64/AE4zVbYpf9IkX6UrZyk5QAx1BEETNUcmueAbgOwBOcM6/nvGlpwHcmfr4TgBPVepM2ZSjxr736Hj645H5oGLXrSTZEbtRp4bdpKWInSAIogapZMR+OYAPAXgbY+xQ6u2dAL4M4FrG2CkA16Y+rwp6jRoGrXKrWxNJjuePTWJtmxkAMNwgwg4AnVYDjbwRBEHUIJpKfSPO+e8AsDxfvrpS5yiGku5zL5+bw6w/ir95+0b8zRNH6lbYF0KxtOucCLnPEQRB1CbkPJeF3ahTTNj3HpmAUavGjdu74DTrMDRXf8KeSHL4IvFlETu5zxEEQdQmJOxZ2IzatIVqKcQTSTx/bBJv29QOk06DXqepLmvsC1nmNCJdNgPcwRiZ1BAEQdQYJOxZWBVKxb98bh5zgShu3NoFAOhzmuoyFZ/tOidCnfEEQRC1CQl7FjajVpFxt71HJ2DSqfHWCwW/nX6nCWOeEGJ1tmDGUyBiB2iWnSAIotYgYc9CieY5MQ1/9aYOGHWCDWuf04REkmOizhrO8kbsKfe5evt5CIIgGh0S9ixsRi0C0URJkfUfz85jPhDFDVs707f1Ok0A6m/kTRR2u4kidoIgiHqgYuNu9YLNKPxKFkIxuCz6FV1j79HxJWl4AOh3CcI+NB/AFWgt/aAVQhR2a1bEbtCq4TBpa67GPjIfxOlpP67aWLWVAxXh9LQfTx8aAy/hGirG8L5Le9P9EgRBNAYk7FnYTToAQm15JcIupuGv2dSR3oYGAB1WA3RqVd1F7Pm64gGhga7WhP07vzuHx/aP4NUvvKPaRykr33rxLH60fwSqfM4QEkhygDHgk9dsUO5gBEFUHRL2LEq1lf3D2Tm4gzHcsK1rye1qFUOPw1h3I2/eUAwGrQp6jXrZ17rtBozVWI19IRxDMJpAOJZY8sKq0ZjyhbGtx4anP3HFiq+x4ws/x5xfuYVHBEHUBlRjz8JaorDvPTIBs06Nt2xYvqimz2WqO5MaTzCaM1oHhL3stVZjD0TiAAC3ghv6apHphQjaVlgqEnFZ9JgLRBQ6EUEQtQIJexa2Eja8xRJJPH98Etds7sgZLfY5TRieC4LzUiqjlSWXT7xIl80ITzCGULR2TGoCEeEs7kB9btKTyow/gnZracLuNOsoYieIBoSEPYtSUvG/PzMHTzCGG7Z25fx6n9MEXyReV+tbvaEY7EZdzq/VYme8LxWxexo4Yk8kOeb8pUfsrRYd5gKN+3siiGaFhD2LtLCvwFb22SMTsOg1uDJHGh4QhB1AXaXjvaH4so54kVp0n1tMxdfPiye5zAUiSHKgzWoo6TpOsw7zJOwE0XCQsGeh06hg0qnTjmtSEdPw1+ZJwwNCjR2or1n2hQKp+G67GLHXorA3rmBNLwh18VIjdqdZD3cwikSyfkpDBEEUh4Q9Bytxn3vp9Cy8ofxpeADoddSfsBeqsXekIsYJT+2k4v1NkIqf8QnCXmqNvdWiA+eN/SKIIJoREvYcrETY9x6ZQItegz/ZkN98xqzXoNWix3CdpOJjiST8OVa2ihi0arjMOozXSMTOOU9H7PMN3DwnCnvpEbvQO0HpeIJoLEjYcyB3w1s0nsTPUmn4XPPemfQ5jXUTsS+a0+S3O+i0GTBZI81zoVgCYla5kSP2aZ/wQqqtpcRxN7Pw+Fk/jbwRRCNBwp4DuRveXjozi4VwfJkpTS7qaX1regGMKXfEDtSW+5yYhgcaO70844vAatCUbMDjslDEThCNCAl7DuSm4vcemUCLQYMr1hf3gO9zmTHhDSEar/31rekFMHnG3QChga5WhF2cYQcauyt+2hdBe4kd8QDgSqXiVzrLnkwulj4IgqgdSNhzYDdq4ZEoDJxz/PLEFK7dVDwNDwgRe5IDYzXUcJaPfAtgMum0GeANxRCMVv8JXhQZo1bd0Kn4GV/pM+yAsBeBMax4lv0/B0fx5i//GpF47RgUEQRBwp4Tm1GLUCwhKaqeXAjDE4xhR59d0rX76mh9a75d7Jl0p2bZx2vAM15Mxfc4jE0QsZcu7GoVg9Okw9wKa+yvT/rgDcUolU8QNQYJew7EmrKUdPzJKT8AYH1Hi6Rr15OwF9rsJtKZcp+brIF0vD8sCHuv04SFcKwh57M554pF7EBpJjVid36j2/cSRL1Bwp4DObayp6Z8AIANEoW9vUUPvUaF4bnAyg9YIcRyhKSIvQY64wPRxYid85Uv8qll/JE4QrGEIhE7UJpfvNid38hlD4KoR0jYcyBnw9vJKR9aLbr0THAxVCqG3jrpjPeGYjBq1dBp8v+ZdNgEgamJiD0jFQ80Zmd8eoa9xFE3kdYSNryJZ5lvwN8zQdQzJOw5kLPh7eSUH+vbpUXrIsLIW/Uj3GJ4QzHYC4y6AYBeo0arRVcTi2ACaWEXyh2NGElOi65zLaV3xQOlpeLFszRyPwNB1CMk7Dmwp4TdEyr8hMc5x+lpPzZ0WGRdX1jfGqj59a2F7GQz6bIZa6R5LgHGgG57KmJvwNqv0hG7y6KDOxhDPCFv/DIcS8CX6mnwUPMcQdQUJOw5kLrhbdwbhj8Sl9w4J9LnNCEQTdR8N7E3FCs46iYiuM/VgLCH4zDrNHCahLJII6biFyN2hYTdLP6u5L0IEl9grOSxBEGUFxL2HCzW2AvPZp+U2TgnUi+d8VIj9m6boTaa5yJxmPVq2M2pjEsDCs6MLwKdWiXp30UKrlR3vdw6+3SGsDdiyYMg6hkS9hxo1SqYdeqizXOnxVG3dnmp+P46Wd8qORVvN8IXji+xdK0G/mgcZr0GLXoNNCrWoBF7GG0tejDGFLleehGMzM74mVRHvFbdmL9ngqhnSNjzIMVWVuiI18MhsSNeRGzuqvUtb9Jr7OIse3Wj9kAkjha9Bowx2E3ahhScGV8ErQql4QFhdSsAzMosC4mp+HVtFkrFE0SNQcKeBykb3k6uoHEOAIw6Ndpb9DUdsccSSQSjiXQjYSG6asR9TkjFC5vo7CZdwzbPKVVfBwBnasPbvEz3uWlfBComCDul4gmitiBhz4MQsed/wuKc4/SUT3Z9XaTfZcJQDQu7lM1uIl014j7nCy8Ku6OBI3alOuIBYQJEtQK/+BlfBE6zPt1VXy1G5oP4+L8eqIldBQRRK5Cw58FuKhyxj3lCCEQTWL+CiB0QbE9H6kHYJUTsHVYDGKu++1wgGoclI2JvtOa5WCKJuUBU0YhdpWKC+5xMYZ9OZQ7sJl1V7XtfPDWD549P4rVJX1W+P0HUIiTseShWYz+VapxbacTe5zRhciGMcKw2N2NJ2ewmotOo0GrRVz1iD0QSMOuFDXuNGLGL1q9KmdOICLay8lLxYubAadJW1b53IlX+kdv8RxCNDAl7HooJe3rUTabrnEi/ywTOgVF39cfEcuGV4BOfSZfNgPEqC7s/EodFL5zXkYrYa90ESA6iN7uSqXgAcJn1sj0Vpn1htLcsNo5W60XUROpvbqW2uATRiJCw58Fm1CIcS+bdNX1yyo/2Fr2kGnQuxFn2Wk3Hy0nFA4KwT1Rxx3w0nkQ0noRFjNjNOkRTDYCNwozC5jQiTou8VHwyyTHrj6ItlYoHqjfLLloZr3SnPEE0IiTseSi24e3U9Mob5wChxg7U7iy7fGE3VjUVL/rEZzbPAY3lPjetsJ2sSKvMDW/zwSgSSS5E7OLvuUoTCOLf3Eo31BFEI0LCngdrAVvZZJLj1JR/xY1zANBm0cOoVWOoRmfZVxKx+yJx+MLVeYL3Zwn7YiTZOA10YsTeqtAudhGnWQ9vKIaYRL/4dObAaoCjiva9nPN0w2at2zMTRCUhYc+DKAy5IvYxTwihWKKkiJ0xltryVrvCbtapoVVL+xPpSi1eqVbULu5it6Qj9sbzi5/2heEwaQuu0V0JrpRJjVuiOGZmDuxVzIx4gjGEY8KLkVmZzX8E0ciQsOehUCr+1LTQOCfXSjabWh55k+o6JyLOslergU5MxVuWpeIbK2JXuiMeWFwEI7VOPb0g/Bu3t+hhSdv3Vv73LDbOMUYRO0FkQsKeh0LCflL0iC8hYgeQjthrsXPbE5S22U1EFPZqNdD5I0KT3PJUfOM84U8rbE4jkl4EI7FOPeNfjNgF+15dVX7PYuPcujYL1dgJIgMS9jwUFnYfOqz6kjds9btMCMUS6SfKWmJBZsQumtRMVCli94eXRuz2Kjd1lQOl7WRFnOmIXdrf4fRCBBa9BiZdhstfFX7P4t/aRd1WzAeiNfkCmSCqAQl7HqwG4UkrZyp+yl9SfV2klkfe5KbitWoV2iz6dBRVaRa74tXp87ToNQ1TY+ecly9iF4VdRsSeeQ6HSVeV3/OENwS1iuHCTiuiiSR8Vd4uSBC1Agl7HjRqFSx6zbKu6mSS4/S0H+tXaEyTiTjyVoud8d5QLB31SqXLbqxexJ5VYwcAu1lblhTx65M+HBvzKn7dQiyE44jGk2URdptRC7WKSa5TzywsFXa7SVuV6YMJTxgdLfp0FoPc5whCgIS9ADajFgtZEfuoW+yIL61xDgB6HEYwVpuz7HIjdgDodRhxdiZQphMVJnuOHRAiyfkyCM5nnzyCj3x/P5IV9EefKZPrHCD4xTtMOsmp+NqJ2MPoshvhtMhr/iOIRoeEvQC5bGVFK9lSG+cAwKBVo9NqqDlhj8aTCMUSsoX94l47xjyhdNd0JfFH49BrVEvG88rR1BWOJXB0zIsJbxgHR9yKXrsQ0wui65zyXfGAsJddaip+eiG8pNYvZEYqb9874Q2h02ZAq1ls/qu9XhWCqAYk7AXIKeziqJsCETtQmyNvcs1pRAb6HQCAweHKCZ6IPxxfkoYHyrMI5vi4F7GEIGDPHJlQ9NqFyOxELwdSN7wFInEEoollEXul7Xs555jwhtFtM6Qjdhp5IwgBEvYC5BL2U1N+dNkMsBpK64gX6Xeaaq7GLu6hlzPuBgBbuq3QqVUYHPaU41gFCUTiS9LwQGoRjMLd2oNDws820GfHs0cnKpaOT0fs1vIIu8sibRHMol/9YuagGva97mAMkXgSXTaj7Dl8gmh0SNgLkFPYp32KpOFF+pwmTPsiCNXQspKVRux6jRoXrbJicKgKEXsksUzY7SYtfJG4ZKtUKRwYcqPPacKdb16NqYVIxbITM/4I9Bqh078cuCSubs2VOaiGfe94yi+hy2aAQauGWaemWXaCSEHCXgCbSQtPhrCLHfEbSnScy6TPJXTGj7prJ2oXhV18wpbDQJ8DR8a8iMaVE1MpBCLxZaInzmcrJTiccwwOuzHQZ8fVmzqg06gqlo6fXgij3SoYwpQDl1mX7rwvfI7lG+acVVjdKloXi1bGwoY6qrETBEDCXhCbUYtoPIlwTIimR9xBhGNJxerrwOIsey2l41casQNCnT0aT+L4eGXHwQLReHqGXURp97kxTwjTvggG+h2w6DW46sK2iqXjZ/wRtCm8/CUTsU5dTJxzdedXw75X9EsQHQ9XslOeIBofqke4AAAgAElEQVQVEvYCZLvPKWUlm0lfDa5vFTfarUTYL0k30FW2zu4P56qxKys44s800Cf8jO/c2oVpXwT7K1B6mF4oj0+8iFinLrZMZdoXgVrF4MzI5lTDvnfCG4ZGxdKb7lxmHWYpFU8QAEjYC7Jc2JVZ/pKJ06yDWaeuLWEPCTPhovueHDqsBqyyGyveGe+P5OqKVzZFPDjkhlGrxsZO4YXd1Zs6oNeo8OzR8qfjs2fHlUb0iy8W9c74Imi16KBSLZYE7MbK2/dOeMPosBqgTp3DZdFhnlLxBAGggsLOGPsXxtg0Y+xYxm27GWNjjLFDqbd3Vuo8UsgW9lNTPnTbDGhRqCMeSK1vdZlrStg9oaiwtUviytZsdvTZcbDCDXS5uuJF5zylIsnBYTe299rSvxchHd+OZ49OIFHGdHwknoAnGCuLT7yIU6Kt7HSODXMatQothsra9457Quk0PCDslCe/eIIQqGTEvgfAO3Lc/g+c84tTb89W8DxFSQt7cDEVr2QaXqTPaawpYV+J61wmA30OjHvDFfONTyY5AtFEgYi99EgyHEvg1fGFdBpe5IZtqXT8+fmSv0c+xBRzOSP2tMmLhIg91zkq7T43uRBON84BgsFOLMGxECa/eIKomLBzzl8EUL5nvzIgipsnFEMiyXFmxq+IlWw2fSmTmkpalBZiISRvZWs2aaOaocrU2YOp5sZsYTfp1NCpVYoIzpFRL+JJvkzY37axHQatCnvLmI5P7z8v0ww7AFiNwl71YiNv03k2zAlmQJVJxYvmNEsjdjKpIQiRWqixf4IxdiSVqncUv3vlEFO53lAMw/NBROLJskTs/S4zIvEkjlZ4sUg+vKFYum66EjZ3WaHXqCpWZxdXtman4oVd4VpFTGrEn2VHn33J7eZ0On6ybOl40RSmzVK+5jnGGJxmXUFhTCQ55gO5hb2SO9nnA1FE48klwr64U57q7ARRbWH/ZwDrAFwMYALA1/LdkTH2UcbYfsbY/pmZmYocTqyle0OxdOOcEutas7lxWxfaWvT4n08cqfj8dy5KTcXrNCps67FVTtizVrZmolSK+MCQG6tdprSAZHLDti7M+iPYV6Z0/LSvvK5zIs4ineVz/giSPHdJoBz2vfkQNwh22RZT8eQ+RxCLVFXYOedTnPME5zwJ4FsALi1w30c55zs55zvb2toqcj61iqHFoMFCKIbT06lRNwU74kXsJh2+9O6teG3Sh//761OKX18upQo7INTZj48tpD0Aykkgx8pWEbsCgsM5x8Fhd7rEkE06HV8ms5oZXwSMLYpXuWi16At2losvMNpyjN3Zy2Dfm49M1zkRqc1/BNEMVFXYGWNdGZ++G8CxfPetFqKt7MkpH1bZjcvSvUpx7eYOvGdgFf7fb87g6Gh1U/KeYAw2mbvYs9nR50A0URmjmlwrW0WEiL00wRmZD2HWH11WXxcx6TS4emMHnjtWnu74aV8ELrNuxVMKUimWik+XBPI0zylt35uPyQXRdS5XjZ1S8QRRyXG3/wDwBwAXMsZGGWN3A/gKY+woY+wIgKsAfKpS55HKorD7FXWcy8X/unELWi06fObxw4jEq+MdH44lEIknS4/Y+4VadCUa6HwFInaHWVty7VcsKeQTdkBMx0fx8rm5kr5XLoTZ8fKm4YHUhrcCEe/iApgcwm4WRwvLH7WPe8LQqlm6kx8QViBb9BoyqSEIVLYr/v2c8y7OuZZz3sM5/w7n/EOc862c822c85s555XbgykRm1GLuUA01RGvfH19yfcyafHQe7bi9SkfHv5VdVLyC6mZ/VK64gFh+1evszJGNYVT8bqSd4UPDrth1qlxYWf+f/+rLmyHUasuSzp+xhdGu7V8jXMirRYh6s73onI6h52sSCXd5ya9IXRYDUtMcgDRpIaEnSCq3TxX89iMWrw2sYBoPFmW+no2b9vYgV2X9OCRF87i8Ejl15+mF8CUKOyAEOEODrvLbhpSOBWvRTzJ01H9Sjgw5MbFffa0y1kujDo13rapHT87Pom4wunoGV95feJFnObC7nMzvgisBg0M2lxNipXzix/3htGd0TgnIuyUp1Q8QZCwF8Fu0iKS6lQvd8Qucv+Nm9Fm0eMzjx+uSPNZJqUsgMlmoM+BqYUIxjzlNarxR3LPsQOLJjUrbewKRuN4bdJXMA0vcuNWIR3/yjnluuM555jxR8reEQ8IES+QvwFtOo85DaC8fW8hJrwhdNqWZzBcZj01zxEESNiLkpmSvqACETsgiOqXb92KU9N+fKPCKXmlhR0o/0KYQCQOFQMM2uV/zqUKzuERLxI5jGly8dYL22HSqfGMgmY1nmAMsQSvSMRebGRsJoedrIjS9r35SCY5pryRJY1zIq4izX8E0SyQsBdBFLgeR/k64nPx1gvbcfvOXnzzhTM4WMGFKp4SNrtls7GrBQatCoNl9o0XF8Dk2lUuNnWtVNjzGdPkwqhT4+pNHXj+mHLp+ErNsAOZi2Byp7MLReyLO9nLm4qfC0QRTSTRlaPnQKyxk1880eyQsBdBFLhKpeEz+fyNm9BpNVQ0Ja9kxK5Vq7C9x172Fya5NruJLDZ1rUxwBofcWNtmTl+nGDds7cR8IIo/nlUmHb/oOleZrnggdyqec56K2HOfw6hVQ6dRxr63EJOiOY09d409nuRYCJFfPNHckLAXQRS4co+65cJq0OLLt27DmZkA/uEXJyvyPb0KdcWLDPQ7cHy8vEY1uTa7iZSSiuec4+CIB5dISMOLiOl4pbzjxU70SnTFWw0aaNUsZyreH4kjFEvkjdgZY3AoZN9biPHUYqFczXNij8AsNdARTQ4JexHSEXt75SN2ALhyQxvef2kfHv3tWRwZXXmt+o9n5/DQsyeK3s8biqHFoCnYAS6HgT4H4kmOIxJMd3zhGP768cM4PxuQ9T38BYTdZtSCsZWliM/PBTEfiOZ1nMuFQavGNZs68PyxCUXS8YVMYZRG9IvP5bc+I6EkUIkNbxOpRsx8zXMALYIhCBL2ImzvteNdF3fjLRdWxsY2F59750aoGcNzxyZXfI1//cMQvvni2bTnfT4WFLCTzUSsTUuZZ/+7vSfw+IFRvHhK3i6AQIFUvFrFYDWszKRG7A2Q0jiXyRXrW+EOxhSZBpj2RWDSqfP+fEoj7jXPdQ6g8CIau0lbdoOaiYUwdGpVTnvdxVICRexEc0PCXgSrQYtvvG9HRZy/8tFi0GJLt7WkJjRRWJ8pYqCihE98Jq0WPfpdpqJn/83r0/jhvhEA8v2+C9XYgZWvFB0cdqNFr5HtX9DvNAEAhueDsr9nNvn2n5eLVosuZypeasQ+X/aIPYwOm36ZOQ2A9P9RWgRDNDsk7HXCjj4HDo96VuTFPe4JYcIbBmPA3iPjBbuGlRZ2ALikz4HBYU/e7+sNxXDvE0exvt0Cq0Ej22QkEEkUnFhY6UpR0Zgml4gUos8lCPvQXOnCPu0L521YKwf5bGWnJTTxVWJ166Q3vGSrWybiBATNshPNDgl7nTDQ70A4lsRrE4VT6bkQo/X3XtKLMzMBnJzy572vpwzCvqPfgVl/BCPzuVPTDz7zKmb8EXztvdvR1pI7FVwIIWJf7oYmspKVov5IHCenpBnTZNPRYoBOo8JIHUbsrjyp+BlfBDq1Kj2vngtHKhVfznGzcW8I3Tnq6wCg16jRYtBQjZ1oekjY64RL+kWzF/np+MEhDwxaFT517QaoUlF7PsoRsQ8UqLP/+rUpPH5gFP/9Leuwrccu2z2Mc16wKx5INXXJ7NY+POJBkkNW45yISsXQ6zAqkoqfLmAKUw5cFh38kfiyKYZpXxhtLfqcXgEiDpOuZPveQiSTHFMLYXTmidgBwaSGUvFEs0PCXid02wzosOpXJuzDbmzrsaPTZsCb1rrwzNGJgmlxpYX9wo4WmHTqZWf3BmP47JNHsbGzBX9x9QUABGGR88QciScRT/KiqXi5EbvYE3Bxb3Fjmlz0OU0lp+LDsQR84XiFI3Zx/enS39eML4LWIudIu8+VaeRtNhBBLMHRncN1TsRl0VPzHNH0kLDXCYyx9FIVOYRjCRwf96ZTyjds68LZmQBem1ye0g/HEojGkyXvYs9GkzKqyT773z5zHLP+KL5623boNUIqvdhO8Gz8qeiwxVC4eS4YTchahXtg2I317ZYVv8jpc5owMh8sKS1dyVE3kXwmNYXMaUTK7Rc/4RFm+jsLzPTL/fshiEZEsrAzxvpYjjwcE+hT9lhELgb6HBiZD6VNS6RwbMyLWIKn0+Hv2NKZSscv745X0nUum0v6HTgx4UMwKgjxL1+dwpODY/jzqy7ARats6fu5LHq4g1EkktIEMb3ZTVcgYjfLc59LJjkODntWVF8X6XOZ4YvESxr/mq6CsLvSneVLo14ptf5S7XuLMZFynevO4Ton0mrR0U52oumRE7GfA5BrmNuZ+hpRZgb6U7XqIelGNWKULNaKXRY9Llvnwt4c6fhyCvtAvx2JJMfhES88wSg+++Oj2NRlxSeuumDJ/VxmHTiXLg7+AitbRRZXikq75tnZALyhWLqvYSX0pUbehkqos8+IrnNVSMVnRuyxRBJzgWjRc5Rq31uMiZTrXFee5jlAiNjdwSiSEl8YEkQjIkfYGYBc/1ssAKSHkMSK2dJtg06tkuW9PjjkQZ/TtGQO/4at3Tg3G8CJrA57JRfAZLOjd7H5b/fTx+EORPHV27ZBp1n6J1hsdWg2gQIrW0WcYopYYu138cXQyurrwKKwl9JAV42I3WlZXmMX/y2KRuzlTsV7w9BpVOlyQS6cZj0SSZ5+kUoQzUhROyvG2MOpDzmAhxhjmc9UagCXAjhUhrMRWRi0amxZZZVcZ+ec48CwG1dc0Lrk9rdv6cD9Tx3D3qPj2NxtTd9ezojdYdZhbasZ//qHIUwuhPGpazZgS7dt2f3SNd5ABEBxG990Kr7AuNtiJClNcA4Ou2E1aLC2deX7AURhL2XkbcYXgYotWqVWgha9Bjq1aonfetqvvkh3fin2vVKY8IbRZTMU7MxvtSyunnUUeAFAEI2MlIh9a+qNAdiU8flWABcAGARwV5nOR2Qx0OfA4VEvovHiRjWj7hBmfJF0fV3EZdHjzetc2HtkaTq+nMIOCCY7kwthbOm24n9ctS7nfdLuYRIjdp+U5rl07Vea4BwYcmNHn0O2MU0mRp0abS16DM3J873PZHohApdFr5hvvxREv/j5jN+/1Ca+Uux7pTDhCRVsnAPIVpYgAAnCzjm/inN+FYDvAbhe/Dz19nbO+cc456fKf1QCEIQ9Gk/i1YmFovdd3CW+vFZ8w9YunJ8L4vj44nVEYbcbyxPpXLmhFUatGl+9bTu06tx/es4841b5CEiqsUtPES+EYzg17S+pcU6kz2kqKRU/4y/eiV4OxL3mIumd8BLOItS4yxexF2qcA2gRDEEAMmrsnPM/5ZwXVxOirCw20BVPxx8c9sCkU2Nj5/KU9nVbOqFWsSXrRb2hGBgrHP2Wws3bu3HwgWuxqcua9z4Okw6MSY+4pAi7QauGQauSFEkeGvaAc5TUOCfS7zTldduTgmgKU2mcZh1mA8sjdin7EoRFMMqLaiJlTlOocQ7IXN1Kwk40L3LG3QyMsf/JGPs5Y+wQY+xI5ls5D0ks0mUzottmkFRnF4xpbNDkiI6dZt2ydPxCKIYWvaakFHQhGGMwaPPXwgEhneswSTep8UsYdwPElaLFI8nBYTcYA7b3Lq//y6XXacK4NyRrfj4TKbPj5aDVosd8Vo3dYdIua3TMRblWt876I4gneVFhF7Mz8zTyRjQxcrri/wnAvQDOA/gJgCey3ogKsaPfgYPDhUfeQtEEXh1fKBh53ritC8PzQRwbExIxnmBUcXOalSDHZCQQicOoVRetQ0tdUDI47MGFHS1oMZT+e+hzmsA5MOaWH7Unkhyz/mjVIvbMHofpBel+9XaTVrZ9rxTEGfZ8C2BEdBoVrAbNkhcmBNFsyMm53gLgNs75L8t1GEIaA30O7D0ygUlvGJ15Ipgjox7Ek7xgrfi6zZ34/I+PYe/RCWztsZXFTnYluPJsGMuFPxKHRULpQMrqVsGYxo0bt3VL+t7F6HctjrytbZPXYS+a9FTSJ17EZdEhGE0gFE3AqFOnav3SzlGuiH0itds+3997Jq0WPaXiiaZGTsQeBDBSroMQ0im0VEVkMBXR52qcE3GYdbj8glbsPSqscq0ZYbfoJK9u9UcSBWfYRaQIzukZP3zhuCL1daC0WfbpBekNa0rjWjJyKC9iX4l9rxSkuM6JZHf1E0SzIUfYvwLg04wx8pevMlu6bdBpVAUb6AaH3VjTai5o5gEI3fEj8yEcHfPCG4qVrSNeDi6zXnKNXdjsVrhuD4gp4sLXFH+f2eOBK6WtRQ+DVoXhFSyDmfFX3pxGJLOznHMuqzu/XO5zE94Q9BpV2kWwEE6z9BeGBNGIyBHpawHcDuAcY+w5xtjTmW9lOh+RA51GhW2rbHkjds6FlPIOCQJ13ZYOaFQMe49MwBuKw1oDEbvTrIMnGEM8UXxW3x+JF22cA4SI3RuKFbQaHRx2w2HSYk2rWdZ588EYW/HI2/SCNFOYcuDMcP9bCMURjSdlROzlcZ8bl2BOI+Ky5N4pTxDNghxhnwXwYwC/BjAJYC7rjaggA/0OHBtbyJnyHJkPYdYflZRStpt0uGJ9K545MoGFGknFi+5h8xLEIRCJS0rF201aJLkwp56PwWEPdvQ5JImHVFYq7NWN2Bfd22b8YVnnSPvyK9xAN+kNF22cE3Glmi/JL55oViQ3z3HO/7ScByHkMdBnx6MvJnFsbHnn+4Hh+dR9pNWKb9jahb9+XZhYrAVhd2akgotFrNKb58RIMpZOF2fiCUZxetqPWy5WpnFOpNdpwu/PzIFzLusFw/RCBC16DYy64mUGpUlvePNHMmr90jIHcu17pTLhCeFNa12S7uuy6JDkgCcUK1qKIohGhOrldYoo2rkWwgwOeWDRa7Cho7jXOiB0x2vVgujUgrDLWQQj1NglCHuRlaIHR4RmQyUc5zLpc5oQjCYk9wyIzPilN6wpjVmnhk6jwnwgKjtzINe+VwqJJMeUL4Iuu7QXF2QrSzQ7kiN2xthR5N7uBgDgnG9T5ESEJNqtBqyyG3PW2QeH3djea5PsMW4zaXHFBa34r9dnakPYM1LBxfBLTMU7ikSSB4fcUDFge68yjXMi4sjb0FxQknObyIyMTnSlYYyh1SyYBKUjdmv1auwzvggSSS45FZ/eNxCIYr1ipyCI+kFOxP6fWGpI8zSAYQC9qY+JCnNJv2PZbvZgNI7XJn2yI88bUrPbtZC6zEwFFyKeSCIcS0pungPy134Hhz3Y2GmVFP3LYaVb3iYXqmMnK+K06DDnj2DGH4Feo0KLxN+LHPteqYxL2MOeiTPHTnmCaCbk1Nj/NtftjLG/BtCv2IkIyQz02fH04XGMe0Lp+d7DI14kkhwDMmexb7m4GxoVw6VrnOU4qizsRi1UrPgij0BUaByUMu5WKJJMJDkOjXhwyw5l6+sA0OOQP8vuCUYxPB/E7W/oVfw8UnGZhc7y6YUw2q16Wf0BUu17pTIp0XVOxJXeKU+peKI5UaLG/iSADypwHUImongfyJhnF1PzA73yhF2jVuGWHasquiI0H6qUX/xskYjLL2Flq0iLQQMVyz1ffXLKB38krnh9HRAi2E6rAUMyZtnLVe+Xg8ss/P5n/BG0ySghANLte6Uy7pEXsYsv4uT2NRBEo6CEsF8JwZWOqDCbuqwwaFVL6uyDQ26sazPXhOd7KQirQwtHXFI2u4moVAz2PO5z6RdDZRLSPqdJVip+cMgNtYopsohmpYirW6cXpNvJikix75XDpDcMg1YFu8S/aa1aBZtRS6l4ommR0zyXXUdnALoA7ACQM01PlBetWoVtq+xp+1jOOQ6OeHD1xvYqn6x0pCyC8csQdkBcKbpccAaHPHCZdelGN6XpdZrw0ulZyfcfHHZjY2cLTBJ6B8qF06xHKJbAiDsoecxMxGHW4cSEchueJ7xhdNuMssoB2TvlCaKZkBOxZxvSTAP4JYDrOedfKMPZCAkM9Dvw6rgX4VgC5+eCmA9EZdfXaxGXRV804hIjdild8UB+v3jBpU9ZY5pM+l0mTC6EEY4V909PJDkODXuqmoYHFicTwrGkbL96R54XUCtl3BuStPwlE6GUQDV2ojkhg5o6Z6DPjkcSHMfGvOk6rlJLTKqJy1x8J3tA4i52EYdJizFPeMlt7kAUZ2cD2LWzZ2UHlYDYGT/qDuKC9sLeAienfAhEExjoV3bsTi5iAxog3/3OkaqxJ5McKgV6Nia9Ybx5Xausx7jMepyZ8Zf8vQmiHpFdY2eMrWWM3cgYu4ExtrYchyKkk9lAd2DYjRaDBhfIXBFai7jMenhDMcQK+MX7wtKb54DcTV0HR8pbXweEVDwgrTNebIS8pK+60wmZY49SZ9hF7CbB+U389ymFeCKJqYWw5MY5Eael+AtDgmhUJAs7Y8zKGHscwGkAPwHwFIBTjLHHGGPSLM4IxWm16NHnNGFw2I3BITcu7rUrEiVVG3ERSaGNbHKa5wCxqWvp9Q6kGtW29ZSvUS3TpKYYg8NutFp06HVKG+0qF5lmOm0W+c1zgDSv/2LM+CNIckh2nRNpNevSO+0JotmQE7F/A8A2AFcBMKberk7d9o/KH42QykCfHa+cm8fJKfnGNLVKaypiLDTyJmeOHRAiyXAsiVB0sdY9OOTB5i5rWRvVXGYdTDq1pIj9YBkW0ayEUiJ2Jd3nxlOlk26JM+wiTrMOnCu/ZY4g6gE5wn4zgA9zzl/gnMdSb78B8FEAt5TldIQkBvodcAdjSHI0ROMcsCgshTqb/ZE4tGoGvUaasGcLTjyRxOFRj2L71/Mhrm8tNvI2H4ji3GygJl6cmXSCgxxji410UhHH0pSYZZ9Iuc7Jbp6zLC4SIohmQ46wG5F7Pes8gMovjSbSiELAGHCxwl7n1SJtK1tgll3qAhiR9ErRlOC8PuVDMJqoyIuhPqepaCr+YHqevvr/howxuMx6uMw6aNTyWnGK2ffKQXSdkxuxu8hWlmhi5OQfXwLwRcbYhzjnQQBgjJkhzLD/vhyHI6QhzDyr0eMw1sQSFyWQ8sTsD0tbACOyuFJUEBxx/r8SEXKf04QXTs4UXN96YMgNjYphW0/1hR0QOuOj8fzNi/mQk4r/59+cweERT96vn5z2wahVw2qUVyqR8sKQIBoVOf9bPg3geQBjjLEjEDa9bYfgOnddGc5GSESjVuHDf7IWHTJrobWMzaiFWsUKPjFL3ewmkr26dXDIjVaLHj2O8jeq9blMiMSTmPFF0G7NneAaHHZjc7e1KjvYc/HuHatW1HxWyL43k2lfGF/52WvoaDHkfUGqValw284e2T0HUko5BNGoyJljP8oYuwDAHQA2QnCe+zcAP+Cch8p0PkIin752Q7WPoCiiX3yhJ+ZAVF4q3pmOJMWI3Y1L+u0VaVQTZ9mH5oM5hT2eSOLwiLeqi1+y+dPL16zocYXsezP52bFJcA58/+5LsaFD2cEasexSbN8AQTQicixl/w7ACOf8kazbP84YW8U5v1/x0xFNjbiIJB/+SEJW6SGdig9EMeuPYGguiA9c2lfyOaUgCvvwXBBvWL18Rv21SR9CsQR21EB9XQny2fdm8syRCaxvtygu6oCQxXKYtLThjWhK5HTFfAjAwRy3DwL4b8ochyAWKeb3HYjEYZE46gYAOo0KZp0a7mAMB8X6eoWmCFY5jGAsv0nNwTIvoqk0+ex7RaYXwnjl/Dxu2NZVtjM4zTpqniOaEjnC3g5gJsftswA6lDkOQSxSbBGM3OY5YNF9bnBYaFTbuqoyG9T0GjW6bca8wn5gyI32lsrU+ytBsQ1vz6XS8DdsLZ+wuyx6cp8jmhI5wj4M4E9y3H4lgFFljkMQi7Ra9AUXecgddwOEBjp3MIoDQ25s6bbCoK1co1qvM7+wD6YWv1TbmEYpiu1k33tkAhd2tGB9GdLwIi6zDnO0CIZoQuQI+zcB/ANj7COMsXWpt48C+BqAR8tzPKKZcZp18IXjOUeuOOcIROVH7A6TULc/MuqpuJlPn9OUU9hn/REMzwervvhFSXLZ94pMLYSxb2ge7yxjtA7Q6laieZHTFf81xlgrgIcBiFZUUQDf4Jx/pRyHI5obccPYfCC6zHksFEsgyaX7xIvYTTr8/swcEkle8Xp2v8uMGd8ogtH4EgvbwaHGqq8Di/a94VhiWVbkuaMTQhp+W2dZz+A06+EOxhBPJGWb7BBEPSPrr51z/lkArQDeBOAyAG2c83ulPJYx9i+MsWnG2LGM25yMsV8wxk6l3jfOMxtRMmmTmhydzX6ZC2BEHCZteja70hG7uOVtZH7pdOjgsAdaNcNFFar3VwJxjjxX1L736AQ2drYUXWFbKq2WpeONBNEsyH4ZyzkPcM73cc5f4ZzLWXi8B8A7sm67F8CvOOfrAfwq9TlBAMhwD8vR2ewXV7auIGIHgA6rHt0y/cdLpS/P+tbBITe2dNsqWu8vN2n73ixb2UlvGPvOu8vaNCdCJjVEs1Kx/BTn/EUIvvKZvAvA91Iffw+0TIbIoNATcyAibnaTH7EDqEqjWr9oUjMXSN8WSyRxZMzTUGl4INO+d+m/3bNHJwAA7yzjmJuIM21LTA10RHNR7cJTB+d8AgBS79urfB6ihmg1CxF7rs74xVS8vChX9DG/pApb8OwmLVr0miVb3k5MLCAcSzZU4xyw+HvO3skupuHXtVnKfobWtF88RexEc1FtYZcMY+yjjLH9jLH9MzO5xumJRsNq1ECjYnkidkHY5XbFr241Q8WAN69rVeSMcmCMoTerM74RG+eAzE16i6n4cU8IB4bcuLEC0TpAETvRvFRb2KcYY10AkHo/ne+OnPNHOec7Oec729raKnZAonowxgjUaBsAACAASURBVODI4x4WiK6see7iXjsO3HctNndbFTmjXPpdJgxlCvuwB51WA7rtjWFMI5Jp3yvy3LFJACj7mJuIw6QDY1RjJ5qPagv70wDuTH18J4CnqngWogZxmXU5U6m+FTbPAYDDrCt+pzLR5zRhdD6EZKoz/8CQuyplgXKTad8rsvfIODZ3WbG2Aml4AFCnFgnNkrATTUbFhJ0x9h8A/gDgQsbYKGPsbgBfBnAtY+wUgGtTnxNEGpdFl3PcLbDCcbdq0+s0IZpIYsoXxvRCGGOeUMMsfskm031uzBPC4LCnrN7wuXCZdZgnv3iiyajYsyLn/P15vnR1pc5A1B8usx6jbs+y2wOROBgDTDWyu1wq/S6xMz6YFr1Kz9NXCtG+FxBMaYDKpeFFnObcLwwJopGpdiqeIArizBNx+SMJmHWauvNWz5xlHxz2QKdWYUuV6v3lRtjwJqTinzkygS3dVqxpNVf0DK20CIZoQkjYiZqm1aKDLxJHJJ5YcruwAKa+onUA6LYboVYxjMwHMTjkxkWrrNBr6u/nkIKYih91B3FopPJpeKD4hkCCaERI2ImaxpmaZc9+cvZH5C+AqQW0ahW67QacnvbjyJi3IRvnRMTVrc8dFbrhK+E2l43LooMnGEMssXyREEE0KiTsRE0jLoLJHnmrV2EHhHT8CydnEI0nG25+PRO7SYeFcAxPHx7HRaus6HdVNg0PLO4byLdpjiAaERJ2oqZZXASz9Il5JbvYa4U+pwnBqFBaaNTGOUCI2DkHjo55ccPW7qqcQcz45PJCIIhGhYSdqGkWF8Es7Wz217WwC5HrKrsRHdbKLqKpJKKtLFCdNDywdPUvQTQLJOxETZNvEUwgWt+peAANO78uYk/Zym7rsaEvNeZXacSMT659AwTRqJCwEzWN1aCBVs0wm11jD9evsIuz7I1cXwcEDwKgzLPrE0eA/7MTCMzlPkMq43NuNpDz61L495eHcdP/+R3cFPUTdQIJO1HTMMZSI0tLI65AJFG3qfgt3VY8eMtFeO8beqt9lLKypduKL75rC+54U3/5vsn53wFzp4CJQzm/7DTrcOWGNnzzhbMYngvmvE8hTk/7sfunx3F0zIvdPz1e6mkJoiKQsBM1j9OsX9L8FI0nEU0kYanDOXZAeLFyx5v66zbjIBWViuFDl60u78/pPr/0fQ6+/J6t0KgY/vo/D6c9+qWQSHJ85vHDMOnUuOvNq/HUoXE8f2yitPMSRAUgYSdqnlbL0kUw9eoTT5QBCcLebTfi/hs34+Vz8/j+H/LfL5tv/fYsDo148IV3XYTP37AJF62y4r6fHKNGPKLmIWEnap5sv28/CTsh4j639H0ebtvZg7de2Ia/f/51nJdQbz815cPXf34S79jSiZu2dUGrVuGrt22HNxTDA08dU+LkBFE2SNiJmsdl1i/xixeFfSUrW4kGIpkE3EPCxwUidkAofzz0nq3QqIun5OOJJD7z+GFYDBo8+O6L0vsINnZa8clrNuCZIxPYe4RS8kTtQsJO1Dwuiw6BaALhmGDqQql4AgDgnwQSEUBrEgSeF66fd9mMeODGzdh33o3v/v583vt988WzODzqxRffdRFaU131Ih+7ci229dhw/1PHaISOqFlI2ImaJ9t9jlLxBABgPpV+778ciCwAwfmiD9l1SQ/etrEd//tnr+HsjH/Z11+f9OEbvzyFG7Z25Vxao1Gr8LXbtsMfjuP+nxwDL/JigiCqAQk7UfOkTWpS6fhARIjcG72rnCiCmH5fd9XSzwsgpuR1ahX++j+PIJGRko+lUvAtBg2+8K4tea+xvqMFn7p2A547NolnKCVP1CAk7ETNI5qMzKYa6BZT8fU57kYohPs8wFTA6itSnxduoBPpsBqw++YtODDkxndfWnzMN184g6NjXjx4y0Xpv7l8fORP1mB7rx33P3UM077wSn8CgigLJOxEzePKith96eY5bdXORNQA7nOArQdwXbD4uUTevWMVrtnUgf/9s9dxZsaPExML+MavTuGm7d24XoJTnpCS34ZgNIH7fkwpeaK2IGEnap706laK2IlM3OcBxxpAZwYsHZJS8SKMMXzpPRfBqFPjM48fxmcePwybUYcv3Jw/BZ/NBe0t+Mx1G/DzV6fw9OFx+ecniDJBwk7UPBa9Bjq1Kt08F4jEodeooFHTn29T4z4POFYLHztWL46+SaS9xYC/vXkLDg57cHx8AX/37ovgMOuKPzCDu69Yi4E+Ox546jhmfNQlT9QG9MxI1DyiX7xoK+uP1O8CGEIhIn4gMLNU2Oelp+JFbt7ejbvevBofu3It3r6lU/bj1SqGB2/ZCm8ohl+emJL9eIIoByTsRF3gsujSVp6BOt7FTiiEmHbPFPaFMSAuL2pmjGH3zVvw2XduWvFRNnW1wGnWYXDIveJrEISSkLATdYEQsQtP2hSxE2lhd64R3jvWAOCAZ6TiR2GMYaDPjsFhEnaiNiBhJ+qCVot+iUENCXuTkytiz7y9wuzoc+DMTACeIC2IIaoPCTtRF2TW2IVd7NQR39S4zwEGG2B0CJ+nhV1+nV0JBvqEcxwc9lTl+xNEJiTsRF3gsugQiiUQiiaoxk4sjrqJtHQCGkPVIvbtvTaoVYzS8URNQMJO1AWLfvERSsUTS0fdAICx1Mjb+aocx6TTYFNXCwk7UROQsBN1gcssWHzO+aMk7M1OMiHMrGcKO7DikTelGOhz4NCwZ4n/PEFUAxJ2oi5wptznZv0RBKMJSsU3MwvjQDK22BEv4lgjROxVsncd6HMgEE3g9UlfVb4/QYiQsBN1QWsqYh+ZDwKgzW5NTXZHvIhjNRALAIHZCh9IQGygo3Q8UW1I2Im6QIzYh1LCThF7EyN2vucS9syvV5hepxGtFj0JO1F1SNibhVioKuYdSmHWqaHTqDA8Jwo7jbs1Le7zgEoDWHuW3i6m5uU20AXmgOB8yccSjWrqdeQtHEtg3BOq9jEIBSBhbxZ++bfAI1cIjUd1CGMMrWZdOmJvMVDE3rS4zwO2XkCd9Tdg71v8uhx+dAfwk/+uxMkw0O/AudlA2v64nnj0xbO47h9eRDSerPZRiBIhYW8Gkgng2BNA2FO1cSAlcFp06Rq7WUfC3rTMn1uehgcArRFo6ZLXGR8LAaOvAFOvKnK0dJ29Dn3jj4x64I/EcW42UO2jECVCwt4MDP0eCEwLH0+fqO5ZSsBl1iOSiiaoxt7EuM8v74gXETvjpTJ+EEjGgYVRIF56lL2txwZNnRrVnJr2p95TV3+9Q8LeDBz/MaAxCh/PvFbds5SAK2NXNnXFNylhLxCazx2xA/JNakZeEd7zJOAtvQfFoFVjS7e17oQ9FE1gOJUNOznlr/JpiFIhYW90EnHgxNPAhe8AbH31LeyWRWGniL1JyTfqJuJYDfjGhRS7FEb3AWCpayvTTb+jz4HDI17EE/VTqz4z40+P/5+aooi93iFhb3SGXgICM8DmW4C2C4Hp+hV2Z2qWHaDmuaYlLex5UvFiit4zXPxanAMjLwNrrlx67RIZ6HcgFEvgtToyqjmZEvO1beb0x0T9QsLe6Lz6E0BrAtZfB7RvBGZP1m1nvBixq1UMeg396TYlaWHvz/11Oetb3edTL3pvBtR65YS9zw6gvoxqTk75oVUzXLu5A+fngojE6/M5ghCgZ8dGJhEHXn0a2PAOQGcC2jYBiUhV/bRLQayxm3VqMMaqfBqiKsyfA4xOYWVrLkRhl/I3PrpPeN/7RuGFgkLCvspuRHuLvq46409N+bC21YLNXVYkkpw64+scEvZGZuh3QHAW2PJu4fP2jcL7mfrsjHdZhFQ8Nc41MYU64gHA3AZozdJEeuQVQGcB2jcLqf15CY+RgGBU48BgHRnVnJz2YX2HBevbW4TPqYGuriFhb2SO/1h4klt/rfB564XC+zqts6cjdhL25iV7XWs2cta3jr4CrBoAVOrFxyi0QOaSfgeG54OY8UUUuV45CUbjGJkPYUNHC9a2maFi1EBX75CwNyqJOHDip0I3vDY16qa3CO5cdRuxC8Juoca55iQRF0bSCgk7kBLpIqn4aACYPAb0XCp87lwDRH2KWMsCwEB//dTZT6fm1zd0WGDQqrHaRQ109Q4Je6Ny/rdAcG4xDS/StqluI3aTTgODVkWp+GZlYVQwk8nXES/ilLC+dWwQ4Amhvg4ovkBmS7cNWnV9GNWIaff1HS2p95a0WQ1Rn5CwNyrHfyzUDy+4Zunt7RuBuVNC9FOHuMx6EvZmpdgMu4hjNRAPA/6p/PcZTRnT9Oxcek2FGugEoxobDg5Vps4+64/gPf/0En5/Rv7K2lNTPujUKvQ7TQCADR0tGKrRznhfOIZd//x7HBhSJrPSqJCwNyKJWCoNf/1iGl6kbROQiALzZ6tzthJ54KbN+MiVa6t9DKIayBF2oHBn/Mg+wLUeMDmFz+2p8TkFV74O9DlwZMyDWJmNajjnuP8nxzA47MHeIxOyH39q2o+1bWZo1IIcXNBuQSLJcXam9jrjf/HqFPYPufHr16arfZSahoS9ETn3omC7mZ2GB+q+M/7tWzrTizaIJmP+HKDWAdbuwvdzFFnfyrkQsfdeunibzgRYOhXrjAeEBrpwLIkTEwuKXTMXzxyZwHPHJqHXqFbUiX9yypdOwwNCxC7eXmuIL1yoa78wJOyNyPEfA7oWYN3Vy7/WukF4P/N6Zc9EEKXiPi80f6rUhe9n7wXA8gv7/Fmh/6TnDUtvl+szXwSxge5AGefZZ3wRPPDUMWzvteMjf7IWr08uwB+RXmYLROIYdYewod2Svm1tmxlqFcOpGhNPbyiGF0/NAKCu/WKQsDcaiRjw2jOpNLxh+dd1ZiHtWMdb3ogmpdiom4hGD1hX5RdpcfGL2Dgn4pS5Ga4IXTYjumyGss2zc85x30+OIhBN4Ku7tmHnageSHDgyIv37iR3xmRG7XqNGv8tUcxH7L1+dQizBceWGNgzNBxGO1V4PQK1Awt5onHsBCLlzp+FF2jfV9TIYoklxnyveES/y/9s78/Cqymv/f97MIQMJkIQpkAAJgciMOCs4MAWVOLW9t7e9tV5tHTrcWq2ttda2VmsHrbW2trX2drQ/lRkBRQUHZFYgzHOAkEASQgJk3r8/1tnkJJyT7DMkJ8P6PA9PknP23ufNJues9/2+37VWn0zv++VH1kF0IqTkNH88OQNOH4X64OWeTxyS3G4V6BZ+eozlBcV864ZsstISmJDu6gXvgxPfDt7ZafHNHs9OTeh0zvglW4sYlBTL7ZMGY1lNkxLlQjSwdzcK5smH1vBrvR+TkgMn98jqXlG6AufKpWWrkxU7tF4itnA9DJoEYS0+/pIzAMtZAxmHTBiSxNFT5yg+XR20awKUnK7msQUFTBiSxF1XiZm0d69IRqTG+yT97ympIioijKF945o9np0Wz6HSM51mVVxxto7395wgb+wARvYXdUH7xntHA3t3oqEOdiyGkbM9y/A2qaOgsa7LOuOVHohTR7xNcoaku9Webf54TSWUFDQ3zp0/pw3TnR9MGupaRQdx1W5ZFt+dt43qugZ+fvs4wsOa+iZMGpLM5sJTWA4r6O0urmR4Snyza4BI842WtHPtDKzYfpy6BovZYwaQ0TeOiDCjBrpW0MDendi/CqpPQe7c1o+zJUjdZ1e6CnbqmuPA7iVIH90IVmNTxblm52Q0f60gkDuwN1ERYUEtVDP/k6O8vaOYB6ePZHhKcwl94tAkTp2tY7/DJi57iqsukOGhyRnfWQx0tgw/brDcz8x+cWqga4VOEdiNMQeNMVuNMZ8YYzaEejxdFicyPLic8Ub32ZWug88rdi+BvdDV0c0uTONOfKq0OA7iij0qIowxg3oHzUBXfLqaHywoYNLQZO688kK/gZ0K6kQhqKqp5+ipc+eDuDuZ/WRV3Bnk7oqzdXyw5yRzxg4439UxO63zeQA6E50isLuYZlnWeMuyPLzjlDapr4WdiyAnT1zBrRHVSz4gdcWudBXKD0rntugLV5ce8VZJ7sg6Uaxiky48x5cGMj4wcUgSW49WUFsfWKEay7J45I2t1DY08sxtYy+QzwGGp8STGBPhaCJhr3izUi+8p1ERYWT0i+sUcvfy7cepb7TIGzvg/GNZafEcLjvLudrO4QHobHSmwK4Ewv73xFzUmhveHXXGK12J8gPOV+sgFeWiEpo74y1LerC3zF93x0kDGR+ZOCSZ2vpGCo5VBHSd1zcd5Z2dJXx7Rg7DUjxPcMLCDBMcOvH3FF+Y6uZOVmp8p5C7l2wpIr1PLGMG9T7/WHZaAlYn8gB0NjpL0W0LWGGMsYDfW5b1UqgHdJ76Glj/R5h854XlWZ1yZANUFsGoG/07v7EB1r0EFUe8H3PoI4juDcOmObtmSg7sWSEr/YgoZ+dsfEXc9K0x9ArIme3seop3diySYiwDxoV6JJ2D8oOQfqnz442BPhnNV9+le8Vd78k4Z5OcIZNky5JrBIGJLgPdr97ew0gP+9lOsCx4dUMhUzL68KXLMzwfsOFlyMlj4pBknl25m9PVdSTGRHq95p6SSqIjwhjiqhHfkqy0BJYVHKe6roGYyDaKArUT5Wdq+XDvSb58VeZ5GR6a0vN2F1dykVvAV4TOEtivsCzrmDEmFXjLGLPTsqzV7gcYY+4G7gYYMmRIx41s2+uw/LsiA469w79rvPUDOLwG7npL0mx8Zc0L8Nb3Zf+PVj5sLrnbeZBOyZFOWWX7ZPXeFqcOw6KvQ3g0hHn5s2mohc1/hQf3Oh+HciGH1sCr/wVx/eDetRDXN9QjCi0NdTKpHZvh23nJGc0rLNqFaTwZ586fkwl1Z+HMCdlzDwJpiTFcNqwvGw6WseGg/81LUhOi+dltYwnzIMFzYhcs+V84uZuJIx7EsuDTwlNclZXi9Xq7i6s8OuJtstPiz+eLhyp4rnDJ8HPGNC8jPLRvHJHh6oz3RqcI7JZlHXN9LTHGzAOmAKtbHPMS8BLA5MmTneVyBIOC+fK1cJ1/gb2hDo65WkTOvxfuXtV6KlpLTuyCd34MOXPgM38L2irifM34kh3OAvv2BfL1vo+hj5cmLLvehH9+VorkZN0QnHH2NGrPwIJ7IWGABJc3vw23vRzqUYWWU4fFye6LFA9y/O4V0NgoOeuFayGmd1NZZW/ngKz0gxTYAf55tw9qgz/Y3eq2L2D81CcwBjYdaj2w7ymuZEpmH6/Pn3fGl4RuVbxk63GG9OnFRYMSmz0eGa7O+NYI+R67MSbOGJNgfw9MB7aFdlQuzpXDvnfk+8K1/l2jeJusAMZ/Xva0Vz3l/NyGepj/VSkDO+dXwQvqIB9uJsz5PnvBPJGFvQV1EDd+dKIcq/jHyiekvsAtL8HUh0UxsidVPRVbTu/jsOqcTXImNNTINhg07a+3LEzT7JwM+RrElLcOwf58qiwioWQT2akJbGwlxa6yuo5jFdVe99eBkOeL2zJ8npsb3p2stAR2dwLXfmck5IEdSAM+MMZ8CqwDlliWtSzEYxJ2LpVCLsOvg+ICWU35ip1eM/U7MPEL8OFzsufuhDXPS97t7GeCunoAxC+QnOEssJcfknG0ZcyLiBZX/o7Fsnev+MbBD2Dt72DKPZB5FVzxTRgwHhb/L5zxvc92t8HXVDcb99V3dYWoU63J8CC+htYayHRWCteLvyUiBgrmMXFoMpsPl9PY6FnctFPFPKW62YQ6X3x5wXEaGi3yxgzw+Hx2agKFZec4W+u86U1PIeSB3bKs/ZZljXP9y7Us6yehHtN5CubJG/2Se0RKP7rJ92scWSeyau/BMP0nkDBQVuF1bZSYLNkB7z4Jo26Ci271b/xtkTIKShwE9u2u7YjRbRS+AQn+NRWw/93AxtbTqD0DC+6TVeb1P5DHwiNg7otQcxqWfCu04wsl5QfF2xHf37fzzgf2AzIxxYL0VhzxINtkiQO7VmA/Vw4nd4lxNusG2L6AiekJVFbXe3WN7/FSI74l2WkJIVuxL9laxNC+vcgdmOjxeXvsWjP+QkIe2DstZ8skOOXmN6XH2PtYvlC4Ts43BmIS4ebn4eRueLeV+YstwUcnQN4vgyvBu5OaI+a5tlbXBfNh4ARnUuiwaeLOVzneN95+XJSRub+VrRebtNGi9myfD9veCNnwQoqd6taahO6JpCGy3VR+0KWcGRjkoExGO6S8tStHNsrX9CnyeVV1nMuj9gLeG8LsLq4iJjKM9GTPjnibrLR4Css7Pl+87EwtH+0rJW+MZxlexta5quN1JjSwe2PXUnGNj54rObH9sptctU6pLIZTh5q3hxx+LUz6Enz0vPfrffgsHNsMeb+AeO/ml4BJGSW/Y+le78eUHxTzn9P8+IgoGDVHtjGC2CWrW3NgtaQzXvIVGHr5hc9f/nUYOFFW7VUlHT++UOO0XWtLwiNFKSs/KHvQqaNkct0WycFt39ruHFknE5hBkyBrBkTEMPDIMpJ6RbLpkOdCNbuLKxmRGu/ZYe9GqPLFz8vwYz3L8AAZfXsRFR6m++we0MDujYJ50rd84AT5efAUMd84bK4ANK3wW+bNTv8R9E53SfLnmj9XXADvPSWB1Gkw9RfbGX+ilQp0dlbA6JudX3f0XJHj96kc3yY1VSLB9xkG1z3m+Rhbkq+tkpQmX/4GuzqWJUqGP4Ed5Lyy/eJraS1/veU5lUUXvjc7K4VrITVXqvJFx0PWdMyOhUxKT/RqoNtTXEV2qvf9dRv3fPGOZMmWIjL7xTF6gPeJWER4GMNS4nTF7gEN7J44WyZFKnLzm2Tw9IvhbKlvHdEK10F41IVFRqITRJIv3SupbDYNdRLsY5Ng9i8C/jXapG+WzPRb22cvmCerRV8+WIdNlbQilePb5q3H4FShBO6oVmTR1ByY9j0pXLPt9Y4bX6g5WyYeA18d8TbJmaJ+1VS0bZw7f06GfC0/5N9rdiSNDSLFu3sHcvOhqpi83gfZW1JFxdnm7ZkrztVx/HTrjnibUOSLl1bV8NG+k63K8DZZaQkdPunoCmhg98TOxSJRu6+Y7Q8FX+T4I+slqHuq3T5sKkz+shSfOfyxPPbBr6DoU9lX74iiJJExslL0tmIvOwBFn/iuHEREQc6Nsp3RlkmwJ7P/PdjwJ7jsPhjiIM/58gfEr7H0Qdnm6Qn464i3Sc6QHHhwvmK3JxFdQY4/sRNqK5tPWrJnQEQsl1W/D8Dmwuar9r0l3mvEtyQU+eLLCo7TaNGqDG+TnRrPkfJznKlRZ7w7Gtg9UTBPZvruK+2UHMnRdmqgq6+VlUJrq4QbnoAklyRfuA5WPQ0X3Qajbwps/L6QkuN9xW674dtqA+uJ3HxZadl1AJTmVJ+GBfdD3xFw7aPOzgkLh5t/Kz3GF3+zZ0jytoktkMAOEJss99qXc7pCYC/0sN0XFQfZM0g7uoII03hBQxhbum4t1c2djs4XX7KliGH94sjp3/b4stQZ75FOUXmuU3G2TPqaX/G15m70sDBp9eh0xX58K9RXt75KiI6XD+q/zIFX8iC2j+SsdyQpOVIxrr7mQmWhYJ64iJP8KOE77BqISZLJQUfWjq+vkVVMa0EvNsn/QBEs3vo+nD4Kd67wrQdBSjZc931Y8Sis/b2zlb43wiMhdbT/WReWJalWvbxXLwsYO7AnDfXvfHv1bWemOKFXX4iK980ZX7oPagIMfr3TfVfqjqyX8bYsHJU7l7Dt87ml72E2H25eA2N3cRWxkeEMTnb2d5edmsCSLUWcra2nV1TrIaPsTC3x0RFERfi3ZjxZVcPH+0u5b9qINmV4aHLG7y6uZFy6h459PRQN7C3ZsUhy1j3Jz4OnwOqfyRs4uo3ZpDfjXEsyr4Ipd4sr+sZn2/dD0hOpo+T3PbkH+l/U9HjpPtkWmO5nWYHwSGl6UzBf5HhfyugGwluPSZGXtrj9lfY3J3pj3zvSUOfyr7WdV+2JS++Vv9NlDwc+lovvkuwLf1j2CGz8M9y1svnfTrCwLJlkJw5q3X/QGsmZEBYpxVucYoxvzvgjG+GP1/o1vGYkZ8LXNvs20SpcJ59LLc/Jmg6Rvbg1eh1fPjyChkbrfE34PSXOHPE27vniYwd7D55Hys8y67n3+czkdB6dM9r57+DGsm3OZXiAoX3EGa+92Zujgb0lBfNk9tt/7IXPpV8s+3VHN8oeeWsUroPEwVLsoi1mPiX77bZLvSNJsZ3xO5t/OG/3ww3fkty50hRm30qpSNfeNNSLsSzzGkkd88bqn0k1t6FXBL+iX1tUV8CCByR9ctr3/LtGWDh8/g04+H5gcvyuJdK5MCdP0jB94cBqWPuifD//K/A/78pkLphsfEV+x7xf+n+N2CT4yvutl0L2RPLQ1tNA3dn6/6SAzm1/AuNnF7Qj6+GDX8pny2AHufYg6mLpHhj/uQufc8nx4/as5lzNrewpqSSnvzjMdxdXcsWIfo6H1rQq9h7YLcvi4de3UFldz4JPj/HI7FFem8u0xpItRQxPiWOkw20C2xmvBrrmaGB358xJ+cC68hueZ812cYvC9W0H9iPrna/GwsJDE9QB+mXJh1FJCwNdwXyRL5PS/b925jWyt1kwr2MC+6EPpXHK5Dtbl//7DIPfXyWpY3f8tf0KAHli+feg8hh8+e3AVIzoeBg5K7CxDJ8Gh9fKROPeNc5yvMEtRW84TH0E3rgL3v+FFNIJFuWHZLth2FT5/wwEJ02OWpKcAXvfbmog443GRpkEZ93gf1tmkPoFa34j7xWngf2Iq1y1Nx9Pbj4xBfOYEraTTYfGk9M/kYpzdRSfrnG8vw5N+eKtGej+vvYwH+4tZerIFN7bdYL1B8u4dJhv2wolldWsPVDK/ddmOZLhbbLTEtjooP98T0LNc+60JsODzP5TRrVtoDt9DCoKnafXhJKIaJcz3s1AV7oPjm8JXKq25fhdb3ZMTnDBPGltmzW99eNClTq2521RMK74Ogz2o31vsImMlTS7jQ5lxgAAH9NJREFUymMSRJ1yPkXvtzD2dhhzB6x+Boq2BGdcjY2w8H7AwE3Pd+zEy6ZPpnhkqtrIPihcKznvgb5XYpOkJ0XBfPn9nVC4TiblgyZ6fn7EDViRcdwave58BTqnpWTdOZ8v7kXuLiw7y5NLd3DliH688B8TiYkMY8mWIsfXt1luy/BeasN7IzstnqOn1BnvjgZ2d7bPl1VIWiv7hekXyxuqtTffeafqJd6P6Uyk5jQP7Hb+eSAyvE1uvhRW2bsy8Gu1RkM97FgI2TOd7cd2dOrYuVOw8AGZGE59pP1fzynpF8te/6a/yAq1Lfa9e2GK3qynxcA1/6vBaf6z8WVRzmb82D/jZjBw6owvmCeNV7JnBP6auflw+ggcddgk6sg62T5zL0HsTlQvTPYMppv1fHJImgjZ+ehZDorTuOMtX7yx0eKh17YQZgxP3TqGuOgIrs1J5c1tUjnOF5ZsLWJEarxPkw57bIDus7uhgd3GluHdi9J4YvAUqD7V+v7bkfXyZu8/JvjjbA9SRknhHTvnvGC+TEp6Dw782hlXi9u/vYvVHPpACgg5Tc3r6NSx5d+V1d/c33quaxBKpj4iXouFXxMPgDeqT8vkpG9W8xS9Xn1gzrPSonh1gFkdZQdgxWOy5z/xi4FdKxCS7Vz2VpzxjY3SUnfE9W2baZ0wcqYUtHLyXmmol6ZUbamCufkkNp4itWwD5Wdq2V1cSa+ocAYl+ZCJgfd88b+vPcSa/aU8mjeKwa6683ljBnKyqoZ1B8ocX19k+DJHRWkuGJubM14RNLDb7Fgoxri2JDXb5d6aHF+4TtptRkQFb3ztSWqO/O6le+DkXijeGjzHeHiE5OW3txxfMA8i42DEDc7PsVPHdi0RA1R7sXs5fPJ3uPKb3mXTUBIZIxOOyuMyAfGGnaI398ULU/RyZsPYz8pe+7FP/BtHY6Pk9oeFh06Ct+md3tRAxhuFH0PV8eC9V2J6yyRh+4K25fiS7aKEtZV1k3UDDRG9mBP2MZsLy9lTUkmWD47485dxBU/3fPHDpWd5culOrs5O4TMXN3lxpuWkEBsZzpKtxxxff9m241g+uOHdGdKnF9ERrXsAehoa2G0K5slKJC239eP6Zskb0Fs+e32NVGvzJ40pVKS4zEUlO2F7EGV4m9x8qDsDe94K3jXdaaiH7QtlxeNrWtSl94o6sfTbEtiCzblyWQmnjoZrHgr+9YPFoEmy97/5b7B7xYXP713pStF7wPvf9qynIC7FJcn70QBo/R9FeZnxZHDUokCIiJKsltYC+3kZfmbwXjc3XyZPtjHOG/bCYnBbbWhjsbJnMTN8HZsPnmR3cZWjUrItaVkzvrHR4sHXPiUizPDULWOarbJ7RUVw7ahUlm07Tn2DM7/A4i1FZKfF+2TqswkPMwxPiQ9Ze9nOiAZ2gKoTcPCDtmV4cBWqmeL9jVe0BRpqu87+OkDf4WLCObHDJcNf6ixNzylDr5Q92PaS4w+uhnNl/q2cbEm+vhoWfSP4kvyyR8SpP/fFzifBt2Tqd2SSt+hrMiGxqa4QCb7fSJjayoo+Nhlu+rWsJlf9zLfXLt0Hb/9AFJcJn/dv/MEmeaj3wN7YICvrrOmSoRAssmdK6lxb75XC9RCX6qjQUsSYW+hjqijbtpITlTWOSsm2ZEiLfPH/W3OQdQfK+P6NoxnoQdbPGzOAk1W1juT4ktPVrD9YxmwfTXPuZKXF64rdDQ3s4CbDO9yfTZ8i6WGe9iML18rXruCIt4mIluC+c4nskwa7cEt4BIy6SSTp2rPBvTbIh2BUvMiY/tBvhHRW2/0mfPqv4I1r51L49J9w1bdg4PjgXbe9iIiG/BelNewytwC+/Lvi/J77YtspetkzYPx/St+DoxudvW5jo6TPhUXKxCCUErw7fTJlz98Th9eIZ8KfcsutEZMoqXPb23DHH1knn0NO7tWI66gJiyX3lJR39mdV7J4vfvDkGZ5atpNpI1O4fZJnZWXayFSXHN+2O/5NW4YPILBnpyVwrKKayuq6tg/uAWhgBwkM/bJFLnXC4IsBS1pBtuTIOnHyJqQFdYjtTortjDftU6veluP3BlmOb6iDHYslp9uX0qwtueQrMOQyePNhSVcMlLNlsPgbkmFx9bcDv15HMXACXPW/8Ok/YNcy2T7Z/DffUvRmPAnxaTD/XmdNgNb+TgLlrKeCqxQFSnIGnCmB2jMXPlcwHyJipf95sMnNl4mUvUhoyZmTYnZtS4a3iYylZOC1zAxfTwT15+ur+0p2WgK7j1fy7dc+JTI8jJ/eMtar0S02KpzrHMrxS7YUMTItwa8tAhtbhVBnvKCBvbJYCps4keFtBk0CjGc5vnB911qt29hFPIZc1j4frkOvkP3XYMvxB1wy/OgAV05h4XDzC7KNsujrgUvybz4sLv25L3YdE6XN1Q/JhGTR18Uf4GuKXmySmN9O7IRVT7V+7Mm9sPIJkaDHeaigFkq8tW+1ZfjsIMvwNtkzWpfj7c8dp93qgJhxt9HHVDEtaqfPjvjzw0qL51hFNesPlvP4jbn07926ejNn7ABKz9SythU5/nhFNesPlfllmms+Npe5T/fZAa0859wN705MoqzuWxroKo5IsQ8f3nCdBru0bHvVT7fl+E//KSsgb7m3vhKoDO9O3+Fw/eNSg/3XE/wvkWo1Sjrk1EdggIfSxJ2diChxyf/hWpngfPbvvvsDsq6HCf8FHz4nWzzeOHNSXm/Os51HgrdxT3lLc1PzDn0kK/n2eq9EJ7jk+AVSbrpl5bvCtRAWIeqKQ/qNn0XVklju6LXB53QyG3tFfV1OKrdMHNTm8VNHptIrKpzFW4q8lrB9c1sRlkVA++sA6S5nvKa8CRrY4/rBRbf6XnYyfQpse6N5yUlPLRS7Clk3wGX3w7jPtN9r5OZLcZM9K4LzodhQJ9XjRs4OXpOZKXeLAnByd2DXycmTvfWuyoBxcMtL8vftb4rejCdle+TMiVYOMjD5S5AY2Ad7u+CtSI3TCoeBkJsPOxdLSt3Qy5s/V7heeln4sPVkImMpH3w9U0vel/eNH5PWK0f0484rMvnK1GGOJgcxkeFcNyqN5QXH+dHNuUSEXygQL9lSRE7/BEb4YehzJzzMMCI1nt0qxQMa2OUN5E+QSZ8ina1O7mqaFBSuk3231irXdVaiE2CGn53cnDL0cnHyFswPTmDfv0qKBQVz5RQWBtNacX73JC66NbDzYxI7vg1xMIlNhujezQO7XeEwa3rwVCdPZM+UVLqCec0De0M9HNsEE7/g8yXTr/wP+NciOLDKL4UrLjqCx270rWtb3pgBLPr0GB/vL+PKrOar9uMV1Ww4VM63bsj2eSyeyE5L4OP9pUG5VldH99j9xd5Hd5fjj6yT1U2wu1x1F8LCxZi3e7lnQ5KvbJ8H0Ym+dyZTFCcYA30ymjvj7UZD7d3yNzpeJg/bF8ievk3xNqg769w4587wa+X90t5VIN2YOjKFuCjPxWqWuhzzswPcX7fJSounqKKa0+qM18DuN32HS6lUu1BEXbXksPvzhutJ5OZD/TkJ7oFQX+tywwdRhleUliRnNF+xb5/f/jK8TW6+pNQdXtP0mB/GufNExsj7Zcfi4NT0d0BMZDjXj05j2bbj1LVwxy/ZWsSoAYkMTwmOATHbVf9+jxroNLD7jTESxAtdb7SiT6Cxrmvur3ckQy6TVKhAVw0HbBk+yHnEiuJOcgacOiReA7vCodNGQ4GSPUO29tzfK4VrIb6/lLz1h9x8ed8cWBWcMTogb8wAys/WsWZfk0x+7NQ5Nh4qZ06QVuvQ5IzXQjUa2AMjfYrssZ8t65qFaUJBWLiUq92zQvp6+0uByvBKB5CcKSmQlcdcjYZOtr8MbxMVJyl12xc2yfGFPhSm8cTwaeIb6EA5/ursFOKjI5q1cj0vwwfohndncHIsMZFhWloWDeyBYa/Oj26UN1xyJsSnhHZMXYHRc6WE6x4/5fj6WnEM5+R1/jKtStfG3RlvNxrK8qHRUKDk5ktq3aGPpCLgqUOBqYIR0dKwZ2cHy/GjUlm+vUmOX7q1iNEDEsnsFzwDYpjLGb+nRFfsGtgDYeBE6QBVuE72vlSGd8aQS0VO9HfVsP89KefbUSsnpediB/bSva7UypmBVTj0lazpsqdfMK/JqBuoKpibL++f/e8FPDyn5I0dyKmzdXy0r5Sjp86x6fCpgIvSeCI7NUH32NHAHhjR8dINruANMbmocc4Z5+X4t6DGj9l1wTyRE4dNC/7YFMWd3unSIGnTX6WSYEdPJqPiZK99x0Ix0YVFSo2BQBjW8XL8VVn9SIiOYMmWY7zpkuEDqQ3vjay0BI6frqbiXM92xmtgD5T0S2Q2D7pi94XcfJHjfXXH19dIJbNRc7peqVal6xEeAUnpcHRD8Coc+kpuvqTYbfqrNBMKNAskIkrePzuX+Nde1w9iIsO5YXQaywuKmf/JUS4alEhGEGV4G7u97N4eLsdrYA8UWxaLjIPUNnq5K02kXwIJA3xfNex7F2oqAq8NryhOseX47A6W4W1G3CByfE1F8My5o+fK9fa9G5zrOWD2mAFUnKtj29HTQTXNuWM743u6gU4De6Cku+T3QRNldq84IyysSY6vPu38vIJ5ENMbhk1tr5EpSnPsmvGh8nRE9ZJJBTR93gTKsKnyPvJlYl1fAxv/4vcq/6pskeMhABm+sRG2vS6ZSB4YlBRLbGR4QDXja8+dYe3rz1FdWe73NUKNBvZASc4UE92odmh12t3JzYeGGudy/PGtsO01yL1FZXil48i4Uto6j7gudGOY9N+icA29MjjXi4iCnBth11LngfqdH8Gir8Gn//LrJaMjwvnMxelMHZnC0L5+yvAb/gSv3Qlv3O2xA2NYmGH0wETWH/TeUa4tDr36IJdsfYxzz18ORzf5fZ1QooE9UIyBu9+FS+4O9Ui6HoOnQMJAZ6uG+lqY/1Wp333t99t/bIpiM+Y2uH99aGR4m2HXwLd2BjedNjcfak7DvnfaPvbwWvjoN/L99vl+v+Sjc0bzypf83E4oOwBv/UAKXO19Cz75u8fDZub2Z9vR0xwq9aNs9cEPyDr4D5Y3TOZsTQ2Nf5oOH78YeBvnDkYDuxI6wsKkctxeB3L8+7+QFfucX0Fc344Zn6J0Z4ZdAzFJbU+sa8/KpLp3unQ/3L8KznRws5XGRlhwv2TU3PU2DL0Clj0irbJbMGtMf0BK1vpETRUN8+7lkJXG2vFPcXevX7LGjIdl34FXPw/nuo40r4FdCS2j50plr11vej+m6FN4/+cw5g4YdWPHjU1RujPhkS53/FLpdeGNd34MZfvg5t/AhP8CqwF2Luq4cQKs/4NU/pvxJCQNgZtfgMZ6WPi1C1bTg5N7MWFIUrNKd454+3HCKg7zYO09zL0ki0dvu4r/PPMNlg96AHYvg99dDUc2BPGXaj80sCuhZfDFkDjI+6qhvhbm3wu9+sKspzt2bIrS3cnNh9pK2LfS8/OH1sDHv4WL75IVfv8x0GdYh+bAU7oP3n5csgMmfF4e65MJNzwh4970fxeckjdmAAXHTnPgpEM5fv8qWP8HViTkczx5AmMG9eay4X3578szuWffZWyb8W857uUZ8NHznV6a18CuhJawMFm171sp1bBasvoZaVU551no1afjx6co3ZnMa8S3UuBh37z2DCy4V1bI1/9QHjNGJgMH3oczJ9t/fI2NsOA+Kcxz06+b18if/GXIuAqWfw9OFTY7bZbLdb/UiRxfUwkL7qcheRjfKr2JvDEDMa7XeWjmSIb27cW9q8I4e+e7kp2w4lH452e9OvM7AxrYldCTm+9Zjj+2WfbWx31O6lsrihJcwiNle2vXUqg71/y5lU9A2X6RvaPdWqvm5oscv6MD5Pi1v5OKe7OegsSBzZ8LC5PtAasRFt7fbBU9KCmWiU7l+Lceg4pC3st5nKrGqGYd53pFRfDMbeMoLD/L0+8dh8/8DWb9TAyHv7tKTIWdEA3sSugZPFmMOe7yXn0NzPsqxKfCzJ+GbmyK0t3JzYfaKtjrJscf/ECC6pS7IfOq5senXQR9R7S/HH9yr0wusmfK5N4TyRkw/UdS937jn5s9lTd2INuLTrP/RCvFava9Cxtehsvu45Uj/Rnatxe5AxObHTIlsw9fujyTv6w5xEf7S+GSe+DLK6RuyZ9nwQe/EmWhE6GBXQk9xkixmr0r4dwpeWzV03BiB9z4nEiFiqK0DxlXQ2yfpkBdUyXyd3IGXP/4hcfbcvzB96HqRPuMqbFBtgEiomQbrrU2tZPvlC2FFd+H8kPnH57tcsd7leOrT8PCB6BvFmWXPMRH+0qZPWbAeRnenW/PGElmvzgeem0LZ2rqYeAEuGe1qB1vPw7/uKNjtiYcooFd6Rzk5kNjnUiCRzfKLHj856UBhqIo7Ud4hEuOf1Pk+Lcflza1N/9WmtB4IjdfJPAdC9tnTB//FgrXwqxnILGNKnXGiCQPIsm7Vs8DescyeWgyi73J8SsehdNHYe6LrNhdQUOj5bUiXmxUOM/cNpajp87x0zd3yIMxveH2VyDvF3BgNfzuSmmv2wnQwK50DgZNEjl+y6vigo/vDzN+EupRKUrPIDcf6s5IAZj1f4BLvgoZV3g/PnU09M1qHzn+xG5JsRuZB2PvcHZO0hD5vDiwGja+fP7h2WMGsPN4JXtLWsjxe1fCpr/A5Q9A+sUs2VpEhgcZ3p3JGX2468pM/vbxYT7Y41qdGyMZA3e9LfX8X8kTw2+IpXkN7ErnwBgpVrP/PTixE256HmKTQj0qRekZZFwlKaXrfi/pbNc91vrxthx/6EOoKgneOGwJPjJWilG1JsG3ZOIXYfi1sOIxyTcvP8iNQ2oZbEr4YL38TPlBmTgsfAD6jYSp36W0qoaP9pWSN9azDO/Ot6aPZFhKHA+/voXKarfWsAPGwj2rpNz1Oz+Gv90S3PviIxrYlc5D7i3ydeIXICsE7TEVpacSHiE+F4xLgu/V9jntIcev+Q0cWQ+zfw4Jab6da4wsCMLC4Y/XwXPjSHl5Ch9Ef4P/Xn8zPDdO/r1wMVQWwdwXITKG5QXFLhl+YJsvERMZzs9vH0dRxTmeXLqz+ZPRCXDrH8UXdHgN/Hk2NNT79jsECW1HpnQeBk2ELyyUlq6KonQs135f3OfpDmu5p46SVW/BfJGjA6VkJ7zzE8iZAxfd6t81eg+Gu1bC0aYKce/vOcH8zcd4aOZI0hJdvez7jYTBkwBYsvUYmf3iGDUgwdFLTBySzP9cPYzfr9rPrIv6c3W2W/1+Y6Rhz+CLpdxtiDp+amBXOhfDrgn1CBSlZ9KrD/TyoUGLLcevehoqi31fYbvTUC/16KPifJfgW5KSLf9cZA+r5o1NKxlSm83Xx2c1O7S0qoY1+0q5d+qINmV4d755fTYrd5Tw8OtbWP7Nq0mMiWx+QFqu/AsRKsUriqIo/pE7F7ACl+M/+jUc2yQO8/jUoAzNJi0xhosz+rBk67ELnltWcJxGC/LG+tYf3pbki09X85PFO4I11KChgV1RFEXxj9RRkJITmDu+eDu891PZ47/oluCNzY28MQPYXVzF7uLKZo8v2VLEsJQ4cvo7k+HdGZ+exD3XDOfVDYW8uyt0RjlPaGBXFEVR/Cc3X/K3K4/7fm5DnUjw0YmQ98vgj83FrIv6YwzNSsyeqKzh4/2lzPFSlMYJ37g+i+y0eL7z+hYqztW1fUIHoYFdURRF8Z/RLjl+ux9y/IfPQtEnMOeXENcv6EOzSU2MYUpGn2Y92m0ZfraPMrw70REiyZ+squVHi7cHY6hBQQO7oiiK4j+pOVKwxlc5/vg2eO9pSXMdfXP7jM2NOWMHsLekSY5fuqWI4SlxjEzzXYZ3Z+zgJL56zXBe23iEd3YWB2OoAaOBXVEURQmM3HzJ3T59oUHNI7YEH5skOesdwIyL+hNmYPGWIkoqq1l7oJS8sQP9luHdeeC6EeT0T+A7r2+l4mzoJflOEdiNMTONMbuMMXuNMd8J9XgURVEUH/BVjn//l3B8i6S2xfVt16HZpCbEcElmX5ZsOcaybSLDzwlAhnfHluRLz9Tyw0UFQblmIIQ8sBtjwoEXgFnAaOBzxpjRoR2VoiiK4piUbEjNdSbHF22B1T+DMbdL85kOZPbYAew7cYbfr9rPiNR4sgOU4d25aFBv7ps2gjc2H+Wt7aGV5EMe2IEpwF7LsvZbllUL/Ato/w0XRVEUJXjk5kPhx1Bx1Psx9bUiwffqC7N+1nFjczEzV+T4o6fOee3kFgj3TxvBqAGJfHfeVsrP1Ab9+k7pDJXnBgGFbj8fAVqtKVpaWsorr7zSnmNSFEVRfCCxrpFbgIrnp1IT7rnWfGRjNcl1x1mZeheF/26nlq9tMDS2NwfORlG3fy2vvBL8NqtXR4fzh6JkvvDcIm4dWNn2Ce1AZwjsnpwL1gUHGXM3cDfAoEGD2ntMiqIoig+cjkxjW+I0kmu99D8H6sJj2Bc/hcJeYzpwZM25pt9ZBlbVkxrd0C7XHxDTwDX9zlJVH0ajBWGBe/N8xljWBTG0YwdgzGXA45ZlzXD9/AiAZVk/9XbO5MmTrQ0bNnh7WlEURVFChmVZQXHbu2OM2WhZ1mQnx3aGPfb1QJYxJtMYEwV8FgiNRqMoiqIoARLsoO4rIZfiLcuqN8bcDywHwoGXLcsKfb6AoiiKonRBQh7YASzLWgosDfU4FEVRFKWr0xmkeEVRFEVRgoQGdkVRFEXpRmhgVxRFUZRuhAZ2RVEURelGaGBXFEVRlG6EBnZFURRF6UZoYFcURVGUboQGdkVRFEXpRmhgVxRFUZRuhAZ2RVEURelGhLy7mz8YY04Ah0I9jhDQDzgZ6kF0A/Q+Bg+9l8FD72Xw6I73cqhlWSlODuySgb2nYozZ4LRtn+IdvY/BQ+9l8NB7GTx6+r1UKV5RFEVRuhEa2BVFURSlG6GBvWvxUqgH0E3Q+xg89F4GD72XwaNH30vdY1cURVGUboSu2BVFURSlG6GBvRNijIkxxqwzxnxqjCkwxvzQ9XimMWatMWaPMeZVY0xUqMfaVTDGhBtjNhtjFrt+1nvpB8aYg8aYrcaYT4wxG1yP9THGvOW6l28ZY5JDPc6ugDEmyRjzmjFmpzFmhzHmMr2XvmGMGen6W7T/nTbGfKOn30cN7J2TGuBay7LGAeOBmcaYS4GngV9ZlpUFlANfDuEYuxpfB3a4/az30n+mWZY13i2d6DvASte9XOn6WWmb54BllmXlAOOQv0+9lz5gWdYu19/ieGAScBaYRw+/jxrYOyGWUOX6MdL1zwKuBV5zPf4XYG4IhtflMMYMBvKAP7p+Nui9DCY3I/cQ9F46whiTCFwN/AnAsqxay7JOofcyEK4D9lmWdYgefh81sHdSXNLxJ0AJ8BawDzhlWVa965AjwKBQja+L8SzwENDo+rkvei/9xQJWGGM2GmPudj2WZllWEYDra2rIRtd1GAacAP7s2iL6ozEmDr2XgfBZ4J+u73v0fdTA3kmxLKvBJS8NBqYAozwd1rGj6noYY+YAJZZlbXR/2MOhei+dcYVlWROBWcB9xpirQz2gLkoEMBF40bKsCcAZephcHExcHpmbgP8X6rF0BjSwd3Jc8tx7wKVAkjEmwvXUYOBYqMbVhbgCuMkYcxD4FyLBP4veS7+wLOuY62sJspc5BSg2xgwAcH0tCd0IuwxHgCOWZa11/fwaEuj1XvrHLGCTZVnFrp979H3UwN4JMcakGGOSXN/HAtcjxpp3gdtch30RWBCaEXYdLMt6xLKswZZlZSBS3TuWZf0nei99xhgTZ4xJsL8HpgPbgIXIPQS9l46wLOs4UGiMGel66DpgO3ov/eVzNMnw0MPvoxao6YQYY8Yiho9wZPL1b8uynjDGDENWnX2AzcDnLcuqCd1IuxbGmKnAg5ZlzdF76TuuezbP9WME8A/Lsn5ijOkL/BsYAhwGbrcsqyxEw+wyGGPGI4bOKGA/8CVc73f0XjrGGNMLKASGWZZV4XqsR/9NamBXFEVRlG6ESvGKoiiK0o3QwK4oiqIo3QgN7IqiKIrSjdDAriiKoijdCA3siqIoitKN0MCuKIqiKN0IDeyKoiiK0o3QwK4oiqIo3QgN7IrSwzHGzDTGvG+MKTfGlBljlhtjRrk9f4kxZpMxptrViWy2McZyVfKzjxltjFlijKk0xpQYY/5pjOkfkl9IUXo4GtgVRYlDGuNMAaYCFcAiY0yUMSYeWAzsBCYh7W+fcT/Z1WRjNVI3fgrS2yAeWGiM0c8YRelgtKSsoijNcDV4OQ1cA+QCPwUGWZZ1zvX8fwB/B6ZZlvWeMeYJpJ3rdW7XSAbKgEssy1rX0b+DovRkdDatKD0cY8xwY8w/jDH7jDGngWLks2EIkANss4O6i7UtLjEJuNoYU2X/Q5pyAAxv7/EritKciLYPURSlm7MIOArc4/paj7QQjQIM0JasFwYsAR708Fyxh8cURWlHNLArSg/G1d5yFHCfZVnvuh6bSNNnww7gC8aYWLdV+5QWl9kE3AEcsiyrrgOGrShKK6gUryg9m3LgJPA/xpgRxphrgN8hq3aQvfQG4A8u5/v1wHddz9kr+ReA3sCrLgf9MGPM9caYl4wxCR33qyiKAhrYFaVHY1lWI/AZYCzian8B+D5Q43q+CrgRMdFtRhzxj7tOr3Ydcwy4AmgElgEFruvU2NdRFKXjUFe8oig+YYy5GZgHpFqWdTLU41EUpTm6x64oSqsYY74I7Eec7hchOe+LNKgrSudEA7uiKG2RBvwQGAAcRxzwD4d0RIqieEWleEVRFEXpRqh5TlEURVG6ERrYFUVRFKUboYFdURRFUboRGtgVRVEUpRuhgV1RFEVRuhEa2BVFURSlG/H/AUfVnO3+K8T1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "ax.axhline(0, c='gray')\n",
    "ax.plot(crosstab.index, crosstab[0], label='non-quitters')\n",
    "ax.plot(crosstab.index, crosstab[1], label='quitters')\n",
    "ax.set_xlabel('age', fontsize=14)\n",
    "ax.set_ylabel('count', fontsize=14)\n",
    "ax.legend(fontsize=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a few combinations of qsmk x age with zero counts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>63</th>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64</th>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65</th>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>6</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>71</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>74</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk  0  1\n",
       "age       \n",
       "63    3  3\n",
       "64    7  1\n",
       "65    3  2\n",
       "66    4  0\n",
       "67    2  0\n",
       "69    6  2\n",
       "70    2  1\n",
       "71    0  1\n",
       "72    2  2\n",
       "74    0  1"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "crosstab.iloc[-10:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a discussion on cells with zero counts, see Fine Point 12.2, pg 155."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 12.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"The effect estimate obtained in the pseudo-population created by weights $0.5 \\, / \\, f(A|L)$\n",
    "is equal to that obtained in the pseudo-population created by weights $1 \\, / \\, f(A|L)$.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "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": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gee = sm.GEE(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    groups=nhefs.seqn,\n",
    "    weights=(0.5 * weights)\n",
    ")\n",
    "res = gee.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Second, we need to estimate Pr[A=1] for the numerator of the weights. We can obtain a nonparametric estimate by the ratio 403/1566 or, equivalently, by fitting a saturated logistic model for Pr[A=1] with an intercept and no covariates.\" pg 154"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "qsmk = (nhefs.qsmk == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.25734355044699875"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# option 1\n",
    "qsmk_mean = qsmk.mean()\n",
    "qsmk_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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": 33,
     "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": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "lgt_pred = res.predict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check for equivalence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_weights = np.zeros(nhefs.shape[0])\n",
    "s_weights[qsmk] = qsmk.mean() * weights[qsmk]    # `qsmk` was defined a few cells ago\n",
    "s_weights[~qsmk] = (1 - qsmk).mean() * weights[~qsmk]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "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": 38,
   "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": 38,
     "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": 39,
   "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": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>qsmk</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sex</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>567.098228</td>\n",
       "      <td>196.513582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>595.423986</td>\n",
       "      <td>205.154456</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "qsmk           0           1\n",
       "sex                         \n",
       "0     567.098228  196.513582\n",
       "1     595.423986  205.154456"
      ]
     },
     "execution_count": 40,
     "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": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1162, 64)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from original dataset\n",
    "\n",
    "intensity25 = nhefs_all.loc[\n",
    "    (nhefs_all.smokeintensity <= 25) & ~nhefs_all.wt82.isnull()\n",
    "]\n",
    "intensity25.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1162, 83)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from reduced dataset\n",
    "\n",
    "intensity25 = nhefs.loc[nhefs.smokeintensity <= 25]\n",
    "intensity25.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the stabilized IP weights $SW^A = f(A) \\, / \\, f(A|L)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"we assumed that the density f(A|L) was normal (Gaussian) with mean $\\mu = E[A|L]$ and variance $\\sigma^2$.  We then used a linear regression model to estimate the mean $E[A|L]$ and variance of residuals $\\sigma^2$ for all combinations of values of L.\" pg 156\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = intensity25.smkintensity82_71\n",
    "X = intensity25[[\n",
    "    'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols = sm.OLS(A, X)\n",
    "res = ols.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_pred = res.predict(X)   # i.e., E[A|L]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The denominator is the distribution, $N(\\mu, \\sigma)$, evaluated at each point of $y = A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "fAL = scipy.stats.norm.pdf(\n",
    "    A,                        # A\n",
    "    A_pred,                   # mu = E[A|L]\n",
    "    np.sqrt(res.mse_resid)    # sigma\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"We also assumed that the density f(A) in the numerator was normal.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFpCAYAAABTSWtMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAF2FJREFUeJzt3XGQnPV93/H3t4hiDZciCHBVhFrRqZox9sU43FBm3D/ujKdg8EQ4E1I81BExrZIZ7HFm1HGFPVO7TZmQaTFtOk5mlEAsN64vjMFBg5Q2ROFCPRPbQZggsOpathUj0IixDarPpmTO+faPfdSsT6e71d3u7Vf7vF8zN9rn9/x29/u9vePD79nnno3MRJIk1fS3hl2AJEk6M4NakqTCDGpJkgozqCVJKsygliSpMINakqTCDGpJkgozqCVJKsygliSpMINakqTC1g27AIBLL700t2zZMuwy+ub73/8+F1544bDLWHNt7Rva27t9t09bex9E3wcPHvx2Zl623LwSQb1lyxaeeuqpYZfRN7Ozs0xNTQ27jDXX1r6hvb3bd/u0tfdB9B0Rf9nLPA99S5JUmEEtSVJhBrUkSYUZ1JIkFWZQS5JUmEEtSVJhBrUkSYUZ1JIkFWZQS5JUmEEtSVJhBrUkSYUZ1JIkFWZQS5JU2LKfnhURbwCeBC5o5n82Mz8aEVcCM8AlwNPAezPzryLiAuBTwDXAd4B/lplHB1S/tCa27Nq35P6dE/Pcscico/fePKiSJLVELyvq14G3Z+ZbgKuBGyPiOuDXgfszcyvwCnBnM/9O4JXM/IfA/c08SZK0AssGdXbMNZvnN18JvB34bDO+B7ilub2t2abZf31ERN8qliSpRXp6jzoizouIZ4CXgceBrwOvZuZ8M+UYsKm5vQl4AaDZfxL48X4WLUlSW0Rm9j45YgPwOeDfAL/bHN4mIjYD+zNzIiKeB27IzGPNvq8D12bmdxY81g5gB8D4+Pg1MzMz/einhLm5OcbGxoZdxpob5b4PvXhyyf3j6+HEa6ePT2y6aEAV1TDKr/lS2to3tLf3QfQ9PT19MDMnl5u37Mlk3TLz1YiYBa4DNkTEumbVfAXwUjPtGLAZOBYR64CLgO8u8li7gd0Ak5OTOTU1dTallDY7O8so9dOrUe57sRPFuu2cmOe+Q6f/Oh29fWpAFdUwyq/5UtraN7S392H2veyh74i4rFlJExHrgXcAh4EngJ9rpm0HHm1u7222afb/SZ7Nsl2SJP1/vayoNwJ7IuI8OsH+UGY+FhFfAWYi4t8DXwYeaOY/APzXiDhCZyV92wDqliSpFZYN6sx8FnjrIuPfAK5dZPz/Arf2pTpJklrOK5NJklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmHLBnVEbI6IJyLicEQ8HxEfbMY/FhEvRsQzzddNXfe5OyKORMRXI+KGQTYgSdIoW9fDnHlgZ2Y+HRE/BhyMiMebffdn5n/snhwRVwG3AW8CfgL444j4R5n5w34WLklSGyy7os7M45n5dHP7e8BhYNMSd9kGzGTm65n5TeAIcG0/ipUkqW3O6j3qiNgCvBX4YjP0/oh4NiIejIiLm7FNwAtddzvG0sEuSZLOIDKzt4kRY8CfAvdk5iMRMQ58G0jgV4GNmfm+iPgE8GeZ+XvN/R4A9mfmwwsebwewA2B8fPyamZmZfvU0dHNzc4yNjQ27jDU3yn0fevHkkvvH18OJ104fn9h00YAqqmGUX/OltLVvaG/vg+h7enr6YGZOLjevl/eoiYjzgYeBT2fmIwCZeaJr/28DjzWbx4DNXXe/Anhp4WNm5m5gN8Dk5GROTU31Uso5YXZ2llHqp1ej3Pcdu/YtuX/nxDz3HTr91+no7VMDqqiGUX7Nl9LWvqG9vQ+z717O+g7gAeBwZn68a3xj17R3A881t/cCt0XEBRFxJbAV+FL/SpYkqT16WVG/DXgvcCginmnGPgy8JyKupnPo+yjwSwCZ+XxEPAR8hc4Z43d5xrckSSuzbFBn5ueBWGTX/iXucw9wzyrqkiRJeGUySZJKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSps2aCOiM0R8UREHI6I5yPig834JRHxeER8rfn34mY8IuI3IuJIRDwbET896CYkSRpVvayo54GdmflG4Drgroi4CtgFHMjMrcCBZhvgncDW5msH8Ft9r1qSpJZYNqgz83hmPt3c/h5wGNgEbAP2NNP2ALc0t7cBn8qOLwAbImJj3yuXJKkFzuo96ojYArwV+CIwnpnHoRPmwOXNtE3AC113O9aMSZKksxSZ2dvEiDHgT4F7MvORiHg1Mzd07X8lMy+OiH3Ar2Xm55vxA8CHMvPggsfbQefQOOPj49fMzMz0p6MC5ubmGBsbG3YZa26U+z704skl94+vhxOvnT4+semiAVVUwyi/5ktpa9/Q3t4H0ff09PTBzJxcbt66Xh4sIs4HHgY+nZmPNMMnImJjZh5vDm2/3IwfAzZ33f0K4KWFj5mZu4HdAJOTkzk1NdVLKeeE2dlZRqmfXo1y33fs2rfk/p0T89x36PRfp6O3Tw2oohpG+TVfSlv7hvb2Psy+eznrO4AHgMOZ+fGuXXuB7c3t7cCjXeO/0Jz9fR1w8tQhckmSdHZ6WVG/DXgvcCginmnGPgzcCzwUEXcC3wJubfbtB24CjgA/AH6xrxVLktQiywZ1815znGH39YvMT+CuVdYlSZLwymSSJJVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYcsGdUQ8GBEvR8RzXWMfi4gXI+KZ5uumrn13R8SRiPhqRNwwqMIlSWqDXlbUnwRuXGT8/sy8uvnaDxARVwG3AW9q7vObEXFev4qVJKltlg3qzHwS+G6Pj7cNmMnM1zPzm8AR4NpV1CdJUqut5j3q90fEs82h8YubsU3AC11zjjVjkiRpBSIzl58UsQV4LDPf3GyPA98GEvhVYGNmvi8iPgH8WWb+XjPvAWB/Zj68yGPuAHYAjI+PXzMzM9OXhiqYm5tjbGxs2GWsuVHu+9CLJ5fcP74eTrx2+vjEposGVFENo/yaL6WtfUN7ex9E39PT0wczc3K5eetW8uCZeeLU7Yj4beCxZvMYsLlr6hXAS2d4jN3AboDJycmcmppaSSklzc7OMkr99GqU+75j174l9++cmOe+Q6f/Oh29fWpAFdUwyq/5UtraN7S392H2vaJD3xGxsWvz3cCpM8L3ArdFxAURcSWwFfjS6kqUJKm9ll1RR8RngCng0og4BnwUmIqIq+kc+j4K/BJAZj4fEQ8BXwHmgbsy84eDKV2SpNG3bFBn5nsWGX5gifn3APespihJktThlckkSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwtYNuwDVtWXXvrOav3NinqnBlCJJreWKWpKkwgxqSZIKM6glSSrMoJYkqTCDWpKkwgxqSZIK88+z+mDhnzHtnJjnjmX+tOnovTcPsiRJ0ohwRS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVtmxQR8SDEfFyRDzXNXZJRDweEV9r/r24GY+I+I2IOBIRz0bETw+yeEmSRl0vK+pPAjcuGNsFHMjMrcCBZhvgncDW5msH8Fv9KVOSpHZaNqgz80nguwuGtwF7mtt7gFu6xj+VHV8ANkTExn4VK0lS26xb4f3GM/M4QGYej4jLm/FNwAtd8441Y8dXXqK0uC279q3ofkfvvbnPlUjS4ERmLj8pYgvwWGa+udl+NTM3dO1/JTMvjoh9wK9l5ueb8QPAhzLz4CKPuYPO4XHGx8evmZmZ6UM7w3HoxZM/sj2+Hk68tvR9JjZdNMCK+mNhX8sZXw+XX7J2fZ1tfaes5Hu/3HOd6TU/F17n1Zibm2NsbGzYZay5tvYN7e19EH1PT08fzMzJ5eatdEV9IiI2NqvpjcDLzfgxYHPXvCuAlxZ7gMzcDewGmJyczKmpqRWWMnx3LFjZ7ZyY575DS39rj94+NcCK+mNhX8vZOTHPz6/h63i29Z2yku/9cs91ptf8XHidV2N2dpZz+Xd3pdraN7S392H2vdI/z9oLbG9ubwce7Rr/hebs7+uAk6cOkUuSpLO37Io6Ij4DTAGXRsQx4KPAvcBDEXEn8C3g1mb6fuAm4AjwA+AXB1CzJEmtsWxQZ+Z7zrDr+kXmJnDXaouSVNupE/l2Tsz3/BaEJ/FJK+OVySRJKsygliSpMINakqTCDGpJkgozqCVJKsygliSpMINakqTCVnoJUUkD4oeNSOrmilqSpMIMakmSCjOoJUkqzKCWJKkwg1qSpMIMakmSCjOoJUkqzKCWJKkwg1qSpMIMakmSCjOoJUkqzKCWJKkwg1qSpMIMakmSCjOoJUkqzKCWJKkwg1qSpMIMakmSCls37AKktbZl175hlyBJPXNFLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFeQlRaYC8XKmk1XJFLUlSYQa1JEmFGdSSJBVmUEuSVNiqTiaLiKPA94AfAvOZORkRlwC/D2wBjgI/n5mvrK5MSZLaqR8r6unMvDozJ5vtXcCBzNwKHGi2JUnSCgzi0Pc2YE9zew9wywCeQ5KkVlhtUCfwRxFxMCJ2NGPjmXkcoPn38lU+hyRJrRWZufI7R/xEZr4UEZcDjwMfAPZm5oauOa9k5sWL3HcHsANgfHz8mpmZmRXXMWyHXjz5I9vj6+HEa0vfZ2LTRQOsqD8W9rWc8fVw+SVr19fZ1jdIvbzmg7aWP1Onvvdn0/e58DPfq7m5OcbGxoZdxlC0tfdB9D09PX2w623jM1pVUP/IA0V8DJgD/iUwlZnHI2IjMJuZP7nUfScnJ/Opp57qSx3DsPDqUzsn5rnv0NLn6R299+ZBltQXZ3tVrZ0T83zg9m0DquZ0la761ctrPmhr+TN16nt/Nn2fCz/zvZqdnWVqamrYZQxFW3sfRN8R0VNQr/jQd0RcGBE/duo28E+B54C9wPZm2nbg0ZU+hyRJbbeaJcA48LmIOPU4/y0z/3tE/DnwUETcCXwLuHX1ZUqS1E4rDurM/AbwlkXGvwNcv5qiJElSh1cmkySpMINakqTCDGpJkgozqCVJKsygliSpMINakqTCDGpJkgozqCVJKsygliSpMINakqTChvtxPzprK/nEqFH61CJpkHr9/do5Mc8dXXP9HdMguaKWJKkwg1qSpMIMakmSCjOoJUkqzKCWJKkwz/pWCSs5m12S2sAVtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQVZlBLklSYQS1JUmEGtSRJhRnUkiQV5sdcSiptpR+BevTem/tciTQcrqglSSrMFbX6aqWrH0nS4lxRS5JUmCtqqcU8AiLVZ1C3gP8xbgdfZ2k0jWRQe5aoVI//IyGtjO9RS5JUmEEtSVJhI3no+1zgYUBJUi8MaknSohZbUOycmOeOJRYa58K5PitZKH3yxgsHUElvBhbUEXEj8J+B84Dfycx7B/VckrTQqB61Wklf50J46swG8h51RJwHfAJ4J3AV8J6IuGoQzyVJ0igb1Mlk1wJHMvMbmflXwAywbUDPJUnSyBrUoe9NwAtd28eAfzyg55IktdCovr2xUGRm/x804lbghsz8F832e4FrM/MDXXN2ADuazZ8Evtr3QobnUuDbwy5iCNraN7S3d/tun7b2Poi+/35mXrbcpEGtqI8Bm7u2rwBe6p6QmbuB3QN6/qGKiKcyc3LYday1tvYN7e3dvtunrb0Ps+9BvUf958DWiLgyIv42cBuwd0DPJUnSyBrIijoz5yPi/cD/oPPnWQ9m5vODeC5JkkbZwP6OOjP3A/sH9fjFjeQh/R60tW9ob+/23T5t7X1ofQ/kZDJJktQffiiHJEmFGdR9EhH/ISL+V0Q8GxGfi4gNXfvujogjEfHViLhhmHUOQkTcGhHPR8RfR8Tkgn2j3vuNTW9HImLXsOsZpIh4MCJejojnusYuiYjHI+Jrzb8XD7PGQYiIzRHxREQcbn7OP9iMj3TvEfGGiPhSRPxF0/e/bcavjIgvNn3/fnPC8MiJiPMi4ssR8VizPbS+Der+eRx4c2b+FPC/gbsBmkun3ga8CbgR+M3mEquj5DngZ4EnuwdHvfcWXir3k3Rex267gAOZuRU40GyPmnlgZ2a+EbgOuKt5nUe999eBt2fmW4CrgRsj4jrg14H7m75fAe4cYo2D9EHgcNf20Po2qPskM/8oM+ebzS/Q+dtx6Fw6dSYzX8/MbwJH6FxidWRk5uHMXOyCNaPee6sulZuZTwLfXTC8DdjT3N4D3LKmRa2BzDyemU83t79H5z/emxjx3rNjrtk8v/lK4O3AZ5vxkesbICKuAG4GfqfZDobYt0E9GO8D/rC5vdjlVDeteUXDMeq9j3p/vRjPzOPQCTTg8iHXM1ARsQV4K/BFWtB7c/j3GeBlOkcNvw682rUoGdWf+f8EfAj462b7xxli334e9VmIiD8G/u4iuz6SmY82cz5C51DZp0/dbZH559yp9r30vtjdFhk753pfwqj3py4RMQY8DPxKZv6fziJrtGXmD4Grm3NuPge8cbFpa1vVYEXEu4CXM/NgREydGl5k6pr1bVCfhcx8x1L7I2I78C7g+vybv3tb9nKq54Llej+Dkeh9CaPeXy9ORMTGzDweERvprLxGTkScTyekP52ZjzTDregdIDNfjYhZOu/Rb4iIdc3qchR/5t8G/ExE3AS8Afg7dFbYQ+vbQ999EhE3Av8a+JnM/EHXrr3AbRFxQURcCWwFvjSMGodg1Hv3Urmdfrc3t7cDZzq6cs5q3p98ADicmR/v2jXSvUfEZaf+eiUi1gPvoPP+/BPAzzXTRq7vzLw7M6/IzC10fqf/JDNvZ4h9e8GTPomII8AFwHeaoS9k5i83+z5C533reTqHzf5w8Uc5N0XEu4H/AlwGvAo8k5k3NPtGvfeb6Pzf9qlL5d4z5JIGJiI+A0zR+RShE8BHgT8AHgL+HvAt4NbMXHjC2TktIv4J8D+BQ/zNe5YfpvM+9cj2HhE/ReekqfPoLOoeysx/FxH/gM6Jk5cAXwb+eWa+PrxKB6c59P2vMvNdw+zboJYkqTAPfUuSVJhBLUlSYQa1JEmFGdSSJBVmUEuSVJhBLUlSYQa1JEmFGdSSJBX2/wBRopJEAU4FtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "A.hist(bins=30, ax=ax);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.057659208261618, 10.467830908151853)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.mean(), A.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "fA = scipy.stats.norm.pdf(A, A.mean(), A.std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the stabilized IP weights are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw = fA / fAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "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": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = intensity25.wt82_71\n",
    "X = pd.DataFrame(OrderedDict((\n",
    "    ('constant', np.ones(y.shape[0])),\n",
    "    ('A', A),\n",
    "    ('A^2', A**2)\n",
    ")))\n",
    "\n",
    "model = sm.GEE(\n",
    "    y,\n",
    "    X,\n",
    "    groups=intensity25.seqn,\n",
    "    weights=sw\n",
    ")\n",
    "res = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 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>    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": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the estimate and confidence interval for \"no change\", you can read off the values in the `constant` row above (because `A` and `A^2` will be zero).\n",
    "\n",
    "Getting Statmodels to calculate the estimate and confidence interval for when smoking increases by 20 cigarettes / day will take a couple extra steps. In Chapter 11, the regression result had a `get_prediction` method. The GEE result doesn't (yet?) have that _method_, so we'll use the hidden `get_prediction` _function_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.regression._prediction import get_prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>mean_ci_lower</th>\n",
       "      <th>mean_ci_upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>2.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.9</td>\n",
       "      <td>-1.7</td>\n",
       "      <td>3.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mean  mean_ci_lower  mean_ci_upper\n",
       "0   2.0            1.4            2.6\n",
       "1   0.9           -1.7            3.5"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred_inputs = [\n",
    "    [1, 0, 0],       # no change in smoking intensity\n",
    "    [1, 20, 20**2],  # plus 20 cigarettes / day\n",
    "]\n",
    "pred = get_prediction(res, exog=pred_inputs)\n",
    "summary = pred.summary_frame().round(1)\n",
    "summary[[\"mean\", \"mean_ci_lower\", \"mean_ci_upper\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can relabel the rows and columns to make this table a little nicer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>estimate</th>\n",
       "      <th>CI lower</th>\n",
       "      <th>CI upper</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>no change</th>\n",
       "      <td>2.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>2.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>+20 per day</th>\n",
       "      <td>0.9</td>\n",
       "      <td>-1.7</td>\n",
       "      <td>3.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             estimate  CI lower  CI upper\n",
       "no change         2.0       1.4       2.6\n",
       "+20 per day       0.9      -1.7       3.5"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = summary[[\"mean\", \"mean_ci_lower\", \"mean_ci_upper\"]]\n",
    "summary.index = [\"no change\", \"+20 per day\"]\n",
    "summary.columns = [\"estimate\", \"CI lower\", \"CI upper\"]\n",
    "summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: since the `get_predictions` function wasn't attached to the GEE regression result, it might not work correctly with other versions of the GEE model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 12.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"if interested in the causal effect of quitting smoking A (1: yes, 0: no) on the risk of death D (1: yes, 0: no) by 1982, one could consider a _marginal structural logistic model_\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "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": 57,
     "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": 58,
   "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": 59,
   "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": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_AV = numer * weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "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": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 83)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "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": 63,
     "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": [
    "The last line gives the estimate $\\hat{\\theta}_2$ = -0.2 and confidence interval (-2.2, 1.9)."
   ]
  },
  {
   "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": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1629, 64)"
      ]
     },
     "execution_count": 64,
     "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": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['constant'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add dummy features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "edu_dummies = pd.get_dummies(nhefs_all.education, prefix='edu')\n",
    "exercise_dummies = pd.get_dummies(nhefs_all.exercise, prefix='exercise')\n",
    "active_dummies = pd.get_dummies(nhefs_all.active, prefix='active')\n",
    "\n",
    "nhefs_all = pd.concat(\n",
    "    [nhefs_all, edu_dummies, exercise_dummies, active_dummies],\n",
    "    axis=1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add squared features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in ['age', 'wt71', 'smokeintensity', 'smokeyrs']:\n",
    "    nhefs_all['{}^2'.format(col)] = nhefs_all[col] * nhefs_all[col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also add a feature to track censored individuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['censored'] = nhefs_all.wt82.isnull().astype('int')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the IP weights for treatment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ip = nhefs_all[[\n",
    "    'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "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": 71,
   "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": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_A = ip_numer / ip_denom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "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": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "# same as previous, but with 'qsmk' added\n",
    "\n",
    "X_ip = nhefs_all[[\n",
    "    'constant', 'sex', 'race', 'edu_2', 'edu_3', 'edu_4', 'edu_5', \n",
    "    'exercise_1', 'exercise_2', 'active_1', 'active_2',\n",
    "    'age', 'age^2', 'wt71', 'wt71^2',\n",
    "    'smokeintensity', 'smokeintensity^2', 'smokeyrs', 'smokeyrs^2',\n",
    "    'qsmk'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "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": 76,
   "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": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_C = ip_numer / ip_denom\n",
    "sw_C[nhefs_all.censored == 1] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "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": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw_AC = sw_A * sw_C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "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": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "wls = sm.WLS(\n",
    "    nhefs.wt82_71,\n",
    "    nhefs[['constant', 'qsmk']],\n",
    "    weights=sw_AC[nhefs_all.censored == 0]\n",
    ")\n",
    "res = wls.fit(cov_type='cluster', cov_kwds={'groups': nhefs.seqn})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "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": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}