chapter15.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 15"
   ]
  },
  {
   "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": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "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": [
    "Subset to uncensored rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "missing = nhefs_all.wt82.isnull()\n",
    "nhefs = nhefs_all.loc[~missing]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, instead of adding dummy variables and squared variables, we'll use formulas to specify Statsmodels models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 15.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 15.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Using the same model as in Section 13.2...\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>Intercept</th>              <td>   -1.5882</td> <td>    4.313</td> <td>   -0.368</td> <td> 0.713</td> <td>  -10.048</td> <td>    6.872</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.2]</th>      <td>    0.7904</td> <td>    0.607</td> <td>    1.302</td> <td> 0.193</td> <td>   -0.400</td> <td>    1.981</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.3]</th>      <td>    0.5563</td> <td>    0.556</td> <td>    1.000</td> <td> 0.317</td> <td>   -0.534</td> <td>    1.647</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.4]</th>      <td>    1.4916</td> <td>    0.832</td> <td>    1.792</td> <td> 0.073</td> <td>   -0.141</td> <td>    3.124</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.5]</th>      <td>   -0.1950</td> <td>    0.741</td> <td>   -0.263</td> <td> 0.793</td> <td>   -1.649</td> <td>    1.259</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(exercise)[T.1]</th>       <td>    0.2960</td> <td>    0.535</td> <td>    0.553</td> <td> 0.580</td> <td>   -0.754</td> <td>    1.346</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(exercise)[T.2]</th>       <td>    0.3539</td> <td>    0.559</td> <td>    0.633</td> <td> 0.527</td> <td>   -0.742</td> <td>    1.450</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(active)[T.1]</th>         <td>   -0.9476</td> <td>    0.410</td> <td>   -2.312</td> <td> 0.021</td> <td>   -1.752</td> <td>   -0.143</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(active)[T.2]</th>         <td>   -0.2614</td> <td>    0.685</td> <td>   -0.382</td> <td> 0.703</td> <td>   -1.604</td> <td>    1.081</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>                   <td>    2.5596</td> <td>    0.809</td> <td>    3.163</td> <td> 0.002</td> <td>    0.972</td> <td>    4.147</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:smokeintensity</th>    <td>    0.0467</td> <td>    0.035</td> <td>    1.328</td> <td> 0.184</td> <td>   -0.022</td> <td>    0.116</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sex</th>                    <td>   -1.4303</td> <td>    0.469</td> <td>   -3.050</td> <td> 0.002</td> <td>   -2.350</td> <td>   -0.510</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>race</th>                   <td>    0.5601</td> <td>    0.582</td> <td>    0.963</td> <td> 0.336</td> <td>   -0.581</td> <td>    1.701</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>                    <td>    0.3596</td> <td>    0.163</td> <td>    2.202</td> <td> 0.028</td> <td>    0.039</td> <td>    0.680</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(age ** 2)</th>            <td>   -0.0061</td> <td>    0.002</td> <td>   -3.534</td> <td> 0.000</td> <td>   -0.009</td> <td>   -0.003</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity</th>         <td>    0.0491</td> <td>    0.052</td> <td>    0.950</td> <td> 0.342</td> <td>   -0.052</td> <td>    0.151</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(smokeintensity ** 2)</th> <td>   -0.0010</td> <td>    0.001</td> <td>   -1.056</td> <td> 0.291</td> <td>   -0.003</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs</th>               <td>    0.1344</td> <td>    0.092</td> <td>    1.465</td> <td> 0.143</td> <td>   -0.046</td> <td>    0.314</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(smokeyrs ** 2)</th>       <td>   -0.0019</td> <td>    0.002</td> <td>   -1.209</td> <td> 0.227</td> <td>   -0.005</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71</th>                   <td>    0.0455</td> <td>    0.083</td> <td>    0.546</td> <td> 0.585</td> <td>   -0.118</td> <td>    0.209</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(wt71 ** 2)</th>           <td>   -0.0010</td> <td>    0.001</td> <td>   -1.840</td> <td> 0.066</td> <td>   -0.002</td> <td> 6.39e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = (\n",
    "    'wt82_71 ~ qsmk + qsmk:smokeintensity + sex + race + age + I(age**2) + C(education)'\n",
    "    '        + smokeintensity + I(smokeintensity**2) + smokeyrs + I(smokeyrs**2)'\n",
    "    '        + C(exercise) + C(active) + wt71 + I(wt71**2)'\n",
    ")\n",
    "\n",
    "ols = sm.OLS.from_formula(formula, data=nhefs) \n",
    "res = ols.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate\n",
      "alpha_1       2.56\n",
      "alpha_2       0.05\n"
     ]
    }
   ],
   "source": [
    "print('           estimate')\n",
    "print('alpha_1     {:>6.2f}'.format(res.params.qsmk))\n",
    "print('alpha_2     {:>6.2f}'.format(res.params['qsmk:smokeintensity']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain the estimates of the effect of quitting smoking, we'll use a `t_test` on the fitted model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll construct a contrast DataFrame for the 5 cigarettes/day example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start with empty DataFrame\n",
    "contrast = pd.DataFrame(\n",
    "    np.zeros((2, res.params.shape[0])),\n",
    "    columns=res.params.index\n",
    ")\n",
    "\n",
    "# modify the entries\n",
    "contrast['Intercept'] = 1\n",
    "contrast['qsmk'] = [1, 0]\n",
    "contrast['smokeintensity'] = [5, 5]\n",
    "contrast['qsmk:smokeintensity'] = [5, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>Intercept</th>\n",
       "      <th>C(education)[T.2]</th>\n",
       "      <th>C(education)[T.3]</th>\n",
       "      <th>C(education)[T.4]</th>\n",
       "      <th>C(education)[T.5]</th>\n",
       "      <th>C(exercise)[T.1]</th>\n",
       "      <th>C(exercise)[T.2]</th>\n",
       "      <th>C(active)[T.1]</th>\n",
       "      <th>C(active)[T.2]</th>\n",
       "      <th>qsmk</th>\n",
       "      <th>...</th>\n",
       "      <th>sex</th>\n",
       "      <th>race</th>\n",
       "      <th>age</th>\n",
       "      <th>I(age ** 2)</th>\n",
       "      <th>smokeintensity</th>\n",
       "      <th>I(smokeintensity ** 2)</th>\n",
       "      <th>smokeyrs</th>\n",
       "      <th>I(smokeyrs ** 2)</th>\n",
       "      <th>wt71</th>\n",
       "      <th>I(wt71 ** 2)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 21 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   Intercept  C(education)[T.2]  C(education)[T.3]  C(education)[T.4]  \\\n",
       "0          1                0.0                0.0                0.0   \n",
       "1          1                0.0                0.0                0.0   \n",
       "\n",
       "   C(education)[T.5]  C(exercise)[T.1]  C(exercise)[T.2]  C(active)[T.1]  \\\n",
       "0                0.0               0.0               0.0             0.0   \n",
       "1                0.0               0.0               0.0             0.0   \n",
       "\n",
       "   C(active)[T.2]  qsmk  ...  sex  race  age  I(age ** 2)  smokeintensity  \\\n",
       "0             0.0     1  ...  0.0   0.0  0.0          0.0               5   \n",
       "1             0.0     0  ...  0.0   0.0  0.0          0.0               5   \n",
       "\n",
       "   I(smokeintensity ** 2)  smokeyrs  I(smokeyrs ** 2)  wt71  I(wt71 ** 2)  \n",
       "0                     0.0       0.0               0.0   0.0           0.0  \n",
       "1                     0.0       0.0               0.0   0.0           0.0  \n",
       "\n",
       "[2 rows x 21 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contrast"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The effect estimate and confidence interval can be calculated with a t-test on row 0 minus row 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'statsmodels.stats.contrast.ContrastResults'>\n",
       "                             Test for Constraints                             \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "c0             2.7929      0.668      4.179      0.000       1.482       4.104\n",
       "=============================================================================="
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t_test(contrast.iloc[0] - contrast.iloc[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the effect estimate with 40 cigarettes/day, we can change a few entries in the `contrast` DataFrame and again use a t-test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'statsmodels.stats.contrast.ContrastResults'>\n",
       "                             Test for Constraints                             \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "c0             4.4261      0.848      5.221      0.000       2.763       6.089\n",
       "=============================================================================="
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contrast['smokeintensity'] = [40, 40]\n",
    "contrast['qsmk:smokeintensity'] = [40, 0]\n",
    "\n",
    "res.t_test(contrast.iloc[0] - contrast.iloc[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we don't use the product term `qsmk:smokeintensity`, we get the following model and effect estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>Intercept</th>              <td>   -1.6586</td> <td>    4.314</td> <td>   -0.384</td> <td> 0.701</td> <td>  -10.120</td> <td>    6.803</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.2]</th>      <td>    0.8185</td> <td>    0.607</td> <td>    1.349</td> <td> 0.178</td> <td>   -0.372</td> <td>    2.009</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.3]</th>      <td>    0.5715</td> <td>    0.556</td> <td>    1.028</td> <td> 0.304</td> <td>   -0.519</td> <td>    1.662</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.4]</th>      <td>    1.5085</td> <td>    0.832</td> <td>    1.812</td> <td> 0.070</td> <td>   -0.124</td> <td>    3.141</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(education)[T.5]</th>      <td>   -0.1708</td> <td>    0.741</td> <td>   -0.230</td> <td> 0.818</td> <td>   -1.625</td> <td>    1.283</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(exercise)[T.1]</th>       <td>    0.3207</td> <td>    0.535</td> <td>    0.599</td> <td> 0.549</td> <td>   -0.729</td> <td>    1.370</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(exercise)[T.2]</th>       <td>    0.3629</td> <td>    0.559</td> <td>    0.649</td> <td> 0.516</td> <td>   -0.734</td> <td>    1.459</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(active)[T.1]</th>         <td>   -0.9430</td> <td>    0.410</td> <td>   -2.300</td> <td> 0.022</td> <td>   -1.747</td> <td>   -0.139</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(active)[T.2]</th>         <td>   -0.2580</td> <td>    0.685</td> <td>   -0.377</td> <td> 0.706</td> <td>   -1.601</td> <td>    1.085</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>                   <td>    3.4626</td> <td>    0.438</td> <td>    7.897</td> <td> 0.000</td> <td>    2.603</td> <td>    4.323</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sex</th>                    <td>   -1.4650</td> <td>    0.468</td> <td>   -3.128</td> <td> 0.002</td> <td>   -2.384</td> <td>   -0.546</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>race</th>                   <td>    0.5864</td> <td>    0.582</td> <td>    1.008</td> <td> 0.314</td> <td>   -0.555</td> <td>    1.727</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>                    <td>    0.3627</td> <td>    0.163</td> <td>    2.220</td> <td> 0.027</td> <td>    0.042</td> <td>    0.683</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(age ** 2)</th>            <td>   -0.0061</td> <td>    0.002</td> <td>   -3.555</td> <td> 0.000</td> <td>   -0.010</td> <td>   -0.003</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeintensity</th>         <td>    0.0652</td> <td>    0.050</td> <td>    1.295</td> <td> 0.196</td> <td>   -0.034</td> <td>    0.164</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(smokeintensity ** 2)</th> <td>   -0.0010</td> <td>    0.001</td> <td>   -1.117</td> <td> 0.264</td> <td>   -0.003</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>smokeyrs</th>               <td>    0.1334</td> <td>    0.092</td> <td>    1.454</td> <td> 0.146</td> <td>   -0.047</td> <td>    0.313</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(smokeyrs ** 2)</th>       <td>   -0.0018</td> <td>    0.002</td> <td>   -1.183</td> <td> 0.237</td> <td>   -0.005</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>wt71</th>                   <td>    0.0374</td> <td>    0.083</td> <td>    0.449</td> <td> 0.653</td> <td>   -0.126</td> <td>    0.200</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>I(wt71 ** 2)</th>           <td>   -0.0009</td> <td>    0.001</td> <td>   -1.749</td> <td> 0.080</td> <td>   -0.002</td> <td>    0.000</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = (\n",
    "    'wt82_71 ~ qsmk + sex + race + age + I(age**2) + C(education)'\n",
    "    '        + smokeintensity + I(smokeintensity**2) + smokeyrs + I(smokeyrs**2)'\n",
    "    '        + C(exercise) + C(active) + wt71 + I(wt71**2)'\n",
    ")  # no qsmk_x_smokeintensity\n",
    "\n",
    "ols = sm.OLS.from_formula(formula, data=nhefs)\n",
    "res = ols.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate   95% C.I.\n",
      "alpha_1       3.5    (2.6, 4.3)\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('alpha_1     {:>5.1f}    ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 15.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 15.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate propensity score, we fit a logistic model for `qsmk` conditional on $L$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "formula = (\n",
    "    'qsmk ~ sex + race + age + I(age**2) + C(education)'\n",
    "    '     + smokeintensity + I(smokeintensity**2) + smokeyrs + I(smokeyrs**2)'\n",
    "    '     + C(exercise) + C(active) + wt71 + I(wt71**2)'\n",
    ")\n",
    "\n",
    "model = sm.Logit.from_formula(formula, data=nhefs_all) \n",
    "res = model.fit(disp=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then propensity is the predicted values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "propensity = res.predict(nhefs_all)\n",
    "nhefs_all['propensity'] = propensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lowest and highest propensity scores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ranked = nhefs_all[['seqn', 'propensity']].sort_values('propensity').reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "seqn          22941.000000\n",
       "propensity        0.052981\n",
       "Name: 0, dtype: float64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ranked.loc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "seqn          24949.000000\n",
       "propensity        0.793205\n",
       "Name: 1628, dtype: float64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ranked.loc[ranked.shape[0] - 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll attempt to recreate Figure 15.1, pg 45\n",
    "\n",
    "First, we'll split the propensities based on whether the subject quit smoking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "propensity0 = propensity[nhefs_all.qsmk == 0]\n",
    "propensity1 = propensity[nhefs_all.qsmk == 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the bins are spaced every 0.05 (except at the right end), with the first bin starting at 0.025."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.arange(0.025, 0.85, 0.05)\n",
    "\n",
    "top0, _ = np.histogram(propensity0, bins=bins)\n",
    "top1, _ = np.histogram(propensity1, bins=bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAF7CAYAAACtnpMfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3X+czXX+///bg0FFkWXaYUiFjJnRMFPqm0TtIFrWj1qyRai2t4pKsWk/Uitq+0Fp+7VEdhu7KpFkS7Kk7CRGTT9EsWuw+bGMX/kx4/n94xzTMQaHOee8zjlzv14ur8uc83y9zuv5eJ4zM6/HeT1fr+fTnHOIiIhIxVbJ6wBERETEe0oIRERERAmBiIiIKCEQERERlBCIiIgISghEREQEJQQiIiKCEgIRERFBCYGIiIighEBERESABK8DiKQ6deq4Ro0aeR2GiIhIxHz22WdbnXN1T7RdhUoIGjVqxLJly7wOQ0REJGLM7N/BbKcuAxEREVFCICIiIkoIREREBCUEIiIighICERERQQmBiIiIoIRAREREUEIgIiIiKCEQERERlBCIR9avX0/79u1JSUkhNTWVCRMmAPDQQw9Rv359MjIyyMjIYO7cuQAcPHiQfv36kZ6eTkpKCmPHjvUyfBGRuFOhhi6W6JGQkMCTTz5Jq1at2LVrF5mZmWRnZwNw9913M2zYsCO2nzFjBvv37+eLL75g7969NG/enD59+qC5KUREQkMJgXgiKSmJpKQkAM4880xSUlLYsGHDMbc3M/bs2UNRURE//vgjVatW5ayzzopUuCIicU9dBuK5devWsWLFClq3bg3AxIkTadGiBQMGDGD79u0A9OrVi+rVq5OUlETDhg0ZNmwYtWvX9jJsEZG4ooRAPLV792569uzJ+PHjOeuss7j99tv57rvvyMvLIykpiXvvvReA3NxcKleuzMaNG1m7di1PPvkk33//vcfRi4jEDyUE4pmDBw/Ss2dP+vbtS48ePQA455xzqFy5MpUqVeKWW24hNzcXgNdee41OnTpRpUoVEhMTufzyyzWVtYhICCkhEE845xg4cCApKSncc889JeWbNm0qeTxz5kzS0tIAaNiwIQsWLMA5x549e1i6dCnNmjWLeNwiIvFKFxWKJ5YsWcK0adNIT08nIyMDgEcffZScnBzy8vIwMxo1asSLL74IwODBg7n55ptJS0vDOcfNN99MixYtvGyCiEhcMeec1zFETFZWltNpZhERqUjM7DPnXNaJtlOXgYiICL5uSjPjm2++Ccn+9u/fz69//WsaN25M69atWbduXUj2Gy5KCERERICcnBzatGnD9OnTQ7K/SZMmcfbZZ7NmzRruvvtuhg8fHpL9hosSAhERqfB2797NkiVLmDRpUsgSglmzZtGvXz/AN5bKBx98QDR30+uiQhERqfDeeustOnXqRNOmTalduzbLly+nVatWR213xRVXsGvXrqPKn3jiCX7xi18cUbZhwwYaNGgA+IZrr1mzJtu2baNOnTrhaUQ5KSEQEZEKLycnh6FDhwLQu3dvcnJyykwIFi9eHPQ+yzobYGanHmSYKSEQEZEKbdu2bSxYsID8/HzMjOLiYsyMxx9//KgD+MmcIUhOTmb9+vUkJydTVFREYWFhVA+5roRAREQqtNdff52bbrqpZNwTgCuvvJKPPvqIK6644ohtT+YMQdeuXZk6dSqXXXYZr7/+OldddVVUnyHQRYUiIlKh5eTk0L179yPKevbsyWuvvVau/Q4cOJBt27bRuHFjnnrqKcaNG1eu/YWbBiYSERGJY8EOTKQuA4m48ePHU1hYGLL91axZs+RiIBEROTVKCCTiCgsLGTVqVMj2N3r06JDtS0SkooqaawjMrIGZfWhmX5vZl2Y2xF/+kJltMLM8/9I54DW/M7M1ZrbKzDp6F72IiEhsi6YzBEXAvc655WZ2JvCZmb3vX/e0c+6JwI3NrDnQG0gF6gHzzaypc644olGLiIjEgag5Q+Cc2+ScW+5/vAv4Gqh/nJd0A6Y75/Y759YCa4BLwh+piIhI/ImahCCQmTUCWgL/8hfdYWafm9lkMzvbX1YfWB/wsgKOn0CIiIjIMURdQmBmNYA3gKHOuZ3A88AFQAawCXjy8KZlvPyoeyjN7FYzW2Zmy7Zs2RKmqEVERGJbVCUEZlYFXzLwV+fcmwDOuR+cc8XOuUPAy/zULVAANAh4eTKwsfQ+nXMvOeeynHNZdevWDW8DREREYlTUJATmG89xEvC1c+6pgPKkgM26A/n+x7OB3mZWzczOA5oAuZGKV0REJJ5E010GlwM3Al+YWZ6/7AGgj5ll4OsOWAfcBuCc+9LM/g58he8OhcG6w0BEROTURE1C4Jz7iLKvC5h7nNeMAcaELSgREZEKImq6DERERMQ7SghERERECYGIiIgoIRARERGUEIiIiAhKCERERAQlBCIiIoISAhEREUEJgYiIiKCEQERERFBCICIiIighEBEREZQQiIiICEoIREREBCUEIiIighICiXPr16+nffv2pKSkkJqayoQJEwD43//+R3Z2Nk2aNCE7O5vt27cDMGvWLFq0aEFGRgZZWVl89NFHXoYvIhIxSggkriUkJPDkk0/y9ddfs3TpUp577jm++uorxo0bx9VXX83q1au5+uqrGTduHABXX301K1euJC8vj8mTJzNo0CCPWyAiEhlKCCSuJSUl0apVKwDOPPNMUlJS2LBhA7NmzaJfv34A9OvXj7feeguAGjVqYGYA7Nmzp+SxiEi8U0IgFca6detYsWIFrVu35ocffiApKQnwJQ2bN28u2W7mzJk0a9aMLl26MHnyZK/CFRGJKCUEUiHs3r2bnj17Mn78eM4666zjbtu9e3e++eYb3nrrLX7/+99HKEIREW8pIZC4d/DgQXr27Enfvn3p0aMHAOeccw6bNm0CYNOmTSQmJh71urZt2/Ldd9+xdevWiMYrIuIFJQQS15xzDBw4kJSUFO65556S8q5duzJ16lQApk6dSrdu3QBYs2YNzjkAli9fzoEDB/jZz34W+cBFRCIswesARMJpyZIlTJs2jfT0dDIyMgB49NFHGTFiBNdffz2TJk2iYcOGzJgxA4A33niDV199lSpVqnD66afzt7/9TRcWikiFoIRA4lqbNm1KvvGX9sEHHxxVNnz4cIYPHx7usEREoo66DOQoAwYMIDExkbS0tJKylStXctlll5Gens4vf/lLdu7cCcD7779PZmYm6enpZGZmsmDBAq/CFhGRclBCIEfp378/8+bNO6Js0KBBjBs3ji+++ILu3bvzxz/+EYA6derw9ttv88UXXzB16lRuvPFGL0IWEZFyUkIgR2nbti21a9c+omzVqlW0bdsWgOzsbN544w0AWrZsSb169QBITU1l37597N+/P7IBi4hIuSkhkKCkpaUxe/ZsAGbMmMH69euP2uaNN96gZcuWVKtWLdLhiYhIOSkhkKBMnjyZ5557jszMTHbt2kXVqlWPWP/ll18yfPhwXnzxRY8iFBGR8tBdBhKUZs2a8d577wHw7bff8s4775SsKygooHv37rz66qtccMEFXoUoIiLloDMEEpTDY/0fOnSIP/zhD/z2t78FYMeOHXTp0oWxY8dy+eWXexmiiIiUQ9QkBGbWwMw+NLOvzexLMxviL69tZu+b2Wr/z7P95WZmz5jZGjP73MxaeduC+NGnTx8uu+wyVq1aRXJyMpMmTSInJ4emTZvSrFkz6tWrx8033wzAxIkTWbNmDY888ggZGRlkZGQcMVGQiIjEhmjqMigC7nXOLTezM4HPzOx9oD/wgXNunJmNAEYAw4FrgCb+pTXwvP+nlFNOTk6Z5UOGDDmq7MEHH+TBBx8Md0giIhJmUZMQOOc2AZv8j3eZ2ddAfaAb0M6/2VRgIb6EoBvwqvMNQ7fUzGqZWZJ/P1LBjR8/nsLCwpDus2bNmgwdOjSk+xQRiRZRkxAEMrNGQEvgX8A5hw/yzrlNZnZ4Wrr6QOC9bwX+MiUEQmFhIaNGjQrpPkePHh3S/YmIRJOouYbgMDOrAbwBDHXO7TzepmWUHTVovZndambLzGzZli1bQhWmiIhIXImqhMDMquBLBv7qnHvTX/yDmSX51ycBh69YKwAaBLw8GdhYep/OuZecc1nOuay6deuGL3gREZEYFjUJgfnmmJ0EfO2ceypg1Wygn/9xP2BWQPlN/rsNLgUKdf2AiIjIqYmmawguB24EvjCzPH/ZA8A44O9mNhD4D3Cdf91coDOwBtgL3BzZcEVEROJH1CQEzrmPKPu6AICry9jeAYPDGpSIiEgFETVdBiKxbMCAASQmJpKWllZS9utf/7pksKZGjRqRkZEBwPvvv09mZibp6elkZmayYMECr8IWESkRNWcIRGJZ//79ueOOO7jppptKyv72t7+VPL733nupWbMmAHXq1OHtt9+mXr165Ofn07FjRzZs2BDxmEVEAikhkBIazOfUtW3blnXr1pW5zjnH3//+95IzAS1btixZl5qayr59+9i/f7+mjRYRTykhkBIazCc8Fi9ezDnnnEOTJk2OWvfGG2/QsmVLJQMi4jklBCJhlpOTQ58+fY4q//LLLxk+fHjJtNIiIl5SQiASRkVFRbz55pt89tlnR5QXFBTQvXt3Xn31VS644AKPohMR+YnuMhAJo/nz59OsWTOSk5NLynbs2EGXLl0YO3Ysl19+uYfRiYj8RAmBSAj06dOHyy67jFWrVpGcnMykSZMAmD59+lHdBRMnTmTNmjU88sgjJbclbt68uazdiohEjLoMREIgJyenzPIpU6YcVfbggw/y4IMPhjkiEZGTozMEIiIiooRARERElBCIiIgISghEREQEJQQiIiKCEgIRERFBCYGIiIighEBERETQwEQi5aIpo0UkXighECkHTRktIvFCXQYiIiKihEBERESUEIiIiAhKCERERAQlBCIiIoISAhEREUEJgYiIiKCEQERERFBCICIiIighEBEREZQQiIiICEoIREREhChLCMxsspltNrP8gLKHzGyDmeX5l84B635nZmvMbJWZdfQmahERkdgXVQkBMAXoVEb50865DP8yF8DMmgO9gVT/a/5kZpUjFqmIiEgciaqEwDm3CPhfkJt3A6Y75/Y759YCa4BLwhaciIhIHIuqhOA47jCzz/1dCmf7y+oD6wO2KfCXiYiIyEmKhYTgeeACIAPYBDzpL7cytnWlC8zsVjNbZmbLtmzZEr4oRUREYljUJwTOuR+cc8XOuUPAy/zULVAANAjYNBnYWMbrX3LOZTnnsurWrRv+gEVERGJQ1CcEZpYU8LQ7cPgOhNlAbzOrZmbnAU2A3EjHJyIiEg8SvA4gkJnlAO2AOmZWAIwC2plZBr7ugHXAbQDOuS/N7O/AV0ARMNg5V+xF3CIiIrEuqs4QOOf6OOeSnHNVnHPJzrlJzrkbnXPpzrkWzrmuzrlNAduPcc5d4Jy70Dn3rpexR8KAAQNITEwkLS2tpOy+++6jWbNmtGjRgu7du7Njxw4AcnNzycjIICMjg4suuoiZM2d6FbaIiMSAU04IzKyxmZ0WymDk+Pr378+8efOOKMvOziY/P5/PP/+cpk2bMnbsWADS0tJYtmwZeXl5zJs3j9tuu42ioiIvwhYRkRgQVEJgZo+aWT//YzOz94FvgU1m1jqcAcpP2rZtS+3atY8o69ChAwkJvp6fSy+9lIKCAgDOOOOMkvJ9+/ZhVtZNGSIiIj7BniHoC6zyP74G3y2AlwKvAuPCEJecgsmTJ3PNNdeUPP/Xv/5Famoq6enpvPDCCyUJgoiISGnBJgTn4LvND6Az8HfnXC7wLNAyHIHJyRkzZgwJCQn07du3pKx169Z8+eWXfPrpp4wdO5Z9+/Z5GKGIiESzYBOCbcC5/scdgAX+xwmUPUCQRNDUqVOZM2cOf/3rX8vsGkhJSaF69erk5+eX8WoREZHgE4I3gNf81w7UBg5f2ZaBbw4B8ci8efN47LHHmD17NmeccUZJ+dq1a0suIvz3v//NqlWraNSokUdRiohItAu2U/ke4N9AQ+B+59wef3kSvqGFJQL69OnDwoUL2bp1K8nJyYwePZqxY8eyf/9+srOzAd+FhS+88AIfffQR48aNo0qVKlSqVIk//elP1KlTx+MWiIhItAo2IaiHbwriQ6XKx3Pk8MESRjk5OUeVDRw4sMxtb7zxRm688cZwhyQiInEi2C6DtUBZXy9r+9eJiIhIDAs2ITDKmEkQqAHo0nUREZEYd9wuAzN7xv/QAWPNbG/A6sr4Zh7MC1NsIiIiEiEnuoYg3f/TgBTgQMC6A8By4IkwxCUiIiIRdNyEwDnXHsDMXgGGOOd2RiQqERERiahgryF4ADirdKGZJZvZOaENSURERCIt2ITgVXxzGJTWEZgWunBERETEC8EmBBcDi8ooXwxkhS4cERER8UKwAxMlANXKKD/tGOUSQuPHj6ewsDCk+6xZsyZDhw4N6T5FRCR2BZsQ/Au43b8EGgx8GtKI5CiFhYWMGjUqpPscPXp0SPcnIiKxLdgug5FAPzP72Mwe8S9LgBvxXXAoIhE2YcIE0tLSSE1NZfz48QDk5eVx6aWXkpGRQVZWFrm5uR5HKSKxIqiEwDm3FLgM3zDFPYCe/seXOec+Dl94IlKW/Px8Xn75ZXJzc1m5ciVz5sxh9erV3H///YwaNYq8vDwefvhh7r//fq9DFZEYEWyXAc65lUDfMMYiIkH6+uuvufTSS0umvL7yyiuZOXMmZsbOnb7hQgoLC6lXr56XYYpIDAk6IfCPN3AjcD7w/5xzW83scmCjc04THIlEUFpaGiNHjmTbtm2cfvrpzJ07l6ysLMaPH0/Hjh0ZNmwYhw4d4uOPdQJPRIITVJeBmWUCq/CdIRjET4MUZQNjwhOaiBxLSkoKw4cPJzs7m06dOnHRRReRkJDA888/z9NPP8369et5+umnjzk9tohIacFeVPgEMME51xLYH1D+D+DykEclIic0cOBAli9fzqJFi6hduzZNmjRh6tSp9OjRA4DrrrtOFxWKSNCCTQgygalllG8CNHSxiAc2b94MwH/+8x/efPNN+vTpQ7169fjnP/8JwIIFC2jSpImXIYpIDAn2GoIfgbPLKG8GbA5dOCISrJ49e7Jt2zaqVKnCc889x9lnn83LL7/MkCFDKCoq4rTTTuOll17yOkwRiRHBJgSzgFFmdp3/uTOzRsBjwBthiEtETmDx4sVHlbVp04bPPvvMg2hEJNYF22UwDKgNbAHOAD4C1gA7gAfDE5qIiIhESlBnCJxzO4E2ZnYV0ApfIrHcOTc/nMGJiIhIZAQ9DgGAc24BsCBMsYiIiIhHjpkQmNk9wJ+cc/v8j49nN5CvYYxFRERi0/HOENyJ71bDff7Hx1MNSDSz8c65YacajJlNBq4FNjvn0vxltYG/AY2AdcD1zrntZmbABKAzsBfo75xbfqp1i4iIVGTHvKjQOXeec25bwOPjLfWAa4B+5YxnCtCpVNkI4APnXBPgA/9z/PU18S+3As+Xs24REZEKK9i7DILxEfCH8uzAObcI+F+p4m78NCjSVOBXAeWvOp+lQC0zSypP/SIiIhVV0AmBmf3KzBaZ2Vb/stjMuh9e75z70Tk3IQwxnuOc2+SvYxOQ6C+vD6wP2K7AXyYiIiInKdjJje7F14+/Crjfv3wDvGZmp3zNQDlZGWXuqI3MbjWzZWa2bMuWLREIS0REJPYEe9vhMOAO59zLAWWTzSwXeBjf5Efh8oOZJTnnNvm7BA4PlVwANAjYLhnYWPrFzrmXgJcAsrKyjkoYRGLB+PHjKSwsDNn+atasydChQ0O2PxGJfcEmBDWAD8so/9C/Lpxm47tYcZz/56yA8jvMbDrQGig83LUgEm8KCwsZNWpUyPY3evTokO1LROJDsNcQvAX0KqO8J74Dc0iYWQ7wCXChmRWY2UB8iUC2ma0Gsv3PAeYC3+MbQvll4P9CFYeIiEhFc6KBiQ5bA4wws/b4DtgAl/qXp0IVjHOuzzFWXV3Gtg4YHKq6RUREKrITDUwUaDvQ1L8ElvXHdx2BiIiIxKhjJgTOufMiGYiIiIh4J5QDE4mIiEiMCuouAzN75njrnXN3hSYcERER8UKwtx2ml3peBWjmf70mFBIREYlxQSUEzrn2pcvM7DRgErA41EGJiIhIZJ3yNQTOuX3AGGBk6MIRERERL5T3osK6hH+kQhEREQmzYC8qvKd0EZAE9MU3YqCIxKkdO3YwaNAg8vPzMTMmT55MQUEBDz30EF9//TW5ublkZWV5HaaIlFOwFxWWHqToELAFeAUYG9KIRCSqDBkyhE6dOvH6669z4MAB9u7dS61atXjzzTe57bbbvA5PREIk2IsKNUiRSAW0c+dOFi1axJQpUwCoWrUqVatWpVatWt4GJiIhd0rXEJhZgpnp2gGROPf9999Tt25dbr75Zlq2bMmgQYPYs2eP12GJSBgcNyEws6vN7PpSZSOA3cAOM5tnZvqqIBKnioqKWL58ObfffjsrVqygevXqjBs37sQvFJGYc6IzBCOA5MNPzOwS4FFgGnA/cBG67VAkbiUnJ5OcnEzr1q0B6NWrF8uXaywykXh0ooQgHfhnwPPrgI+dc7c4554C7gK6his4EfHWz3/+cxo0aMCqVasA+OCDD2jevLnHUYlIOJwoIagFbA54fjkwL+D5p0D9UAclItHj2WefpW/fvrRo0YK8vDweeOABZs6cSXJyMp988gldunShY8eOXocpIuV0orsMNgEXAOvNrBrQEvh9wPozgf1hik1EokBGRgbLli07oqx79+50797do4hEJBxOdIbgXeBxM7sKeAzYw5FzF7QA1oQpNhEREYmQE50h+H/Am8B8fHcW9HPOHQhYPwB4P0yxiYiISIQcNyFwzm0F2ppZTWC3c6641CbX4UsUREREJIYFO1Jh4THK/xfacERERMQL5Z3tUEREROKAEgIRERFRQiAiIiJKCEQkihQXF9OyZUuuvfZaACZOnEjjxo0xM7Zu3epxdCLxTQmBiESNCRMmkJKSUvL88ssvZ/78+Zx77rkeRiVSMQR1l8GxmNnXQBPnXLn2IyLRYfz48RQWlnlT0SmpWbMmQ4cODWrbgoIC3nnnHUaOHMlTTz0FQMuWLUMWi4gcX3kP5M8BPwtFICLivcLCQkaNGhWy/Y0ePTrobYcOHcrjjz/Orl27Qla/iASvXF0GzrmJzrng/+JFRMowZ84cEhMTyczM9DoUkQrrpM4QmNlpQGPAAd855/aFJSoRqVCWLFnC7NmzmTt3Lvv27WPnzp385je/4S9/+YvXoYlUGEGdITCzBDP7I7AdWAl8AWw3s8fNrEo4AxSR+Dd27FgKCgpYt24d06dP56qrrlIyIBJhwXYZPA78Bvgt0BRoAtwO3AiMDU9oRzKzdWb2hZnlmdkyf1ltM3vfzFb7f54diVhEJDKeeeYZkpOTKSgooEWLFgwaNMjrkETiVrBdBjcAA5xzcwPKvjOzLcCfgWEhj6xs7f0TLh02AvjAOTfOzEb4nw+PUCwiEgbt2rWjXbt2ANx1113cdddd3gYkUkEEe4agJvBdGeXfAbVCF85J6wZM9T+eCvzKw1hERERiVrAJwUqgrDR9CJAXunCOywHvmdlnZnarv+wc59wmAP/PxAjFIiIiEleC7TK4H5hrZtnAJ/gOzpcB9YBrwhRbaZc75zaaWSLwvpl9E8yL/MnDrQANGzYMZ3wiIiIxK6gzBM65RfguJpwB1ADO8j++0Dn3UfjCOyKGjf6fm4GZwCXAD2aWBOD/ubmM173knMtyzmXVrVs3EqGKiIjEnKDHIfAfkEeGMZZjMrPqQCXn3C7/4w7Aw8BsoB8wzv9zlhfxiYiIxLrjniHw39Z3wiUCcZ4DfGRmK4Fc4B3n3Dx8iUC2ma0Gsv3PRUTKtH79etq3b09KSgqpqalMmDABgN///ve0aNGCjIwMOnTowMaNGz2OVCTyTnSGYCu+6wWOxwWxn3Jxzn0PXFRG+Tbg6nDWLSLxIyEhgSeffJJWrVqxa9cuMjMzyc7O5r777uORRx4BfGMfPPzww7zwwgseRysSWSc6kLc/zrpO+O4yKApdOCIi4ZOUlERSUhIAZ555JikpKWzYsIHmzZuXbLNnzx7MzKsQRTxz3ITAOffP0mVm1gp4DGgLvAg8Ep7QRETCZ926daxYsYLWrVsDMHLkSF599VVq1qzJhx9+6HF0IpEX9Kl+MzsPGANcB7wJNHfOlTVYkYjIMY0fP57CwsKQ7rNmzZoMHTo06O13795Nz549GT9+PGeddRYAY8aMYcyYMYwdO5aJEyee1NTNIvHghAmBmf0M+H/45jFYAlzmnFsW7sBEJD4VFhYyatSokO7zZA7eBw8epGfPnvTt25cePXoctf6GG26gS5cuSgikwjnRXQYP4Bue+Eqgm3PuKiUDIhKrnHMMHDiQlJQU7rnnnpLy1atXlzyePXs2zZo18yI8EU+d6AzBH4AfgQLg/8zs/8rayDnXNdSBiYiE2pIlS5g2bRrp6elkZGQA8OijjzJp0iRWrVpFpUqVOPfcc3WHgVRIJ0oIXuXEtx1KKQMGDGDOnDkkJiaSn58P+O5znjVrFpUqVSIxMZEpU6ZQr149jyMVqVjatGmDc0f/S+vcubMH0YhElxPdZdA/QnHElf79+3PHHXdw0003lZTpPmcREYlmwc52KCehbdu21K595ACOh69kBt3nLCIi0UcJQQSNHDmSBg0a8Ne//pWHH37Y63BEJAImTJhAWloaqampjB8/3utwRI5JCUEEjRkzhvXr19O3b18mTpzodTgiEmb5+fm8/PLL5ObmsnLlSubMmXPEHQ0i0UQJgQduuOEG3njjDa/DEJEw+/rrr7n00ks544wzSEhI4Morr2TmzJlehyVSJiUEEaL7nEUqnrS0NBYtWsS2bdvYu3cvc+fOZf369WGpq1GjRiW3U2ZlZYWlDolvYZ2lsKLq06cPCxcuZOvWrSQnJzN69Gjmzp2r+5xFKpiUlBSGDx9OdnY2NWrU4KKLLiIhIXz/dj/88EPq1KkTtv1LfFNCEAY5OTlHlQ0cONCDSETEawMHDiyIs4O0AAAgAElEQVT5+3/ggQdITk72OKJTt2/fPtq2bcv+/fspKiqiV69eGuI5jqjLQEQkjDZv3gzAf/7zH95880369OkTlnrMjA4dOpCZmclLL70UljqqVavGggULWLlyJXl5ecybN4+lS5eGpS6JPJ0hEJG4FOpZFU92RsXDevbsybZt26hSpQrPPfccZ599dshiCrRkyRLq1avH5s2byc7OplmzZrRt2zakdZgZNWrUAHyTRB08eDAsY6rMmzePIUOGUFxczKBBgxgxYkTI64jHespLCYGIxKVQz6p4qqfGFy9eHLIYjufwUOiJiYl0796d3NzckCcEAMXFxWRmZrJmzRoGDx5M69atQ77/wYMH8/7775OcnMzFF19M165dad68ueoJM3UZiIjEuD179rBr166Sx++99x5paWlhqaty5crk5eVRUFBAbm5uyXwtoZKbm0vjxo05//zzqVq1Kr1792bWrFkhrSMe6wkFJQQiIjHuhx9+oE2bNlx00UVccskldOnShU6dOoW1zlq1atGuXTvmzZsX0v1u2LCBBg0alDxPTk5mw4YNIa0jHusJBXUZiIjEuPPPP5+VK1eGvZ4tW7ZQpUoVatWqxY8//sj8+fMZPnx4SOsoazbKcFynEG/1hIISAhERCcqmTZvo168fxcXFHDp0iOuvv55rr702pHUkJycfMXhTQUFBWKaKj7d6QkEJgYiIBKVFixasWLEirHVcfPHFrF69mrVr11K/fn2mT5/Oa6+9pnoiQAmBiIhEjYSEBCZOnEjHjh0pLi5mwIABpKamqp4IUEJQTtFyr7OISLzo3LkznTt3Vj0RpoSgnKLlXmcREZHy0G2HIiIiojMEIiLloW5DiRdKCEREykHdhhIv1GUgIiIisX+GwMw6AROAysCfnXPjPA5JRCSkQt0tAeqakKPFdEJgZpWB54BsoAD41MxmO+e+8jYyEZHQCXW3BKhrQo4W0wkBcAmwxjn3PYCZTQe6AUoIREROks5EVGyxnhDUB9YHPC8AQjs5t4hIBaEzERWblTUTU6wws+uAjs65Qf7nNwKXOOfuDNjmVuBWgKSkpMzbbrvNk1hFRES88NBDD33mnMs60XaxnhBcBjzknOvof/47AOfc2LK2z8rKcsuWLQtpDOHIfktn6OHKsOOpnnhqS6TqKeuboD6bk6sjUvXE8ntWVj1PPPEEe/bsCdn+q1evzrBhw0K2v5MR6rZA6NtjZkElBLHeZfAp0MTMzgM2AL2BGyIZQPXq1UP+iy0iEs+8OniHQzy1JaYTAudckZndAfwD322Hk51zX0Yyhnj6ZRCRii3UX3AO71NiQ0wnBADOubnAXK/jEBGJdfqCU7FppEIRERFRQiAiIiJKCGJCOPrg1K8nIiKBYv4agopA/XoiIqE3b948hgwZQnFxMYMGDWLEiBFeh+QpnSEQEZGoMmDAABITE0lLSwtbHcXFxQwePJh3332Xr776ipycHL76KvSj3keiLaGihCDM1q9fT/v27UlJSSE1NZUJEyZ4HZKIxJiK1m3Yv39/5s2bF9Y6cnNzady4Meeffz5Vq1ald+/ezJo1K+T1RKItoaIugzBLSEjgySefpFWrVuzatYvMzEyys7Np3ry516GJSAhEYnCyitZt2LZtW9atWxfWOjZs2ECDBg1KnicnJ/Ovf/0r5PVEoi2hooQgzJKSkkhKSgLgzDPPJCUlhQ0bNighEIkTFe1gDb7T7VlZWdSvX585c+Z4Hc4pKWvYfjPzIJLooS6DCFq3bh0rVqygdWtNyCgisWvChAmkpKR4HUa5JCcns379T5PlFhQUUK9ePQ8j8p4SggjZvXs3PXv2ZPz48Zx11llehyMickoKCgp45513GDRokNehlMvFF1/M6tWrWbt2LQcOHGD69Ol07drV67A8pYQgAg4ePEjPnj3p27cvPXr0CEsdTz/9NKmpqaSlpdGnTx/27dsXlnpEYkWoL5qL5ovwImno0KE8/vjjVKoU24ePhIQEJk6cSMeOHUlJSeH6668nNTXV67A8pWsIwsw5x8CBA0lJSeGee+4JSx0bNmzgmWee4auvvuL000/n+uuvZ/r06fTv3z8s9YnEgorYtx9uc+bMITExkczMTBYuXBi2evr06cPChQvZunUrycnJjB49moEDB4a8ns6dO9O5c+eQ7zdQpNoSCkoIwmzJkiVMmzaN9PR0MjIyAHj00UdD/ktYVFTEjz/+SJUqVdi7d2+F7wuT6KUZ9WLXkiVLmD17NnPnzmXfvn3s3LmT3/zmN/zlL38JaT05OTkh3Z+XYqktSgjCrE2bNmVezRpK9evXZ9iwYTRs2JDTTz+dDh060KFDh7DWKXKq9M09NAYMGFDyjT0/Px+A++67j7fffpuqVatywQUX8Morr1CrVq2Q1Tl27FjGjh0LwMKFC3niiSdCngyId2K7E0gA2L59O7NmzWLt2rVs3LiRPXv26I9UxENljU73+9//nhYtWpCRkUGHDh3YuHFjueooa8Cb7Oxs8vPz+fzzz2natGnJwVskGEoI4sD8+fM577zzqFu3LlWqVKFHjx58/PHHXoclUmGVdbC+7777+Pzzz8nLy+Paa6/l4YcfLlcdbdu2pXbt2keUdejQgYQE34nfSy+9lIKCgnLVcTzt2rWL2TEIpGxKCOJAw4YNWbp0KXv37sU5xwcffBDz9wiLxLKyDtaBtxvv2bMn7IPgTJ48mWuuuSasdUh80TUEcaB169b06tWLVq1akZCQQMuWLbn11lu9DktEShk5ciSvvvoqNWvW5MMPPwxbPWPGjCEhIYG+ffuGrQ6JPzpDECdGjx7NN998Q35+PtOmTaNatWpehyQSdcrq21+5ciWXXXYZ6enp/PKXv2Tnzp1hq3/MmDGsX7+evn37MnHixLDUMXXqVObMmcNf//rXCj8Ur5wcJQQiUmGU1bc/aNAgxo0bxxdffEH37t354x//GPY4brjhBt54442Q73fevHk89thjzJ49mzPOOCPk+5f4poRARDx3rGnCZ8yYQWpqKpUqVWLZsmXlrqesvv1Vq1bRtm1bwHeVfjgO1ACrV68ueTx79myaNWtWrv316dOHyy67jFWrVpGcnMykSZO444472LVrF9nZ2WRkZPDb3/62vGFLBaJrCETEc8eaJjwtLY0333yT2267LWx1p6WlMXv2bLp168aMGTOOmPDmVJU1Ot3cuXNZtWoVlSpV4txzz+WFF14oVx1lDXgTrSPgSWxQQiAinjvWNOHZ2dlhr3vy5MncddddPPzww3Tt2pWqVauWe586WEssUkIgIlEl0tOEN2vWjPfeew+Ab7/9lnfeeSci9YpEG11DICJRw4tpwjdv3gzAoUOH+MMf/qB+d6mwlBBIiXBMEKNJZyRYkZgmvKwL8XJycmjatCnNmjWjXr163HzzzWGpWyTaqctASmjSGfFKJKYJh2PPPDdkyJCw1SkSK5QQiIjnjjVN+P79+7nzzjvZsmULXbp0ISMjg3/84x8eRysSn5QQiIjnjjdNePfu3SMcjUjFpGsIRERERAmBiPwk1BeB6qJSkdgR9V0GZvYQcAuwxV/0gHNurn/d74CBQDFwl3NOnYsi5aALS0UqrqhPCPyeds49EVhgZs2B3kAqUA+Yb2ZNnXPFXgQoIiISy2IlIShLN2C6c24/sNbM1gCXAJ94G1Z82bdvH23btmX//v0UFRXRq1cvRo8ezRVXXMGuXbsA38Aul1xyCW+99ZbH0YqIyKmKlYTgDjO7CVgG3Ouc2w7UB5YGbFPgL5MQqlatGgsWLKBGjRocPHiQNm3acM0117B48eKSbXr27Em3bt08jFJERMorKi4qNLP5ZpZfxtINeB64AMgANgFPHn5ZGbs66r4lM7vVzJaZ2bItW7aU8RI5HjOjRo0agG8kuYMHD2L201u/a9cuFixYwK9+9SuvQhQRkRCIijMEzrlfBLOdmb0MzPE/LQAaBKxOBjaWse+XgJcAsrKyyr7RWY6ruLiYzMxM1qxZw+DBg4+YdGbmzJlcffXVERt3XkREwiMqzhAcj5klBTztDuT7H88GeptZNTM7D2gC5EY6voqgcuXK5OXlUVBQQG5uLvn5+SXrcnJy6NOnj4fRiYhIKER9QgA8bmZfmNnnQHvgbgDn3JfA34GvgHnAYN1hEF61atWiXbt2zJs3D4Bt27aRm5tLly5dPI5Mwq24uJiWLVty7bXXHlF+5513lnQpiUhsi/qEwDl3o3Mu3TnXwjnX1Tm3KWDdGOfcBc65C51z73oZZ7zasmULO3bsAODHH39k/vz5NGvWDIAZM2Zw7bXXctppp3kZYoXg9YBBEyZMICUl5YiyZcuWlfxuiEjsi4prCCR6bdq0iX79+lFcXMyhQ4e4/vrrS74lTp8+nREjRngcYcXg5YBBBQUFvPPOO4wcOZKnnnoK8J0xuO+++3jttdeYOXOmZ7GJSOgoIZDjatGiBStWrChz3cKFCyMbjHhi6NChPP744yXjTgBMnDiRrl27kpSUdJxXikgsifouAxHxzpw5c0hMTCQzM7OkbOPGjcyYMYM777zTw8hEJNR0hkBEjmnJkiXMnj2buXPnsm/fPnbu3ElqairVqlWjcePGAOzdu5fGjRuzZs0aj6MVkfLQGQIROaaxY8dSUFDAunXrmD59OldddRXbt2/nv//9L+vWrWPdunWcccYZSgZE4oASAhEREVGXgYgEp127drRr1+6o8t27d0c+GBEJOZ0hECmHUI8PcLL7bNSoEenp6WRkZJCVlQXAQw89RP369cnIyCAjI4O5c+eGPEYRiT86QyARV716dfbs2RPS/XnFy/EBDvvwww+pU6fOEWV33313VMQmIrFDCYFEnA5UIiLRR10GIjHMzOjQoQOZmZm89NJLJeUTJ06kRYsWDBgwgO3bt3sYoYjECiUEIjFsyZIlLF++nHfffZfnnnuORYsWcfvtt/Pdd9+Rl5dHUlIS9957r9dhikgMUEIgEsPq1asHQGJiIt27dyc3N5dzzjmHypUrU6lSJW655RZyczUruIicmBICkRi1Z8+ekvkF9uzZw3vvvUdaWhqbNpVMCMrMmTNJS0vzKkQRiSG6qFAkRv3www90794dgKKiIm644QY6derEjTfeSF5eHmZGo0aNePHFFz2OVERigRICkTDYsWMHgwYNIj8/HzNj8uTJ/OMf/+Dll1+mbt26ADz66KN07tz5lOs4//zzWbly5VHl06ZNO+V9ikjFpS4DiQo7duygV69eNGvWjJSUFD755BMAnn32WS688EJSU1O5//77PY4yeEOGDKFTp0588803rFy5kpSUFMA3PkBeXh55eXnlSgZEREJNZwgkKhw+gL7++uscOHCAvXv38uGHHzJr1iw+//xzqlWrxubNm70OMyg7d+5k0aJFTJkyBYCqVatStWpVb4MSETkBnSEQzx0+gA4cOBDwHUBr1arF888/z4gRI6hWrRrgu5I+Fnz//ffUrVuXm2++mZYtWzJo0KCSkRk1PoCIRCslBOK5Yx1Av/32WxYvXkzr1q258sor+fTTT4Pep5dzDBQVFbF8+XJuv/12VqxYQfXq1Rk3bpzGBxCJcjNnzsTM+Oabb0Kyv0WLFtGqVSsSEhJ4/fXXQ7LPcFKXgXju8AH02WefpXXr1gwZMoRx48ZRVFTE9u3bWbp0KZ9++inXX38933//PWZ2wn16OTxycnIyycnJtG7dGoBevXoxbtw4zjnnnJJtbrnlFq699lqvQhSRMuTk5NCmTRumT5/OQw89VO79NWzYkClTpvDEE0+UP7gI0BkC8VxZB9Dly5eTnJxMjx49MDMuueQSKlWqxNatWz2O9sR+/vOf06BBA1atWgXABx98QPPmzTU+gEgU2717N0uWLGHSpElMnz49JPts1KgRLVq0oFKl2DjU6gyBeC7wAHrhhReWHEAvuOACFixYQLt27fj22285cODAUbP6Ratnn32Wvn37cuDAAc4//3xeeeUV7rrrLo0PIBKl3nrrLTp16kTTpk2pXbs2y5cvp1WrVkdtd8UVV5QMCBboiSee4Be/+EUkQg0bJQQSFco6gFavXp0BAwaQlpZG1apVmTp1alDdBdEgIyODZcuWHVGm8QFEoldOTg5Dhw4FoHfv3uTk5JSZECxevDjSoUWMEgKJCmUdQAH+8pe/eBCNiFQk27ZtY8GCBSUDiRUXF2NmPP7440d9CdEZAhERkTj1+uuvc9NNNx3RjXfllVfy0UcfccUVVxyxbTyfIYiNKx1ERETCJCcnp2RekMN69uzJa6+9Vq79fvrppyQnJzNjxgxuu+02UlNTy7W/cDPnnNcxRExWVpYr67S0VBxPP/00f/7znzEz0tPTeeWVV9i0aRO9e/fmf//7H61atWLatGkaWVBE4oaZfeacyzrRdjpDIBXGhg0beOaZZ1i2bBn5+fkUFxczffp0hg8fzt13383q1as5++yzmTRpktehiohEnBICqVCKior48ccfKSoqYu/evSQlJbFgwQJ69eoFQL9+/Xjrrbc8jlJEJPKUEEiFUb9+fYYNG0bDhg1JSkqiZs2aZGZmUqtWLRISfNfXJicns2HDBo8jFRGJvKhICMzsOjP70swOmVlWqXW/M7M1ZrbKzDoGlHfyl60xsxGRj1pizfbt25k1axZr165l48aN7Nmzh3ffffeo7WJlrAMRkVCKioQAyAd6AIsCC82sOdAbSAU6AX8ys8pmVhl4DrgGaA708W8rckzz58/nvPPOo27dulSpUoUePXrw8ccfs2PHDoqKigAoKCigXr16HkcqIhJ5UZEQOOe+ds6tKmNVN2C6c26/c24tsAa4xL+scc5975w7AEz3bytyTA0bNmTp0qXs3bsX51zJEMnt27cvmYls6tSpdOumXyURqXiiIiE4jvrA+oDnBf6yY5WLHFPr1q3p1asXrVq1Ij09nUOHDnHrrbfy2GOP8dRTT9G4cWO2bdvGwIEDvQ5VRCTiIjZSoZnNB35exqqRzrlZx3pZGWWOshOZMgdUMLNbgVvB9w1RKrbRo0czevToI8rOP/98cnNzPYpIRCQ6RCwhcM6dyiDPBUCDgOfJwEb/42OVl673JeAl8A1MdAoxiIiIxL1o7zKYDfQ2s2pmdh7QBMgFPgWamNl5ZlYV34WHsz2MU0REJKZFxeRGZtYdeBaoC7xjZnnOuY7OuS/N7O/AV0ARMNg5V+x/zR3AP4DKwGTn3JcehS8iIhLzNJeBiIhIHNNcBiIiIhI0JQQiIiKihEBERESUEIiIiAhKCERERAQlBCIiIoISAhEREaGCjUNgZluAf5cqrgNs9SCccImn9qgt0Sue2hNPbYH4ao/aEhrnOufqnmijCpUQlMXMlgUzYEOsiKf2qC3RK57aE09tgfhqj9oSWeoyEBERESUEIiIiooQA/FMjx5F4ao/aEr3iqT3x1BaIr/aoLRFU4a8hEBEREZ0hEBERESpQQmBmncxslZmtMbMRZayvZmZ/86//l5k1inyUwQmiLW3NbLmZFZlZLy9iPBlBtOceM/vKzD43sw/M7Fwv4gxGEG35rZl9YWZ5ZvaRmTX3Is5gnag9Adv1MjNnZlF7FXUQn01/M9vi/2zyzGyQF3EGI5jPxcyu9//dfGlmr0U6xpMRxGfzdMDn8q2Z7fAizmAE0ZaGZvahma3w/0/r7EWcZXLOxf0CVAa+A84HqgIrgealtvk/4AX/497A37yOuxxtaQS0AF4Fenkdcwja0x44w//49hj/bM4KeNwVmOd13OVpj3+7M4FFwFIgy+u4y/HZ9Acmeh1riNrSBFgBnO1/nuh13OX9PQvY/k5gstdxl+OzeQm43f+4ObDO67gPLxXlDMElwBrn3PfOuQPAdKBbqW26AVP9j18HrjYzi2CMwTphW5xz65xznwOHvAjwJAXTng+dc3v9T5cCyRGOMVjBtGVnwNPqQDRfxBPM3w3AI8DjwL5IBneSgm1LLAimLbcAzznntgM45zZHOMaTcbKfTR8gJyKRnbxg2uKAs/yPawIbIxjfcVWUhKA+sD7geYG/rMxtnHNFQCHws4hEd3KCaUssOdn2DATeDWtEpy6otpjZYDP7Dt9B9K4IxXYqTtgeM2sJNHDOzYlkYKcg2N+znv7TuK+bWYPIhHbSgmlLU6CpmS0xs6Vm1ili0Z28oP8H+LsLzwMWRCCuUxFMWx4CfmNmBcBcfGc8okJFSQjK+qZf+ptZMNtEg1iJM1hBt8fMfgNkAX8Ma0SnLqi2OOeec85dAAwHHgx7VKfuuO0xs0rA08C9EYvo1AXz2bwNNHLOtQDm89MZw2gTTFsS8HUbtMP3jfrPZlYrzHGdqpP5n9YbeN05VxzGeMojmLb0AaY455KBzsA0/9+S56IiiAgoAAKz/WSOPk1Tso2ZJeA7lfO/iER3coJpSywJqj1m9gtgJNDVObc/QrGdrJP9bKYDvwprROVzovacCaQBC81sHXApMDtKLyw84WfjnNsW8Lv1MpAZodhOVrD/z2Y55w4659YCq/AlCNHoZP5uehO93QUQXFsGAn8HcM59ApyGb54Dz1WUhOBToImZnWdmVfH9Us0utc1soJ//cS9ggfNf9RFlgmlLLDlhe/ynpV/ElwxEc19oMG0J/KfcBVgdwfhO1nHb45wrdM7Vcc41cs41wnd9R1fn3DJvwj2uYD6bpICnXYGvIxjfyQjmf8Bb+C7Gxczq4OtC+D6iUQYvqP9pZnYhcDbwSYTjOxnBtOU/wNUAZpaCLyHYEtEoj8XrqxojteA7NfMtvitAR/rLHsb3Dwx8H8oMYA2QC5zvdczlaMvF+DLVPcA24EuvYy5ne+YDPwB5/mW21zGXoy0TgC/97fgQSPU65vK0p9S2C4nSuwyC/GzG+j+blf7PppnXMZejLQY8BXwFfAH09jrm8v6e4et7H+d1rCH4bJoDS/y/Z3lAB69jPrxopEIRERGpMF0GIiIichxKCEREREQJgYiIiCghEBEREZQQiIiICEoIRMrNP8tf1M8qGcjM5pjZFK/j8JKZLTSziV7HIRItlBBIhWFmU/wH79LL0pN4fVlj9ifhG/Y2rMxsnZkNC3c9/rralXqPtpjZu2Z2USTqj5AewO8OPwnV+2tmZ5jZo/7pb/eZ2Vb/nAJ9yrtvkXBK8DoAkQibD9xYquxAeXbonPtveV4f5VLxDeHdEHgGmGdmzZxzhaU3NLOqzjfDW0xwzoVraPIXgMuBIUA+UBto7f8ZFrH23kuU8npkJC1aIrUAU4A5J9jmNnyjjO3DN5zoP/Alzg/hm6QkcGnnf40DevkfN/I/7w38E/gR37z0LfCN+/8xvhEkPwLOC6j3AmAW8F//+uXAtQHrF5auP2Dd/+evay+wAXgeOCtg/Rn+tu/GN+LjA8AcfBOsHOt9aOevp05A2eX+so7+5+v878tkYAcww1+eji/x+hFfMjEFqFn6c8A3sdMP/rheAU4P2MaA+/GN9vYjvtH2fhOw/vD73BN439/2r4DsgG2q4EtiNgL78c1CN67UezrxWO8vvumpdx7+bANelw0cBM45xnu3Axh0gt8zwzcp1Gp/bAXA2ID1wb6Hw/2v3ewvrwo8xk8jlX56+PPSouVEi7oMRPz8k/I8B4wGLgR+Aczzr34C34Qk8/F1ESThO7gfy2h8/5hb4jtAvAY8i2+CpkvwDZX9TMD2NfBN65wNXAS8AbxpZs3863vg+yf/cED9mFk68B6+8dIv8m+Xge8gfdgT/v32xDeGekugbTDvSSk/+n9WCSi7B/gG3yyUD5jZGfjes93+dnbHl7AExgNwpT/eq/1xdcD3fh32B3yTwAzGN9TrWOBFM+tSaj9j8L2PF+E7+E03sxr+dXf56++Nb2KfX+Ob5KcsR72/zrk9+CbSGVBq2wH4EssfjrGv/wKdzKzmMdYDPAr83t+uVOA6/NPmnuR72ALohH9sfHyJ1ZXADfiSiqnA23HW1SPh4nVGokVLpBZ836qK8P2jDVwe86/vARQCZx7n9UedYaDsMwS3Bay/1l/WI6CsP7D7BPEuBR4MeL4OGFZqm1eBSaXKMvz1JeJLNPYDfQPW18CXpEw5Tt3tCDhDAPwM3xmMnUBiQDxvl3rdLaXfw4B9NQ54H3cANQK2+Y0/zur+5UfgilL7Hg/MPc77XN9f1sb//BngA/AN0V5GGxfiP0NwnPc3y/87U9///Gx/bNeWtU//Nm3xHdwP4jvTM5Ejz1zUwHcG6rfHeH2w7+EWoFrANhcAh4CGpfb3FvAnr//+tET/omsIpKJZBNxaqmyH/+f7wL+BtWb2D3zfvN90zu06hXo+D3h8+JvkF6XKqpvZGc65vWZWHRiFL3lIwvct/LRS+ylLJtDYzH4dUHZ4TvYL8J1Kr0rADHHOud1mFhjL8awzM/AdpFcD17kjZ5wsPbNhCvB5qffsY3wHqub4Jg/Dv83ugG0+8cd5AVANX9vnmVngZCtV8B20AwW+P4enmU30/5yC7zP91szeA+YC7zrnDh2rsaU555b536t++L7V3wBsx3c251ivWWRm5+ObDvpy4CrgPTN7yTl3G773oRq+ZKUswb6H+e7IqcBb4fvsv/J/ZodVAxYE0Vyp4JQQSEWz1zm3pqwVzrldZtYK3ze8bHxXoD9qZhc75441P/uxHAzc9XHKDnfbPYHv1O8wfAfevfi+/Vc9QT2VgD8DT5exbgO+ro/yaI+vD3uLc25nGev3lHpu/NS20oKdSe3we/L/t3c/oVFdURzHv78oWGrdtC5ECnUhVARXIohoQbGUgEK7cyUogpsuRG2hUCMo4v+FuGiTEBHRVVBx6cY/cSMiFgWVpppmI7X1Hyq0xJieLs4NPpOZZKIJUvP7QJjMzJuZO3cec86997w3q8mfiq3qr3c9IqIEwqZy/ZqkOWS/riCnz69L+nIsSQHZv5vIhGA9ObMyMCw6ZRIAAAL5SURBVNIDIqIfuFT+9kj6EdgpaTevErZ6Gu3DoX3fVO5fxPB++gezUTghMKuIiJfkaOqcpO3AX+SovY08GmHKBL30UuBYRJwEkPQBOVrurmxT6/WvkT+hXDPJkXSHDA6LgZ5y23SywPFuA+36PSIejuF93ALWS5pRGeEuIYPV7cp2CyRNj1ynp7TvRWlTE7l88FlEvNXItrShE+gs5124DMzl9X4dVO/zPQ7sl/QtOQpf8wZNuVUuPyr/95Hr/r/V2baRPhzqFzKZmBUR59+gjTbJOSGwyWaapFlDbhuIiAeSVpFBuIscFS8HZvDqS7gXaJb0OfAIeFpGguOhG/hG0hkygG8np82reoFlko4DfSVQ7wUuS/oZaAWeA/OA1RGxsSwPdAB7JT0gp9VbmLjE5gRZUHlMUgu55t5KLr1Uk5apwBFJO4DZwB6gfTBBkHQAOKAc8neRgXQx8G9EtDXSEEmbgT/I35zvJ6f7n5HFg7X0Mrx/iYinkjqBg0BXRNQK4tXXvUAWI14l95P55OzCr8DtiBiQdAjYLamvvL9PgIUR8RON9+FrIqJb0gngqKQtZLL4MVl/0BMRp0Zqt5kTAptsVpJBouoe8ClZS/A1GTA/JEerGyLiUtmunfxyvUoGqOVkYdp42Ax0kFPMT8gCuqEJQQsZGO6S68KKiBuSviCr8i+Sgb4HOF153FayBuA0uRRxuFwfd6Ue4qvS/itk8dwZ8pj8qovATeA82dcnycMMB20j6yy2kodRPiMD+74xNOc58B15hEGQI+jmiPi7zvbD+rdyXwewtlyO5ix5rotd5H5yn6xl2FFZaviB/Jy3kfven+QS0Vj6sJZ15JEs+8rzPi7P4RkDG5UiGl3WMzN7e2XqfmZErHrXbWlUKdpsBWaPkFCY/a95hsDMrI5yToA55Mmc2p0M2PvMJyYyM6vve+A6OfW+8x23xWxCecnAzMzMPENgZmZmTgjMzMwMJwRmZmaGEwIzMzPDCYGZmZnhhMDMzMyA/wD4ejbEsZXg3AAAAABJRU5ErkJggg==\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.set_ylim(-115, 295)\n",
    "\n",
    "ax.axhline(0, c='gray', linewidth=1)\n",
    "\n",
    "bars0 = ax.bar(bins[:-1] + 0.025, top0, width=0.04, facecolor='white')\n",
    "bars1 = ax.bar(bins[:-1] + 0.025, -top1, width=0.04, facecolor='gray')\n",
    "\n",
    "for bars in (bars0, bars1):\n",
    "    for bar in bars:\n",
    "        bar.set_edgecolor(\"gray\")\n",
    "\n",
    "for x, y in zip(bins, top0):\n",
    "    ax.text(x + 0.025, y + 10, str(y), ha='center', va='bottom')\n",
    "\n",
    "for x, y in zip(bins, top1):\n",
    "    ax.text(x + 0.025, -y - 10, str(y), ha='center', va='top')\n",
    "\n",
    "ax.text(0.75, 260, \"A = 0\")\n",
    "ax.text(0.75, -90, \"A = 1\")\n",
    "\n",
    "ax.set_ylabel(\"No. Subjects\", fontsize=14)\n",
    "ax.set_xlabel(\"Estimated Propensity Score\", fontsize=14);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                  mean propensity\n",
      "    non-quitters: 0.245\n",
      "        quitters: 0.312\n"
     ]
    }
   ],
   "source": [
    "print('                  mean propensity')\n",
    "print('    non-quitters: {:>0.3f}'.format(propensity0.mean()))\n",
    "print('        quitters: {:>0.3f}'.format(propensity1.mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 15.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 15.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"only individual 22005 had an estimated $\\pi(L)$ of 0.6563\", pg 186"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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>propensity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1088</th>\n",
       "      <td>22005</td>\n",
       "      <td>0.656281</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       seqn  propensity\n",
       "1088  22005    0.656281"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all[['seqn', 'propensity']].loc[abs(propensity - 0.6563) < 1e-4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the deciles and check their counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs_all['decile'] = pd.qcut(nhefs_all.propensity, 10, labels=list(range(10)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    163\n",
       "1    163\n",
       "2    163\n",
       "3    163\n",
       "4    163\n",
       "5    162\n",
       "6    163\n",
       "7    163\n",
       "8    163\n",
       "9    163\n",
       "Name: decile, dtype: int64"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nhefs_all.decile.value_counts(sort=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create a model with interaction between `qsmk` and deciles of $\\pi(L)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sm.OLS.from_formula('wt82_71 ~ qsmk * C(decile)', data=nhefs_all)\n",
    "res = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "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>Intercept</th>           <td>    3.9952</td> <td>    0.630</td> <td>    6.338</td> <td> 0.000</td> <td>    2.759</td> <td>    5.232</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.1]</th>      <td>   -1.0905</td> <td>    0.916</td> <td>   -1.190</td> <td> 0.234</td> <td>   -2.888</td> <td>    0.707</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.2]</th>      <td>   -1.3831</td> <td>    0.918</td> <td>   -1.506</td> <td> 0.132</td> <td>   -3.184</td> <td>    0.418</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.3]</th>      <td>   -0.5205</td> <td>    0.926</td> <td>   -0.562</td> <td> 0.574</td> <td>   -2.337</td> <td>    1.295</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.4]</th>      <td>   -1.8964</td> <td>    0.940</td> <td>   -2.017</td> <td> 0.044</td> <td>   -3.741</td> <td>   -0.052</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.5]</th>      <td>   -2.1482</td> <td>    0.954</td> <td>   -2.252</td> <td> 0.024</td> <td>   -4.019</td> <td>   -0.277</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.6]</th>      <td>   -2.4352</td> <td>    0.949</td> <td>   -2.566</td> <td> 0.010</td> <td>   -4.297</td> <td>   -0.574</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.7]</th>      <td>   -3.7105</td> <td>    0.974</td> <td>   -3.810</td> <td> 0.000</td> <td>   -5.621</td> <td>   -1.800</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.8]</th>      <td>   -4.8907</td> <td>    1.034</td> <td>   -4.731</td> <td> 0.000</td> <td>   -6.918</td> <td>   -2.863</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.9]</th>      <td>   -4.4996</td> <td>    1.049</td> <td>   -4.288</td> <td> 0.000</td> <td>   -6.558</td> <td>   -2.441</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>                <td>   -0.0147</td> <td>    2.389</td> <td>   -0.006</td> <td> 0.995</td> <td>   -4.700</td> <td>    4.671</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.1]</th> <td>    4.1634</td> <td>    2.897</td> <td>    1.437</td> <td> 0.151</td> <td>   -1.520</td> <td>    9.847</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.2]</th> <td>    6.5591</td> <td>    2.884</td> <td>    2.275</td> <td> 0.023</td> <td>    0.903</td> <td>   12.215</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.3]</th> <td>    2.3570</td> <td>    2.808</td> <td>    0.839</td> <td> 0.401</td> <td>   -3.151</td> <td>    7.865</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.4]</th> <td>    4.1450</td> <td>    2.754</td> <td>    1.505</td> <td> 0.132</td> <td>   -1.257</td> <td>    9.547</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.5]</th> <td>    4.5580</td> <td>    2.759</td> <td>    1.652</td> <td> 0.099</td> <td>   -0.853</td> <td>    9.969</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.6]</th> <td>    4.3306</td> <td>    2.776</td> <td>    1.560</td> <td> 0.119</td> <td>   -1.115</td> <td>    9.776</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.7]</th> <td>    3.5847</td> <td>    2.744</td> <td>    1.307</td> <td> 0.192</td> <td>   -1.797</td> <td>    8.966</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.8]</th> <td>    2.3155</td> <td>    2.700</td> <td>    0.858</td> <td> 0.391</td> <td>   -2.981</td> <td>    7.612</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk:C(decile)[T.9]</th> <td>    2.2549</td> <td>    2.687</td> <td>    0.839</td> <td> 0.402</td> <td>   -3.016</td> <td>    7.526</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the effect estimates, we'll use t-tests with contrast DataFrames, like we did in Program 15.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           estimate    95% C.I.\n",
      "\n",
      "decile 0     -0.0    (-4.7, 4.7)\n",
      "decile 1      4.1    ( 0.9, 7.4)\n",
      "decile 2      6.5    ( 3.4, 9.7)\n",
      "decile 3      2.3    (-0.6, 5.2)\n",
      "decile 4      4.1    ( 1.4, 6.8)\n",
      "decile 5      4.5    ( 1.8, 7.2)\n",
      "decile 6      4.3    ( 1.5, 7.1)\n",
      "decile 7      3.6    ( 0.9, 6.2)\n",
      "decile 8      2.3    (-0.2, 4.8)\n",
      "decile 9      2.2    (-0.2, 4.7)\n"
     ]
    }
   ],
   "source": [
    "# start with empty DataFrame\n",
    "contrast = pd.DataFrame(\n",
    "    np.zeros((2, res.params.shape[0])),\n",
    "    columns=res.params.index\n",
    ")\n",
    "\n",
    "# modify the constant entries\n",
    "contrast['Intercept'] = 1\n",
    "contrast['qsmk'] = [1, 0]\n",
    "\n",
    "# loop through t-tests, modify the DataFrame for each decile,\n",
    "# and print out effect estimate and confidence intervals\n",
    "print('           estimate    95% C.I.\\n')\n",
    "for i in range(10):\n",
    "    if i != 0:\n",
    "        # set the decile number\n",
    "        contrast['C(decile)[T.{}]'.format(i)] = [1, 1]\n",
    "        contrast['qsmk:C(decile)[T.{}]'.format(i)] = [1, 0]\n",
    "    \n",
    "    ttest = res.t_test(contrast.iloc[0] - contrast.iloc[1])\n",
    "    est = ttest.effect[0]\n",
    "    conf_ints = ttest.conf_int(alpha=0.05)\n",
    "    lo, hi = conf_ints[0, 0], conf_ints[0, 1]\n",
    "\n",
    "    print('decile {}    {:>5.1f}    ({:>4.1f},{:>4.1f})'.format(i, est, lo, hi))\n",
    "    \n",
    "    if i != 0:\n",
    "        # reset to zero\n",
    "        contrast['C(decile)[T.{}]'.format(i)] = [0, 0]\n",
    "        contrast['qsmk:C(decile)[T.{}]'.format(i)] = [0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the estimates above to the estimate we get from a model without interaction between `qsmk` and deciles of $\\pi(L)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sm.OLS.from_formula(\n",
    "    'wt82_71 ~ qsmk + C(decile)', \n",
    "    data=nhefs_all\n",
    ")\n",
    "res = model.fit()"
   ]
  },
  {
   "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>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>      <td>    3.7505</td> <td>    0.609</td> <td>    6.159</td> <td> 0.000</td> <td>    2.556</td> <td>    4.945</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.1]</th> <td>   -0.7391</td> <td>    0.861</td> <td>   -0.858</td> <td> 0.391</td> <td>   -2.428</td> <td>    0.950</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.2]</th> <td>   -0.6182</td> <td>    0.861</td> <td>   -0.718</td> <td> 0.473</td> <td>   -2.307</td> <td>    1.071</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.3]</th> <td>   -0.5204</td> <td>    0.858</td> <td>   -0.606</td> <td> 0.544</td> <td>   -2.204</td> <td>    1.163</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.4]</th> <td>   -1.4884</td> <td>    0.859</td> <td>   -1.733</td> <td> 0.083</td> <td>   -3.173</td> <td>    0.197</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.5]</th> <td>   -1.6227</td> <td>    0.868</td> <td>   -1.871</td> <td> 0.062</td> <td>   -3.324</td> <td>    0.079</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.6]</th> <td>   -1.9853</td> <td>    0.868</td> <td>   -2.287</td> <td> 0.022</td> <td>   -3.688</td> <td>   -0.282</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.7]</th> <td>   -3.4447</td> <td>    0.875</td> <td>   -3.937</td> <td> 0.000</td> <td>   -5.161</td> <td>   -1.729</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.8]</th> <td>   -5.1544</td> <td>    0.885</td> <td>   -5.825</td> <td> 0.000</td> <td>   -6.890</td> <td>   -3.419</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(decile)[T.9]</th> <td>   -4.8403</td> <td>    0.883</td> <td>   -5.483</td> <td> 0.000</td> <td>   -6.572</td> <td>   -3.109</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>           <td>    3.5005</td> <td>    0.457</td> <td>    7.659</td> <td> 0.000</td> <td>    2.604</td> <td>    4.397</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         estimate   95% C.I.\n",
      "effect      3.5    (2.6, 4.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('effect    {:>5.1f}    ({:>0.1f}, {:>0.1f})'.format(est, lo, hi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 15.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will do \"outcome regression $\\mathrm{E}[Y|A, C=0, p(L)]$ with the estimated propensity score $\\pi(L)$ as a continuous covariate rather than as a set of indicators\" (pg 187)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhefs['propensity'] = propensity[~nhefs_all.wt82_71.isnull()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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>Intercept</th>  <td>    5.5945</td> <td>    0.483</td> <td>   11.581</td> <td> 0.000</td> <td>    4.647</td> <td>    6.542</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>qsmk</th>       <td>    3.5506</td> <td>    0.457</td> <td>    7.765</td> <td> 0.000</td> <td>    2.654</td> <td>    4.448</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>propensity</th> <td>  -14.8218</td> <td>    1.758</td> <td>   -8.433</td> <td> 0.000</td> <td>  -18.269</td> <td>  -11.374</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sm.OLS.from_formula('wt82_71 ~ qsmk + propensity', data=nhefs)\n",
    "res = model.fit()\n",
    "res.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the coefficient on `qsmk` we can see the effect estimate is 3.6.\n",
    "\n",
    "We'll use bootstrap to get confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def outcome_regress_effect(data):\n",
    "    model = sm.OLS.from_formula('wt82_71 ~ qsmk + propensity', data=data)\n",
    "    res = model.fit()\n",
    "    \n",
    "    data_qsmk_1 = data.copy()\n",
    "    data_qsmk_1['qsmk'] = 1\n",
    "    \n",
    "    data_qsmk_0 = data.copy()\n",
    "    data_qsmk_0['qsmk'] = 0\n",
    "    \n",
    "    mean_qsmk_1 = res.predict(data_qsmk_1).mean()\n",
    "    mean_qsmk_0 = res.predict(data_qsmk_0).mean()\n",
    "    \n",
    "    return mean_qsmk_1 - mean_qsmk_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonparametric_bootstrap(data, func, n=1000):\n",
    "    estimate = func(data)\n",
    "    \n",
    "    n_rows = data.shape[0]\n",
    "    indices = list(range(n_rows))\n",
    "    \n",
    "    b_values = []\n",
    "    for _ in tqdm(range(n)):\n",
    "        data_b = data.sample(n=n_rows, replace=True)\n",
    "        b_values.append(func(data_b))\n",
    "    \n",
    "    std = np.std(b_values)\n",
    "    \n",
    "    return estimate, (estimate - 1.96 * std, estimate + 1.96 * std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2000/2000 [00:19<00:00, 100.81it/s]\n"
     ]
    }
   ],
   "source": [
    "data = nhefs[['wt82_71', 'qsmk', 'propensity']]\n",
    "\n",
    "info = nonparametric_bootstrap(\n",
    "    data, outcome_regress_effect, n=2000\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         estimate   95% C.I.\n",
      "effect      3.6    (2.6, 4.5)\n"
     ]
    }
   ],
   "source": [
    "print('         estimate   95% C.I.')\n",
    "print('effect    {:>5.1f}    ({:>0.1f}, {:>0.1f})'.format(info[0], info[1][0], info[1][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
}