chapter14.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 14"
   ]
  },
  {
   "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 namedtuple\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats\n",
    "import scipy.optimize\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": "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": "markdown",
   "metadata": {},
   "source": [
    "Add a variable for censored weight, then add`'constant'`, dummy variables, and squared variables, as done in previous chapters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['censored'] = nhefs_all.wt82.isnull().astype('int')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['constant'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": "code",
   "execution_count": 9,
   "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": [
    "Subset the data as described in the margin, pg 149."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1566, 81)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 14.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 14.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"In our smoking cessation example, we will use the nonstabilized IP weights $W^C = 1 \\, / \\, \\Pr[C = 0|L,A]$ that we estimated in Chapter 12. Again we assume that the vector of variables $L$ is sufficient to adjust for both confounding and selection bias.\" pg 174"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ip = nhefs_all[[\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",
    "    'qsmk'\n",
    "]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can reuse a function from chapter 12 to help us create IP weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.142836\n",
      "         Iterations 8\n"
     ]
    }
   ],
   "source": [
    "weights = 1 / logit_ip_f(nhefs_all.censored, X_ip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ip_censor = weights[nhefs_all.censored == 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   min     mean      max\n",
      "------------------------\n",
      "  1.00     1.04     1.82\n"
     ]
    }
   ],
   "source": [
    "print('   min     mean      max')\n",
    "print('------------------------')\n",
    "print('{:>6.2f}   {:>6.2f}   {:>6.2f}'.format(\n",
    "    ip_censor.min(),\n",
    "    ip_censor.mean(),\n",
    "    ip_censor.max()\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 14.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still Program 14.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"all individuals can be ranked according to the value of their observed outcome Y\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "ranked = nhefs.sort_values('wt82_71', ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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>seqn</th>\n",
       "      <th>wt82_71</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1366</th>\n",
       "      <td>23522</td>\n",
       "      <td>48.538386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>259</th>\n",
       "      <td>6928</td>\n",
       "      <td>47.511303</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       seqn    wt82_71\n",
       "1366  23522  48.538386\n",
       "259    6928  47.511303"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ranked[['seqn', 'wt82_71']][:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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>seqn</th>\n",
       "      <th>wt82_71</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1328</th>\n",
       "      <td>23321</td>\n",
       "      <td>-41.28047</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       seqn   wt82_71\n",
       "1328  23321 -41.28047"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ranked[['seqn', 'wt82_71']][-1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 14.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 14.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"In our smoking cessation example, we first computed each individual’s value of the 31 candidates ...\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to need a few different things from the regressions, wo we'll first create a container that holds that information together.\n",
    "\n",
    "We'll want\n",
    "1. absolute value of the coefficient (the basis of comparison),\n",
    "2. the value of the coefficient\n",
    "3. the value of psi that produced the coefficient, and\n",
    "4. the p-value (for finding the 95% confidence interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "GInfo = namedtuple('GInfo', ['abs_alpha', 'alpha', 'psi', 'pvalue'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll create a function that will perform the regression and return the info we need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logit_g_info(psi, data, y, X_cols, weights):\n",
    "    \"\"\"\n",
    "    Return logistic regression parameters to identify best `psi`\n",
    "    \n",
    "    Note: this is written specifically for the problem in this program\n",
    "    \n",
    "    Paramters\n",
    "    ---------\n",
    "    psi : float\n",
    "    data : Pandas DataFrame\n",
    "        needs to contain the given `X_cols`\n",
    "    y : Pandas Series or Numpy array\n",
    "    X_cols : list of strings\n",
    "        column names for `X`\n",
    "    weights : Pandas Series or Numpy array\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    GInfo namedtuple, containing\n",
    "    - absolute value of H_of_psi coefficient\n",
    "    - H_of_psi coefficient\n",
    "    - psi value\n",
    "    - p-value for H_of_psi coefficient\n",
    "    \n",
    "    \"\"\"\n",
    "    data['H_of_psi'] = data.wt82_71 - psi * data.qsmk\n",
    "    X = data[X_cols]\n",
    "    \n",
    "    gee = sm.GEE(y, X, groups=data.seqn, weights=weights, family=sm.families.Binomial())\n",
    "    res = gee.fit()\n",
    "    \n",
    "    alpha = res.params.H_of_psi\n",
    "    pvalue = res.pvalues.H_of_psi\n",
    "    return GInfo(abs(alpha), alpha, psi, pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all uses here, `y` and the `X` columns are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = nhefs.qsmk\n",
    "X_cols = [\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",
    "    'H_of_psi'\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll run the regression once for the known right answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_info = logit_g_info(3.446, nhefs, y, X_cols, weights=ip_censor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "psi: 3.446  regression coefficient: -1.9e-06\n"
     ]
    }
   ],
   "source": [
    "print('psi: {}  regression coefficient: {:>0.2g}'.format(g_info.psi, g_info.alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll do the course-grained search for $\\psi$ with $H(2.0), H(2.1), H(2.2), \\ldots , H(4.9)$ and $H(5.0)$ (pg 178)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi_vals = np.arange(2.0, 5.0, 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_info = [\n",
    "    logit_g_info(psi, nhefs, nhefs.qsmk, X_cols, ip_censor)\n",
    "    for psi in psi_vals\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# by default, `min` will minimize the first value,\n",
    "# which in this case is the absolute value of the coefficient\n",
    "best = min(g_info) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best psi: 3.4000  best alpha: 0.00086\n"
     ]
    }
   ],
   "source": [
    "print('best psi: {:>0.4f}  best alpha: {:>0.5f}'.format(best.psi, best.alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below shows $p$-value as a function of $\\psi$, with a red line at 0.05.\n",
    "\n",
    "To find the 95% confidence interval, we find the last values that are below the line, coming from the left and the right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGPCAYAAABBBJmLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXyU5b3//9cn+wYJWSAQAoGwg7KjLALVqqh1aV3b2mOrrdWentPlfE9Pe5Zqbfs7red72p62frW22u20tVar4i7WBVlUNhEIWyAsCSFkIwGyZ67fHzPRGBOYkGTuzMz7+XjMw8w9d+555/ZmPnNf931dlznnEBERkcgR43UAERER6V8q7iIiIhFGxV1ERCTCqLiLiIhEGBV3ERGRCKPiLiIiEmFU3EVERCKMiruIiEiEUXEX6QdmdreZhXxEKDO71cz2mlmLmR0P9fsPpI59amZxXmfpC6+ODYluKu4iYcrMRgEPAuuAC4GPeptIRAaLsP5GLBLlJgKxwG+dc2v6Y4Nmluica+6PbYmId3TmLmHNzN42s0fN7B4z22dmTWb2rplddJrfuSHQ3HtuN689b2bvdHo+wcx+b2YlZtZoZvvN7H4zG3aGXL8xswPdLH/NzF7rsmymma00s9rAe6w1swvOtH2gYzt/C/w9v+n0+gozWx/YXp2ZPWlmk7tso6PZe4aZvWhmJ4FHT/OeZ8zZm/0V2N4TZlYdWHe3mX2rm7ceZ2bPmtlJMztoZt82s9N+dpnZUDPzmdntgeefCzxPCzz/m5k918PvBnV8hOLYCCzv9fEhouIuYStwLfYc4ApgCfA14JP4W6T+ambZPfzqSqAOuLnL9kbgb9r+fafFo4BS4KvApcA9wEVAt4XhLP6GOfib1TOBLwDXAtXAy2Y29zS/+l3gHwM//z2wMLAMM1sBPAucBG4E7gRmAGvMLK+bbT0FvA5cBfy4jzmD2l9mtgBYDxTi//92BfAjYHQ3b/8E8ApwDfAk8B3glu5ydjIbMGBL4PkcYK9z7mSn1zf38LvBHh8DemwE3vNsjw+Jds45PfQIywdwLuDwF6bYTsuXBZZ//DS/+0v8H8wxnZZ9FWgDRp7m9+Lwf5FwwOxOy+/2/3N67/lvgAPd/P5rwGudnv8N2AkkdFoWG1j25Bn+/o8GcizvsnwjsBeI67RsHNAK/KhrZuArQezrs8p5mv21GjgMpJzmdzvyfa7L8m3AS2fI+7XA35sUeL4W+FOnfeGAT/Tn8dHfx0Zfjw89ovuhM3cJZx1nLv/qnGvvtHxX4L9Z5hfX+RF47fdAHv4b0Tp8BnjZOVfescDMEszsX81sl5k14i8YbwRe/kAzd2+ZWTL+LyJ/AXyd8hnwMrD0LLaZiv8s9c/OubaO5c65EvwFblk3v/ZEf+UMZn+ZWQqwGPiDc64hiD/r2S7PtwNjzvA7c4Ai51xToAl/Ju+fqc8J/LenM3cI4vgYyGMjsP1+Pz4keqi4SzibAxxxzq3tsnxU4L+l+D8cW7s8wP8hfAD/BzZmNjWwvc5N8gD/if/M63/xNx0vAD4ReC2pj/kz8Z+F/Uc3Gb8MDDvTteVuDMP/4V/ezWtHA+/ZVXfrnm3OYPbXMPyfPaVn/Gv8aro8b+bM+34O7xfvyUBqp+ezgRrn3IHT/H4wx8dAHhswMMeHRAndLS/hbA5Q1s3yG4EG/B/QMcD8ris455yZ/S/wVTO7E/+H+Ek+fBZ7E/A759z3OhZ03JR1Bk1AQjfLs/BfMwU4DviA+4DfdbcR55wviPfqrBZ/s3BuN6/ldnrvD7zNGbbZm5zB7K/awPa6u/7fZ4GCN4n3bw7seqZ+Mf7r2D0K8vgYyGMDBub4kCih4i5hqVNT6ykzi+togjZ/3+8vAT93zp0KrL6xh838Hvh3/GdbnwYe76aZOIX3z/Y7fC6IiAeBEWaW7ZyrCmQrxH8WuQ7AOXfKzN4I/B2b++ODOrDNTcD1ZnZ3x+UKMxsLLAJ+dpbbDDbnGfeXc67BzNYAN5vZPc65xt5mOgML/LejBWcO/mvctWZ2Of4z7GuD2M6Zjo8BOzZgYI4PiR4q7hKupuBvaq0BfmNmv8Z/p/W38V9zv+tMG3DO7TGzt4Af4D+L7NokD/ACcIuZbQOK8X/QLwoi31/w373+BzP7EZANfAuo6rLe1/HfXPaimT2Ev4k8G39BinXOfTOI9+rqP/Bfp37GzP4fkIb/DvM64L/PYnu9yRns/vo/+G+EXG9m/42/iX48MMs59w9nmREA51y7mT0OfN7MYoHzgUYz+yXwd8ADzrm/BrGdMx0fA31swMAcHxIFdL1GwlVHU+vlQAbwNHAv/m5IFznnmoLcTseNU2XAq928/g/4u0Z9H/gzMAR/d7vTcs4VA9cFtv0k8A38H9R7uqy3Gf9lg2rgp8BLwP/g7+K3Osi/oet7v4D/GnAG/qbpB/DfXb3EOXfkLLcZbM6g9pdzbgP+m+oO429NeA74Z4K/Dn8mtwDfDOSbjv+LIPh7FtzZi+2c7vgY0GMjsG6/Hx8SHcw5DXks4cfMfgxc55zL9zqLDF6B5u5i4LLAlx6RqKAzdwlXc4BNXoeQQS+Ybm8iEUfFXcKOmRkwCxV3ObM5QJlz7pjXQURCSc3yIiIiEUZn7iIiIhFGxV1ERCTChH0/9+zsbFdQUOB1DBERkZDYtGlTlXMu53TrhH1xLygoYOPGngYgExERiSxmdvBM66hZXkREJMKouIuIiEQYFXcREZEIo+IuIiISYVTcRUREIoyKu4iISIRRcRcREYkwKu4iIiIRRsVdREQkwqi4i4iIRBgVdxERkQij4i4i3WpqbWdbaR1t7T6vo4hIL4X9xDEiMjB+8fp+fvzyHjJS4rloyggumT6CpRNzSE6I9TqaiJyBiruIdGv13krGZ6cya0wGL++s4PHNpSTGxXDBxBwumT6Ci6YMJyst0euYItINFXcR+ZCTzW1sPXycLy4bzz9fOoXWdh8bDtTw0o4KVhVV8PLOCmIM5hVkcsm0EVwyLZcxWSlexxaRABV3EfmQt0uqafM5FhVmAxAfG8OiwmwWFWZz15XT2HGknpeKKnhpx1G+9+xOvvfsTqbkDuGS6blcMm0E00cNxcw8/itEopeKu4h8yLriahLiYpg7dtiHXjMzZuSlMyMvna9fPIlD1Q28VHSUl4oq+Pkre/np3/YyMz+Dx+9YSFys7tkV8YKKu4h8yNp91cwdM4yk+DPfPDcmK4XPXzCez18wnuqTzTy8toT7Xt3H1tK6br8ciMjA09dqEfmAmlMt7CyvZ/GErF7/blZaIrctGY8ZrCuuGoB0IhIMFXcR+YD1+6oBWDQh+6x+PzM1gemjhrJGxV3EMyruIvIBa/dVkZYYx7l56We9jcWF2Ww+VEtDS1s/JhORYKm4i8gHrN9XzXnjMvt0M9ziCdm0tjs2HKjtx2QiEiwVdxF5z5HjjZRUnTrrJvkO8wsySYiNYa2a5kU8oeIuIu9Z13G9vbD3N9N1lpwQy9yxw1izV8VdxAsq7iLynnXFVWSlJjB5xJA+b2vxhCyKyuupPtncD8lEpDdU3EUEAOcca/dVcX5hFjExfR9dbnGgaX/9/uo+b0tEekfFXUQA2F91ior6ZhYX9u16e4dz8tIZkhSn6+4iHlBxFxHg/UFnzmbwmu7ExcZw/vgs9XcX8YCKu4gAsLa4mryMZMZk9t/sbksmZHO4ppFD1Q39tk0ROTMVdxHB53Os31/NosKsfp3NreO6+9p9OnsXCSUVdxGhqLyeusZWFvVTk3yHwpxUcocmqWleJMRU3EWEdYEz60X9dDNdBzNj0YQs1hVX4fO5ft22iPRMxV1EWFtczYThaYwYmtTv214yIZvahlaKyuv7fdsi0j0Vd5Eo19Lm4+2Smj6PSteTjuvu63TdXSRkVNxFotzW0uM0trb3e5N8hxFDk5g4PI01xRrMRiRUVNxFotza4irMYOH4gTlzB//Z+9sl1TS3tQ/Ye4jI+1TcRaLcuuJqZoxKJz0lfsDeY/GEbJpafWw5dHzA3kNE3qfiLhLFGlra2HK4tt+7wHV13vhMYmNMQ9GKhIiKu0gU23CgltZ2N2DX2zsMTYpn5uh09XcXCREVd5Eotm5fFfGxxvyCYQP+XosnZLP18HHqm1oH/L1Eop2Ku0gUW1dczewxw0hJiBvw91o8IRufg7f21wz4e4lEOxV3kSh1vKGF7UfqBqx/e1ezx2SQHB+r6+4iIaDiLhKl3txfg3PvDzIz0BLjYpk/LlPX3UVCQMVdJEqt21dFSkIsM0dnhOw9l0zIovjYSY7WNYXsPUWikYq7SJRaW1zF/IJMEuJC9zGgoWhFQkPFXSQKVdQ3sa/yFIsHuH97V1Nzh5KZmqCmeZEBpuIuEoUGaorXM4mJMRYWZrG2uArnNAWsyEBRcReJQuuKq8lIiWfayKEhf+8lE7KpqG9mX+WpkL+3SLQIaXE3sxVmttvMis3sm928PsbMXjWzLWb2rpldHsp8ItHAOce6fdUsHJ9FTIyF/P2XBK67q0ucyMAJWXE3s1jgPuAyYBrwSTOb1mW1fwcedc7NBm4C/l+o8olEi4PVDZQdbwxZ//au8jNTGJOZouvuIgMolGfuC4Bi59x+51wL8AhwdZd1HNDRTpgOHAlhPpGosG6ff171RSHq396dxROyeHNfNW3tPs8yiESyUBb3POBwp+elgWWd3Q3cbGalwHPAP4Qmmkj0WLuvityhSYzPTvUsw+IJ2ZxobmNbWZ1nGUQiWSiLe3cX97reLvtJ4DfOudHA5cDvzexDGc3sdjPbaGYbKysrByCqSGTy+Rzr91WzqDALs9Bfb+/QcZe+rruLDIxQFvdSIL/T89F8uNn9NuBRAOfceiAJ+FDboXPuQefcPOfcvJycnAGKKxJ5dlecoOZUi6dN8gCZqQlMHzVU191FBkgoi/sGYKKZjTOzBPw3zK3sss4h4CIAM5uKv7jr1Fykn3ScKXt1M11niydks/ngcRpb2r2OIhJxQlbcnXNtwJeBF4Gd+O+K32Fm95jZVYHV/gn4gpltBf4EfNZppAuRfrN+XzXjslMZlZHsdRQWT8impd3HhgOaAlakvw38JM6dOOeew3+jXOdl3+70cxGwOJSZRKJFW7uPt0pquHrWKK+jADC/YBgJsTGsLa5i6SRdXhPpTxqhTiRKbC2t42RzW8iHnO1JSkIcs8dk6Lq7yABQcReJEusD48kvHATX2zssmZBNUXk9NadavI4iElFU3EWixNriaqaN9M/KNlgsnpiNc/57AUSk/6i4i0SBptZ2Nh2qHRR3yXd2bl46QxLj1DQv0s9U3EWiwKaDtbS0+Vjscf/2ruJiYzhvfJYGsxHpZyruIlFgbXEVcTHG/HGZXkf5kCUTsjhU08Dhmgavo4hEDBV3kSiwbl81M/MzSEsMae/XoCyZqKFoRfqbirtIhKtvauXd0uMsHmTX2zsU5qQxYmiirruL9CMVd5EI99b+GnwOFg6S/u1dmRmLC7NZt68an08DUor0BxV3kQi3triKpPgY5ozN8DpKjxZPyKbmVAu7jp7wOopIRFBxF4lwq/dUcv74LBLjYr2O0qOOu/h13V2kf6i4i0SwwzUN7K86xdKJg3vs9tz0JCYMT9N1d5F+ouIuEsFW7/XPmBwOE7MsLszi7ZIamts0BaxIX6m4i0Sw1XsqGZWeRGFOqtdRzmjppBwaW9vZUFLrdRSRsKfiLhKh2tp9rCuuZumkHMzM6zhntLAwi4TYGF7bfczrKCJhT8VdJEK9c/g4J5rbwqJJHvxTwJ43PpPX9lR6HUUk7Km4i0So1XsqiTFYPEj7t3dn2aQcio+dpLRWQ9GK9IWKu0iEen1vFbPyM0hPifc6StCWTx4OwGu7dfYu0hcq7iIRqPZUC++WHueCQd4FrqvCnFRGD0tWcRfpIxV3kQi0prgK58KjC1xnZsbyyTms21elLnEifaDiLhKB3thbydCkOGaOTvc6Sq8tnzSchpZ2Nh5QlziRs6XiLhJhnHOs3lPFkonZxMWG3z/xRRPUJU6kr8LvX76InNbeYyc5Wt806Iec7UlKQhwLxmXqurtIH6i4i0SY1YF+4heE2fX2zpZPzmGvusSJnDUVd5EI8/qeSgpzUsnLSPY6yllTlziRvlFxF4kgTa3tvF1SE3Z3yXelLnEifaPiLhJB/LOq+cK+uKtLnEjfqLiLRJDVeypJiIvh/HFZXkfpM3WJEzl7Ku4iEWT13koWFGSSnBDrdZQ+U5c4kbOn4i4SIcrrGtlTcZKlk8JnopjTUZc4kbOn4i4SId7YUwUQduPJn05Hl7iy441eRxEJKyruIhHi9b2VDB+SyJTcIV5H6TfLJ/u/qKhpXqR3VNxFIkC7z7G2uIoLJuZgZl7H6TeFOWnkZahLnEhvqbiLRIBtZXUcb2iNmOvtHd7rEldcRUubz+s4ImFDxV0kAqzeU4lZZF1v77B88nBOtbSz8UCN11FEwoaKu0gEWL2nknPy0slMTfA6Sr9bVBjoErdHTfMiwVJxFwlz9U2tbDl8nAsmRlaTfIfUxI4ucbqpTiRYKu4iYW5dcRXtPhe2U7wGY/nkHPZUqEucSLBU3EXC3Oq9VaQlxjFn7DCvowwYdYkT6R0Vd5Ew5pxj9Z5KFhZmER8buf+c1SVOpHci99NAJAqUVJ2itLYx7GeBOxN1iRPpHRV3kTC2OnAH+dIIvZmuM3WJEwmeirtIGFu9t4qxWSmMzUr1OsqAU5c4keCpuIuEqea2dtbvq47ou+Q7S02MY/64YbqpTiQIKu4iYWrTwVoaW9sj/np7Z8snDWdPxUmOqEucyGmpuIuEqdV7qoiLMRYWZnkdJWTe7xKnpnmR01FxFwlTq/dUMnfsMNIS47yOEjIThnd0iVPTvMjpqLiLhKHKE80UlddHVZM8+LvELZucw1p1iRM5LRV3kTD0xt6OLnDRVdwBlk/K8XeJO6gucSI9UXEXCUOr91SSlZrA9FFDvY4ScosnZBMfa7yu6+4iPVJxFwkzPp9jTXEVSyZmExNjXscJufdniVNxF+mJirtImCkqr6fqZEtUNsl3WD5pOLsrTqhLnEgPVNxFwszqwPX2SJ2/PRgdXeJe12h1It1ScRcJM6v3VDIldwjDhyZ5HcUzHV3iXt2lLnEi3VFxFwkjp5rb2HSwlmVR1gWuK3WJEzk9FXeRMLJ+XzWt7S7q+rd3R13iRHoW0uJuZivMbLeZFZvZN3tY5wYzKzKzHWb2x1DmExns3thbSXJ8LPMKhnkdxXOL1CVOpEchK+5mFgvcB1wGTAM+aWbTuqwzEfgWsNg5Nx34aqjyiYSD1XurOH98JolxsV5H8VxaYhzzC9QlTqQ7oTxzXwAUO+f2O+dagEeAq7us8wXgPudcLYBzTnfLiAQcrmmgpOoUF0RxF7iulk/OUZc4kW6EsrjnAYc7PS8NLOtsEjDJzNaa2ZtmtqK7DZnZ7Wa20cw2VlbqW7tEh45uX7re/r7lk4cD6hIn0lUoi3t3Q2m5Ls/jgInAcuCTwK/MLONDv+Tcg865ec65eTk5+qCT6PDyzgrGZKZQmJPqdZRBY+LwNEalJ2mWOJEuQlncS4H8Ts9HA0e6Wecp51yrc64E2I2/2ItEtZPNbawrrubiaSMwi74hZ3tiZiyfMpy1xdXqEifSSSiL+wZgopmNM7ME4CZgZZd1ngQ+AmBm2fib6feHMKPIoLR6TyUt7T4unjbC6yiDzsVTR3CyuY21+6q8jiIyaISsuDvn2oAvAy8CO4FHnXM7zOweM7sqsNqLQLWZFQGvAv/snKsOVUaRwWpVUQUZKfHMG6sucF0tmpBFWmIcL24/6nUUkUEjLpRv5px7Dniuy7Jvd/rZAV8PPEQEaG338cquY1w0dThxsRp3qqvEuFgunDKcl4oq+N41Pu0jETRCncigt+FADXWNrVyiJvkeXTYjl5pTLWw4UOt1FJFBQcVdZJBbVVRBQlyM+refxrLJOSTGxfDiDjXNi4CKu8ig5pxjVVEFSyZkk5oY0qtoYSUlIY5lk3J4YftRfL6uPWxFoo+Ku8ggtuvoCUprG3WXfBAuOyeXo/VNbC097nUUEc+puIsMYquKKjCDi6YO9zrKoHfhlBHExRgvqGleRMVdZDBbVVTBrPwMhg9J8jrKoJeeHM+iCdm8sP0o/o43ItFLxV1kkCqva2RbWZ2a5Hvhshm5HKxuYNfRE15HEfGUirvIIPVyUQWAusD1gn94XnhBA9pIlFNxFxmkXiqqYFx2KoU5aV5HCRvZaYnML8hUcZeop+IuMgjVN7Xy5n5NFHM2LpuRy+6KE+yvPOl1FBHPqLiLDEKv766ktd3pevtZuHR6LgAv7qjwOImId1TcRQahVUUVZKUmMGeMJorprVEZycwcnc4L28u9jiLiGRV3kUGmtd3Hq7uPceGU4cTGqEn+bKyYMZKtpXWUHW/0OoqIJ1TcRQaZt/bXcKKpTU3yfXDpdP++0zSwEq16XdzNbISZ6UuByABZVXSUpHhNFNMX43PSmDxiiEark6gVVJE2s3gzu9fMTgBlQEFg+Q/N7EsDmE8kqrw/UUwOyQmxXscJaytm5LLhQA2VJ5q9jiIScsGegd8FXAncDHT+l/I28Nl+ziQStXYcqedIXZMGrukHK2bk4pz/5kSRaBNscf8kcIdz7inA12n5dmBSv6cSiVIv7/RPFHOhJorpsym5QxiblaKmeYlKwRb3UcDBbpbHBR4i0g9WFVUwd8wwstMSvY4S9syMFTNyWVdcRV1jq9dxREIq2OK+A1jazfIbgE39F0ckepUdb2THkXrdJd+PVkzPpc3n+NtONc1LdAn2rPs7wP+aWT4QC1xvZlOATwFXDFQ4kWjSMVHMR1Xc+83M0RnkDk3ihe1H+cSc0V7HEQmZoM7cnXNP4z9LvwT/Nfe7gInAlc65lwcunkj0WFVUwfgcTRTTn2Ji/E3zr++ppKGlzes4IiETdH9159yLzrllzrk051yKc26Jc+6lgQwnEi3qGt+fKEb616XTc2lu8/Ha7kqvo4iEjAajERkEXtt9jDafUxe4ATC/YBiZqQmaBlaiSlDX3AOD17ieXnfODe23RCJRaFVRBdlpCczK10Qx/S0uNoZLpo3gmXfLaW5rJzFOgwNJ5Av2hrovd3keD8wGrgW+36+JRKJMS5uP13dXcvk5IzVRzAC5dEYuj2w4zNriKi6cotYRiXxBFXfn3G+7W25mm4GLgJ/1ZyiRaPLm/mpONGuimIG0qDCLIYlxvLD9qIq7RIW+XnN/Ff+wtCJyllYVVZAcH8uSidleR4lYiXGxXDR1OKuKKmhr9535F0TCXF+L+01AVX8EEYlGzjle3lnBBROzSYrXteCBtGJGLrUNrbxdUuN1FJEBF+wNddv44A11BowAMoE7ByCXSFTYXlZPeV0TX79YUzQMtKWTckiKj+GFHUdZNEGtJBLZgr2h7rEuz31AJfCac25X/0YSiR6rio4SY3DRVF0HHmgpCXEsnzScF7Yf5e4rpxOjmxclggV7Q913BjqISDR6qaiCeWMzyUxN8DpKVFgxI5cXdhxly+HjzB2rbocSuTSIjYhHDtc0sOvoCd0lH0IfmTKc+FjjRU0DKxGux+JuZifMrD6YRygDi0SKVYGJYlTcQyc9OZ7FE7J5fns5zvU4LpdI2Dtds3zXgWtEpB+tKqpg4vA0CrJTvY4SVVZMz+Wbf91GUXk900elex1HZED0WNx7GrhGRPrueEMLbx+o4YtLx3sdJepcPG0E//rENl7cflTFXSKWrrmLeODV3cdo9zk1yXsgKy2RBeMyeV4TyUgEC6q4m1mCmX3HzPaYWZOZtXd+DHRIkUizqqiC4UMSmTk6w+soUWnF9Fz2HjtJ8bGTXkcRGRDBnrl/F7gF+G/8fdz/GbgPqAa+NDDRRCJTc1s7r++u5KKpI9TX2iOXzsgF0F3zErGCLe43AHc4534BtANPOef+EbgLuHigwolEonX7qjnV0q652z00Mj2ZWfkZmuNdIlawxX0EUBT4+STQ0Zb4AnBJf4cSiWTPvlvOkMQ4FhZmeR0lql1+Ti7byurYV6mmeYk8wRb3Q8CowM/FwKWBnxcCjf0dSiRSNba088L2o1x2Tq4mivHYNbPyiDF4fFOp11FE+l2wxf0J/PO2A/wP8B0zKwF+A/xqAHKJRKRVOys42dzGNbPzvI4S9YYPTWLZpBz+urmMdp8GtJHIElRxd859yzn3/cDPjwFLgJ8Bn3DO/dsA5hOJKE9uKWNkehLnj1OT/GBw/bx8jtY38cbeSq+jiPSrYLvCfWB+ROfcW865HznnnhmYWCKRp+pkM6/vqeTqWXm6S36QuGjqcDJS4vmLmuYlwgTbLH/EzJ42sxvMLGlAE4lEqGe2HqHd5/i4muQHjcS4WK6ZlceqHRUcb2jxOo5Ivwm2uH8Mf5/2XwIVZvZrM7vQzHT6IRKkJ945wrSRQ5mcO8TrKNLJdXNH09LuY+XWI15HEek3wV5zf8k591n8XeJuBzKB54HDZnbvwMUTiQz7Kk+y9fBxnbUPQjPy0pk6cih/2aimeYkcvRpb3jnX5Jz7s3PuamAWUAn804AkE4kgT20pI8bgqlmjzryyhNz1c0ezrayOXUc1g7VEhl4VdzNLNbObzex5YCswBPjegCQTiRDOOZ54p4zFE7IZMVS3rAxGV88aRVyM8ZjO3iVCBHu3/BVm9kegAvgxUAIsd85NcM7dNZABRcLdpoO1HK5p5JpZapIfrLLSErlo6nCefKeM1naf13FE+izYM/e/4J/7/VPASOfcl5xz6wYulkjkeGJLGcnxsawITFYig9P1c/OpOtnCq7uOeR1FpM/iglwv1zmni1EivdTc1s4z75ZzyfQRpCYG+89NvLB8cg7ZaYn8ZVMpl0zXFzEJb8HeLa/CLnIWXttdSV1jq4abDQNxsTF8Yk4er+46RtXJZq/jiPRJr26oE5HeeXJLGdlpCVwwIfvMK4vnrp87mjaf48ktZV5HEemTkBZ3M1thZrvNrNjMvnma9a4zM2dm80KZT6Q/1TW08redx7hy5ijiYvU9OhxMHDGEmfkZPLapFOc0mYyEr7nlhkgAACAASURBVJB94phZLHAfcBkwDfikmU3rZr0hwD8Cb4Uqm8hAeG57OS3tPg1cE2aunzuaXUdPsL1MVyMlfPW6uJtZmpmlncV7LQCKnXP7nXMtwCPA1d2s913gXqDpLN5DZNB4YksZhTmpnJOX7nUU6YUrzx1FQlwMf9l02OsoImct6OJuZl81s0NAHVBnZofN7Gu9GF8+D+j8r6U0sKzze8wG8jXbnIS7wzUNvF1Sw8dn56EpGMJLeko8l07P5al3jtDU2u51HJGzEuwgNvcCdwO/AC4OPB4Avg38MMj36u4T7r2LWmYWg3+AnDMOZ2tmt5vZRjPbWFmpeZhl8OmYhORqDVwTlq6fO5q6xlZe3lnhdRSRsxLsmfvngc87577vnHsl8Pg+8AXgtiC3UQrkd3o+Gug8DdMQYAbwmpkdAM4HVnZ3U51z7kHn3Dzn3LycnJwg314kNJxz/HVzKQsKMsnPTPE6jpyFxROyGZmepMlkJGz15pr7uz0sC3YbG4CJZjbOzBKAm4CVHS865+qcc9nOuQLnXAHwJnCVc25jLzKKeG57WT37Kk+pb3sYi40xrp0zmjf2VnK0Trf/SPgJtjD/Dvj7bpbfCfw+mA0459qALwMvAjuBR51zO8zsHjO7KsgcIoPeE1vKSIiN4YpzRnodRfrgurmj8Tn46xadvUv4CXY8zETgU2Z2Kf4zaoDzgFHAH8zspx0rOuf+saeNOOeeA57rsuzbPay7PMhsIoNGW7uPlVuPcOGU4aSnxHsdR/qgIDuVBQWZPLaxlDuXFerGSAkrwZ65TwE2A+XA2MDjaGDZVOCcwGPGAGQUCRtriquoOtmsJvkIcd3c0eyvOsXmQ7VeRxHplaDO3J1zHxnoICKR4IktZaQnx/ORKbrRMxJcfu5I7lq5g79sLGXu2Eyv44gETWNiivSTk81tvLjjKFecO5LEuFiv40g/SEuM4/JzRvLMu+U0tLR5HUckaCruIv3kpR1HaWrVcLOR5vp5o9/74iYSLlTcRfrJE1vKGD0smXljh3kdRfrReeMyGZOZoj7vElZU3EX6QUV9E2uLqzTcbAQyM66bO5p1+6o5XNPgdRyRoKi4i/SDp7cewefQXfIR6tq5ozGDxzfr7F3Cg4q7SD/46+YyZo5OpzDnbCZMlMEuLyOZRYVZPLapFJ9P87zL4KfiLtJHu4+eoKi8XmftEe76ufmU1jbyZkm111FEzkjFXaSPnthSRmyMceXMUV5HkQF06fRchiTG8ZhurJMwoOIu0gc+n+Opd8pYOjGb7LREr+PIAEpOiOVjM0fx3PZyTjS1eh1H5LRU3EX64K2SGsrrmtQkHyWunzeaplYfz20r9zqKyGmpuIv0wRNbSklNiOWSableR5EQmJ2fQWFOqvq8y6Cn4i5ylppa23l+21FWzBhJcoKGm40GZsaN8/PZeLCW7WV1XscR6ZGKu8hZWlVUwYnmNg03G2VuWjCGtMQ4frF6v9dRRHqk4i5yln63/gCjhyWzsDDL6ygSQkOT4vn0eWN49t0jHKrWiHUyOKm4i5yFrYePs+FALZ9bPI7YGA03G21uXTKOuJgYfvmGzt5lcFJxFzkLD60pIS0xjhvmjfY6inhgxNAkPj47j0c3HqbqZLPXcUQ+RMVdpJeOHG/kuW3l3Dg/nyFJ8V7HEY98Yel4Wtp9/G7dAa+jiHyIirtIL/12/QF8zvHZRQVeRxEPTRiexsVTR/Db9Qc51dzmdRyRD1BxF+mFU81t/OmtQ6yYkUt+ZorXccRjdywvpK6xlUc2HPY6isgHqLiL9MLjm0upb2rjtiXjvY4ig8CcMcNYUJDJQ2/sp7Xd53UckfeouIsEyedzPLymhFn5GcwdO8zrODJI3LF8PEfqmnh66xGvo4i8R8VdJEh/23WMA9UN3LZknNdRZBD5yOThTB4xhF+8vh/nNNe7DA4q7iJBemjNfvIykrlshsaRl/eZGbcvHc/uihO8uvuY13FEABV3kaBsL6vjzf013LJoLHGx+mcjH3TVrFGMSk/igdc1qI0MDvqUEgnCw2tKSEmI5cb5Y7yOIoNQfGwMt10wnrdLath8qNbrOCIq7iJnUlHfxNPvHuGGefmkJ2vQGuneTfP9x8cDr+3zOoqIirvImfxu/QHafI7PLS7wOooMYqmJcfzdwrGs2llB8bGTXseRKKfiLnIajS3t/OGtQ1wybQRjs1K9jiOD3C2LCkiIjeGXmg5WPKbiLnIaj28u5XhDqwatkaBkpyVyw7x8nthSRkV9k9dxJIqpuIv0wOdzPLy2hHPy0plfoEFrJDhfuGA8bT4fD68t8TqKRDEVd5EevL6nkv2Vp/j8BeMw05ztEpwxWSlcfs5I/vjmIeqbWr2OI1FKxV2kB79as5/coUlcfs5Ir6NImLljWSEnmtv4w5uHvI4iUUrFXaQbO8vrWVtczd8tGku8Bq2RXpqRl86SCdk8vLaEptZ2r+NIFNKnlkg3Hl5TQnJ8LJ9aoEFr5OzcsayQyhPNPLmlzOsoEoVU3EW6OHaiiafeOcJ1c0eTkZLgdRwJU4snZDEjbygPrt5Pu08TykhoqbiLdPG/bx6i1efToDXSJ2bGF5cWsr/qFKuKjnodR6KMirtIJ02t7fzhzYNcNGU443PSvI4jYe6yGbmMyUzhfk0HKyGm4i7SyZNbyqg+1cKtmrNd+kFcbAxfWDqerYeP81ZJjddxJIqouIsEOOcftGbayKEsHJ/ldRyJENfPHU1WagIPvK4JZSR0VNxFAt7YW8WeipPctkSD1kj/SYqP5bOLCnhtdyU7y+u9jiNRQsVdJOBXa0rIGZLIlTNHeR1FIsxnFo4lJSGWBzWhjISIirsIsLfiBKv3VHLLwrEkxOmfhfSvjJQEPrlgDCu3HmHXUZ29y8DTp5gI8PDaEhLjYvjUeWO9jiIR6ssfmcCQpDjuXrlDd87LgFNxl6hXfbKZxzeXce3c0WSmatAaGRjDUhP4p0sm8+b+Gp7dVu51HIlwKu4S9X61poSWNh+3atAaGWCfWjCGaSOH8v1nd9LQ0uZ1HIlgKu4S1Q7XNPDQmhI+MTuPCcOHeB1HIlxsjHHP1dMpr2vivleLvY4jEUzFXaLaD57fRawZ31gxxesoEiXmFWTyidl5/HJ1CQeqTnkdRyKUirtErbf2V/PstnLuXF5IbnqS13EkinzzsikkxMVwzzNFXkeRCKXiLlGp3ee455kiRqUn8YULxnsdR6LM8KFJfOWiibyy6xh/21nhdRyJQCruEpUe31zKjiP1/MtlU0hOiPU6jkShWxYVUJiTyj3PFNHU2u51HIkwKu4SdU42t/FfL+5m9pgMrtJodOKRhLgY7r5qOger/Td1ivQnFXeJOve/VkzliWa+/bFpGkNePHXBxBxWTM/l568Uc+R4o9dxJIKouEtUOVzTwC/fKOHjs/OYPWaY13FE+LcrpuJzju8/t9PrKBJBVNwlqvzghV3EGHxjxWSvo4gAkJ+ZwpeWT+DZd8tZV1zldRyJECEt7ma2wsx2m1mxmX2zm9e/bmZFZvaumf3NzDTQt/Sbt0tqePbdcu5YVsjI9GSv44i854vLxjN6WDJ3P72D1naf13EkAoSsuJtZLHAfcBkwDfikmU3rstoWYJ5z7lzgMeDeUOWTyObzOe55Zgcj05P44tJCr+OIfEBSfCzf/tg09lSc5PfrD3odRyJAKM/cFwDFzrn9zrkW4BHg6s4rOOdedc41BJ6+CYwOYT6JYI9vLmV7WT3fVNc3GaQunjaCpZNy+PGqPVSeaPY6joS5UBb3POBwp+elgWU9uQ14vrsXzOx2M9toZhsrKyv7MaJEolPNbdyrrm8yyJkZd105jaa2du59YZfXcSTMhbK4d9fnqNtJjc3sZmAe8F/dve6ce9A5N885Ny8nJ6cfI0okuv+1fVSeaOY/1PVNBrnCnDRuXTKOv2wqZfOhWq/jSBgLZXEvBfI7PR8NHOm6kpl9FPg34CrnnNqmpE9Kaxt48I39XDNrFHPU9U3CwD9cOJHhQxK5e+UOfL5uz39EziiUxX0DMNHMxplZAnATsLLzCmY2G/gF/sJ+LITZJEL94PmOrm+a9U3CQ1piHP92xVTeLa3j0Y2Hz/wLIt0IWXF3zrUBXwZeBHYCjzrndpjZPWZ2VWC1/wLSgL+Y2TtmtrKHzYmc0cYDNTzzbjlfXFrIqAx1fZPwcdXMUSwoyOTeF3dT19DqdRwJQ+ZceDf7zJs3z23cuNHrGDLI+HyOq+9bS+WJZl75P8tISYjzOpJIrxQdqedjP3uDz5w/lu9cPcPrODKImNkm59y8062jEeokIv11Sxnbyur4l8smq7BLWJo2aig3nz+W3795kKIj9V7HkTCj4i4R51RzG/e+sIuZ+RlcPfN0vS1FBrevXzyJjJQE7lq5nXbdXCe9oOIuEeeB1/dx7EQzd105jZgYdX2T8JWRksC/Xj6VDQdq+cnLe7yOI2FExV0iSmltAw+u3s/V6vomEeLaOXncMG80P3ulmJeLKryOI2FCxV0iyg9f2I0Z/Iu6vkmEMDPuuXoGM/KG8rVH3+FA1SmvI0kYUHGXiPF2SQ1Pbz3C7er6JhEmKT6W+z89l9gY44u/30RDS5vXkWSQU3GXiFB9spmvPLKF/Mxk7lg23us4Iv0uPzOFn940mz3HTvCtv24j3Lsxy8BScZew1+5zfOWRd6g+1cL9n56rrm8SsZZOyuGfLp7EU+8c4TfrDngdRwYxFXcJez95eQ9riqv47tXTmZGX7nUckQH1peUT+OjUEXz/2Z1sOFDjdRwZpFTcJay9squCn71SzA3zRnPj/DFexxEZcDExxo9unEl+Zgpf+sNmjtU3eR1JBiEVdwlbh2sa+NqftzJt5FDu0fCcEkWGJsXzwM1zOdnUxpf+sJnWdp/XkWSQUXGXsNTU2s6df9iEzzkeuHkuSfGxXkcSCanJuUP44XXnsvFgLd9/dqfXcWSQUXGXsPSdp4vYXlbPj26YxZisFK/jiHjiqpmjuHXxOH6z7gBPvVPmdRwZRFTcJew8tqmUP719iC8tL+TiaSO8jiPiqW9dPoUFBZn8y+PvsrNcE8yIn4q7hJWiI/X82xPbWDg+i69fPMnrOCKei4+N4eefns3QpHju+N9N1DVq/ndRcZcwUtfYyp1/2ERGSjw//eRs4mJ1+IoADB+SxP03z6GstpGv//kdfJpBLurp01HCgnOOf/7LVspqG7nvU3PIGZLodSSRQWXu2Ez+42PT+NuuY/z81WKv44jHVNwlLDy4ej8vFVXwrcunMq8g0+s4IoPS3y0cyzWzRvHjl/fw2u5jXscRD6m4y6C3fl81P3xhF1ecM5JbFxd4HUdk0DIz/vMT5zJ5xBC+8sg7HK5p8DqSeETFXQa1Y/VN/MOftlCQncoPrzsXM/M6ksiglpwQyy8+MxfnHF/43UZqT7V4HUk8oOIug1Zru4+//+NmTjW38cDNc0lL1IQwIsEYm5XKzz81h5KqU1z/i/WU1zV6HUlCTMVdBq17X9jFhgO1/ODac5g0YojXcUTCytJJOfz21gUcrWviuvvXU1J1yutIEkIq7jIoPb+tnF++UcItC8dy9aw8r+OIhKXzx2fxyO3n09TazvUPrGN7WZ3XkSREVNxl0NlbcYJ/fuxdZuVn8G9XTPM6jkhYm5GXzqN3LCQhNoZPPvgmb5domthooOIug8o7h49zwy/WkxQfy//79BwS4nSIivRVYU4aj925iOFDE/nMQ2/xyq4KryPJANMnpwwaq/dU8qlfvsmQpHgev3MhozKSvY4kEjFGZSTz6BcXMmnEEG7/3Sae3KKJZiKZirsMCk+9U8atv9nA2KxUHrtzIWOzUr2OJBJxstIS+eMXzmN+QSZf/fM7/HbdAa8jyQBRcRfP/XptCV955B3mjh3Gn794PsOHJHkdSSRiDUmK59efm8/F00Zw18od/M/Le3FOY9FHGhV38Yxzjv/74m6+83QRl04fwW9vXcDQpHivY4lEvKT4WO7/9ByumzuaH7+8h+88XaTJZiKMRgURT7S1+/iPp7bzp7cP88kF+XzvmnOIjdHocyKhEhcbw73Xnkt6cjwPrSmhrrGVe687l3jNthgRVNwl5Jpa2/nHP23hpaIK/uHCCXz94kkaVlbEAzExxr9fMZVhKfH835f2cKKplZ9/ag5J8bFeR5M+0lc0Can6plZuefhtXiqq4O4rp/FPl0xWYRfxkJnx5Qsn8t1rZvC3Xcf4u4ffpr6p1etY0kcq7hIyx040ceMv3mTzoVr+56ZZfHbxOK8jiUjAZ84fy09unMXmg7Xc8MB6dpbXex1J+kDFXULiYPUprrt/PQerT/HQLfM1pKzIIHT1rDwe+ux8qk42c+XP1vCjl3bT3NbudSw5CyruMuC2l9Vx7f3rONHUyh+/cD5LJ+V4HUlEerBsUg6rvraMq2aO4qevFPOxn65h86Far2NJL6m4y4Bat6+Kmx58k8S4WP5yxyJm5Wd4HUlEzmBYagI/unEWv/7cfE41t3Ht/eu45+kiGlravI4mQVJxlwHR7nP8Zm0Jn314A6MyknjszoVMGJ7mdSwR6YWPTB7Oi19bys3njeXhtSVc+pPVrC2u8jqWBEHFXfrd1sPHufq+Ndz9dBHnF2bx6BcXMjJd48SLhKMhSfF895oZ/Pn284mLieHTv3qLf3nsXeoadUf9YGbhPuzgvHnz3MaNG72OIUBdQyv/9dIu/vDWIXLSEvn2ldO44pyR6uomEiGaWtv5yct7+eUb+8lKTeB718zgkum5XseKOma2yTk377TrqLhLXznneGJLGf/fczupOdXCLYsK+PrFkxiioWRFItK20jq+8fi77Cyv54pzR3L3ldPJGZLodayoEUxx1wh10id7K07w709u562SGmblZ/Cbzy1gRl6617FEZACdMzqdlV9ezAOv7eNnrxSztriKu66cxjWz8tRSN0jozF3OSmNLOz99ZS+/XL2f1MQ4/mXFFG6an0+MxocXiSp7K07wjcffZcuh41wwMZs7lxeycHyWivwAio5meTOn0h5aLxcu4K6Lv0hZ+giu3fYy33rt12Q31HkdS0Q80m4x/G7OFfx84Y1Up2YwtWI/t258iqt2vk5iu7rP9TeDKCjuOnMPmdLaBu5eWcTLOyuYODyN710zg/PGZ3kdS0QGiabWdp56p4yH1pSwp+Ik2WmJfOb8sXz6/DFkp+mafH+JjjN3FfcB19Tazq/XHuB//rYHw/jKRydy25JxmhpSRLrlnGNNcRUPrynh1d2VJMTF8PFZedy6ZByTc4d4HS/s6YY66ZOd5fU88vYhnthSRn1TG5dMG8G3r5zG6GEpXkcTkUHMzLhgYg4XTMyh+NhJfr22hMc3l/LnjYdZMiGb25aMY9mkHN2jM4B05i4fcLK5jae3HuGRtw+xtbSOhNgYLp2Ry6cWjGFhoZrgReTs1J5q4U8bDvG7dQc5Wt/E+JxUbl08jk/MySMlQeeZvaFmeQmKc44th4/zyNuHeObdchpa2pk0Io2b5o/h47PzGJaa4HVEEYkQre0+nttWzkNrSni3tI705HhumDeaj04dwZyxw3S5Lwgq7nJatada+OuWMv684RB7Kk6SkhDLleeO4sYF+czOz1BXFhEZMM45Nh2s5aE1JawqqqDN5xiSFMfSiTksn5zDssk5DB+S5HXMQUnX3OVDfD7H+v3VPLLhMC9uP0pLu49Z+Rn84BPn8LGZo0hL1CEhIgPPzJhXkMm8gkxONLWytriK13ZX8uruYzy7rRyAc/LS+cjkHJZNHs6s/AxidY0+aDpzjwLHTjSxoaSWDQdqeGXXMQ7VNJCeHM/HZ+dx4/x8po4c6nVEERHAf0a/s/wEr+4+xmu7j7HpYC0+B8NS4lk6KYePTB7O0kk5ZEbx5UI1y0ch5xyHahp4u6SGDQdqeLukhgPVDQAkx8cyf1wm187J49LpuSTFx3qcVkTk9OoaWlm9139Gv3pPJVUnWzCDWfkZLJmQzfRRQ5k6cij5w1Ki5u57Ffco4PM5dh094S/kB2rYUFLDsRPNAGSkxDNvbCbnjctk/rhMpo8aqptVRCRs+XyO7UfqeHVXJa/sPsa20uP4AiUsNSGWKSOHMnXkEKaO9Bf8KblDIvJOfBX3CNPc1s7hmkYO1Zxi99GTbDhQw8YDNdQ3+Yd3HJmexPyCTBaM8z8m5KRFzTdZEYk+jS3t7Kk4wc7y+sDD//OJZv9nohkUZKUyJff9gj915BDyMpLD+oZh3VAXhhpa2jhU08CBqgYOVp/iQHUDh2pOcaCqgSN1jXT+LjY+J5XLzxnJgnGZzC/IZPSw8D5gRUR6Izkhlpn5GczMz3hvmXOO0trGDxT7ovJ6nt9+9L11UhNiGZmRzMj0JEalJzMyI4mR6UmMTE9mVIb/v6lhfnNxSNOb2Qrgf4BY4FfOuR90eT0R+B0wF6gGbnTOHQhlxoHi8znqGlupaWih9lQLNadaqG1oofJEMwerGzhY3cCB6lPvNal3yExNYExmCvMLhjE2azQF2SmMyUxlfHaq+p+LiHRhZuRnppCfmcIl03PfW36yuY3dR/0Ff1/lScqPN1Fe18iuoyeo7PK5CzA0KY5RgS8AIzOSGZWexPChSWQkx5OeHE9GSgLpgZ+T4mMG3YlVyIq7mcUC9wEXA6XABjNb6Zwr6rTabUCtc26Cmd0E/BC4MVQZe+Kco7nNR2NLOw2t7TS2tNHY4qOhpY3G1nYaW9o52dzG8YYPF2//f1s53tDy3rWhrkYMTWRsZirLJuVQkJ3K2KwUxmamMiYrhfTk+ND+sSIiESgtMY65YzOZOzbzQ6+1tPmoqG/iyPFGyuuaOFLX+F7xP3K8iXcOH6e2obXHbSfExfgLfqDYpyfHk57y/s9ZqQl8ZmHBAP51HxbKM/cFQLFzbj+AmT0CXA10Lu5XA3cHfn4M+LmZmQvRjQGriir4+avFNLW009DaRmOLv3A3trb3WJi7io81hqUkkJmawLCUBKbkDmVYajyZKQkMS31/eWaq/3lWaoLuWhcR8VBCXMx7Z/s9aWxpp/JEM3WNre89jje2vP+84f3l5XVN7Dp6grrGVk42tzEsJT6ii3secLjT81LgvJ7Wcc61mVkdkAVUhSJgx7evkUOTSEmIJSkhlpT4WJIT/I/3f4774PKEWFLi4xiWGk9aYtyga54REZG+SU6IZUxW7yfNamv3caq5fQASnV4oi3t3Fa/r+XAw62BmtwO3A4wZM6bvyQKWTcph2aScftueiIhEt7jYGNJTQt8FOZTvWArkd3o+GjjS0zpmFgekAzVdN+Sce9A5N885Ny8nR8VYRESks1AW9w3ARDMbZ2YJwE3Ayi7rrARuCfx8HfBKqK63i4iIRIqQNcsHrqF/GXgRf1e4h51zO8zsHmCjc24l8BDwezMrxn/GflOo8omIiESKkPZzd849BzzXZdm3O/3cBFwfykwiIiKRRgONi4iIRBgVdxERkQij4i4iIhJhVNxFREQijIq7iIhIhFFxFxERiTAq7iIiIhFGxV1ERCTCqLiLiIhEGAv3odvNrBI42I+bzCZEU8yGGe2X7mm/dE/7pXvaL93TfuleT/tlrHPutLOmhX1x729mttE5N8/rHION9kv3tF+6p/3SPe2X7mm/dK8v+0XN8iIiIhFGxV1ERCTCqLh/2INeBxiktF+6p/3SPe2X7mm/dE/7pXtnvV90zV1ERCTC6MxdREQkwkRlcTezfDN71cx2mtkOM/tKN+uYmf3UzIrN7F0zm+NF1lAKcr8sN7M6M3sn8Pi2F1lDycySzOxtM9sa2C/f6WadRDP7c+B4ecvMCkKfNLSC3C+fNbPKTsfL573I6gUzizWzLWb2TDevRd3x0uEM+yUqjxczO2Bm2wJ/88ZuXu91PYobmKiDXhvwT865zWY2BNhkZqucc0Wd1rkMmBh4nAfcH/hvJAtmvwC84Zz7mAf5vNIMXOicO2lm8cAaM3veOfdmp3VuA2qdcxPM7Cbgh8CNXoQNoWD2C8CfnXNf9iCf174C7ASGdvNaNB4vHU63XyB6j5ePOOd66uvf63oUlWfuzrly59zmwM8n8B9oeV1Wuxr4nfN7E8gws5EhjhpSQe6XqBM4Bk4GnsYHHl1vVrka+G3g58eAi8zMQhTRE0Hul6hkZqOBK4Bf9bBK1B0vENR+ke71uh5FZXHvLNAcNht4q8tLecDhTs9LiaJCd5r9ArAw0BT7vJlND2kwjwSaEt8BjgGrnHM9Hi/OuTagDsgKbcrQC2K/AFwbaEp8zMzyQxzRKz8BvgH4eng9Ko8XzrxfIDqPFwe8ZGabzOz2bl7vdT2K6uJuZmnA48BXnXP1XV/u5lei4qzkDPtlM/6hD2cCPwOeDHU+Lzjn2p1zs4DRwAIzm9Fllag8XoLYL08DBc65c4GXef9sNWKZ2ceAY865TadbrZtlEX28BLlfou54CVjsnJuDv/n9781saZfXe328RG1xD1wjfBz4g3Pur92sUgp0/tY4GjgSimxeOtN+cc7VdzTFOueeA+LNLDvEMT3jnDsOvAas6PLSe8eLmcUB6UBNSMN5qKf94pyrds41B57+Epgb4mheWAxcZWYHgEeAC83sf7usE43Hyxn3S5QeLzjnjgT+ewx4AljQZZVe16OoLO6Ba1sPATudcz/qYbWVwN8F7lI8H6hzzpWHLKQHgtkvZpbbcW3QzBbgP4aqQ5cy9Mwsx8wyAj8nAx8FdnVZbSVwS+Dn64BXXIQPIhHMfulyXfAq/PdxRDTn3Lecc6OdcwXATfiPhZu7rBZ1x0sw+yUajxczSw3cwIyZpQKXANu7rNbrehStd8svBj4DbAtcLwT4V2AMgHPuAeA54HKgGGgAPudBzlALZr9cB9xpZm1AI3BTpH8oASOB35pZLP4vM486554xs3uAjc65lfi/FP3eTfK08gAAA7VJREFUzIrxn4Hd5F3ckAlmv/yjmV2FvydGDfBZz9J6TMdL93S8MAJ4InDOFAf80Tn3gpndAWdfjzRCnYiISISJymZ5ERGRSKbiLiIiEmFU3EVERCKMiruIiEiEUXEXERGJMCruIiIiEUbFXUREJMKouItIvzKz33Q3V7eIhI6Ku4iISIRRcRcREYkwKu4iEcrMPmZmp8wsxsxGmJkzs8mB146a2ae7rP9FM6sIzFLWefkfzeypwM8rzOwNM6s1sxoze9HMpp4hx2tm9vMuyz7QdB+YEOMbZrbPzBrNbJuZdZ1sRUSCpOIuErlmAe8653zAbPwTTuw1sxH4J6vY2mX9R4EM/LO7Ae/NUnU10DE1ZyrwE/xTUi4H6oCnzSyhj1m/B9wG/D0wDfhP4BdmdkUftysSlaJ1VjiR/7+9OwapKorjOP79BUJJtFRDRUQIRUIpLkVjg00VUtDQ1BANbtUqBI05BNHSIA1Si9VgmVAZWS0tFZrQZBRRUCCURZL4bzgnvLz0WTZ13u8Dj3c4757Lvcv7vf+553EaQRvwvNIei4g5Se3ADDXbs0bElKQh4BgwnLu7SDt0DeZjrlfHSDoOfCaF/ePlXGT+AXEK6IyIR7l7Mm8p3A3cXs55zRqZw92sXO1Ab6X9otIej4jZBcb0A1ckNUfEN1LQD0TEdwBJLcA5YDewnjT7t4K8LfAytQIrgWFJ1W0qm4DX/3Bes4blcDcrkKRmoAUYy11twMXc7gCeLTL0FqlSPyTpPmmKvrPy+SDwDjiZ32eBCaDetPwcoJq+pkr71+PBA8CbmuN+1DmvmS3C4W5Wpo2kQJ2WtArYBjyXtAbYTwrn30TEjKQBUsW+DvgAPASQtBbYAXRHxIPc18HS3yMfgQ01fW3MV+UTpMcEWyJi5C/u0cwW4XA3K9N70gK608A1UnUcwHVSqN6oM7YfuAdsBa7mBXkAU8An4ISkt8Am4Dypeq9nBLgg6SDwivTDYnO+DiLii6ReoFeSgFFgNbAHmIuIy39812YGeLW8WZEi4itwBNgL3CFV8UOkRXT7IqLedPcoacq9lflV8uSQPwrsAsaBS0APqequp6/yegJMAzdrjukBzgJngJfAXeAwMLnEuc1sAYqIpY8ys/+WpD5gLdBVqcLNrGCu3M3Ktx146mA3axwOd7OC5WfYO5n/v7uZNQBPy5uZmRXGlbuZmVlhHO5mZmaFcbibmZkVxuFuZmZWGIe7mZlZYRzuZmZmhXG4m5mZFcbhbmZmVpifZq6SBIPLik0AAAAASUVORK5CYII=\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.05, c='r', linewidth=1)\n",
    "ax.plot(psi_vals, [g.pvalue for g in g_info])\n",
    "ax.set_xlabel('$\\psi$ value', fontsize=14)\n",
    "ax.set_ylabel('p value', fontsize=14)\n",
    "ax.set_title('$p$-value for each $\\psi$ value', fontsize=16);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval: (2.50, 4.50)\n"
     ]
    }
   ],
   "source": [
    "cutoff = 0.05\n",
    "ci_lo = max([g.psi for g in g_info[:len(g_info)//2] if g.pvalue < cutoff])\n",
    "ci_hi = min([g.psi for g in g_info[len(g_info)//2:] if g.pvalue < cutoff])\n",
    "\n",
    "print('95% confidence interval: ({:>0.2f}, {:>0.2f})'.format(ci_lo, ci_hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can run a finer search between 3.4 and 3.5, with steps of 0.001. That can be done using the same steps as above, but using\n",
    "\n",
    "```\n",
    "psi_vals = np.arange(3.4, 3.5, 0.001)\n",
    "```\n",
    "\n",
    "Here, we'll use automatic function optimization to find a more exact $\\psi$.\n",
    "\n",
    "The second cell below can take a while to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def just_abs_alpha(psi):\n",
    "    g_info = logit_g_info(psi, nhefs, nhefs.qsmk, X_cols, ip_censor)\n",
    "    return g_info.abs_alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 1.4009161459933774e-10\n",
       " hess_inv: array([[19.82605991]])\n",
       "      jac: array([-1.2353978e-09])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 87\n",
       "      nit: 2\n",
       "     njev: 29\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([3.4458988])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.optimize.minimize(\n",
    "    fun=just_abs_alpha,\n",
    "    x0=4.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The x value in the output is the estimated $\\psi$, which rounds to 3.446, as expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 14.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 14.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can solve for $\\psi$ directly. From Technical Point 14.2, we have\n",
    "\n",
    "$$\n",
    "\\hat\\psi = \\frac{\n",
    "    \\sum_{i=1}^{N} W_i^C Y_i\\left( A_i - \\mathrm{E}[A|L_i]\\right)\n",
    "}{\n",
    "    \\sum_{i=1}^{N} W_i^C A_i\\left( A_i - \\mathrm{E}[A|L_i]\\right)\n",
    "},\n",
    "$$\n",
    "\n",
    "where the sum is over the uncensored observations, $W^C$ is the IP weights, $Y$ is `wt82_71`, $A$ is `qsmk`, and $\\mathrm{E}[A|L_i]$ is the predicted `qsmk` from the model below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = nhefs.qsmk\n",
    "X = 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",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "glm = sm.GLM(\n",
    "    A,\n",
    "    X,\n",
    "    freq_weights=ip_censor,\n",
    "    family=sm.families.Binomial()\n",
    ")\n",
    "res = glm.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "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.4029</td> <td>    1.314</td> <td>   -1.829</td> <td> 0.067</td> <td>   -4.978</td> <td>    0.172</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sex</th>              <td>   -0.5137</td> <td>    0.150</td> <td>   -3.422</td> <td> 0.001</td> <td>   -0.808</td> <td>   -0.219</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>race</th>             <td>   -0.8609</td> <td>    0.206</td> <td>   -4.178</td> <td> 0.000</td> <td>   -1.265</td> <td>   -0.457</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>              <td>    0.1152</td> <td>    0.049</td> <td>    2.328</td> <td> 0.020</td> <td>    0.018</td> <td>    0.212</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age^2</th>            <td>   -0.0008</td> <td>    0.001</td> <td>   -1.478</td> <td> 0.140</td> <td>   -0.002</td> <td>    0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_2</th>            <td>   -0.0289</td> <td>    0.193</td> <td>   -0.150</td> <td> 0.881</td> <td>   -0.407</td> <td>    0.349</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_3</th>            <td>    0.0877</td> <td>    0.173</td> <td>    0.507</td> <td> 0.612</td> <td>   -0.252</td> <td>    0.427</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_4</th>            <td>    0.0664</td> <td>    0.266</td> <td>    0.249</td> <td> 0.803</td> <td>   -0.455</td> <td>    0.588</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>edu_5</th>            <td>    0.4711</td> <td>    0.221</td> <td>    2.130</td> <td> 0.033</td> <td>    0.038</td> <td>    0.904</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity</th>   <td>   -0.0783</td> <td>    0.015</td> <td>   -5.271</td> <td> 0.000</td> <td>   -0.107</td> <td>   -0.049</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity^2</th> <td>    0.0011</td> <td>    0.000</td> <td>    3.854</td> <td> 0.000</td> <td>    0.001</td> <td>    0.002</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs</th>         <td>   -0.0711</td> <td>    0.027</td> <td>   -2.627</td> <td> 0.009</td> <td>   -0.124</td> <td>   -0.018</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs^2</th>       <td>    0.0008</td> <td>    0.000</td> <td>    1.830</td> <td> 0.067</td> <td>-5.78e-05</td> <td>    0.002</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exercise_1</th>       <td>    0.3363</td> <td>    0.175</td> <td>    1.921</td> <td> 0.055</td> <td>   -0.007</td> <td>    0.679</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exercise_2</th>       <td>    0.3800</td> <td>    0.182</td> <td>    2.093</td> <td> 0.036</td> <td>    0.024</td> <td>    0.736</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>active_1</th>         <td>    0.0341</td> <td>    0.130</td> <td>    0.262</td> <td> 0.793</td> <td>   -0.221</td> <td>    0.289</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>active_2</th>         <td>    0.2135</td> <td>    0.206</td> <td>    1.036</td> <td> 0.300</td> <td>   -0.190</td> <td>    0.617</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71</th>             <td>   -0.0077</td> <td>    0.025</td> <td>   -0.312</td> <td> 0.755</td> <td>   -0.056</td> <td>    0.041</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71^2</th>           <td> 8.655e-05</td> <td>    0.000</td> <td>    0.574</td> <td> 0.566</td> <td>   -0.000</td> <td>    0.000</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the equation at the top of this section, $\\hat{\\psi}$ is calculated as below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_pred = res.predict(X)\n",
    "Y = nhefs.wt82_71\n",
    "\n",
    "estimate = (\n",
    "    (ip_censor * Y * (A - A_pred)).sum() /\n",
    "    (ip_censor * A * (A - A_pred)).sum()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4458988033708127"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"If $\\psi$ is D-dimensional...\"\n",
    "\n",
    "The following is a direct translation of what's in the R and Stata code examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = A - A_pred\n",
    "diff2 = ip_censor * diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   292.07618742,   5701.54685801],\n",
       "       [  5701.54685801, 153044.85432634]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lhs = np.array([\n",
    "    [\n",
    "        (A * diff2).sum(),\n",
    "        (A * nhefs.smokeintensity  * diff2).sum()\n",
    "    ],\n",
    "    [\n",
    "        (A * nhefs.smokeintensity * diff2).sum(), \n",
    "        (A * nhefs.smokeintensity**2 * diff2).sum()\n",
    "    ]\n",
    "])\n",
    "\n",
    "lhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1006.46498472],\n",
       "       [20901.0679999 ]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rhs = np.array([\n",
    "    [(Y * diff2).sum()],\n",
    "    [(Y * nhefs.smokeintensity * diff2).sum()]\n",
    "])\n",
    "\n",
    "rhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.85947039],\n",
       "       [0.03004128]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi = np.linalg.solve(lhs,rhs)\n",
    "psi"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}