11-Propensity-Score.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Psychology of Growth\n",
    "\n",
    "The field of positive psychology studies what are the human behaviours that lead to a great life. You can think of it as the intersection between self help books with the academic rigor of statistics. One of the famous findings of positive psychology is the **Growth Mindset**. The idea is that people can have a fixed or a growth mindset. If you have a fixed mindset, you believe that abilities are given at birth or in early childhood. As such, intelligence is fixed and can't change throughout life. If you don't have it by now, you can't acquire it. The corollary of this though is that you should not waste time on areas where you don't excel, since you will never learn how to handle them. On the other hand, if you have a growth mindset, you believe that intelligence can be developed. The direct consequence of this is you see failure not as lack of intelligence, put as part of a learning process. \n",
    "\n",
    "I don't want to debate which of these mindsets is the correct one (its probably somewhere in the middle). For our purpose, it doesn't matter much. What does matter is that psychologists found out that people who have a growth mindset tend to do better in life. They are more likely to achieve what they've set to.\n",
    "\n",
    "As versed as we are with causal inference, we've learned to see those statements with skepticism. Is it that a growth mindset causes people to achieve more? Or is simply the case that people who achieve more are prone to develop a growth mindset as a result of their success? Who came first, the egg or the chicken? In potential outcome notation, we have reasons to believe that there is bias in these statements. \\\\(Y_0|T=1\\\\) is probably larger than \\\\(Y_0|T=0\\\\), which means that those with a growth mindset would have achieved more even if they had a fixed mindset. \n",
    "\n",
    "To settle things, researches designed the [The National Study of Learning Mindsets](https://mindsetscholarsnetwork.org/about-the-network/current-initatives/national-mindset-study/#). It is a randomised study conducted in U.S. public high schools which aims at finding the impact of a growth mindset. The way it works is that students receive from the school a seminary to instil in them a growth mindset. Then, they follow up the students in their college years to measure how well they've performed academically. This measurement was compiled into an achievement score and standardised. The real data on this study is not publicly available in order to preserve students' privacy. However, we have a simulated dataset with the same statistical properties provided by [Athey and Wager](https://arxiv.org/pdf/1902.07409.pdf), so we will use that instead. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import style\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "import graphviz as gr\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "style.use(\"fivethirtyeight\")\n",
    "pd.set_option(\"display.max_columns\", 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Besides the treated and outcome variables, the study also recorded some other features:\n",
    "\n",
    "* schoolid: identifier of the student's school;\n",
    "* success_expect:  self-reported expectations for success in the future, a proxy for prior achievement, measured prior to random assignment;\n",
    "* ethnicity: categorical variable for student race/ethnicity;\n",
    "* gender: categorical variable for student identified gender;\n",
    "* frst_in_family: categorical variable for student first-generation status, i.e. first in family to go to college;\n",
    "* school_urbanicity: school-level categorical variable for urbanicity of the school, i.e. rural, suburban, etc;\n",
    "* school_mindset: school-level mean of students’ fixed mindsets, reported prior to random assignment, standardised;\n",
    "* school_achievement: school achievement level, as measured by test scores and college preparation for the previous 4 cohorts of students, standardised; \n",
    "* school_ethnic_minority: school racial/ethnic minority composition, i.e., percentage of student body that is Black, Latino, or Native American, standardised;\n",
    "* school_poverty: school poverty concentration, i.e., percentage of students who are from families whose incomes fall below the federal poverty line, standardised;\n",
    "* school_size: total number of students in all four grade levels in the school, standardised."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>schoolid</th>\n",
       "      <th>intervention</th>\n",
       "      <th>achievement_score</th>\n",
       "      <th>...</th>\n",
       "      <th>school_ethnic_minority</th>\n",
       "      <th>school_poverty</th>\n",
       "      <th>school_size</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>259</th>\n",
       "      <td>73</td>\n",
       "      <td>1</td>\n",
       "      <td>1.480828</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.515202</td>\n",
       "      <td>-0.169849</td>\n",
       "      <td>0.173954</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3435</th>\n",
       "      <td>76</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.987277</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.310927</td>\n",
       "      <td>0.224077</td>\n",
       "      <td>-0.426757</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9963</th>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.152340</td>\n",
       "      <td>...</td>\n",
       "      <td>0.875012</td>\n",
       "      <td>-0.724801</td>\n",
       "      <td>0.761781</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4488</th>\n",
       "      <td>67</td>\n",
       "      <td>0</td>\n",
       "      <td>0.358336</td>\n",
       "      <td>...</td>\n",
       "      <td>0.315755</td>\n",
       "      <td>0.054586</td>\n",
       "      <td>1.862187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2637</th>\n",
       "      <td>16</td>\n",
       "      <td>1</td>\n",
       "      <td>1.360920</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.033161</td>\n",
       "      <td>-0.982274</td>\n",
       "      <td>1.591641</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 13 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      schoolid  intervention  achievement_score  ...  school_ethnic_minority  \\\n",
       "259         73             1           1.480828  ...               -0.515202   \n",
       "3435        76             0          -0.987277  ...               -1.310927   \n",
       "9963         4             0          -0.152340  ...                0.875012   \n",
       "4488        67             0           0.358336  ...                0.315755   \n",
       "2637        16             1           1.360920  ...               -0.033161   \n",
       "\n",
       "      school_poverty  school_size  \n",
       "259        -0.169849     0.173954  \n",
       "3435        0.224077    -0.426757  \n",
       "9963       -0.724801     0.761781  \n",
       "4488        0.054586     1.862187  \n",
       "2637       -0.982274     1.591641  \n",
       "\n",
       "[5 rows x 13 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"./data/learning_mindset.csv\")\n",
    "data.sample(5, random_state=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the study was randomised, it doesn't seem to be the case that this data is free from confounding. If we look at the additional features, we will notice that they vary systematically between treatment and control. One possible reason for this is that the treatment variable is measured by the student's receipt of the seminar. So, although the opportunity to participate was random, participation itself is not. We are here dealing with a case of non-compliance here. One evidence of this is how the student's success expectation is correlated with the participation in the seminar. Students with higher self-reported participation are more likely to have joined the growth mindset seminar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "success_expect\n",
       "1    0.271739\n",
       "2    0.265957\n",
       "3    0.294118\n",
       "4    0.271617\n",
       "5    0.311070\n",
       "6    0.354287\n",
       "7    0.362319\n",
       "Name: intervention, dtype: float64"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"success_expect\")[\"intervention\"].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still, let's see what the difference in means \\\\(E[Y|T=1] - E[Y|T=0]\\\\) looks like. This will be a useful baseline to compare against. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>   -0.1538</td> <td>    0.012</td> <td>  -13.201</td> <td> 0.000</td> <td>   -0.177</td> <td>   -0.131</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>intervention</th> <td>    0.4723</td> <td>    0.020</td> <td>   23.133</td> <td> 0.000</td> <td>    0.432</td> <td>    0.512</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smf.ols(\"achievement_score ~ intervention\", data=data).fit().summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simply comparing those with and without the intervention, we can see that the treated have an achievement score that is, on average, 0.3185 (0.4723 - 0.1538) higher than the untreated. But is this big or small? I know that interpreting standardised outcomes can be challenging, but bear with me for a moment. I think it is worth going through this because it won't be the last time you will encounter standardized scores.\n",
    "\n",
    "The outcome variable being standardised means that it is measured in standard deviations. So, the treated are 0.3185 deviations above the untreated. That is what this means. As for if this is small or big, let's remember some stuff about the normal distribution. We know that 95% of its mass is between 2 standard deviations, leaving 2.5% on one tail and 2.5% on another. This also means that if someone is 2 standard deviations above the mean, 97.5% (95% plus the left 2.5% tail) of all the individuals are below that person. By looking at the normal CDF, we also know that about 0.85% of its mass is below 1 standard deviation and that 70% of its mass is below 0.5 standard deviations. So, this means that the average of the treated is above 70% of the individual achievements. The untreated mean, on the other hand, have only 44% of the individuals below. \n",
    "\n",
    "Here is what this looks like in a picture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAD1CAYAAAAS9/GUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAeDUlEQVR4nO3de3SV1bnv8e9DkIRyD+46QPCy25zuXS+tlyKgokOsGFTw1Dpq9z67ORaH1WLr1rPlsjvG8RwcVK0dXuqwum3BiqXaVtxF6w22lIqDotZKVRBdObZCDDU1gRAiIYQ854/1Jl2EtXJZt7nWyu8zRkbWO9+55vtMspKH+b7zna+5OyIiIqEMCR2AiIgMbkpEIiISlBKRiIgEpUQkIiJBKRGJiEhQQ0MHkKi5uVlT+EREStyYMWMscVsjIhERCUqJSEREghpUiSgWi4UOIWvUl8KkvhSuUupPKfUFBlkiEhGRwqNEJCIiQRXUrDkRkXxzd/bu3UtnZ2foUPqtoqKC5ubm0GGkNGTIEEaOHImZ9V0ZJSIRGeT27t1LeXk5w4YNCx1Kv5WXl1NRURE6jJTa29vZu3cvo0aN6lf9Pk/NmdlyM2sws7eS7Ps3M3MzOzLaNjP7gZnVmtkbZnZqQt0aM4tFXzUD6JOISM50dnYWVRIqBsOGDRvQCLM/14h+AlzYs9DMJgNfBLYnFFcDVdHX1cD9Ud1K4GbgDGAKcLOZjet3lCIiUrL6TETu/iLQlGTXXcACIHE1hLnACo/bBIw1swnALGCtuze5+y5gLUmSm4jIYPXUU08xduxY3n33XQDef/99pk2bBsCGDRv4yle+EjK8nErrGpGZzQE+cPc/9rgYdTSwI2G7LipLVS4SxLPb92XcRvUxw7MQiRSabHw2EvX3c7Jq1SqmTZvGqlWrWLx4cVZjKHQDTkRm9gngO8AFyXYnKfNeylPK1Q1bpXQjmPqSvvrGsozbiO0/mLxcP5eClaw/FRUVlJeXd28fOHAgq8dsa+t75lhrayubNm1i1apV1NTUcMMNN7B//346Oztpa2ujvb2dgwcP0tbWltBuWy8thrdnzx4aGhq6t6uqqlLWTWdE9CngeKBrNDQJ+IOZTSE+0pmcUHcSUB+Vn9ujfH1vB+kt6HTFYrGctBuC+pKZ2vLM/9dbleR/uvq5FK5U/Wlubj5kBtoRR2R37eX+zG5bvXo1559/PieccAKVlZVs27aNcePGMWTIECoqKhg2bBhlZWXdbbW1tRX0rDmA0aNHM3ny5L4rksYNre7+prt/0t2Pc/fjiCeZU939L8CTwNei2XNTgWZ33wk8D1xgZuOiSQoXRGUiIoPeqlWruOyyywD40pe+xKpVqwJHlF99jojM7FHio5kjzawOuNndl6Wo/gwwG6gFPgauBHD3JjO7BXg1qrfE3ZNNgBARGVSampp48cUX2bp1K2ZGZ2cnZsa8efNCh5Y3fSYid/9qH/uPS3jtwPwU9ZYDywcYn4hISVu9ejVXXHEFd999d3fZ7Nmzqa+vDxhVfmmtORGRgB5//HEuvvjiQ8rmzJnDnXfeGSii/NMSPyIiCfI9Lf/pp58+rOyaa67hmmuu6d4+++yzOfvss/MZVl5pRCQiIkEpEYmISFBKRCIiEpQSkYiIBKXJCiJpSrYmWX1j2YBWbdB6dSIaEYmISGBKRCIigSU+8qHLrbfeyr333pvyPW+88QZr1qzJWgzptnfRRRfx+uuvZ3RsJSIRkSL05ptvsnbt2qT7Ojo6stperukakYhIgtZ12X1m54jznsvo/RdddBGnn346GzZsoLm5mXvvvZcTTzyRW2+9lX379vG73/2OG2+8kXfeeYe//OUvbN++nfHjx3Pbbbdxww03UFdXB8RHWFOnTuW1115j8eLF7Nu3j+HDh3Pfffdx7LHHHtberFmzWLBgAVu3bqWjo4NFixZx0UUXsW/fPubPn8+2bdv4zGc+k5XHUSgRiYgUuI6ODtatW8eaNWu4/fbb+fnPf87ixYvZvHkzd9xxBxBPNJs3b+a5555j+PDhXHXVVXzzm99k2rRp7Nixg8suu4xXXnmFqqoqnnnmGYYOHcr69etZsmQJjzzyyGHtLVmyhBkzZnDfffexe/duZs6cybnnnstDDz3E8OHD2bhxI2+99RbnnHNOxv1TIhIRCazHk64Pc8kllwDw+c9/nu3bt6esV11dzfDh8ZmY69evZ9u2bd37WlpaaGlpYc+ePVx77bW89957mFnKBwGuW7eOZ599tvs61f79+6mrq2Pjxo184xvfAODEE0/khBNO6H9HU1AiEhEJrLKykt27dx9StmvXLo499liA7ifIlpWV9Xr9Z8SIEd2vOzs7Wbt2bXdi6rJgwQLOPvtsVq5cyfvvv3/Ygqtd3J0VK1YkfZhgX4lzoJSIREQSZHpNJx0jR47kqKOO4re//S3nnHMOu3bt4oUXXuDaa69l5cqVKd/T0tKSss3zzjuPH/3oR3z7298G4rPiTj75ZPbs2cOECRMA+NnPfpayvZkzZ/Lggw/yve99DzPjj3/8I5/73OeYPn06v/zlL5kxYwZbt25ly5YtGfdfs+ZERArAAw88wB133MFZZ53FJZdcwsKFCzn++ONT1p8xYwbvvPMOZ511Fk888cRh+2+//XZef/11pk+fzhlnnMFDDz0EwPXXX8+SJUuYNWsWBw8eTNneTTfdxIEDBzjzzDOZNm0a3/3udwGYN28era2tTJ8+nXvuuYfTTjst475b/Fl2haG5uTmnwaR6Zn0xUl8yk2xVhGyor9/JxIkT+l2/kFdWKKXPGKTuT3NzM2PGjAkQUfra2tqoqKgIHUavevt3HTNmzCHn9jQiEhGRoJSIREQkKCUiEREJSolIRESC6jMRmdlyM2sws7cSyu4ws21m9oaZ/aeZjU3Yt9jMas3sHTOblVB+YVRWa2aLst8VEREpRv0ZEf0E6Ln40lrgRHc/GXgXWAxgZp8FrgBOiN7zQzMrM7My4D6gGvgs8NWoroiIDHJ93tDq7i+a2XE9yhLXCt8EfDl6PRd4zN33A38ys1pgSrSv1t3fAzCzx6K6WzOKXqTIZWMaeSFPAZe+NTU1MWfOHAAaGhooKytj/PjxQHyZnWHDhg24zQ0bNjBs2DDOOOOMAb3vpJNOYv369d3Hz5dsrKzwdeDn0eujiSemLnVRGcCOHuUD+xcSESlBlZWVvPTSS0B84dKRI0fyrW99q3t/R0cHQ4cO7E/1Sy+9xMiRIweciELJKBGZ2XeADqBrDYpkCxA5yU8B9nrzaiwWyyS0vLcbgvqSvvrGsty1Xb8zZ20nE9t/sO9K6bZdQp8xSN6fioqK7rXcAJY+9o2sHvM7V/xHv+t2dHRw4MABrr76asaNG8ebb77JSSedxJVXXsnixYtpbGxk+PDhfP/736eqqoo1a9Zw1113ceDAAcaNG8cPf/hD9u3bx/LlyykrK+Oxxx5j6dKlfPrTn2bhwoXdj4S45ZZbmDJlCk1NTVx77bU0NjZyyimn0NnZSVtbW1Ye7bBnzx4aGhq6t3u7OTrtRGRmNcDFwEz/2/IMdcDkhGqTgProdarypHJxR3cp3SmuvmSmtrwwVlbIhqocnZorpc8Y9L6yQi5XKRhI20OHDuWII46grKyMP//5zzz11FOUlZUxZ84c7rrrLj71qU/x+9//nkWLFvH0008zY8YMLrnkEsyMFStW8MADD7B06VK+/vWvHzKyuuqqq7juuusOeyTE3XffzfTp01m4cCHPP/88jzzyCBUVFVn59xg9ejSTJ0/uuyJpJiIzuxBYCJzj7h8n7HoS+JmZ3QlMBKqAV4iPlKrM7HjgA+ITGv4pnWOLiAwGc+fOpaysjL179/LKK69QU1PTvW///v0AfPDBB1x55ZV8+OGHtLe3d6/W3VOqR0Js3LiRn/70pwDMmjWLsWPHJn1/rvWZiMzsUeBc4EgzqwNuJj5LrhxYGy0Hvsndr3H3LWb2C+KTEDqA+e5+MGrnOuB5oAxY7u6ZL9kqIlKiuh7p0NnZyZgxY7qvIwHdp84WLFjA/PnzmT17Nhs2bOC2225L2laqR0IUiv7MmvtqkuJlvdRfCixNUv4M8MyAohMRybNb/ufDoUM4xOjRozn22GP51a9+xaWXXoq7s2XLFk477TT27NnDxIkTAXj00Ue739PzkQ6pHgnR9UiHm266ibVr1x72TKR80coKIiIF7sEHH+SRRx7hzDPPZOrUqTz3XPyZSYsWLaKmpobq6upDplxXV1fz61//mrPOOouNGzemfCTEokWL2LhxIzNmzGDdunVMmjQpSP/0GIgipb5kplAeA5ENubqPqJQ+Y6DHQOSbHgMhIiJFQ4lIRESCUiISEZGglIhERCQoJSIRGdSGDBlCe3t76DBKSnt7O0OG9D+9ZGPRUxGRojVy5Ej27t3Lvn25mUmZC3v27GH06NGhw0hpyJAhjBw5st/1lYik6ORq6rUMTmbGqFGjQocxIA0NDf1ex60Y6NSciIgEpUQkIiJBKRGJiEhQSkQiIhKUEpGIiASlRCQiIkEpEYmISFBKRCIiEpQSkYiIBKVEJCIiQSkRiYhIUH0mIjNbbmYNZvZWQlmlma01s1j0fVxUbmb2AzOrNbM3zOzUhPfURPVjZlaTm+6IiEix6c+I6CfAhT3KFgEvuHsV8EK0DVANVEVfVwP3QzxxATcDZwBTgJu7kpeIiAxufSYid38RaOpRPBd4OHr9MHBpQvkKj9sEjDWzCcAsYK27N7n7LmAthyc3EREZhNK9RnSUu+8EiL5/Mio/GtiRUK8uKktVLiIig1y2n0dkScq8l/KUYrFYVgLKV7shDNa+1DeW5TCSzNXX78zr8WL7D+au7RL6jEFp9afY+lJVVZVyX7qJ6EMzm+DuO6NTbw1ReR2Q+LSmSUB9VH5uj/L1vR2gt6DTFYvFctJuCIO5L7XlhftgvPr6nUycOCGvx6zNQhvVxww/rKyUPmNQWv0ppb5A+qfmngS6Zr7VAKsTyr8WzZ6bCjRHp+6eBy4ws3HRJIULojIRERnk+hwRmdmjxEczR5pZHfHZb7cBvzCzecB24PKo+jPAbOL/SfsYuBLA3ZvM7Bbg1ajeEnfvOQFCREQGoT4Tkbt/NcWumUnqOjA/RTvLgeUDik5EREqeVlYQEZGglIhERCQoJSIREQlKiUhERILK9g2tIr16dvvh9wDVN5YV9L1BIpJbGhGJiEhQSkQiIhKUEpGIiASlRCQiIkEpEYmISFBKRCIiEpQSkYiIBKX7iEQila2v9l2pF00jvpClSEQGF42IREQkKCUiEREJSolIRESC0jUiKRmZXuMRkTA0IhIRkaCUiEREJCglIhERCUqJSEREgsooEZnZDWa2xczeMrNHzazCzI43s5fNLGZmPzezYVHd8mi7Ntp/XDY6ICIixS3tRGRmRwPfBk539xOBMuAK4HbgLnevAnYB86K3zAN2ufungbuieiIiMshlempuKDDczIYCnwB2AucBj0f7HwYujV7PjbaJ9s80M8vw+CIiUuTM3dN/s9n1wFJgH7AGuB7YFI16MLPJwLPufqKZvQVc6O510b7/B5zh7h91tdfc3NwdTCwWSzsuKVwvNpblrO2JHW/krO3+qB96ctDjZ2LG+IOhQ5ASV1VV1f16zJgxhwxC0r6h1czGER/lHA/sBn4JVCep2pVcko1+UmbBxKCzJRaL5aTdEIq1L7Xl+w4rq6/fycSJEzJuu7K1LuM2MjJiQtb6km9Vxww/rKxYP2OplFJ/SqkvkNnKCucDf3L3vwKY2RPAdGCsmQ119w5gElAf1a8DJgN10am8MUBTBscXKSiVra9CR1NaCVErd8tglsk1ou3AVDP7RHStZyawFfgN8OWoTg2wOnr9ZLRNtH+dZ3JeUERESkLaicjdXyY+6eAPwJtRWw8CC4EbzawWGA8si96yDBgfld8ILMogbhERKREZLXrq7jcDN/cofg+YkqRuG3B5JscTkdx4dnuSa3eNZUmv6aVSneQ6k0h/aGUFEREJSolIRESCUiISEZGglIhERCQoJSIREQlKiUhERILKaPq2iGRHZeurGb1fKzNIMdOISEREglIiEhGRoJSIREQkKF0jkoKR6XUSESlOGhGJiEhQSkQiIhKUEpGIiASla0TSL8keEyAikg0aEYmISFBKRCIiEpQSkYiIBKVEJCIiQSkRiYhIUBklIjMba2aPm9k2M3vbzKaZWaWZrTWzWPR9XFTXzOwHZlZrZm+Y2anZ6YKIiBSzTEdE9wDPufs/AJ8D3gYWAS+4exXwQrQNUA1URV9XA/dneGwRESkBaSciMxsNzACWAbh7u7vvBuYCD0fVHgYujV7PBVZ43CZgrJlNSDtyEREpCZmMiP4e+CvwkJm9bmY/NrMRwFHuvhMg+v7JqP7RwI6E99dFZSIiMohlsrLCUOBU4Fvu/rKZ3cPfTsMlY0nKPFXlWCyWQWip5ardEPLZl/rGsty2X78TOppyeox8aWrMfz/qm3fmpt36/rcb238wJzFkk37/w6mqqkq5L5NEVAfUufvL0fbjxBPRh2Y2wd13RqfeGhLqT054/ySgPp2g0xWLxXLSbgj57kttee6W+Kmv38nEiROobK3L2THypamxicrxlfk/8Ijsn+Xu+rn0V9Uxw7MeQzbp979wpX1qzt3/Auwws89ERTOBrcCTQE1UVgOsjl4/CXwtmj03FWjuOoUnIiKDV6aLnn4LWGlmw4D3gCuJJ7dfmNk8YDtweVT3GWA2UAt8HNWVEpPWw+06mkpiNCQi6ckoEbn7ZuD0JLtmJqnrwPxMjiciIqVHKyuIiEhQSkQiIhKUEpGIiASlRCQiIkEpEYmISFCZTt8WkQKQ1rT5BE0jvpClSEQGTolIRLLi2e2Zr75RXeCrM0hu6NSciIgEpUQkIiJBKRGJiEhQSkQiIhKUEpGIiASlRCQiIkEpEYmISFBKRCIiEpQSkYiIBKVEJCIiQSkRiYhIUEpEIiISlBY9lUN0fLQpaXlla3ueIxGAD3dtz7iNo8Ydk4VIRHJHIyIREQkq4xGRmZUBvwc+cPeLzex44DGgEvgD8C/u3m5m5cAK4DSgEfiKu/850+OLSOaSPs+oo4nK1rp+vV/PM5JMZOPU3PXA28DoaPt24C53f8zMHgDmAfdH33e5+6fN7Iqo3leycHzpw0CeE6NTcIca6KmxltYWDgzZ272t02Iifcvo1JyZTQIuAn4cbRtwHvB4VOVh4NLo9dxom2j/zKi+iIgMYpmOiO4GFgCjou3xwG5374i264Cjo9dHAzsA3L3DzJqj+h8lazgWi2UYWnK5ajeE/valvrGs/412NKUZTWaaGsMcty8trS0Df0/L395zRGdm/Urn+D1lEkN/fy71zTvTPkai2P6DWWknZfuD8Pe/UFRVVaXcl3YiMrOLgQZ3f83Mzu0qTlLV+7HvML0Fna5YLJaTdkMYSF9qywdyaq5/1wSyqamxicrxlXk/bn8knmbrj5aWFkaNGtW9XTkus34N9PjJpBvDgH4uIyakdYyeqnL4qPDB+vtfDDIZEZ0JzDGz2UAF8WtEdwNjzWxoNCqaBNRH9euAyUCdmQ0FxgCF+d9gkSzJxvRrkVKX9jUid1/s7pPc/TjgCmCdu/8z8Bvgy1G1GmB19PrJaJto/zp3TzkiEhGRwSEX9xEtBG40s1ri14CWReXLgPFR+Y3AohwcW0REikxWVlZw9/XA+uj1e8CUJHXagMuzcTwRESkdWllBRESC0lpzIlIwBnLzdSrVOZx5J7mhEZGIiASlEZGUNE2fFil8SkQiJS7dZNy1bp7Wy5Nc06k5EREJSolIRESCUiISEZGglIhERCQoTVYQ6YeZbT89ZPuFiv+R9WO8+u7aQ7a/8N++mPVjiBQiJSIR6VV/Zt217CtPua/yyJOzGY6UIJ2aExGRoJSIREQkKCUiEREJSolIRESC0mQFEcnYqLZYyn2Vrft7fW/TiC9kOxwpMhoRiYhIUBoRlZiOjzYdVlbZ2h4gkuzQ6tkipU8jIhERCUqJSEREgko7EZnZZDP7jZm9bWZbzOz6qLzSzNaaWSz6Pi4qNzP7gZnVmtkbZnZqtjohIiLFK5MRUQfwv9z9H4GpwHwz+yywCHjB3auAF6JtgGqgKvq6Grg/g2OLiEiJSHuygrvvBHZGr1vM7G3gaGAucG5U7WFgPbAwKl/h7g5sMrOxZjYhakdSeHb7vqTl9Y1l1JYfvq+YJyaIyOCUlVlzZnYccArwMnBUV3Jx951m9smo2tHAjoS31UVlSkQiJayvmY+9LZgKWjR1MMg4EZnZSGAV8K/uvsfMUlZNUuapKsdiqW+Qy0Su2s2V+say1Pvqk+TwjqYcRpM7TY3J425pbclzJCkccehmS0vquHrbNxDZaqfQY2hq7f0z29Y+sP+rLqtPtaeMFxvf61cbM8YfHNAxQyi2v2VVVVUp92WUiMzsCOJJaKW7PxEVf9h1ys3MJgANUXkdMDnh7ZOAlB+Z3oJOVywWy0m7uZTs9BvEk9DEiRMOK69srct1SFnX1NhE5fjKpPsODNmb52hSaDt0c9SoUUmrtbS0pNzXpx5/b9NuJ0sy6stAVCT/2XepPPLwz3k6Uv3OJFN1zPCsHDNXivFvWW/STkQWH/osA9529zsTdj0J1AC3Rd9XJ5RfZ2aPAWcAzbo+VNr6ezNqS2tL4SQcEcm7TEZEZwL/ArxpZpujsn8nnoB+YWbzgO3A5dG+Z4DZQC3wMXBlBscWEZESkcmsuZdIft0HYGaS+g7MT/d4IlKaelswFbRo6mCglRVERCQoJSIREQlKq2+LSEHL9D4k0L1IhU4jIhERCUojogLT83lCKZfs6WgqynuGRER60ohIRESCUiISEZGgdGpOUtJjukUkHzQiEhGRoDQiyqFUzxLqjZ4nJDIwfa3MANHqDCkm+GhlhvA0IhIRkaA0IhKRkvfhru0pV3lPdkPsNht2yPY/TD4lZ7GJRkQiIhKYRkQiIj289tdDr9X+yQd+vbe6wB+uV0g0IhIRkaA0IsqBrmV6Qs+A031AIlIMNCISEZGglIhERCQonZoTkUGtPzfEHqhLXaeloippeeIUcE3/7p0SkYhIDiTOvEtn1l2XwTD7TqfmREQkqLyPiMzsQuAeoAz4sbvflu8Y+iOddeK6ZGu2XG+z3lLdJS4iUmzymojMrAy4D/giUAe8amZPuvvWfMbRl46PNgWfei0ixSGTa0ypri8lWvkRVB558iFl9Y1l1Jb3/z/LhX56z9w9fwczmwb8H3efFW0vBnD3WwGam5vzF4yIiAQxZswYS9zO9zWio4EdCdt1UZmIiAxS+U5ElqRMoyARkUEs35MV6oDJCduTgPqujZ7DNRERKX35HhG9ClSZ2fFmNgy4AngyzzGIiEgByWsicvcO4DrgeeBt4BfuviWfMZjZLWb2hpltNrM1ZjYxn8fPJjO7w8y2Rf35TzMbGzqmdJnZ5Wa2xcw6zez00PGkw8wuNLN3zKzWzBaFjiddZrbczBrM7K3QsWTKzCab2W/M7O3o83V96JjSZWYVZvaKmf0x6sv/DR1TtuR11lwhMLPR7r4nev1t4LPufk3gsNJiZhcA69y9w8xuB3D3hYHDSouZ/SPQCfwH8G/u/vvAIQ1IdGvCuyTcmgB8tdBuTegPM5sB7AVWuPuJoePJhJlNACa4+x/MbBTwGnBpkf5cDBjh7nvN7AjgJeB6d98UOLSMDbqVFbqSUGQERTxZwt3XRKNMgE3Er7kVJXd/293fCR1HBqYAte7+nru3A48BcwPHlBZ3fxFoCh1HNrj7Tnf/Q/S6hfiZmKKcqetxXXexHxF9Fe3fr0SDLhEBmNlSM9sB/DPwv0PHkyVfB54NHcQgplsTCpyZHQecArwcNpL0mVmZmW0GGoC17l60fUlUkonIzP7LzN5K8jUXwN2/4+6TgZXEr1kVrL76EtX5DtBBvD8Fqz99KWK6NaGAmdlIYBXwrz3OihQVdz/o7p8nfvZjipkV9anTLiW5+ra7n9/Pqj8DngZuzmE4GemrL2ZWA1wMzPQCv+A3gJ9LMer11gQJJ7qesgpY6e5PhI4nG9x9t5mtBy4Ein5SSUmOiHpjZomLO80BtoWKJVPRArILgTnu/nHoeAY53ZpQgKIL/MuAt939ztDxZMLM/q5rZqyZDQfOp4j/fiUajLPmVgGfIT5D633gGnf/IGxU6TGzWqAcaIyKNhXxDMD/DtwL/B2wG9jctSZhsTCz2cDdxFeWX+7uSwOHlBYzexQ4FzgS+BC42d2XBQ0qTWZ2FrABeJP47zzAv7v7M+GiSo+ZnQw8TPzzNYT47S9LwkaVHYMuEYmISGEZdKfmRESksCgRiYhIUEpEIiISlBKRiIgEpUQkIiJBKRGJiEhQSkQiIhKUEpGIiAT1/wFC5+/pW55thAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(data[\"achievement_score\"], bins=20, alpha=0.3, label=\"All\")\n",
    "plt.hist(data.query(\"intervention==0\")[\"achievement_score\"], bins=20, alpha=0.3, color=\"C2\")\n",
    "plt.hist(data.query(\"intervention==1\")[\"achievement_score\"], bins=20, alpha=0.3, color=\"C3\")\n",
    "plt.vlines(-0.1538, 0, 300, label=\"Untreated\", color=\"C2\")\n",
    "plt.vlines(-0.1538+0.4723, 0, 300, label=\"Treated\", color=\"C3\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we still think this result is biased. The difference between treated and untreated is probably smaller than this, because we think the bias is positive. We've already seen that more ambitious people are more willing to go to the seminar, so they probably would have achieved more even if they had attended it. To control for this bias, we could use regression or matching, but it's time to learn about a new technique.\n",
    "\n",
    "## Propensity Score\n",
    "\n",
    "Propensity score comes from the realisation that you don't need to directly control for confounders X to achieve conditional independence \\\\((Y_1, Y_0) \\perp T | X\\\\). Instead, it is sufficient to control for a balancing score \\\\(E[T|X]\\\\). This balancing score is often the conditional probability of the treatment, \\\\(P(T|X)\\\\), also called the propensity score \\\\(P(x)\\\\). The propensity score makes it so that you don't have to condition on the entirety of X to achieve independence of the potential outcomes on the treatment. It is sufficient to condition on this single variable with is the propensity score:\n",
    "\n",
    "$\n",
    "(Y_1, Y_0) \\perp T | P(x)\n",
    "$\n",
    "\n",
    "There is a formal proof for why this is, but we can forget it for now and approach the matter in a more intuitive way. The propensity score is the conditional probability of receiving the treatment right? So we can think of it as some sort of function that converts X into the treatment T. The propensity score makes this middle ground between the variable X and the treatment T. If we show this in a causal graph, this is what it would look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"90pt\" height=\"260pt\"\n",
       " viewBox=\"0.00 0.00 90.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-256 86,-256 86,4 -4,4\"/>\n",
       "<!-- T -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"55\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"55\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">T</text>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Y</text>\n",
       "</g>\n",
       "<!-- T&#45;&gt;Y -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>T&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M48.222,-72.5708C45.0128,-64.3187 41.1065,-54.2738 37.5242,-45.0623\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"40.6992,-43.5697 33.8126,-35.5182 34.1752,-46.1069 40.6992,-43.5697\"/>\n",
       "</g>\n",
       "<!-- X -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"27\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-229.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">X</text>\n",
       "</g>\n",
       "<!-- X&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>X&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M23.6685,-216.0736C21.8913,-205.5982 19.8913,-192.0982 19,-180 15.4733,-132.1297 15.4733,-119.8703 19,-72 19.6267,-63.4935 20.8015,-54.2939 22.0616,-45.9399\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"25.5398,-46.3546 23.6685,-35.9264 18.6282,-45.2455 25.5398,-46.3546\"/>\n",
       "</g>\n",
       "<!-- P(x) -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>P(x)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"55\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"55\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">P(x)</text>\n",
       "</g>\n",
       "<!-- X&#45;&gt;P(x) -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>X&#45;&gt;P(x)</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M33.778,-216.5708C36.9872,-208.3187 40.8935,-198.2738 44.4758,-189.0623\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"47.8248,-190.1069 48.1874,-179.5182 41.3008,-187.5697 47.8248,-190.1069\"/>\n",
       "</g>\n",
       "<!-- P(x)&#45;&gt;T -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>P(x)&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M55,-143.8314C55,-136.131 55,-126.9743 55,-118.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"58.5001,-118.4132 55,-108.4133 51.5001,-118.4133 58.5001,-118.4132\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a1bae39b0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "g.edge(\"T\", \"Y\")\n",
    "g.edge(\"X\", \"Y\")\n",
    "g.edge(\"X\", \"P(x)\")\n",
    "g.edge(\"P(x)\", \"T\")\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If I know what P(x) is, X alone tells me nothing more that can help me learn what T would be. Which means that controlling for P(x) acts the same way as controlling for X directly. Think of it in terms of our mindset program. Treated and non treated are initially not comparable because the more ambitious are both more likely to take the treatment and of achieving more in life. However, if I take 2 individuals, one from the treated and one from the control, but with the same probability of receiving the treatment, they are comparable. Think about it. If they have the exact same probability of receiving the treatment, the only reason one of them received it and the other did not is pure chance. Holding the propensity score constant acts in a way of making the data look as good as randomly assigned. \n",
    "\n",
    "Now that we got the intuition, let's look at the mathematical proof. We want to show that \\\\((Y_1, Y_0) \\perp T | P(x)\\\\) is equivalent to saying that \n",
    "\n",
    "$\n",
    "E[T|P(x), X] = E[T|P(x)] \n",
    "$\n",
    "\n",
    "This simply says that once I condition on P(x), X can give me no extra information about T. The proof of this is quite weird. We will show that the equation above is true by converting it to a trivial statement. First take a look at the left hand side \\\\(E[T|P(x), X]\\\\).\n",
    "\n",
    "$\n",
    "E[T|P(x), X] = E[T|X] = P(x)\n",
    "$\n",
    "\n",
    "We use the fact that P(x) is just a function of X, so conditioning on it gives no further information after we've conditioned on X itself. Then, we use the definition of the propensity score \\\\(E[T|X]\\\\). For the right hand side, we will use the law of iterated expectations \\\\(E[A] = E[E[A|B]]\\\\) This law says that we can compute the expected value of A by looking at the value of A broken down by B and then averaging that. \n",
    "\n",
    "$\n",
    "E[T|P(x)] = E[E[T|P(x),X]|P(x)] = E[P(x)|P(x)] = P(x)\n",
    "$\n",
    "\n",
    "The first equality comes from the law of iterated expectations. The second comes from what we've figured out when dealing with the left hand side. Since both the left and right hand side equals, \\\\(P(x)\\\\), this equation is trivially true.\n",
    "\n",
    "## Propensity Weighting\n",
    "\n",
    "![img](./data/img/ps/balance.png)\n",
    "\n",
    "OK, we got the propensity score. Now what? Like I've said, all we need to do is condition on it. For example, we could run a linear regression that conditions only on the propensity score, instead of all the Xs. For now, let's look at a technique that just uses the propensity score and nothing else. The idea is to write the conditional difference in means with the propensity score\n",
    "\n",
    "$\n",
    "E[Y|X,T=1]−E[Y|X,T=0] = E\\bigg[\\dfrac{Y}{P(x)}|X,T=1\\bigg]P(x) - E\\bigg[\\dfrac{Y}{(1-P(x))}|X,T=0\\bigg](1-P(x))\n",
    "$\n",
    "\n",
    "We can simplify this further, but let's take a look at it like this because it gives us some nice intuition of what the propensity score is doing. The first term is estimating \\\\(Y_1\\\\). It is taking all those that are treated and scaling them by the inverse probability of treatment. What this does is makes those with very low probability of treatment have a high weight. This makes sense right? If someone has a low probability of treatment, that individual looks like the untreated. However, that same individual was treated. This must be interesting. We have a treatment that looks like the untreated, so we will give that entity a high weight. What this does is create a population with the same size as the original, but where everyone is treated. By the same reasoning, the other term looks at the untreated and gives a high weight to those that look like the treated. This estimator is called the Inverse Probability of Treatment Weighting, since it scales each unit by the probability of receiving some treatment other than the one it received. \n",
    "\n",
    "In a picture, here is what this weightin does.\n",
    "\n",
    "![img](./data/img/ps/iptw.png)\n",
    "\n",
    "The upper left plot shows the original data. The blue dots are the untreated and the red dots are the treated. The bottom plot shows the propensity score P(x). Notice how it is between 0 and 1 and it grows as X increases. Finally, the upper right plot is the data after weighting. Notice how the red (treated) that are more to the left (lower propensity score) have a higher weight. Similarly, the blue plots that are to the right have also a higher weight. \n",
    "\n",
    "Not that we got the intuition, we can simplify the terms above to\n",
    "\n",
    "$\n",
    "E\\bigg[Y \\dfrac{T-P(x)}{P(x)(1-P(x))}\\bigg|X\\bigg]\n",
    "$\n",
    "\n",
    "which if we integrate over X becomes our propensity score weighting estimator.\n",
    "\n",
    "$\n",
    "E\\bigg[Y \\dfrac{T-P(x)}{P(x)(1-P(x))}\\bigg]\n",
    "$\n",
    "\n",
    "Notice that this estimator requires that neither \\\\(P(x)\\\\) nor \\\\(1-P(x)\\\\) be strictly positive. In words, this means that everyone needs to have at least some chance of receiving the treatment and of not receiving it. Another way of stating this is that the treated and untreated distributions overlap. This is the **positivity assumption** of causal inference. It also makes intuitive sense. If treated and untreated don't overlap, it means they are very different and I won't be able to extrapolate the effect of one group to the other. This extrapolation is not impossible (regression does it), but it is very dangerous. It is like testing a new drug in an experiment where only men receive the treatment and then assume women will respond to it equally well.\n",
    "\n",
    "\n",
    "## Propensity Score Estimation\n",
    "\n",
    "In an ideal world, we would have the true propensity score \\\\(P(x)\\\\). However, in practice, the mechanism that assigns the treatment is unknown and we need to replace the true propensity score by an estimation of it \\\\(\\hat{P}(x)\\\\). One common way of doing so is using logistic regression but other machine learning methods, like gradient boosting, can be used as well (although it requires some additional steps to avoid overfitting). \n",
    "\n",
    "Here, I'll stick to logistic regression. This means that I'll have to convert the categorical features in the dataset to dummies. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10391, 32)\n"
     ]
    }
   ],
   "source": [
    "categ = [\"ethnicity\", \"gender\", \"school_urbanicity\"]\n",
    "cont = [\"school_mindset\", \"school_achievement\", \"school_ethnic_minority\", \"school_poverty\", \"school_size\"]\n",
    "\n",
    "data_with_categ = pd.concat([\n",
    "    data.drop(columns=categ), # dataset without the categorical features\n",
    "    pd.get_dummies(data[categ], columns=categ, drop_first=False)# dataset without categorical converted to dummies\n",
    "], axis=1)\n",
    "\n",
    "print(data_with_categ.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's estimate the propensity score using logistic regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>intervention</th>\n",
       "      <th>achievement_score</th>\n",
       "      <th>propensity_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.277359</td>\n",
       "      <td>0.315560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>-0.449646</td>\n",
       "      <td>0.263827</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0.769703</td>\n",
       "      <td>0.344060</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>-0.121763</td>\n",
       "      <td>0.344060</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1.526147</td>\n",
       "      <td>0.367804</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   intervention  achievement_score  propensity_score\n",
       "0             1           0.277359          0.315560\n",
       "1             1          -0.449646          0.263827\n",
       "2             1           0.769703          0.344060\n",
       "3             1          -0.121763          0.344060\n",
       "4             1           1.526147          0.367804"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "T = 'intervention'\n",
    "Y = 'achievement_score'\n",
    "X = data_with_categ.columns.drop(['schoolid', T, Y])\n",
    "\n",
    "ps_model = LogisticRegression(C=1e6).fit(data_with_categ[X], data_with_categ[T])\n",
    "\n",
    "data_ps = data.assign(propensity_score=ps_model.predict_proba(data_with_categ[X])[:, 1])\n",
    "\n",
    "data_ps[[\"intervention\", \"achievement_score\", \"propensity_score\"]].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we can make sure that the propensity score weight indeed reconstructs a population where everyone is treated. By producing weights \\\\(1/P(x)\\\\), it creates the population where everyone is treated and by providing the weights  \\\\(1/(1-P(x))\\\\) it creates the population where everyone is untreated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original Sample Size 10391\n",
      "Treated Population Sample Size 10388.538689767902\n",
      "Untreated Population Sample Size 10391.464663106926\n"
     ]
    }
   ],
   "source": [
    "weight_t = 1/data_ps.query(\"intervention==1\")[\"propensity_score\"]\n",
    "weight_nt = 1/(1-data_ps.query(\"intervention==0\")[\"propensity_score\"])\n",
    "print(\"Original Sample Size\", data.shape[0])\n",
    "print(\"Treated Population Sample Size\", sum(weight_t))\n",
    "print(\"Untreated Population Sample Size\", sum(weight_nt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the propensity score to find evidence of confounding. If a segmentation of the population has a higher propensity score than another, it means that something that is not random is causing the treatment. If that same thing is also causing the outcome, we have confounding. In our case, we can see that students that reported to be more ambitious also have a higher probability of attending the growth mindset seminar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEXCAYAAADMVxF8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfZyVdZ3/8dcbCDUzHSNBRRNrrLVSQQVJ0NHIdG2xX7NtylqitK1bZnaz3fozNNufua2JrltuZLqWUjm7xbZu6ixMQKEComt4E6goI4SiI95x0wyf3x/XNXTmeGbmnOGcc805834+Hucx53vdfq5zYD7z/V7f6/tVRGBmZlZrhmUdgJmZ2UA4gZmZWU1yAjMzs5rkBGZmZjXJCczMzGqSE5iZmdUkJzAbEiSNkHSDpOckhaSmQRDTbElreivXgvSzPLufbdokza1WTDZ0OIFZZiS9SdKVkh6VtFXSM5IWSfqYpBFlPl0zMAP4C2B/4LdlPn45fBs4rhonSpNl9PIaVcKh9gduq1ScZn0p9y8Js6JIGgv8BugELgFWAn8E3gN8Afhf4P4ynrIReDoiBmPiAiAiXgZeruIp1wKTCyx/rtgDRMQfyhaNWYlcA7OsfBfYDZgQET+OiIciYnVE3AQcDawGkPQ6SVdIelrSdkkPSZqRe6C01vBJSTdLeknSOklfzFnfBnwDODTddm2Jxz47b1mrpBtzymslXSZpjqTnJW2U9G1Jw3O22U3SdyVtltQhqfv6c49bsElR0hmSHpH0iqSFkt6at99Zkh5La7G/lfSBNO4p/XwHXRHxhwKvkPQ+SV2SDso710fS8+xT6POR9BZJv5K0RdJTkj6df9K0OXe2pCfSY62S9LcFPvdev9Oc41ySXvu29Hu8Nmf9G9Lv5GlJr0paKelD/XwmVksiwi+/qvoC9gW6gIuL2PYfSWoEHwYOA74K7ADem7NNABuBvwHeClyYLjsp53zfBp4AxgBvLvHYZ+fF1ArcmFNeC3QAXyap6X2EpGZ5bs423wGeAc4A3pHG8yKwJmeb2QXKrwC/IknqR5LUVH+ds83RacyXA28HPgisSeOe0sfn2uNcBdYPA9qBr+Qt/yXwk0KfDyDgPmAZMAk4Crgrvc65OfvcSFLDPgUYl35eLwCziv1O021uSj/Tj6bbHAd8NieWhUAbMAU4FPgEsD33+/Wrtl+ZB+DX0HsBE9NfRh/qZ7vXA9uAT+Yt/w9gQU45gGvytnkE+H855fzkUMqxi0lg8/O2+RVwa/p+T2Ar8Dd52ywvIoF1kibcdNmZacLaPS3/GFicd9zzi0xgO0iaLHNfD+RscwXwUE55P5Jm3tMLfT7AtLR8WM76NwNbSBNYmrB2AO/Ii+cS4P5iv1Pgbek2f9nL9TWln/neectvAH6e9f8Bv8rz8j0wy4LSn/2NJP02YCSwKG/5r4Gv5C3Lv1/2NDC6TMcuRqHzj0vfv5WkuTD//tsS4AP9HHd9RDybd1yRJJOngMNJEmqupUXGvA54b96y7TnvbwK+JOnYiFgGnEVSY72jl+MdDmyKiN93L4iIZyU9mrPNMWn8yyXl7juCpFaeq6/vdEL6885eYjmW5Pt9Ou88I0mbp632OYFZFlaT/BX+TpIaT3/yE50KLNueVw6Ku8fb37GDPyXcbq8rcJy+zl9swi6k0HGh57UNdEqJP0ZEr932I+JhScuBj5E0C34MuCUiOnvZpdD3kq877vcAr+afMq880O+0+zybSRJZvvzjWo1yJw6ruoh4Hvhv4AJJe+evTztX7ElyL2cbcGLeJicAq3YxjGKP/QxwQE5su5HUNEo913bg+Lzl7ynxOIU8xGt7EpazK/6/AWdKOpKk1nNTH9uuAt4sqbF7Qdol/7CcbVakPw+OiDV5r8dKiOu+9OcpvaxfDuxD0tSaf56nSjiPDWKugVlWPknSjX6FpEtImou2k/zy/XvgnIi4X9I1wDckPZtu82GSjhDv25WTR8SrRR67FThf0iLgJeBrJM1QpZzrFUnfAy6XtBF4FJhF0pnjmV25DuAqYJmky4Afpcf8fPep+9l3uKQxBZZvyqll3Qr8E2nHi4h4oI/j/Q/wAPCjtPfhduBbJPfxkoAi1ki6Afh+2qtwKck9wqNJ7vV9q5+Yc4/zY+BfJO2eHmdf4D0RMQdYQPLd/bukL6VxNZD80bA1Ir5fzHlscHMCs0xExFOSJpD03JsNHEzSW+1hkt6Bv0s3/RpJc+PVJB0C1pB0GvifMoRRzLG/AHyf5L7PZuAf0m1L9WVgd+DmtPwT4DqSpDlgEbFC0l+T9EL8EknN5OL0+Fv72f0QYEOB5ceS1GCIiE2S/oukd+MX+oklJH0Q+FeSe4ubSL7L3fI2/QRJkv0aSe/AF0lqb//cT7z5ziXp/HE5SS35GdKHqtNYpgNfJ0nyBwLPk/yhcmWJ57FBShGekdmsnkj6GPBD4E0R8ULW8ZhVimtgZjVO0hdInnl6nqT29C3gZ05eVu+cwMxq3xEkTXL7knSN/xFJ05lZXXMTopmZ1aSqdaOXdKqSUcfXSPpygfUzJT0r6f709fGcdV05y+dXK2YzMxu8qlIDSwc1/T1J9+R2kociz4qIh3K2mQkcExEXFNj/5Yh4Q2/H37x5s6uRZmZ1bu+99+4xqEC1amATScZ4ezwitgPzSJ63MTMzG5BqJbADSW4ud2tPl+VrlvS/km7Lm8Zhd0nLJd2dPmdiZmZDXLV6IeaPJQevHSXgP0lG794m6XySIWtOTtcdHBHrJR0KLJD0YG/Dzqxe7XE6zczqRWNjY6/rqpXA2oHcGtVYYH3uBhGROwvs90meZeletz79+biSyQnHAwUTWF8Xa2Zm9aNaTYjLgEZJ4ySNJJnTqEdvQkn75xSnkwwphKSGdADV7oFBjycZwNTMzIawqtTAIqJT0gUk48kNB26IiFXpAKTLI2I+cGE6dlknyYgCM9Pd/wy4XtIOkoR7RW7vRTMzG5rq4kFmd6M3s8Goo6ODOXPmcNFFF7HPPvtkHU7Ny6obvZnZkNPS0sIjjzxCS0tL1qHUJScwM7MK6OjooK2tjYigra2NF17w2Mrl5gRmZlYBLS0tdN+i2bFjh2thFeAEZmZWAUuWLKGzM5mMurOzk8WLF2ccUf1xAjMzq4ApU6YwYkTS0XvEiBFMnTo144jqjxOYmVkFNDc3IyWd5oYNG0Zzc3PGEZVPR0cHs2fPzvy+nhOYmVkFNDQ00NTUhCSamprqqhv9YOld6QRmZlYh06ZNY/fdd2fatGlZh1I2g6l3pROYmVmFtLa2snXrVlpbW7MOpWwGU+9KJzAzswro6Ohg4cKFRAQLFy7M/H5RuQym3pVOYGZmFdDS0kJXVxeQ/KLP+n5RuQym3pVOYGZmFbB48eKdTW0RwaJFizKOqDwGU+9KJzAzswoYNWpUn+VaNZh6V1ZrQkszsyFl06ZNfZZrWXNzM+3t7Zk/2+YamJlZBUycOLFHedKkSRlFUn4NDQ3Mnj0782fbnMDMzKwmVS2BSTpV0qOS1kj6coH1MyU9K+n+9PXxnHXnSFqdvs6pVsxmZgO1bNmyHuV77703o0jqV1USmKThwHXAacDhwFmSDi+w6U8i4qj0NTfdd1/g68AkYCLwdUkN1YjbzGygpkyZwvDhwwEYPnx4XQ3mO9TGQpwIrImIxyNiOzAPOKPIfd8P3BURz0dEB3AXcGqF4jQzK4vm5maGDUt+xQ4fPjzzDg/lNFjGQqxWL8QDgXU55XaSGlW+ZkknAL8HPhsR63rZ98DeTrR69epdj9bMrAzGjx/PsmXLOOqoo3j22Wd59tlnsw5pl7300ks7RxhZsGAB48ePZ6+99qrY+RobG3tdV60EpgLLIq/8n8CtEbFN0vnATcDJRe67U18Xa2ZWTeeddx4vvfQSs2bNyrzHXrnMnTu3xwPaK1euZNasWZnEUq0mxHbgoJzyWGB97gYR8VxEbEuL3weOLnZfM7PBaLB0Ny+nJUuW7Bwiq6ura0iMhbgMaJQ0TtJI4Exgfu4GkvbPKU4HHk7f3wGcIqkh7bxxSrrMzMyq7Nhjj+1Rzn/erZqq0oQYEZ2SLiBJPMOBGyJilaTLgOURMR+4UNJ0oBN4HpiZ7vu8pG+QJEGAyyLi+WrEbWZmg5e62zJr2ebNm2v/IszMasDMmTPZsmXLzvIee+zBjTfeWJVz77333j36RHgkDjOzClm7di0zZ87kySefzDqUshlMTYhOYGZmFXLttdeyZcsWrrnmmqxDqUtOYGZmFbB27Vra29sBaG9vr5ta2GAaIssJzMysAq699toe5XqphbkJ0cysznXXvnor265zAjOzQWGwDBBbLmPHju2zXKvuueeeHuW77747o0icwMxskBgsA8SWyznn9Jz56dxzz80okvIaMWJEn+VqcgIzs8x1dHTQ1tZGRLBw4cK6qIXlD7G0aNGijCIpr1deeaXPcjU5gZlZ5lpaWujs7ASgs7OzLmphv/nNb3qUlyxZklEk5bXffvv1Wa4mJzAzy9zixYt7jHBeL7WVerRjx44+y9XkBGZmmWto6DnJ+r777ptRJOVz/PHH9yhPmTIlo0jKa9OmTX2Wq8kJzMwyt3Hjxh7lP/zhDxlFUj4zZsxASobuGzZsGDNmzMg4ovrjBGZmVgENDQ1MnToVgKlTp9bNnGCTJk3qUZ48eXJGkVRvRmYzs16NHj2aDRs27CyPGTMmw2jK5/TTT2fZsmWcfvrpWYdSlI985CMl77N06VKWLl3a5zY/+clPBhpSn1wDM7PMdXR09Cg//3x9TPnX2trK1q1baW1tzTqUulS1BCbpVEmPSloj6ct9bPeXkkLSMWn5EElbJN2fvr5XrZjNrDq6m9q6nXDCCRlFUj4dHR0sXLiwrp5tG2yqksAkDQeuA04DDgfOknR4ge32Ai4E7slb9VhEHJW+zq94wGZWVc3NzTs7PEiiubk544h2XUtLC11dXUD9PNs22FTrHthEYE1EPA4gaR5wBvBQ3nbfAK4EvlCluMxskJBEROxMZLWu0LNts2bNyjiqvhVzryr/Plml7m8Vo1pNiAcC63LK7emynSSNBw6KiF8W2H+cpJWSfi1paoH1ZlbDWlpaetTA6qG2MmrUqD7LtuuqVQMr9CdV7FwpDQO+A8wssN0G4OCIeE7S0cDPJb0zIl4sdKLVq1eXIVwzq6ZFixbtbG7r6uri17/+dc3fB8t/tm3jxo11+fup0tfU2NjY67pqJbB24KCc8lhgfU55L+BdQFv6V9gYYL6k6RGxHNgGEBErJD0GHAYsL3Sivi7WzAanSZMm9Rg+6rjjjqv5/8tvfOMbee6553qUa/2aCsnymqrVhLgMaJQ0TtJI4ExgfvfKiNgcEaMi4pCIOAS4G5geEcslvTntBIKkQ4FG4PEqxW1mNiC5yatQ2XZdVWpgEdEp6QLgDmA4cENErJJ0GbA8Iub3sfsJwGWSOoEu4PyIqI+HRMwMgHvvvbdH+Z577uGTn/xkRtEUZyAP/fa3T5YdImpR1UbiiIjbgdvzll3Sy7ZNOe9bgNq/o2tmvRo1ahTt7e09ymb98UgcZjWmo6OD2bNn19WDsYNphHOrHU5gZjWmpaWFRx55pC66mnebOnVqj270td4D0arDg/ma1ZCOjg7a2tqICNra2mhubq6LUc6bm5u56667gOSh31oYiaO/+1X597uuv/76uviuBhPXwMxqSEtLy87RHXbs2FE3tbDNmzf3Wa4HTl7l5wRmVkOWLFlCZ2cnkIyvt3jx4owjKo+rr766z7JZISUlMElvkvRRSV9MywdIGluZ0Mws35QpUxgxImn5HzFixGtGca9VuXOBAaxfv76XLc3+pOgEJulE4FHgr4H/my5uBL5bgbjMrIDcUduHDRtWE/eKzCqllBrY1cBHIuJUoDNddg/JSPNmVgUNDQ00NTUhiaamprq5r7LHHnv0WTYrpJReiIdExP+k77sH4t1e4jHMbBc1NzfT3t5eV7WvLVu29Fk2K6SUGthDkt6ft2wa8GAZ4zGzfjQ0NDB79uy6qX2ZDVQptafPA7+U9F/AHpKuB/6CZGJKMzOzqiqlBnYvcASwCrgBeAKYGBHLKhGY2a6qxyGXzOxPikpg6XQmrwDPRcSVEfGpiLgiItr729csK/U45FK9GjZsWJ9ls0KK+lcSEV3A74E3VTYcs/LIH3LJtbDBbceOHX2WzQop5c+cH5PcAztH0nslndz9qlRwZgNVr0Mu1avuh7N7K5sVUsq/kr9Lf87OWx7AoWWJxqxMCg25NGvWrIyjst64CdEGouh/JRExrpdXUclL0qmSHpW0RtKX+9juLyWFpGNyln0l3e/RAl35zV6jXodcqlcnnnhij3JTU1M2gVhNKameLmkE8B7gQKAdWBoRnX3vtbMTyHXA+9L9lkmaHxEP5W23F3AhyQgf3csOB84E3gkcALRKOiy9L2dWUHNzMwsXLgSS+aXq6aHfjo4O5syZw0UXXVQTz4LlTytSjDvvvJM777yz1/X9TWViQ0MpYyG+A3gYuIUkydwKPCLpz4rYfSKwJiIej4jtwDwKPz/2DeBKYGvOsjOAeRGxLSKeANbg4ausHw0NDYwZMwaA0aNH18Qv+mL98Ic/5OGHH+bGG2/MOhSzTJXS0PwvwL8CB0XE5IgYC3wvXd6fA4F1OeX2dNlOksanx/5lqfua5evo6GDjxo0AbNy4sW56IXZ0dHDPPUkDxdKlS+vmuswGopQmxKOA90V3167E1cDXithXBZbtPI6kYcB3gJml7ptv9erVRYRj9W7+/Pk7u2Lv2LGDH/zgB0yfPj3jqHbdrbfe2qN8zTXXcNZZZ2UUTXZq8f95LcZcjEpfV2NjY6/rSklg64ETgQU5y6amy/vTDhyUUx6bt99ewLuAtnSqiDHAfEnTi9i3h74u1oaOBx98kK6u5DZpV1cXDz74IJ///OczjmrXrVq16jXlwf5vvpj7Vfn3yerxHle1v6eB3HsciIsvvrjsxyz2+y+lCfGrJEllnqRvSZoHzE+X92cZ0ChpnKSRJJ0y5nevjIjNETEqIg6JiEOAu4HpEbE83e5MSbtJGkcyB9m9JcRtQ5B7IZrVv1K60c8HJgC/I6kx/Q44OiJ+UcS+ncAFwB0kHUF+GhGrJF2W1rL62ncV8FPgIeBXwKfcA9H6U68TP+633349yqNHj84oErPsFd2EKGk34ImIuDxn2esk7RYR2/rbPyJuB27PW3ZJL9s25ZW/CXyz2FjNuid+bG1trauJHw899FCeeeaZHmWzoaqUe2B3AV8kad7rdjRwBdBUxpjMyqIeJ3584IEHepTvv//+jCKxWnPje6/POoRezfyfvx3QfqXcA3s3OQ8Yp+4FjhzQmc0qrB4nfpwyZcrOplFJvrdnQ1opCWwzkN/gPppkmhUzq4Jp06btHKQ4Ipg2bVrGEZllp5QE1gLcIuldkl4v6d3Av5F0sDAbdOpxQsvW1tYeNbDW1taMIzLLTin3wL4G/BNJs+FuwDaSmZmL6UZvVnW5E1rWwkj0pT63ExH9jhkI9flMlRmU1o1+a0R8CtiT5EHjPSPigojY2s+uZlXnCS3N6l8p3egPB56LiI2StgCzJXUB346IVysWodkAFJrQshZqYWaVMtCefoNZKffAbgG6u3N9GzgBmAwM3r6ZNmQVmtDSzOpLKffADomIR5XcQf4/JPNzbQGeqEhkZrtgypQpLFy4kM7OzpoZSspjBmanGuMGVuocQ/nfQCk1sG3phJMTgXURsYmkI8fuFYnMbBfU61BSZvYnpdTAbiEZiX4v4J/TZRNwDcwGoYaGBiZPnsyiRYuYPHlyXT3MbDYQ9TgSR9EJLCI+K+kU4I8RsTBdvAP47IDObGZmtgtKqYEREXfmlZfnliW9GBFvLEdgZruio6ODpUuXAvDb3/6WGTNmuBZmRTvrEzdnHUKvbv3Xj2YdwqBRyj2wYhSaPdms6lpaWnr0Qmxpack4IjMrt3InsCjz8cwGZPHixT3GDFy0aFHGEZlZuZU7gZkNCqNGjeqzbGa1r2oJTNKpkh6VtEbSlwusP1/Sg5Lul7QkHfkDSYdI2pIuv1/S96oVs9WujRs39lk2s9pXlXtgkoYD1wGnAYcDZ3UnqBy3RMS7I+Io4Ergqpx1j0XEUenr/DLHPOTV46jtO3bs6LNsZrWvlLEQpwO3R0RnH5ud1svyicCaiHg8PdY84Azgoe4NIuLFnO33xPfTqqbWRm0vRldXV59ls6GmHsdCLKUb/TeAH0j6CXBzROTPzkxELOll3wOBdTnldmBS/kaSPgV8DhgJnJyzapyklcCLwMUR0evAdqtXr+7vOizHSy+9xMKFC4kIFixYwPjx49lrr72yDqsi6vHfRj1eE9TvdVVCPX5WudfU2NjY63alPMh8pKQjgbOBFkmvADcDP4qItf3sXqhp8TU1rIi4DrhO0gzgYuAcYANwcEQ8J+lo4OeS3plXY9upr4u115o7d26P8sqVK+uiFjZ+/HhWrly5s3z00UfX5b+NerwmqN/rqoR6/KyKvaZSH2R+AHhA0heB95JMcHmppN+QjEp/a0QUutnQDhyUUx4LrO/jVPOA76bn3EYy5iIRsULSY8BhwPLed7diFRq1fbAnsIEMirpixYp+9xvKg6JaT35YuDaUlMAAJL2VpBZ2NslQUpcATwEXAM3AhwrstgxolDQOeBo4E5iRd9zGiOiuN54OrE6Xvxl4PiK6JB0KNAKPlxq3FVaLo7abWf8q9QfZYJoRoZROHJ8CPgq8Dfgp8NGIuDtnfQvwTKF9I6JT0gXAHcBw4IaIWCXpMmB5RMwHLpA0Dfgj0EHSfAjJvGOXSeoEuoDzI+L5Eq/TetHc3ExbWxvgUdutf9WYdqRS53ENu/6UUgM7jaTJ8BcRsT1/ZUS8KqlQ7at7/e3A7XnLLsl5/5le9msBPA5QhTQ0NNDU1ERraytNTU0eL9DMakYpCawtIn6Wv1DS5yLiKnjtYL9WG5qbm2lvb6+Z2pcnfrRK82C+taGUB5kv6WX5xeUIxLLT0NDA7NmzXfsys5rSbw1MUvfzWCMknUTPLvGHAi9VIjAzG/x2+/jgnQ5w29zvZB2CVVgxNbAfpK/dgBtyynOB84BPVyw6q4q1a9cyc+ZMnnzyyaxDMTMrWr8JLCLGRcQ44Mfd79PXoRHxnrQHodWwq6++mi1btnD11VdnHYqZWdGKvgcWER+rZCCWjbVr17JhwwYA1q9f71qYmdWMPhOYpIdz3q+T9FShV+XDtErJr3W5FmZmtaK/Thx/k/P+7EoGYtnorn11W7++rxG+zIYGd1WvDX0msNzR5SPi15UPx8zMrDilDCX1OWBBRNwv6TiS4aQ6gb+OiKWVCtAqa9KkSdxzz59mxpk8eXKG0dSPoTLkkruqW5ZKGYnjsyTd5wH+H8mMyS8BV1Ngbi8bHEr9Bbd06VKWLu377xGPamFmg0EpCWzviNgsaS/gSGBaOkL8P1UotkGpo6ODOXPmcNFFF3nkCrM6UYk/yjycWeWVMpTUOknvIZkKZVGavN5IMkL8kNHS0sIjjzxCS4vHFzYzy1IpNbC/B24DtpPM+wXwAeDecgc1WHV0dNDW1kZE0NbWRnNzs2thVpThnzg16xB61fWvvxrwvh5KyrJUdAJLp0M5IG/xz9LXkNDS0kJEALBjxw5aWloG/ezF/TVbuJnDzGpVKU2ISNpb0kRJJ6eD/E5NX8Xse6qkRyWtkfTlAuvPl/SgpPslLZF0eM66r6T7PSrp/aXEXE5Lliyhs7MTgM7OThYvXpxVKGZmQ14p3ehnAtcBLwOv5qwKklHp+9p3eLrv+4B2YJmk+RHxUM5mt0TE99Ltp5P0cjw1TWRnAu8kqQG2SjosIsp6720g3ZG3bNnS736u0ZiZVUYp98C+CfxlRPz3AM4zEVgTEY8DSJoHnAHsTGAR8WLO9nuSJEbS7eZFxDbgCUlr0uP52TMzsyGslAQ2AhjojMsHAutyyu0UeHZM0qeAzwEjge55yA4E7s7b98ABxmFmZnWilAT2LeBiSd+IiB0lnkcFlsVrFkRcB1wnaQbJTM/nFLtvt9WrV5cYWmUNtnj6U2vxFqter6sS6vWzyvq6sj5/pVT6uhobG3tdV+pIHGOAL0p6LndFRBzcz77twEE55bFAX6PGzgO+O5B9+7rYvhRzr2oo9Ngb6Oc32NXrdVVCKZ9VLXVVz/rfQNbnr5Qsr6uUBLYro9EvAxoljQOeJumUMSN3A0mNEdGdyk8Hut/PB26RdBVJJ45GhtCzZ0NBNcYNrNQ5iv0jZleetTKzwkp5DmzAo9FHRKekC4A7gOHADRGxStJlwPJ0VucLJE0D/gh0kDQfkm73U5IOH53Ap8rdA9HMzGpPKd3odwMuAc4C3hQRe0s6BTgsIv65v/3TB6Fvz1t2Sc77z/Sx7zdJekGamZkBpTUhfoek999fA91d6Vely/tNYGZW+yp133co3F+28islgf0f4G0R8YqkHQAR8bQkd2m3sjppZq+dTDO38MZCnWL7V69jIVp9Gci94mL2qdQfJKUMJbWdvIQn6c3Ac4U3NzMzq5xSEtjPgJvSnoRI2p+k6XBeJQIzMzPrSykJ7KvAWuBBYB+Sbu7rgUvLH5aZmVnfSulGvx24CLgobTrcFN1zi5iZWc0r5l7Vueeey6uv/mk899e//vX88Ic/rGRYvSqlEweSGoG/InmgeL2kn+Y8fGxmvXBHCasXucmrULmaim5CTMcnXAkcAbwCvBu4L11uZmZWVaXUwC4H/jwiFnUvkDQVuBm4pdyB2WsNhSGXYOBd1c2s8iSRe/dIyu7/aymdOPbitXNw3U0yd5eZmQ0Bkyb1nAnruOOOyyiS0mpgVwH/IOn/RsRWSXuQ9EC8qjKhmdUuj1hh9WrkyJF9lquplBrYJ0l6Ib4oaSOwmWSKlb+T9FT3qxJBmpnZ4LBs2bIe5XvvzW5ykGpNp2IV8LXmkjqRVtU3WzoHvG89DiVlVi+OPfZYFi3a2RWCiRMnZhZLVaZTMTMzK7dSutG/TtKlkh6XtDX9eamk7BpAzcysqgZTE2Ip98CuBKYB5wNHpj9PBr5VgWfW6UQAABI5SURBVLjMzGwQOuKII3qUjzrqqIwiKe0e2IeBIyOie/T5RyXdBzxA0pmjT5JOBeaQzMg8NyKuyFv/OeDjJLMuPwucFxFPpuu6SMZgBHgqIqaXELfVGN9nMhu8nnqqZ1+9J598MqNISktgvf1W6fe3jaThwHXA+4B2YJmk+RHxUM5mK4FjIuJVSX9HUuPr7jO8JSKKTvPVeOC3UucppVv0rnSUMDMbiA0bNvQor1+/PqNISktgPwP+U9KlwFPAW4CLgZ8Wse9EYE1EPA4gaR5wBrAzgUXEwpzt78a9Hs1sECv1D9gsJ34spz333JNXXnmlRzkrpdwD+yLQSlKTWgFcCywEvlTEvgcC63LK7emy3swC/junvLuk5ZLulvTBEmI2M7My6uzs7LNcTUXVwNImwO8Dn4iISwZwnkLNjAUf9pF0NnAMcGLO4oMjYr2kQ4EFkh6MiMcK7b96dW0Pjl/r8femr+u6/PLLy36+iy++uOLngOy/r6zPXyn1el39qYXrPvLII3v0PDzyyCMrGndjY2Ov64pKYBHRJekUYMcAY2gHDsopjyWZDLMHSdOArwEnRsS2nPOvT38+LqkNGA8UTGCFLvbmMwZva+RHf/GjHuW+vqx8tfQgcynXVQlZn79SfF31pRau+7zzzuO+++6js7OT173udcyaNYt99tknk1hK+Q34HeBSSbPTyS1LsQxolDQOeBo4E+gxDYuk8cD1wKkR8UzO8gbg1YjYJmkUcDxJBw8zs8z0d7+q0D2vWrjH1Z+GhgZOOukkWltbOemkkzJLXlBaAvs0MAb4nKRnSZoABUREHNzXjhHRKekC4A6SbvQ3RMQqSZcByyNiPvCPwBuAn6XD83d3l/8z4HpJO0ju2V2R13vRzMyqqLm5mfb2dpqbmzONo2pjIUbE7cDtecsuyXk/rZf9fksyeaaZmQ0CDQ0NzJ49O+swSkpgS0m6zZ8FHEByD2se8M0KxGVF8HNgZjaUlZLAvgu8HbgQeJLkObCvkHSHP6/8oZVPfkcJMzOrfaUksA8Cb42IF9LyQ5LuAdYwyBOYmZnVn1IS2B+A1wMv5CzbA9hQeHMzs4ENudbfPvXQm892XSkJ7GbgV5Ku5U/PdX0K+DdJJ3dvFBELyhuidavEf1pPUW9mtaqUBPa36c+v5i0/P31B0rX+0F0Nqtxq6UFmMzMrTikzMo+rZCBmZvVk0qRJ3HPPPTvLkydPzjCa+jR4xyIys7pQTLN0PY5ace655/ZIYDNnzswumDpVymj0ZmZWgnRUIYYN86/aSvCnamZWAS0tLTsTlyRaWloyjqj+DIkmRHeUMLNqW7JkCV1dXQB0dXWxePFiZs2alXFU9WVIJDCzWlCJ56Wg9u8l1aopU6awcOFCOjs7GTFiBFOnTs06pLrjJkQzswpobm7ucQ8s65Hb65ETmJllrvsXfW/lWtTQ0EBTUxOSaGpqynTerHpVl02IlWoy8agVZpUxbNiwnfeLusv1YLDMm1WvqpbAJJ0KzCGZ0HJuRFyRt/5zwMeBTuBZ4LyIeDJddw7JVC4Al0fETdWK26xahurzUgDHHHNMj2emJk6cmGE05TNY5s2qV1X5M0fScOA64DTgcOAsSYfnbbYSOCYijgBuA65M990X+DowCZgIfF1SQzXiNhts9ttvvx7l0aNHZxRJee222249yiNHjswoEqsl1aqnTwTWRMTjEbGdZCLMM3I3iIiFEfFqWrwbGJu+fz9wV0Q8HxEdwF3AqVWK22xQGTeu54huhx466IYeHZBly5b1KN97770ZRWK1pFpNiAcC63LK7SQ1qt7MAv67j30PLGt0dazUrtm10C17KHc3f+CBB3qUV65cmVEk5TVlyhQWLFhAV1cXw4cPd5dzK0q1ElihLkVRcEPpbOAY4MRS9wVYvXp1ycENVDXPNZjU63XXwnXtueeebN26tUe5FuLuz4QJE1i4cCGQ9EAcP358XVyX7brGxsZe11UrgXXPH9ZtLLA+fyNJ04CvASdGxLacfZvy9m3r7UR9XWy5VfNcg0m9XnctXNcLL7zwmnItxF2Mk046idbWVk4++WQmTJiQdThWA6p1D2wZ0ChpnKSRwJnA/NwNJI0HrgemR8QzOavuAE6R1JB23jglXWY25EREn+Va1tzczDve8Q53ObeiVaUGFhGdki4gSTzDgRsiYpWky4DlETEf+EfgDcDP0ocYn4qI6RHxvKRvkCRBgMsi4vlqxF0PauG+TqmKuaYZM2b0eK5o+PDh3HLLLZUMqypGjx7Nhg0bdpbHjBmTYTTl5S7nVqqqPQcWEbcDt+ctuyTn/bQ+9r0BuKFy0Vm9Of7441m0aNHO8pQpUzKMpnw6Ojp6lJ9/3n/L2dBVH4+7m+WZMWPGzuGIJDFjxoyMIyqPqVOn9riuE044IeOIzLLjBGZ1qaGhYWfz2pgxY+pmHLrm5mZGjEgaTkaMGOH7RTakOYFZXero6GDTpk0AbNq06TW992pV7gCxJ510Ut0kZrOBcAKzutTS0rKzh15E1NVsuO6tZ5ZwArO6tGTJEjo7OwHo7Oxk8eLFGUdUPt299Vz7sqHOCcxYu3YtM2fO5Mknn8w6lLKZMmVKj84OHprIrP44gRnXXnstW7Zs4Zprrsk6lLKZNm1ajybEadN6fUrDzGqUE9gQt3btWtrb2wFob2+vm1pYa2trjxpYa2trxhGZWbk5gQ1x1157bY9yvdTClixZ0qMGVk/3wMws4QQ2xHXXvnor16opU6b0eF7K98DM6o8T2BA3duzYPsu1qrm5eWcT4rBhw9zl3KwOOYENcZ/+9Kd7lC+88MKMIimv3Ad+m5qa3OXcrA5VbTDfwW6ozvJ7yCGHsP/++7Nhwwb2339/3vKWt2QdUtk0NzfT3t7u2pdZnXINzHYmrUMOOSTbQMrMD/ya1TcnsCGuo6ODFStWALBixYq6GTPQzOqfE9gQlztm4I4dO+pqzEAzq2+q1pTkkk4F5pDMyDw3Iq7IW38CcDVwBHBmRNyWs64LeDAtPhUR03P33bx5c/3Mq15lM2fOZMuWLTvLe+yxBzfeeGN2AZmZ9WLvvfdWbrkqNTBJw4HrgNOAw4GzJB2et9lTwEyg0LzvWyLiqPQ1vcB6GyA/L2VmtapaTYgTgTUR8XhEbAfmAWfkbhARayPif4EdVYrJ8PNSZla7qpXADgTW5ZTb02XF2l3Sckl3S/pgeUMb2vy8lJnVqmo9B6YCy0q5b3VwRKyXdCiwQNKDEfFYoQ1Xr149oACHsgkTJrB69WrGjx/vz8/MBpXGxsZe11UrgbUDB+WUxwLri905ItanPx+X1AaMBwomsL4u1no3YcKErEMwMytJtZoQlwGNksZJGgmcCcwvZkdJDZJ2S9+PAo4HHqpYpGZmVhOqksAiohO4ALgDeBj4aUSsknSZpOkAko6V1A58GLhe0qp09z8Dlkt6AFgIXBERTmBmZkNc1Z4Dq6RqPgfW0dHBnDlzuOiii9zhwcysijJ5DqyetLS08Mgjj3jECjOzjDmBlaCjo4O2tjYigra2No8baGaWISewEnjcQDOzwcMJrARLliyhs7MTgM7OThYvXpxxRGZmQ5cTWAk8bqCZ2eDhBFYCjxtoZjZ4OIGVwOMGmpkNHtUaSqpuNDc3097e7tqXmVnG/CCzmZnVBD/IbGZmdcEJzMzMapITmJmZ1SQnMDMzq0lOYGZmVpOcwMzMrCbVRTd6MzMbelwDMzOzmuQEVgJJN0h6RtLvso6lXCQdJGmhpIclrZL0maxjKgdJu0u6V9ID6XVdmnVM5SJpuKSVkn6ZdSzlImmtpAcl3S9pedbxlIukfSTdJumR9P/Y5Kxj2lWS3p5+T92vFyVdlEksbkIsnqQTgJeBf4uId2UdTzlI2h/YPyLuk7QXsAL4YEQ8lHFou0TJqMt7RsTLkl4HLAE+ExF3ZxzaLpP0OeAY4I0R8YGs4ykHSWuBYyJiU9axlJOkm4DFETFX0kjg9RFRNzPhShoOPA1Miognq31+18BKEBGLgOezjqOcImJDRNyXvn8JeBg4MNuodl0kXk6Lr0tfNf/XmqSxwOnA3Kxjsb5JeiNwAvADgIjYXk/JK/Ve4LEskhc4gVkOSYcA44F7so2kPNKmtvuBZ4C7IqIerutq4IvAjqwDKbMA7pS0QtInsg6mTA4FngV+mDb5zpW0Z9ZBldmZwK1ZndwJzACQ9AagBbgoIl7MOp5yiIiuiDgKGAtMlFTTzb6SPgA8ExErso6lAo6PiAnAacCn0ub6WjcCmAB8NyLGA68AX842pPJJm0SnAz/LKgYnMCO9R9QC/Dgi/j3reMotbbZpA07NOJRddTwwPb1fNA84WdKPsg2pPCJiffrzGeA/gInZRlQW7UB7Ts3/NpKEVi9OA+6LiI1ZBeAENsSlnR1+ADwcEVdlHU+5SHqzpH3S93sA04BHso1q10TEVyJibEQcQtJ0syAizs44rF0mac+0AxFpE9spQM339I2IPwDrJL09XfReoKY7R+U5iwybD8ETWpZE0q1AEzBKUjvw9Yj4QbZR7bLjgY8CD6b3iwC+GhG3ZxhTOewP3JT2khoG/DQi6qbbeZ0ZDfxH8rcUI4BbIuJX2YZUNp8Gfpw2tz0OnJtxPGUh6fXA+4C/zTQOd6M3M7Na5CZEMzOrSU5gZmZWk5zAzMysJjmBmZlZTXICMzOzmuQEZmZmNckJzMwqSlJIelvWcVj9cQIzM7Oa5ARmQ56kL0l6WtJLkh6V9F5JN0q6PGebpnT0le7yQZL+XdKzkp6T9M856/4mnbzwJUkPSZqQLj9AUku6zxOSLszZZ6Kk5enkgBslXZUu313Sj9JzvCBpmaTR/VzP3pJ+IGlDel2XpyPzj0wnIPx0ut1wSb+RdElanp1OvviTNPb7JB2Zc9y+4h8u6auSHkv3XZF+RovSTR6Q9LKkjwzwazJ7DQ8lZUNaOk7dBcCxEbE+nVJmeD/7DAd+CSwgGYari2SCSSR9GJgNfBBYDrwV+KOkYcB/Ar8gGUNuLNAq6dGIuAOYA8yJiJvTmQG6R84/B9gbOAjYBhwFbOnnsm4CNgJvA/ZMY10XEddLOhtYLKkV+FB6rd/M2feMNL6zgc8AP5d0WHqNfcX/uXT5nwO/B44AXo2IEyQFcGRErOknbrOSOIHZUNcF7AYcLunZiFgLkI7L15uJwAHA30dEZ7psSfrz48CVEbEsLa9JjzcJeHNEXJYuf1zS90kG5b0D+CPwNkmj0lmJu2eO/iPwJuBtEfG/JDNm9yqtnZ0G7BMRW4BXJH0H+ARwfUT8Lq1Z/gfJGIQTI6Ir5xArIuK29FhXAZ8HjgO29xP/x4EvRsSj6foH+orTrBzchGhDWloruIik1vSMpHmSDuhnt4OAJ3OSV/66xwosfwtwQNoM+IKkF4CvkiQRgFnAYcAjaTPhB9LlN5MkiHmS1ku6Mp3+pjdvIZl9ekPOea4H9svZ5ibgEOD2iFidt/+67jcRsYNkSpADioi/t+s2qxjXwGzIi4hbgFuUTAF/PfAt4EXg9Tmbjcl5vw44WNKIAklsHUmzYb51wBMR0dhLDKuBs9Kmxg8Bt0l6U0S8AlwKXJo2b94OPEo6TX0v59kGjOolwQL8C0mz4vslTYmIJTnrDup+k8YyFlgPdPYVP3+67pqfBsVqh2tgNqRJerukkyXtBmwlub/UBdwP/LmkfSWNIamldbsX2ABcoWQuq90lHZ+umwt8QdLRSrxN0lvSfV5MO4zskXZ6eJekY9M4zpb05rTW80J6rC5JJ0l6d3rf7UWSJsXcJr8eImIDcCfwT5LeKGmYpLdKOjE9z0eBo4GZwIUkU868IecQR0v6kKQR6TVvI2nO7DP+9Lq/Iakxve4jJL0pXbcROLSoL8SsBE5gNtTtBlwBbAL+QNLU9lWSprsHgLUkCeEn3Tuk94z+gqSTxFMkzWwfSdf9jKRTxC3AS8DPgX1z9jkKeCI931ySDhqQzBa9StLLJB06zoyIrSQ1v9tIktfDwK+B/mZh/hgwkmTyxI50//0lHQxcDXwsIl5Oa57Lge/k7PuL9Fo6SDqofCgi/lhE/FcBP00/qxdJaoh7pOtmkyTKFyT9VT+xmxXN84GZGZB0oyfpLFLzszzb0OAamJmZ1SR34jCrQWlTYyGnRcTiqgZjlhE3IZqZWU1yE6KZmdUkJzAzM6tJTmBmZlaTnMDMzKwmOYGZmVlNcgIzM7Oa9P8BexzvcwnP+8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.boxplot(x=\"success_expect\", y=\"propensity_score\", data=data_ps)\n",
    "plt.title(\"Confounding Evidence\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have to check that there is overlap between the treated and untreated population. To do so, we can see the empirical distribution of the propensity score on the untreated and on the treated. Looking at the image below, we can see that no one has a propensity score of zero and that even in lower regions of the propensity score we can find both treated and untreated individuals. This is what we called a nicely balanced treated and untreated population. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEXCAYAAADiEjDuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZxU1Z338c+XZlPZXYgKER/TGrc4RiOIxBidiaIjJBrnZcaMwGAyI5o4URMwPplkzPiASwZJJiYuiXsUshiXcQkuJDiMW6LBfRrUKGK6RbCBKLL9nj/u6bZsq7ur6equW833/XrVq27de+69v7pN1Y9z7qlzFBGYmZlVWq9KB2BmZgZOSGZmlhNOSGZmlgtOSGZmlgtOSGZmlgtOSGZmlgtOSNbjSRolKSSNa6fcZEkbO3jsBZKu7lyEnZPe2xe7+BwlXUOzznBCsoqTdG36sgtJGyX9SdKPJW1fplO8CuwMPJLONyKd64gW5eYCu3bw2CcAZze9kHSfpGu3PNT3k7SvpBskLZP0bro2t0r6dLnOYZYXTkiWFwvJksYo4KvAicD15ThwRGyKiD9HxIZ2yr0TEfUdPPbKiFjduQiLk3Q08DiwC/AlYB/geOBh4IquOKdZJTkhWV6sT0ljWUTcBlwGHCNpG2XOlfSipPWSlkr6l8KdJU2U9ISktyW9JelRSQembS2bm15Nzw+m9S+ncs1NdpIGpWP9fYvz7Cxpk6Rj0uvmJrtUMzoKmFRQ4ztC0m8lXdniOErv4zvFLoakbckS8oKIOCoi7o6IpRGxOCIuAka32GVQqkmtkfSqpG+0OF5vSd+R9JKkdZKekfRPLcoMkHRZ2v9dSS9L+max+FL56ZJWFalpmm0RJyTLq3fI/n32BqYB3wVmAfsClwCzJE0FkPQh4OfAzWn7oWQJrbX7QR9PzyeS1co+0bJAqvXcBkxqsekUoB6YX+S4Z5HV9Oal4+4MLAJ+DHxB0oCCskeS1QZ/2kqMnwF2Ai4stjEiVrVY9W3gd8BfkV2fi1o0611N1rz4T8DewAWpTNM1FHAnMAH4SipzKvBGy3NL6iXpB+n9fioiFrTyHsw6pHelAzBrSdI+wBnAIxGxRtIM4AcR0VTLqJO0F3A+8BOyL/4+wLyIeDmVea6NUzR9ya6MiD+3Ue464E5Ju0TE8rTuH4CbImJTy8IR0ShpPfBO4XEl/RL4PnAyWWIAOA24JyJeaeXce6bnZ9uIr9DciLgqLX9f0jSypPagpN3Jkss+EfF8KvNSuoZfIbuGRwKfAj4REY+nMi+SJblC/cjute0PHBoRfyoxPrN2uYZkeXGEpLWS3gGeJvsy/HtJg4ARfPCL8bfAqNS0tRi4F3g63fA/S9LIMsQ0H2ggqxUh6QDgY3Tw3lZErAeuJbsPROqs8TngqjZ2UwdjfbLF69eA4Wn54HS8x9M1XitpLfBNoDaVOQhYVZCMWnMNWTI6zMnIys0JyfLiEbLmpr2BbSLibyLixYLtLYelb/7CTrWV8WT/y3+MrCnufyX9bWcCSse9iax2QXp+IiKe2oLDXQF8QtLHyGpZK8mayFrzQnrep8Tjr2/xOnjv8930PJbsGjc99iNLsIX7tOe/gN2BY0qMy6xkTkiWF+9ExJKIeDki3m1ame7lLCNrTip0OPBSRLydykVEPBoR/y8iDierQU1p5VxNX941JcR1HbCfpIOBL6TXbVlf7LgRsQR4gKyWdBpwTUS09Zun35DVzs4vtlHS0PZDb/b79PzhdI0LH0sLygxL77MtN5HdV/uppJb318w6xfeQrBrMBL4nqQ5YQFYTOp3sPhOSxpL1bvsN8DpZM9THyO6NFLMCWAt8RtIzwLtFOgkAEBFPS3qCrHltR7KOE215Cfi0pD2ARqCxoLv5FcCNZPe7jm/rIBHxtqTJwK2S7gO+R1Zr2g44miyx7dVOLE3HWiLpp8BVqffd/6TjHATsmHrtPUDWIWOupLPJmkF3AfaOiKtbHO8WSRuAmyT1Lbh3ZdYpriFZNfgR8K9k9zyeBaYDMyKiKeE0kvWsuw2oI+u5dhNZz7wPiIjNZMns78i6gD/RzvmvI2viuiciGtop+z2yhPdHss4ThxVs+3WKdX5EvNTOcYiIu8l6ANaTJdfnyZrMxpLuR3XAl4HZZDWuZ4H7yWo6L6ZzBXAccBdZr8AXyJLnDq3E9kuy69fUgcKs0+QZY826h6RhZJ0Nvpi+0M2sgJvszLqYpD5kPd6+BSwnqymZWQtOSGZd7zDgQbL7S6cW+w2TmbnJzszMcsKdGszMLBcq0mTX2NjoapmZWQ82ePDgjo424hqSmZnlgxOSmZnlwlabkOrq6iodQkmqJU6onlirJU6onlirJU6onlirJU4oX6xbbUIyM7N8cUIyM7Nc8A9jzaxqRQRr165l8+bNJe/Tv39/GhsbuzCq8sh7nL169WLAgAFkkw2XhxOSmVWttWvX0q9fP/r27VvyPv369aN///5dGFV55D3O9evXs3btWgYOHFi2Y7rJzsyq1ubNmzuUjKx8+vbt26GaaSmckMzMLBeckMzMOmHIkCGcf/57E/v+4Ac/YObMmZ065sqVKznqqKMYN24ce+65J3vvvTfjxo1j3LhxrF/fcrb60i1cuJBHHnmkw/vtv//+vPnmm1t83lL5HpJZB1z7wl9a3TZ5r+26MRIrpq2/T5MNGzbQp09pA66X8jft168fd9xxB2effTbbb799Scdtz7Bhw7j//vvp378/M2fOZMCAAXzlK195X5mNGzfSu3fHvsIfeughBgwYwOjRo8sSZ7k5IZmZdULv3r2ZPHkyl19+Od/61rfet+2VV17hzDPPZMWKFeywww788Ic/ZOTIkZx++ukMHDiQJ598kvr6ei644AImTpzY7rlOP/10hg4dyuLFiznggAM47bTTOPfcc1mxYgXbbrstc+bMYc899+Tuu+/m0ksvZf369QwbNoyrrrqKd955h2uuuYaamhrmzp3LxRdfzJ577snXvvY1li1bBsDMmTMZM2YMK1euZOrUqaxYsYKDDjqI7poVwk12ZmaddNpppzFv3rwPdNP++te/zsknn8yiRYs46aSTmD59evO2+vp67rnnHubOnct3vvOdks+1dOlSbrvtNi688ELOOussLr74Yn7729/y3e9+l3POOQeAQw89lPvuu4+FCxdy4oknMmfOHHbbbTemTJnCtGnTeOihhxg7diwzZsxg2rRpPPjgg1x//fV89atfBWDWrFmMGTOGhQsXMn78+OaE1dVcQzIz66RBgwZx8sknc8UVV7DNNts0r3/ssce48cYbATj55JP59re/3bztuOOOo1evXnz0ox/ljTfeKPlcEydOpKamhrVr1/Loo48yadKk5m1N95dee+01pkyZQn19PevXr2e33XYreqwFCxbw/PPPN79es2YNa9asYdGiRc1xH3300QwZMqTk+DrDCcnMrAymTZvG4YcfzimnnNJqmcIfkfbr1695uSNNYtttl93X2rx5M4MHD+ahhx76QJlvfOMbnHHGGRx77LEsXLiQWbNmFT3W5s2bmT9//vuSaCW5yc7MrAyGDh3K5z73ueaaBcAhhxzCL3/5SwDmzZvHmDFjyna+QYMGsdtuu/HrX/8ayJLaU089BcDq1avZZZddALj55pub9xkwYABr1qxpfn3kkUdy1VVXNb9evHgxAGPHjuXnP/85APPnz+ett94qW9xtcUIyMyuTM888833doy+66CJuuukmxo4dy9y5c1utqWypK6+8khtuuIHDDjuMMWPGcNdddwEwY8YMJk2axPjx49/X82/8+PHceeedjBs3jkWLFnHRRRfxxBNPMHbsWEaPHs0111zTvP+iRYs4/PDDeeCBBxgxYkRZ426Nuqv3RKE8zBhbV1dHbW1tpcNoV7XECdUTa2fi7O5u31vDNe2MxsZGBg8e3KF91q1bl+sheZpUQ5xN17/Y398zxpqZWdVyQjIzs1xwQjIzs1xwQjIzs1xwQjIzs1xwQjIzs1zwSA1mZlto5cqVTJgwAYCGhgZqamqaf/fzwAMPbNHkgQsXLqRv374ccMABHdpv//33Z8GCBWUbcbwSnJBsq+RpJHqm3g/e0W6Zfhs20LtPn5KOt/HTx7e5fdiwYc1D9xSbJqIzU0R0NCH1BCVdKUkvA2uATcDGiDhY0jBgLjAKeBn4u4hYpWywpjnAscDbwOSI+EP5Qzczy59yTRFxyy23cMkll+Ruioiu1JHU/emIWFHwegZwf0TMkjQjvZ4OjAdq02M08KP0bGa2VWiaIqKmpoYJEyYwe/Zs9thjDx5//HHOOecc7rjjjuYpIiRx/fXXM2fOHC688EKmTJnCgAED+NKXvkT//v057bTTmDZtGoceeiivvvoqJ554Io8++mjzFBHTp0/n3nvv5dprr6302+60zjTZTQSOSMvXAQvIEtJE4PrI0vXDkoZI2jkiXu9MoGZm1aInTxHRlUpNSAH8RlIAV0TElcDwpiQTEa9L2imV3RV4tWDfZWmdE5JtlXy/auvTk6eI6EqlJqTDImJ5SjrzJT3fRtliA+q12rhZV1dXYgjlV8lzd0S1xAnVE2t9Q32r2+p6bWpjv5oO77cl+7yvTJVc00rE2b9///fNK9Rvw4aS9ttQYrl3160rOZaNGzeyYcMGNm3axPr161m3bh19+/Zl5MiRzJs3jwkTJhARPPvss+y77740Njay/fbbs27dOm688UY2b97cPKDqqlWrgGyA1U996lNcfvnlnHHGGQA8/fTT7LfffowePZqbb76Zr33ta9x///289dZbrFu3jnUdiLmzVq9eTUNDA1CeAXZLSkgRsTw9N0i6FTgEqG9qipO0M9CQii8DRhbsPgJY3tqxKzWSsUdRLr9qibWuro7hOw1vdXttbeu1luGbW6/ttLbfluzTpJquaaVG+y4cEbuU3nMbNmygT4m97NSB0bZ79+5Nnz59qKmpoW/fvs1xXX311ZxzzjnMmTOHjRs3csIJJ3DQQQdx3nnn8eUvf5lddtmFgw8+mGXLltG/f3+OP/54Tj31VO69914uueQSLr30Us4991yOPPJINm3axNixY5k9ezbnn38+U6dO5TOf+QyHHXYYI0aMoH///t06QvigQYMYOXJk2f7+7SYkSdsBvSJiTVr+DHABcDswCZiVnm9Lu9wOnCnpFrLODI2+f2Rm3aG9btqQ1Xo6kmhKdd555xVdP2rUqOZJ+godd9xxHHfccR9Y/5GPfIRFixa9b/qJpnmKCg0bNoxbb721+fXMmTO3NPTcKKWGNBy4NU292xv4WUTcI+kxYJ6kqcArwEmp/F1kXb6XkHX7nlL2qM3MrMdpNyFFxIvAB36hFRFvAkcVWR/AGWWJzszMthoey87MzHLBCcnMzHLBCcnMqlavXr2af2hq3Wv9+vX06lXeFOLBVc2sag0YMIC1a9fyzjvvlLzP6tWrGTRoUBdGVR55j7NXr14MGDCgrMd0QjKzqiWJgQMHdmifhoYGRo4c2X7BCquWOMvJTXZmZpYLriFZj9XaGHL1DTUM36noJjOrINeQzMwsF5yQzMwsF9xkZ1YmbU0zYWbtc0Iya8GJxawy3GRnZma54IRkZma54IRkZma54IRkZma54IRkZma54F52ZhXUVo++yXtt142RmFWea0hmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLTkhmZpYLJSckSTWSnpB0Z3q9u6RHJNVJmiupb1rfL71ekraP6prQzcysJ+lIDeks4LmC1xcBsyOiFlgFTE3rpwKrIuIjwOxUzszMrE0lJSRJI4DjgKvTawFHAr9IRa4DPpuWJ6bXpO1HpfJmZmatUkS0X0j6BTATGAicC0wGHk61ICSNBO6OiP0kPQ0cExHL0ralwOiIWNF0vMbGxuaT1tXVle/dmBX41Z9rKh1Cp5zwoU2VDsGsQ2pra5uXBw8e3OGKSLsT9En6W6AhIn4v6Yim1UWKRgnbPqDwDXSnurq6ip27I6olTshfrMM3F5/8rr6hnuE7De/maDruvze3HmveJu/L29++LdUSa7XECeWLtZQZYw8DJkg6FugPDAIuA4ZI6h0RG4ERwPJUfhkwElgmqTcwGFjZ6UjNzKxHa/ceUkScFxEjImIUcDLwQEScAjwIfD4VmwTclpZvT69J2x+IUtoFzcxsq9aZ3yFNB86WtATYHvhJWv8TYPu0/mxgRudCNDOzrUEpTXbNImIBsCAtvwgcUqTMOuCkMsRmZmZbEY/UYGZmueCEZGZmueCEZGZmueCEZGZmueCEZGZmueCEZGZmudChbt9meXPtC8WHBzKz6uMakpmZ5YITkpmZ5YITkpmZ5YITkpmZ5YITkpmZ5YITkpmZ5YITkpmZ5YITkpmZ5YJ/GGtVwT+ANev5XEMyM7NccEIyM7NccEIyM7NccEIyM7NccEIyM7NccEIyM7NccEIyM7NccEIyM7NccEIyM7Nc8EgNZlWorZErJu+1XTdGYlY+riGZmVkuOCGZmVkutJuQJPWX9KikP0p6RtK/pfW7S3pEUp2kuZL6pvX90uslafuorn0LZmbWE5RSQ3oXODIiDgD+CjhG0hjgImB2RNQCq4CpqfxUYFVEfASYncqZmZm1qd2EFJm16WWf9AjgSOAXaf11wGfT8sT0mrT9KEkqW8RmZtYjlXQPSVKNpCeBBmA+sBR4KyI2piLLgF3T8q7AqwBpeyOwfTmDNjOznkcRUXphaQhwK/CvwDWpWQ5JI4G7ImJ/Sc8AR0fEsrRtKXBIRLzZdJzGxsbmk9bV1ZXljVjP9qs/11Q6hKpxwoc2VToE20rV1tY2Lw8ePLjDLWMd+h1SRLwlaQEwBhgiqXeqBY0Alqdiy4CRwDJJvYHBwMrWjln4BrpTXV1dxc7dEdUSJ3RtrMM3l2/G2PqGeobvNLxsx+tKWxJrbW33/w7J/07Lr1rihPLFWkovux1TzQhJ2wB/DTwHPAh8PhWbBNyWlm9Pr0nbH4iOVMPMzGyrVEoNaWfgOkk1ZAlsXkTcKelZ4BZJ/w48Afwklf8JcIOkJWQ1o5O7IG4zM+th2k1IEbEYOLDI+heBQ4qsXwecVJbozMxsq+GRGszMLBeckMzMLBeckMzMLBeckMzMLBeckMzMLBeckMzMLBeckMzMLBc8hbnlRlvTcle7fZ+4u+j6Zw4c382RmOWXa0hmZpYLTkhmZpYLbrIz62HaavqcvFf3jwRuVionJLOtSGvJyonK8sBNdmZmlgtOSGZmlgtOSGZmlgtOSGZmlgtOSGZmlgtOSGZmlgvu9m3WAa0NAQQeBsiss1xDMjOzXHANybpVngZQdW3HLF9cQzIzs1xwDcmsiLZqT+Xcx8ze4xqSmZnlgmtIZhXk+1hm73ENyczMcsEJyczMcsEJyczMcsEJyczMcqHdhCRppKQHJT0n6RlJZ6X1wyTNl1SXnoem9ZL0fUlLJC2W9PGufhNmZlb9SqkhbQTOiYi9gTHAGZL2AWYA90dELXB/eg0wHqhNjy8DPyp71GZm1uO0m5Ai4vWI+ENaXgM8B+wKTASuS8WuAz6blicC10fmYWCIpJ3LHrmZmfUoiojSC0ujgN8B+wGvRMSQgm2rImKopDuBWRHxUFp/PzA9Ih5vKtvY2Nh80rq6us6+B6siv/pzTaVDaHbw8wsqHUKbHv/oEZUOAYATPrSp0iFYlaitrW1eHjx4sDq6f8k/jJU0APgl8C8RsVpq9VzFNrSa9QrfQHeqq6ur2Lk7olrihNJiHb65ewdXLfbD09VrVjNo4CAYOKhbY+mo4TsNp76hnuE7Da9oHLW127Vbpqf9O82DaokTyhdrSb3sJPUhS0Y3RcSv0ur6pqa49NyQ1i8DRhbsPgJY3ulIzcysR2u3hqSsKvQT4LmI+I+CTbcDk4BZ6fm2gvVnSroFGA00RsTrZY3abCvX2pBDHm7IqlkpTXaHAf8APCXpybTum2SJaJ6kqcArwElp213AscAS4G1gSlkjNjOzHqndhJQ6J7R2w+ioIuUDOKOTcZmZ2VbGIzWYmVkuOCGZmVkueD4ks5za94m7GblmNYNey3f3dLNycQ3JzMxywQnJzMxywQnJzMxywQnJzMxywQnJzMxywQnJzMxywd2+rUdobWw3M6seriGZmVkuOCGZmVkuuMnOrAdpq+nSU1NY3jkhWdXwfaLOcbKyvHNCMjNP+Ge54HtIZmaWC05IZmaWC26ys7K79oW/VDoEM6tCriGZmVkuOCGZmVkuOCGZmVku+B6SmW2xpvuF9Q01DN/8/nuHk/farhIhWRVzDcnMzHLBCcnMzHLBTXZm1iZ347fu4hqSmZnlghOSmZnlghOSmZnlQrsJSdJPJTVIerpg3TBJ8yXVpeehab0kfV/SEkmLJX28K4M3M7Oeo5RODdcC/wlcX7BuBnB/RMySNCO9ng6MB2rTYzTwo/RsPUyxG93FfotiZlaqdmtIEfE7YGWL1ROB69LydcBnC9ZfH5mHgSGSdi5XsGZm1nNt6T2k4RHxOkB63imt3xV4taDcsrTOzMysTeX+HZKKrIu2dqirqytzCKWr5Lk7Io9x1jfUtLK+vsvOOXLN6rIda3UZj9XVKhlrR/6eLcvW9dpU7nDKJo+fqWKqJU7IYq2tre3UMbY0IdVL2jkiXk9Ncg1p/TJgZEG5EcDytg7U2Tewpcpx8bpDXuMsdq+ovqGe4TsN77JzDnptUFmOs3rNagYNLM+xulqlYy3171nsb19bm8+x7PL6mWqpWuKE8sW6pU12twOT0vIk4LaC9aem3nZjgMampj0zM7O2tFtDknQzcASwg6RlwLeBWcA8SVOBV4CTUvG7gGOBJcDbwJQuiNnMzHqgdhNSRHyhlU1HFSkbwBmdDcq2Xvs+cXelQzCzCvFIDWZmlgtOSGZmlgtOSGZmlgtOSGZmlgtOSGZmlgueMdbMWtVWr8dnDhzf5r5tzTQ7ea98/mjWKssJyVrlqavNrDu5yc7MzHLBNSSrCP8AtvoV/g1Hrln9vrEG22vOMyvGCcncNGdmueAmOzMzywXXkMys7FprknVTnrXFCcnMcqW1JmR3Fe/53GRnZma54BqSdRn3pDOzjnANyczMcsE1pK2Eu3abWd65hmRmZrngGpKZdTvX2K0Y15DMzCwXXEOykviHjmbW1ZyQrFPaGmDTrKUt+SmA/9Oz9XBC6kHcLm9m1cz3kMzMLBdcQ6oyrgXZ1qapma/38n4f2Lbx08d3dzjWhZyQtjIezsfM8spNdmZmlguuIeVQYbNcfUMNwze7mc7Mej4nJDOreh25t1r4nzzPsZQvXZKQJB0DzAFqgKsjYlZXnMda53tF1tP8T/27H1j3TBuJqLXPQP2uHy9bTFZeZU9IkmqAHwJ/AywDHpN0e0Q8W+5zVTP3ljPrPP/Hq2dRRJT3gNKhwHci4uj0+jyAiJjZVKaxsbG8JzUzs1wZPHiwOrpPV/Sy2xV4teD1srTOzMysVV2RkIplRdeIzMysTV3RqWEZMLLg9QhgeWGBLanKmZlZz9YVNaTHgFpJu0vqC5wM3N4F5zEzsx6k7AkpIjYCZwL3As8B8yLimXKfpzWSjpH0gqQlkmYU2X64pD9I2ijp8y22bZL0ZHp0eRItIdazJT0rabGk+yXtVrBtkqS69JiU4zjzdk3/WdJTKZ6HJO1TsO28tN8Lko7OY5ySRkl6p+Ca/rgr4ywl1oJyn5cUkg4uWJeba9panHm8ppImS3qjIKbTCrbl6bPfVpwd/+xHRI95kP3uaSnwf4C+wB+BfVqUGQV8DLge+HyLbWtzFuungW3T8unA3LQ8DHgxPQ9Ny0PzFmdOr+mgguUJwD1peZ9Uvh+wezpOTQ7jHAU8nadrmsoNBH4HPAwcnMdr2kacubumwGTgP4vsm7fPftE407YOf/Z72lh2hwBLIuLFiFgP3AJMLCwQES9HxGJgcyUCLFBKrA9GxNvp5cNk9+MAjgbmR8TKiFgFzAeOyWGc3a2UWFcXvNyO9zrcTARuiYh3I+IlYEk6Xt7i7G7txpp8F7gYWFewLlfXtI04u1upsRaTq89+ufW0hNTZLuf9JT0u6WFJny1vaB/Q0VinAk2/AuzOrvWdiRNyeE0lnSFpKdkX01c7sm8O4gTYXdITkn4r6ZNdFGOTdmOVdCAwMiLu7Oi+ZdSZOCFn1zQ5MTWD/0JSU0exXF3TpFicsAWf/Z6WkDrb5fzDEXEw8PfAZZL2KE9YRZUcq6QvAgcDl3R03zLoTJyQw2saET+MiD2A6cD/7ci+ZdKZOF8nu6YHAmcDP5PUlfPGtxmrpF7AbOCcju5bZp2JM1fXNLkDGBURHwPuA67rwL7l0pk4YQs++z0tIbXb5bwtEbE8Pb8ILAAOLGdwLZQUq6S/Bs4HJkTEux3ZNwdx5vKaFrgFaPqfW+6uaYHmOFPz15tp+fdkbfx7dlGc0H6sA4H9gAWSXgbGALenDgN5uqatxpnDa0pEvFnwOboKOKjUfXMS55Z99rviZlilHmS/q3qR7AZq0024fVspey0FnRrIbhD2S8s7AHUUuSnanbGmP+BSoLbF+mHASynmoWl5WA7jzOM1rS1YPh54PC3vy/tvwL9I192A70ycOzbFRXaz+bWu+tuXGmuL8gt4r7NArq5pG3Hm7poCOxcsfw54OC3n7bPfWpxb9NnvkgteyQdwLPC/6Qvy/LTuArL/uQN8gizz/wV4E3gmrR8LPJUu+lPA1BzEeh9QDzyZHrcX7PuPZDeJlwBT8hhnTq/pHOCZFOeDhR8wshreUuAFYHwe4wROTOv/CPwBOL7S17RF2QWkL/q8XdPW4szjNQVmFsT0IPDRgn3z9NkvGueWfvbLPriqmZnZluhp95DMzKxKOSGZmVkuOCGZmVkuOCGZmVkuOCGZmVkuOCGZdRFJn5T0QqXjMKsW7vZt1k3SCAGnRcR9lY7FLI9cQ7IeQVJXzH681fF1tEpyQrJck/RymuTtWUmrJF0jqb+kIyQtkzRd0p+BayT1k3SZpOXpcZmkfuk4TeW/KWlFOu4pBefpJzh73z0AAARSSURBVOlSSa9Iqpf0Y0nbtNj3HEkNkl6XNKVg32NTfGskvSbp3ML90vINwIeBOyStlfQNSf8l6Sst3u/itkZGVmZ2iqMxld8vbdtG0vck/Slte6jgPUyQ9IyktyQtkLR3i2s8XdJi4C+SeksaI2lRKv9HSUd08k9p1r6uHiLDDz868wBeBp4mG+RxGPDfwL8DRwAbgYvIxkrbhmxIk4eBncjGJ1sEfDcdp6n8f6TynyIbPmqvtP0y4PZ0joFkoxjPbLHvBUAfsuFU3iZNjEY2WvQn0/JQ4OMF+y1r8V7+uuD13wGPFLw+gGw4q75tXI+jgd8DQ8hGY96bNJ4Y8EOyIXF2JZtcbWx6r3um9/o3Kf5vkA0707cgrifTNd4m7f9mep+90n5vAjtW+t+DHz37UfEA/PCjrUf6svzngtfHko2rdQSwHuhfsG0pcGzB66OBl9NyU1LZrmD7POBb6Yv9L8AeBdsOBV4q2PcdoHfB9gZgTFp+BfgnCmZ5LdivrYTUD1hJGkgVuBS4vJ3rcSTZ2GJjgF4F63ulGA8oss+3gHktyr4GHFEQ1z8WbJ8O3NDiGPcCkyr978GPnv1wk51Vg8JJwv4E7JKW34iIwpk/d0nbi5UFWBURfymyfUdgW+D3qYnqLeCetL7JmxGxseD128CAtHwiWaL8U5rg7dBS3lRkw/bPA76Y5uv5AnBDO/s8APwnWW2oXtKVae6eHYD+ZEm5pfddl4jYTHZNCydbK7zGuwEnNV2LdD3GATuX8r7MtpQTklWDwjlZPsx7c7K07CK6nOzLtFhZgKGStiuyfQVZ7WLfiBiSHoMjYgAliIjHImIiWVPhr8mSTNGiRdZdB5wCHAW8HRH/U8L5vh8RB5FN77An8PX0HtYBxSZBe991kSSya/paK7G9SlZDGlLw2C4iZrUXm1lnOCFZNThD0ghJw4BvAnNbKXcz8H8l7ShpB+BfgRtblPk3SX2VTVP9t8DPU43hKmC2pJ0AJO0q6ej2AkvHOkXS4IjYAKwGNrVSvJ5svp1mKQFtBr5HO7WjdL5PSBotqQ9ZM+M6YFN6Dz8F/kPSLpJqJB2aOnXMA46TdFTa7xzgXbJ7bMXcCBwv6eh0nKZOJCPai8+sM5yQrBr8DPgN2WRhL5J1aijm34HHgcVkc7D8oUXZPwOryGoMN5Hdm3o+bZtOdqP/YUmryeZ42qvE+P4BeDnt98/AF1spN5MsYb7V1BMvuR7Ynw8mz2IGkSXPVWTNcG+S3XsCOJfsfT9Gdm/qIrL7TC+kmH5AVpM6nmzOn/XFThARrwITyZL/G2Q1pq/j7wvrYv5hrOVauX5Mmrot3xgRuftfvqRTgS9HxLhKx2JWSf4fj1kFSdoWmAZcWelYzCrNCcmsQtI9qjfI7i39rGD9J9OPZz/wqFiwZt3ATXZmZpYLriGZmVkuOCGZmVkuOCGZmVkuOCGZmVkuOCGZmVkuOCGZmVku/H8ew3XCaj9CUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(data_ps.query(\"intervention==0\")[\"propensity_score\"], kde=False, label=\"Non Treated\")\n",
    "sns.distplot(data_ps.query(\"intervention==1\")[\"propensity_score\"], kde=False, label=\"Treated\")\n",
    "plt.title(\"Positivity Check\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use our propensity score weighting estimator to estimate the average treatment effect. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Y1: 0.259590031959361\n",
      "Y0: -0.12892940555758584\n",
      "ATE 0.3885194375169472\n"
     ]
    }
   ],
   "source": [
    "weight = ((data_ps[\"intervention\"]-data_ps[\"propensity_score\"]) /\n",
    "          (data_ps[\"propensity_score\"]*(1-data_ps[\"propensity_score\"])))\n",
    "\n",
    "y1 = sum(data_ps.query(\"intervention==1\")[\"achievement_score\"]*weight_t) / len(data)\n",
    "y0 = sum(data_ps.query(\"intervention==0\")[\"achievement_score\"]*weight_nt) / len(data)\n",
    "\n",
    "ate = np.mean(weight * data_ps[\"achievement_score\"])\n",
    "\n",
    "print(\"Y1:\", y1)\n",
    "print(\"Y0:\", y0)\n",
    "print(\"ATE\", np.mean(weight * data_ps[\"achievement_score\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Propensity score weighting is saying that we should expect treated individuals to be 0.38 standard deviations above their untreated fellows, in terms of achievements. We can also see that if no one got the treatment, we should expect the general level of achievements to be 0.12 standard deviation lower than what it is now. By the same reasoning, we should expect the general level of achievement to be 0.25 standards deviation higher if we've given everyone the seminar. Compare this to the 0.47 ATE estimate we've got by simply comparing treated and untreated. This is evidence that the bias we have is indeed positive and that controlling for X gives us a more modest estimate of the impact of the growth mindset.\n",
    "\n",
    "## Standard Error\n",
    "\n",
    "![img](./data/img/ps/bootstrap.png)\n",
    "\n",
    "To compute the standard error for the IPTW estimator, we can use the formula of the variance of a weighted average.\n",
    "\n",
    "$\n",
    "\\sigma^2_w = \\dfrac{\\sum_{i=1}^{n}w_i(y_i-\\hat{\\mu})^2}{\\sum_{i=1}^{n}w_i}\n",
    "$\n",
    "\n",
    "However, we can only use this if we have the true propensity score. If we are using the estimated version of it, \\\\(\\hat{P}(x)\\\\), we need to account for the errors in the estimation process. The easiest way of doing this is by bootstrapping the whole procedure. This is achieved by sampling with replacement from the original data and computing the ATE like we did above. We then repeat this many times to get the distribution of the ATE estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed # for parallel processing\n",
    "\n",
    "# define function that computes the IPTW estimator\n",
    "def run_ps(df, X, T, y):\n",
    "    # estimate the propensity score\n",
    "    ps = LogisticRegression(C=1e6).fit(df[X], df[T]).predict_proba(df[X])[:, 1]\n",
    "    \n",
    "    weight = (df[T]-ps) / (ps*(1-ps)) # define the weights\n",
    "    return np.mean(weight * df[y]) # compute the ATE\n",
    "\n",
    "np.random.seed(88)\n",
    "# run 1000 bootstrap samples\n",
    "bootstrap_sample = 1000\n",
    "ates = Parallel(n_jobs=4)(delayed(run_ps)(data_with_categ.sample(frac=1, replace=True), X, T, Y)\n",
    "                          for _ in range(bootstrap_sample))\n",
    "ates = np.array(ates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ATE is then the mean of the bootstrap samples and the standard error is the standard deviation of these samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATE 95% CI: 0.387743504939349 +- 0.0338185831687405\n"
     ]
    }
   ],
   "source": [
    "print(f\"ATE 95% CI: {ates.mean()} +- {1.96*ates.std()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also have a visual on what the bootstrap samples look like, along with the confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEGCAYAAABRvCMcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAdJElEQVR4nO3deZwdVZ338c+3O5sJJMEliwQeQNsYQFSeyJpxIjiAywACKiqEIOr4GBEEhk1G4yCyiALPsDzjEElwQRBhwAhKyBh5gcCoGJFFuDFkQiBJo5DGmIQQ8nv+qOpw+6aXm7uc27f7+3696tW3Tp2q+p271K9P1bl1FRGYmZnVW0ujAzAzs8HBCcfMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMWtykhZJurZO254taUlP83XY31xJd9dr+9ZYTjiDlKQJkjZIWiVpaF42XVL0Mc3N6/a0/Mpe9jmzpO56SX+UdHod2jcp38f0bVxviaTZtY5nW0laVvQ8vSRppaS7JH1K0pCS6kcBp5W53W19Xi4F9tuG0Msi6ThJ3X0J8BTgw7Xen/UPpW9cGzw+CfwUeCtwBHAz8CtgYlGdU4FPAO8qKltf9PjzwI9Ltvu3Pvb7CjApf/wa4BDgSkntEfHdbWlAI0kaFhEb67ybi4HLyT6nE4D3kiWAEyQdGhHrACLi+VrvWFILoIhYC6yt9fZ7EhEdqfZlDRARngbZRNazfQo4HDgTuKuHeucBy3pYFsBx27jfmcCmbsp/C1xRNC/gDGApsBH4E3BqyTrbA/8OPAdsAH4DHFISX/G0LC+fRJYk/0yWPJcC/5wvW9TNersA0/PHHwDuzff3eWAH4HvA8nxbTwCnkx2oO+OYC9xN1gN5BliX7//1fTxXy4Dzuil/B/Ay8OWiskXAtUXz04D7gL/m0++BQ/t4XmYDS4CPAn8ENgF7dpYXbbuz3sfz525D3r5dS+uUxD2tm+ezeJpb/Hxt43thGfCvwBXA88BqssTc2ujPmqeuk0+pDU6HAKOAO4DvAtMl7ZY6CGWmA1PIeledPgecD1wE7AF8A7hI0klFdb4DHAocB7yT7AA7X9Jb8+V753+PJuu1dfbSrgbGkPUWpgAnASvyZUeRHby+ma8zEXi6aJ/fBC7J1/tPYDjwB+BIYPc85q+SJdZi+5AdZA8D3g/slce/zSJiMfBz4CPdLZfUCtwOPEj2HOxNlgDW5VV6el4A3kj23M/M2/M/PYQxMa/3UeDvyJL/f0pSmc34FVnC7tzWRLJTad0p570AcDKwEtgX+AJZ73xGmfFYKo3OeJ7ST8CtwGVF83cAX++mXl89nA1kp1uKp4/2st+Z+XqddV/O5y8tqfc0cElJ2WXA0vzxm/P13l9S5yHgO/njSXmd6SV1fg/M7iXGJaXLefU/8uPLeG6vABYUzc/N2zqmqOyQfHttvWxnGd30cPJlFwHriuYXkfdwyHpdW7W7qG5Pz8tsYDOwczflpT2cAN5cVPaWvOy93a2Tl23p4eTzx2WHn63im0vXHk6v74Wi5+r2kjo/A25I+bny1PfkHs4gI2ki8EFgXlHxXODEbi5G9+VLZKd4iqef9rHOK0V13052amZG54V6SaPJDor3lKz3S2AXSSPJ/vummzr3kP0X3JvLgXMlPSjpYknv7qN+sf8unpHUIulsSYsl/VnSWuCzwP8qWe+x6Hpt4r7875Rt2HeXXZMdvLcSES8A1wI/l3RnHt/kMre7OiKWl1HvuYjYMlItIp4kO0W5e8+rbLsy3wudFpfUeQYYX8t4rHpOOIPPSWQXoX8jaZOkTcAPyC5KH76N21odEUtKpj4vMBfVfSwibiD7j/VcScOLq5WsVs7pmh4PxEX7vo4sIfw/slM5d0r6Xhnbhq0HRJwOnAP8G/APZEn0WmBYmdur1J5k1zK6FRGfBv43sAD4e+ARSf9Uxnb7GvDRm+LXZzNbv15Dq9h2Oe+F0gEcgY9v/Y5fkEEkH3n0KeDrbN0z+R7wmQaFtoksCQ6PiBfJrqn8fUmddwNPRTYy69GismJ/V7Ss8wDUWrqziFgZEddFxAyyBPyJ/L/pzvW2WqcH7wZ+FhFzIuJ3+X/9bd3Um1K0fYAD8r+Pl7mfLSS9g+za1Y291YuIRyLiWxHxPmAOr762PT4v2+ANkt5UFNNbgNfxanvagXH59aROe9PVxnzdHuMo871gTcTDogeXw4CdgX8vPXUi6TpggaRdImJZmdsbI2lCSdmGiFjT20pF6wwjO612KrAwP8AAXAh8U1KB7PrEQcD/AWYBRMSfJP0IuDr/z/1/8uV7kp2ig+wUz1rgEEmPAi9FxAv594TuIBtRNoJsoMDTZKO5IBu9d6CknckutPc25PgJ4HhJ7yE7hTOD7KL1CyX1Arhe0nnAa4GrgJ9GRKG35wnYLn+uiodFn0U2Uu5b3a0g6c3Ap4Gf5O16I1kifqi356WPOEqtA66T9EWy3sa/kQ2e6PzC5i+AkcD5kuaQJZtZJdt4Kv97uKR7gfU99I57fS9Yk2n0RSRP6SbgNuD+Hpa1AquArxWV9TVooLtpfi/7n1lS92WyA8+VwOuK6gn453zZy2RDYkuHwo7m1WHRL1EyLDqvM6NoG8vysquAJ8mGMf+F7JrTHkXrTCUbpr2erYfxTirZ/hjgJuDFfFtXkY2oWlZUZy7ZgfgMslFU68kGbbyhj9dqWdHztDF/be4i66G2ltRdxKuDBiYCt5D1DF4CngX+g66DFrp7XmZTcqG/u3JeHRZ9XB7jS8B/AW8qWe+T+eu2HrgTOJaiQQN5ncvJhjBvpvdh0X29F5ZRMsCC7NTmokZ/5jx1nZS/OGZWB/mdGSZFxHsbHYtZo/kajpmZJeGEY2ZmSfiUmpmZJeEejpmZJdGQYdEdHR3uVpmZDXBjxozp8iVd93DMzCwJJxwzM0vCCacGCoW+vjA+sAy29sLga/Ngay8MvjY3or1OOGZmloQTjpmZJeGbd5qZ9SEiWLt2LZs3b250KDUzYsQIOjo6+q7Yi5aWFrbbbjvK/bFXJxwzsz6sXbuW4cOHM2xYvX/qKJ3hw4czYsSIqraxceNG1q5dy/bbb19WfZ9SMzPrw+bNmwdUsqmVYcOGbVOvzwnHzMyScMIxM2sC11xzDfvvvz/77bcfV1999ZbyCy+8kClTpjBt2jSmTZvGXXfdBcADDzzAAQccwHve8x6WLl0KwJo1azjqqKPo6R6aL7/8MrNnz2bvvfdm//3356CDDmLBggUAvO1tb+Mvf/lLVW3wNRyzBpr7xN8qWm/m5FE1jsT6s8cee4zrr7+ehQsXMmzYMI4++mgOPfRQ3vSm7Je+P/e5z3HyySd3WefKK6/ku9/9LsuXL2fOnDlccMEFfOMb3+C0007r8SL/BRdcwKpVq7j//vsZPnw47e3t3HfffTVrh3s4ZmbbaOzYsV2mnsydO7dLvVNOOaWi/T355JNMnTqVkSNHMmTIEA488EDmz5/f6zpDhw5l/fr1rFu3jqFDh/LUU0+xcuVKpk2b1m39devWMW/ePC655BKGDx8OwLhx4/jQhz5UUczdccIxM+vnpkyZwq9+9Suef/551q1bx4IFC1ixYsWW5d/+9rc54IADmDVrFmvWrAHgi1/8IqeeeirXXHMNn/70pzn//PP50pe+1OM+li5dyqRJkxg9enTd2tFnwpH0HUntkh4pKvuGpD9KeljSrZLGFi07R9ISSU9IOrRegZuZDRaTJ0/mlFNO4cgjj+Too49mzz33ZMiQ7IrISSedxOLFi7n33nuZMGHClqSy1157cffddzN//nyWLVvGhAkTiAhOPPFEPvOZz/Dcc88lb0c5PZy5wGElZQuAPSNiL+BJ4BwASbsDxwJ75OtcLam1ZtGamQ1SM2bM4J577uHOO+9khx122HL9Zty4cbS2ttLS0sKMGTN46KGHuqwXEVx66aWceeaZXHzxxZxzzjl85CMf4dprr+1Sb7fddmPFihX89a9/rVsb+kw4EXEP8HxJ2V0RsSmffQCYlD8+AvhhRLwUEU8BS4B9ahivmVnDrVmzpsvUk5kzZ3apd8UVV1S8z84eydNPP81PfvITjjnmGABWrVq1pc78+fOZMmVKl/V+8IMfcMghhzB27FjWrVtHS0sLLS0trF+/vku9kSNHcvzxx3PWWWexcePGLdu+8cYbK465VC1GqX0S6IxoR7IE1GlFXmZmZlWYMWMGzz//PEOGDOHSSy/dMljhy1/+Mo88kl3x2Hnnnbn88su3rLNu3TpuuOEGbr31VgBmzZrFjBkzGDp0aJeh1Z3OO+88vva1r7HvvvsyYsQIRo4cybnnnluzNqin8dhdKkm7APMjYs+S8i8BU4GjIiIkXQXcHxHfy5fPAe6IiB8Xr1f8i5+D7ZbgZsVuWVXZGeejJrxS40isNyNGjOANb3hDo8Pol5577jk2bNiwZb6trW3L49Jf/Ky4hyPpBOCDwMHxatZaAexUVG0S8Gxv2ykOrlkVCoUB0Y5yDbb2Qv3aPH5zZd/DaWur7/dw/Bp31dHRUfV9x/qbDRs21KRNo0ePZqedduq7IhUOi5Z0GHAWcHhErCtadDtwrKThknYF2oD/rmQfZmY2sPTZw5F0AzAdeL2kFcBXyEalDQcW5N9YfSAiPhsRj0q6CXgM2ATMigj3/c3MrO+EExEf66Z4Ti/1LwAuqCYoMzMbeHynATOzPrS0tGwZKmyv2rhxIy0t5acR37zTzKwP2223HWvXrt3quyvN7MUXX6z6Njadv/hZLiccM7M+SCr7Vy2bRXt7e9mjy2rFp9TMzCwJJxwzM0vCCcfMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMTOzJJxwzMwsCSccMzNLwgnHzMyScMIxM7MknHDMzCwJJxwzM0vCCcfMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMTOzJJxwzMwsCSccMzNLos+EI+k7ktolPVJU9lpJCyQV8r875OWS9H8lLZH0sKS96xm8mZk1j3J6OHOBw0rKzgYWRkQbsDCfB3gf0JZPnwGuqU2YZmbW7Ib0VSEi7pG0S0nxEcD0/PE8YBFwVl5+fUQE8ICksZImRsTKWgVsZjD3ib9VtN7MyaNqHIlZ+fpMOD0Y35lEImKlpHF5+Y7A00X1VuRlPSacQqFQYQj9y0BpR7kGW3uhPm1e3d5a8232ptDySvl1/RoPePVob1tbW4/LKk04PVE3ZdHbCr0F1ywKhcKAaEe5Blt7oX5tHr+5sp5Kpdrayuvh+DUe+BrR3kpHqa2WNBEg/9uel68AdiqqNwl4tvLwzMxsoKg04dwOnJA/PgG4rah8Rj5abT+gw9dvzMwMyjilJukGsgECr5e0AvgKcBFwk6STgOXAh/PqdwDvB5YA64AT6xCzmZk1oXJGqX2sh0UHd1M3gFnVBmVmZgOP7zRgZmZJ1HqUmtmgVel3Y8wGC/dwzMwsCSccMzNLwgnHzMyScMIxM7MknHDMzCwJJxwzM0vCCcfMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMTOzJJxwzMwsCSccMzNLwj9PYFaiu58ZWN3eyvjN/vkBs2q4h2NmZkk44ZiZWRJOOGZmloQTjpmZJeGEY2ZmSTjhmJlZEk44ZmaWRFXfw5H0ReBTQAB/AE4EJgI/BF4LPAQcHxEbq4zTzGqgu+8Ydaf4e0czJ4+qZ0g2iFTcw5G0I/AFYGpE7Am0AscCFwOXRUQb8AJwUi0CNTOz5lbtKbUhwGskDQFGAiuBg4Cb8+XzgCOr3IeZmQ0AFSeciHgGuBRYTpZoOoDfAmsiYlNebQWwY7VBmplZ86v4Go6kHYAjgF2BNcCPgPd1UzV6206hUKg0hH5loLSjXAO5vavbW3soX504ksbqbG+h5ZUGR5LOQH5fd6ce7W1ra+txWTWDBt4LPBURzwFIugU4ABgraUjey5kEPFtpcM2iUCgMiHaUa6C3t7ubdK5uX834ceMbEE1jFLe3rW1wDBoY6O/rUo1obzXXcJYD+0kaKUnAwcBjwC+AY/I6JwC3VReimZkNBNVcw3mQbHDAQ2RDoluAbwNnAadJWgK8DphTgzjNzKzJVfU9nIj4CvCVkuKlwD7VbNfMzAYe32nAzMyS8C9+2oBV7rfqzSwN93DMzCwJ93DMrFeV9hR9DzYr5R6OmZkl4YRjZmZJOOGYmVkSTjhmZpaEE46ZmSXhhGNmZkk44ZiZWRJOOGZmloQTjpmZJeGEY2ZmSTjhmJlZEk44ZmaWhBOOmZkl4YRjZmZJOOGYmVkSTjhmZpaEE46ZmSXhhGNmZkk44ZiZWRJOOGZmlkRVCUfSWEk3S/qjpMcl7S/ptZIWSCrkf3eoVbBmZta8qu3hXAH8LCLeCrwdeBw4G1gYEW3AwnzezMwGuYoTjqTRwLuBOQARsTEi1gBHAPPyavOAI6sN0szMml81PZzdgOeA6yT9TtK1kkYB4yNiJUD+d1wN4jQzsyaniKhsRWkq8ABwYEQ8KOkK4EXg5IgYW1TvhYjoch2no6Njy04LhUJF+zfryy2rWhsdwqB21IRXGh2CNUBbW9uWx2PGjFHxsiFVbHcFsCIiHsznbya7XrNa0sSIWClpItBebnDNqlAoDIh2lKtZ2jt+899qtq3V7asZP258zbbX39WivW1to2oUTRrN8r6ulUa0t+KEExGrJD0taXJEPAEcDDyWTycAF+V/b6tJpDaozX2idsnDzBqjmh4OwMnA9yUNA5YCJ5JdF7pJ0knAcuDDVe7DzMwGgKoSTkQsBqZ2s+jgarZrZmYDT7U9HLNt4lNjZoOXb21jZmZJOOGYmVkSTjhmZpaEE46ZmSXhhGNmZkk44ZiZWRJOOGZmloQTjpmZJeGEY2ZmSTjhmJlZEk44ZmaWhBOOmZkl4YRjZmZJOOGYmVkSTjhmZpaEE46ZmSXhH2Azs7qo9Mf2Zk4eVeNIrL9wD8fMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMTOzJJxwzMwsiaoTjqRWSb+TND+f31XSg5IKkm6UNKz6MM3MrNnVoodzCvB40fzFwGUR0Qa8AJxUg32YmVmTqyrhSJoEfAC4Np8XcBBwc15lHnBkNfswM7OBQRFR+crSzcCFwPbAGcBM4IGIeHO+fCfgzojYs3i9jo6OLTstFAoV798a55ZVrY0OwQaooya80ugQrAptbW1bHo8ZM0bFyyq+l5qkDwLtEfFbSdM7i7up2mtGKw6uWRUKhQHRjnIVCgXGjxvf6DCSWt2+elC1uZHtbWtrzL3UBuPnOHV7q7l554HA4ZLeD4wARgOXA2MlDYmITcAk4NnqwzQzs2ZX8TWciDgnIiZFxC7AscB/RcQngF8Ax+TVTgBuqzpKMzNrevX4Hs5ZwGmSlgCvA+bUYR9mZtZkavJ7OBGxCFiUP14K7FOL7ZqZ2cDhOw2YmVkS/sVPM+tXKvmlUP9KaHNwD8fMzJJwwjEzsySccMzMLAknHDMzS8KDBsys6VUy0AA82CA193DMzCwJJxwzM0vCCcfMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMTOzJJxwzMwsCSccMzNLwgnHzMyS8L3UbJvvQ7W6vZXx4+oUjJkNWO7hmJlZEk44ZmaWhBOOmZkl4YRjZmZJOOGYmVkSFSccSTtJ+oWkxyU9KumUvPy1khZIKuR/d6hduGZm1qyq6eFsAk6PiCnAfsAsSbsDZwMLI6INWJjPm5nZIFdxwomIlRHxUP74r8DjwI7AEcC8vNo84MhqgzQzs+ZXk2s4knYB3gk8CIyPiJWQJSXAXxE0MzMUEdVtQNoO+CVwQUTcImlNRIwtWv5CRHS5jtPR0bFlp4VCoar9W9/e9a53dZn/9a9/3WX+llWtKcMx6zeOmvDKVmV9fV6sd21tbVsejxkzRsXLqrq1jaShwI+B70fELXnxakkTI2KlpIlAe7nBNatCodCv27Fo0aIu86Wxjt+8rbe2Wc34ceOrDaupDLY2D5b2trWN2vK483Pc1+dloGjEcavihCNJwBzg8Yj4VtGi24ETgIvyv7dVFaFV7R3veEejQzBrGv681E81PZwDgeOBP0hanJedS5ZobpJ0ErAc+HB1IZqZ2UBQccKJiHsB9bD44Eq3a2ZmA5PvNGBmZkk44ZiZWRL+AbYBpKcfUnv68Ye7zO80Za8U4Zg1pcWLF3eZ9yCC2nHCGQS+OfN9XeYvf/CZBkVi1v9Nnz69y/yaNWsaE8gA5FNqZmaWhHs4ZjZoFZ+GXt3e2u2XoLs7VT1z8qityqxvTjiDwKTJb2t0CGZNw5+X+nHCGQTOuP5njQ7BrGn481I/voZjZmZJOOGYmVkSTjhmZpaEE46ZmSXhhGNmZkl4lNogcOPXz+wy/9FzL2lQJGb9nz8v9eOEMwjcf9v3u8z7A2TWM39e6sen1MzMLAknHDMzS8Kn1AaBj5x9caNDMGsa/rzUjxNOnfX0GzUpHfCh4xodglnTKOfzUsnn2jf89Ck1MzNLxAnHzMyScMIxM7MkfA3HzKwfq/Q6cH+8ZuSEU6beXvSefinQzKxTfxhA1GhOOIPAqfvu2GX+8gefaVAkZv2fPy/1U7eEI+kw4AqgFbg2Ii6q9T48NNHMrHnUZdCApFbgKuB9wO7AxyTtXo99mZlZc1BE1H6j0v7A7Ig4NJ8/ByAiLgTo6Oio/U7NzKxfGTNmjIrn6zUsekfg6aL5FXmZmZkNUvVKOOqmzL0aM7NBrF6DBlYAOxXNTwKe7Zwp7WaZmdnAV68ezq+BNkm7ShoGHAvcXqd9mZlZE6hLwomITcDngZ8DjwM3RcSj9dhXPUk6TNITkpZIOrub5Z+V9AdJiyXdWzoST9LOktZKOiNd1NWpps2S9pJ0v6RH8zoj0ka/7Sptr6Shkublyx7vHBjTDPpqc1G9YySFpKlFZefk6z0h6dA0EVen0vZK+gdJv81f499KOihd1NWp5jXOy+tz7IoIT91MZN8f+hOwGzAM+D2we0md0UWPDwd+VrL8x8CPgDMa3Z56t5ns9OzDwNvz+dcBrY1uUx3b+3Hgh/njkcAyYJdGt6kWbc7rbQ/cAzwATM3Lds/rDwd2zbfT9K9xL+19J/DG/PGewDONbk+921y0rC7HLt+8s2f7AEsiYmlEbAR+CBxRXCEiXiyaHUXRwAhJRwJLgWbq2VXT5kOAhyPi93m9v0TEKwlirkY17Q1glKQhwGuAjUBx3f6qzzbnzgcuATYUlR1BlmRfioingCX59vqzitsbEb+LiM5rz48CIyQNr3fANVDNa1zXY5cTTs/KGtotaZakP5G9cF/Iy0YBZwFfTRBnLVXcZuAtQEj6uaSHJJ1Z92irV017bwb+BqwElgOXRsTz9Q23Jvpss6R3AjtFxPxtXbcfqqa9xY4GfhcRL9U+xJqruM31PnY54fSsrKHdEXFVRLyJ7EU6Ly/+KnBZRKytY3z1UE2bhwDTgE/kfz8k6eB6BVoj1bR3H+AV4I1kp5dOl7RbvQKtoV7bLKkFuAw4fVvX7aeqaW9nnT2Ai4F/qnl09VFNm+t67PLNO3vW69DubvwQuCZ/vC9wjKRLgLHAZkkbIuLKukRaO9W0eQXwy4j4M4CkO4C9gYV1iLNWqmnvx8mu57wMtEu6D5hKdiqiP+urzduTXa9YJAlgAnC7pMPLWLc/qri9EfEbSZOAW4EZEfGnRDFXq5rXuL7HrkZf4OqvE1kyXkr232vnhbc9Suq0FT3+R+A33WxnNs0zaKDiNgM7AA+RXUAfAtwNfKDRbapje88CriP7b3IU8BiwV6PbVIs2l9RfxKsX0feg66CBpfT/QQPVtHdsXv/oRrcjVZtLymt+7HIPpwcRsUlS59DuVuA7EfGopH8lO+jcDnxe0nuBl4EXgBMaF3H1qmlzRLwg6Vtk38EK4I6I+GlDGlKmKl/jq8gSziNkSee6iHg4eSO2UZlt7mndRyXdRJZcNwGzop8PDKmmvWRf7Xgz8C+S/iUvOyQi2usbdXWqbHNd1eXmnWZmZqU8aMDMzJJwwjEzsySccMzMLAknHDMzS8IJx8zMknDCMTOzJJxwzMwsCSccMzNL4v8DuVGI3Mo9m0MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(ates, kde=False)\n",
    "plt.vlines(ates.mean()-1.96*ates.std(), 0, 20, linestyles=\"dotted\")\n",
    "plt.vlines(ates.mean()+1.96*ates.std(), 0, 20, linestyles=\"dotted\", label=\"95% CI\")\n",
    "plt.title(\"ATE Bootstrap Distribution\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common Issues with Propensity Score\n",
    "\n",
    "As a data scientist, I know it can be tempting to use all the power of the machine learning toolkit to make propensity score estimation as precise as possible. You can quickly get taken away by the all AUC optimisation, cross validation and bayesian hyper-parameter tuning. Now, I'm not saying you shouldn't do that. In fact, all of the theory about propensity score and machine learning is very recent, so there are lots we don't know yet. But it pays to understand something first. \n",
    "\n",
    "The first thing is that the predictive quality of the propensity score does not translate into its balancing properties. Coming from the field of machine learning, one of the most challenging aspects of getting acquainted with causal inference is letting go of treating everything as a prediction problem. In fact, maximising the prediction power of the propensity score can even hurt the causal inference goal. **Propensity score doesn't need to predict the treatment very well. It just needs to include all the confounding variables**. If we include variables that are very good in predicting the treatment but have no bearing on the outcome this will actually increase the variance of the propensity score estimator. This is similar to the problem linear regression faces when we include variables correlated with the treatment but not with the outcome. \n",
    "\n",
    "![img](./data/img/ps/ml-trap.png)\n",
    "\n",
    "To see this, consider the following example (adapted from Hernán's Book). You have 2 schools, one of them apply the growth mindset seminar to 99% of its students and the other to 1%. Suppose that the schools have no impact on the treatment effect (except through the treatment), so it's not necessary to control for it. If you add the school variable to the propensity score model, it's going to have a very high predictive power. However, by chance, we could end up with a sample where everyone in school A got the treatment, leading to a propensity score of 1 for that school, which would lead to an infinite variance. This is an extreme example, but let's see how it would work with simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>T</th>\n",
       "      <th>school</th>\n",
       "      <th>intercept</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.309526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1.571468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2.982024</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2.445420</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2.693187</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   T  school  intercept         y\n",
       "0  1       0          1  0.309526\n",
       "1  1       0          1  1.571468\n",
       "2  1       0          1  2.982024\n",
       "3  1       0          1  2.445420\n",
       "4  1       0          1  2.693187"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "school_a = pd.DataFrame(dict(T=np.random.binomial(1, .99, 400), school=0, intercept=1))\n",
    "school_b = pd.DataFrame(dict(T=np.random.binomial(1, .01, 400), school=1, intercept=1))\n",
    "ex_data = pd.concat([school_a, school_b]).assign(y = lambda d: np.random.normal(1 + 0.1 * d[\"T\"]))\n",
    "ex_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having simulated this data, we run bootstrap with the Propensity Score algorithm twice. The first including school as a feature to the propensity score model. The second time, we don't include school in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "ate_w_f = np.array([run_ps(ex_data.sample(frac=1, replace=True), [\"school\"], \"T\", \"y\") for _ in range(500)])\n",
    "ate_wo_f = np.array([run_ps(ex_data.sample(frac=1, replace=True), [\"intercept\"], \"T\", \"y\") for _ in range(500)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD1CAYAAABz79PWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAaUklEQVR4nO3dfXBV9b3v8fc3CU8VggRPYgp0ELvHKpYnI4WIlYKca4+3oKfSg20xXKm9DBfH01aN9taZ1rlTy9GpHh3HKz21MK1WKdbROlUuoGiVQ7AiCBTsBir4ADuHBxPAByR87x9Z7CbkaSf57afk85rJ7LXWXnutb35Ze3/yW0/b3B0REZFQCrJdgIiI9CwKFhERCUrBIiIiQSlYREQkKAWLiIgEVZSuBdfV1el0MxGRHm7w4MF2+jT1WEREJCgFi4iIBKVgaUc8Hs92CXlB7ZQ6tVXq1FapycV2UrCIiEhQChYREQkqbWeFtcXdOXr0KCdPnsz0qjutf//+1NXVZbuMrCooKGDgwIGYtTjxQ0SkVRkPlqNHj9KvXz/69u2b6VV3Wr9+/ejfv3+2y8iq48ePc/ToUQYNGpTtUkQkT3S4K8zMzjOzTU1+6s3sX82sxMxWmVk8ehySygpPnjyZF6Eijfr27ZsXvUsRyR0dBou7v+Xu49x9HHAR8CHwFHAbsMbdY8CaaFxERHq5zh68nw7scvc9wCxgWTR9GXBVyMLSqaSkhClTpjB58mSqqqr48MMPAbjnnnuYNGkSlZWVTJkyhY0bNzZ73ZYtW5gyZUpyfMWKFZSXl/Ppp58CsG3bNiorK1us77XXXmP69OlMmTKFiRMnctddd7Vb37Bhw7r7KwJw11138cADDwRZlohIqjp7jGUO8NtouMzd9wG4+z4zK23rRU3Ps+7fvz/9+vVLjv9658edLKF9cz/f8TGR/v37s3r1agAWLlzIkiVLqKio4LnnnmPlypX069ePgwcP8umnn/Lxx3+v79xzz+Wdd97hwIEDDBw4kHXr1vH5z3+e1157jQkTJvDqq69SUVHR7DUACxYsYMmSJYwePZqGhgZ27tzZYp6m3L3d51N14sSJFr9DV9TX11NbW9vuPLl4Ln13/H5/YZqWXAj7dyfH/vnshjStp2foadtVumS6nWKxWLvPpxwsZtYXmAnc3p0i6urqmh0Q79Mn7BsrlYPtZpacb8qUKWzbto1Ro0Zx1llnMXjwYKCx1/Dxxx+3WN6ECRPYunUrU6dOZevWrdxwww1s2rSJyspKNm7cyNSpU1u85sCBA3zuc59LTh87dizQeCLDrbfeyqZNmwCorq5m1qxZmBl33303zz//PAMGDOCxxx6jtLSUvXv3smjRIg4cOMBZZ53Fgw8+yIgRI9qcXlRURJ8+fbp9AkJxcTEjRoxo8/l4PN7hhpZvyk4eS8tyE7UJykrLkuOx2BlpWU9P0BO3q3TIxXbqzK6wrwIb3T0RjSfMrBwgemz/X9ocdOLECVatWsUFF1zAtGnTeO+997jooov4wQ9+wCuvvNLqayZOnEhNTQ3Hjh2joKCASy+9lA0bNgCwYcMGJk2a1OI1CxcupKKigm9961v86le/SvYg7r77boqLi1m3bh3r1q3jsssuA+DYsWNUVFTw6quvUllZybJljXscb7nlFubMmcO6deuYPXs21dXV7U4XEcmGzgTLtfx9NxjAM0BVNFwFPB2qqHT76KOPmDJlClOnTmX48OHMnTuXgQMH8tJLL3HfffcxdOhQrr/+eh5//PEWr500aRIbNmzg9ddfZ/z48Zxzzjns3r2bAwcOcOzYMUaOHNniNdXV1axdu5Zp06bxu9/9jq9//esArF27lhtuuCE535lnngk0nol1xRVXADBu3Dj27t0LNB6rmT17NgBz5sxh/fr17U4XEcmGlHaFmdlngBnA/2wy+WfAcjObD+wFZocvLz0GDBjQao+ksLCQSy+9lEsvvZTRo0fzm9/8hnnz5jWbp6Kigo0bN7J+/XomTpwIwGc/+1mefPLJ5HhrzjnnHObPn09VVRXnnnsuhw4dwt1bvfCwT58+yemFhYWcOHGi1WW2ddGiLmYUkWxKqcfi7h+6+1B3r2sy7aC7T3f3WPR4KH1lpl88HmfXrl3J8S1btjB8+PAW8w0aNIhhw4bx2GOPcfHFFwONu8ceeuihNoNl5cqVuDd+Pc2uXbsoLCxk8ODBTJs2jSVLliTn++CDD9qtceLEiTz55JMALF++PLnbra3pIiLZkPEr73PVsWPHuPXWW6mrq6OwsJBRo0axePHiVuedNGkSf/zjH5PBc/HFF3PnnXfypS99qdX5n3jiCX74wx8yYMAAioqK+MUvfkFhYSE333wzN998M5MnT6agoIDq6mpmzpzZZo2LFy9m0aJF3H///cmD9O1N74mWvpWeg+oiEo6d+k86tLa+QbKuri555lWua+2ssN6oo79ZJs9KyfdgOf2ssHnn6aywtuTi2U65KNvtpG+QFBGRtNOuMJFeIlO9PfXCRD0WEREJSsEiIiJBaVeYiPQIRS/+ITl84itfy2Iloh6LiIgE1SuDJVO3zV+xYgX33HMPAM8++yyVlZVcfPHFVFZW8uyzz7ZaWzwe58orr0zeYv+mm25q93f54he/yMGDBzvfCKd59NFHueWWW7q9HBGRrO8Ka9p9DSGVLnDTW7rccMMNPPLII0ycOJGVK1fy0ksvJW+bf+TIkWavGz16NO+88w5Hjhxh0KBB1NTUEIvFePPNN7nooota3IRy9erVLFiwgC1btnDHHXfw1FNPMXLkSN5++22uvvpqRo4cyYUXXthsHdXV1SxcuJArr7wSaAwrEZF80it7LE1NnjyZv/3tb+zfv5+SkpLkd8UMHTqUs88+u9m8BQUFTJgwgddffx2AzZs3853vfIeamhoAampqkrd1cXe2bNnC2LFjeeCBB/j+97+fvEHlyJEj+d73vsf999/fop79+/c3+6Kv0aNHA9DQ0MCPfvQjKisrqays5OGHH07O8/DDD/PlL3+ZyspK/vrXvwJw+PBhvvnNb1JZWcnll1/O1q1b250uIhJKrw6WdN42f/PmzVx44YWYGTt27GDcuHHNljN+/Hh27NjRYvkLFy5k5syZXHPNNTz44IPJ+4ctXbqUPXv28PLLL7Nu3Tq+8Y1vJF8zdOhQXn75Za6//vrkN0b+9Kc/ZcyYMaxbt4477riDBQsWtDtdRCSUXhksmbht/po1a5gxYwZAq3cxbuvOxt/+9repqalh1qxZvPLKK8yYMYNPPvmEtWvXcv3111NU1Lj3csiQIcnXfO1rjbv/mt5if/369cyZMweAyy67jMOHD1NXV9fmdJF8U/TiH5r9SO7I+jGWbMjEbfNfeOEFfv3rXwNw/vnn88YbbzQ7nrJ582bOO++8VusrLy9n7ty5zJ07l8mTJ7N9+/Y2gwhI7r5reov91u4BZ2ZtThcRCaVX9lhaE/K2+XV1dTQ0NFBSUgLAokWLuPfee9mzZw8Ae/bs4ec//zmLFi1qsfzVq1cnzzJLJBIcOnSI8vJypk2bxiOPPJIMjsOHD7f7+1RWVrJ8+XIA/vSnP1FSUkJxcXGb00VEQumVPZbWhLxt/tq1a5NfMwwwZswYfvzjH3Pttdfy6aef0qdPH37yk58wZsyYFst+4YUXuO2225K9kDvvvJOysjKuu+46du7cySWXXEJRURFVVVV897vfbfP3uf3221m4cCGVlZV85jOf4aGHHmp3uohIKLptfju6etv8G2+8keuuuy7Zo8l3um1+ONm8bX6+3YSyo+2qveMqvenK+1y8bb56LGlw6swsEZHeSMdYREQkKAWLiIgElVKwmNmZZrbCzHaY2XYzm2xmJWa2yszi0eOQjpckIiI9Xao9ln8Hnnf3LwBjge3AbcAad48Ba6LxjldYUMDx48e7UqtkwfHjxykoUMdWRFLX4cF7MysGvgzMA3D348BxM5sFTI1mWwasBao7Wt7AgQM5evQoH330UdcqzqD6+vpef41HQUEBAwcOzHYZIpJHUjkrbBTwX8CvzGws8DpwE1Dm7vsA3H2fmZWmskIzY9CgQV2tN6Nqa2sZMWJEtssQEckrqQRLETABuNHda8zs30lxt9cp8Xi8K7XlhHyuPZMy1U6J2sKMrCedErWJ5HC8oCGD681M24X8ndrbrobW1rb53MFe9r7N9OdUR9fNpBIs7wLvuntNNL6CxmBJmFl51FspB9r8K2fz4p3uyPaFR/kik+1UdrJnXSAZi2XuAslMtV2o36nDCyTfbXl38FNKetH7Nhc/pzo8Kuvu+4F3zOzUHROnA38BngGqomlVwNNpqVBERPJKqlfe3wg8amZ9gd3A/6AxlJab2XxgLzA7PSWKiEg+SSlY3H0TUNHKU9PDliMiIvlOFyiIiEhQugmliOQFfUtk/lCPRUREglKwiIhIUNoVJiI93um70XrTF4Flg3osIiISlIJFRESCUrCIiEhQChYREQlKwSIiIkEpWEREJCgFi4iIBKVgERGRoHSBpIgEtfStMF8olqgtbPblZKMTnzR7fnJZvyDrkfDUYxERkaAULCIiEpR2hYlIj6Nb7GeXeiwiIhKUgkVERIJSsIiISFAKFhERCSqlg/dm9jZwBGgATrh7hZmVAE8AI4G3gW+4++H0lCkiIvmiMz2Wr7j7OHeviMZvA9a4ewxYE42LiEgv151dYbOAZdHwMuCq7pcjIiL5zty945nM/gYcBhx42N2XmNkH7n5mk3kOu/uQU+N1dXXJBcfj8bBVS6/1+/2F2S4hqH8+uyFj68r3tqvYsbbZ+Pjijj+72nJwwpe7WU3vFovFksODBw+2059P9QLJS9z9fTMrBVaZ2Y6uFpFP4vF43taeSZlsp6b3jspHidoEZaVlyfFY7IyMrTvf2u70tip+r7jZ86WlXb9XWEkPel/n4udUSrvC3P396LEWeAqYCCTMrBwgeqxNV5EiIpI/OgwWMzvDzAadGgb+EdgKPANURbNVAU+nq0gREckfqewKKwOeMrNT8z/m7s+b2WvAcjObD+wFZqevTBERyRcdBou77wbGtjL9IDA9HUWJiEj+0pX3IiISlIJFRESCUrCIiEhQChYREQlKwSIiIkEpWEREJCh9572I5KzRbzyX7RKkC9RjERGRoBQsIiISlIJFRESCUrCIiEhQChYREQlKwSIiIkEpWEREJCgFi4iIBKVgERGRoBQsIiISlIJFRESC0r3CRKTXKXrxD83GT3zla1mqpGdSj0VERIJKOVjMrNDM3jCzZ6Pxc8ysxsziZvaEmfVNX5kiIpIvOrMr7CZgO1AcjS8G7nX3x83s/wLzgYcC1yfdsPStYxlZT6K2kOpYRlYlInkgpR6LmQ0HrgT+Ixo3YBqwIpplGXBVOgoUEZH8kuqusPuAW4GT0fhQ4AN3PxGNvwsMC1ybiIjkoQ53hZnZfwdq3f11M5t6anIrs3pby4jH412rLgfkc+2J2sKMrWvxn3ZnbF35LlGbSA7HCxoyuN7MbQ+h1B+pb/O5WmvzI6fTDubx+xwy/zkVi7W/7zuVYyyXADPN7J+A/jQeY7kPONPMiqJey3Dg/a4Wkavi8Xje1g5QdjJTx1gSlJWWZWRd+e70torFzsjYujO1PYSSqE1QPKi4zedLS/sFW1dJHr/Pc/FzqsNdYe5+u7sPd/eRwBzgBXf/FvAicE00WxXwdNqqFBGRvNGd61iqge+b2U4aj7n8MkxJIiKSzzp15b27rwXWRsO7gYnhSxIRkXymK+9FRCQo3StMJIsydRGrSCapxyIiIkEpWEREJCjtChORvPSfiU+CLWtbB7sk552XueuNegL1WEREJCgFi4iIBKVgERGRoBQsIiISlIJFRESCUrCIiEhQChYREQlKwSIiIkEpWEREJCgFi4iIBKVgERGRoBQsIiISlIJFRESCUrCIiEhQChYREQlKwSIiIkF1GCxm1t/MNpjZZjPbZmY/iaafY2Y1ZhY3syfMrG/6yxURkVyXSo/lE2Cau48FxgFXmNkkYDFwr7vHgMPA/PSVKSIi+aLDYPFGR6PRPtGPA9OAFdH0ZcBVaalQRETySkrHWMys0Mw2AbXAKmAX8IG7n4hmeRcYlp4SRUQknxSlMpO7NwDjzOxM4Cng/NZma+v18Xi8a9XlgHyuPVFbmMF1JTK2rnyntkpd/ZH6jKxnxMtPJIf//IWpLZ6PFzRkpI6uyvTnVCwWa/f5lILlFHf/wMzWApOAM82sKOq1DAfe72oRuSoej+dt7QBlJ49lZD2J2gRlpWUZWVe+U1ulLlGboHhQccbX29rfJxY7I+N1pCoXP6dSOSvsH6KeCmY2ALgc2A68CFwTzVYFPJ2uIkVEJH+k0mMpB5aZWSGNQbTc3Z81s78Aj5vZ/wHeAH6Zxjp7lKVvZaYnISKSDR0Gi7u/CYxvZfpuYGI6ihIRkfylK+9FRCQoBYuIiASlYBERkaA6dbqxiEhvlKkTbuadl7unNXeGeiwiIhKUgkVERIJSsIiISFAKFhERCUrBIiIiQemsMBHJGaPfeC45POJIPWThJpTSfeqxiIhIUAoWEREJSsEiIiJBKVhERCQoBYuIiASlYBERkaAULCIiEpSCRUREglKwiIhIUAoWEREJqsNgMbMRZvaimW03s21mdlM0vcTMVplZPHockv5yRUQk16XSYzkB/MDdzwcmAf/LzC4AbgPWuHsMWBONi4hIL9dhsLj7PnffGA0fAbYDw4BZwLJotmXAVekqUkRE8kenjrGY2UhgPFADlLn7PmgMH6A0dHEiIpJ/Ur5tvpkNBJ4E/tXd680s5ZXE4/EulJYb0lF7orYw+DKzLVGbyHYJeUNt1bYRR+qbjdefNp4J2fz7xAsauva6DH/GxmKxdp9PKVjMrA+NofKou/8+mpwws3J332dm5UBtV4vIVfF4PC21l508FnyZ2ZSoTVBWWpbtMvKC2qp9xe/9/ftX6o/UU5yF72PJ5t8nFjuj069J1+dUd6RyVpgBvwS2u/vPmzz1DFAVDVcBT4cvT0RE8k0qPZZLgLnAFjPbFE37IfAzYLmZzQf2ArPTU6KI9FRNvzFSeo4Og8XdXwHaOqAyPWw5IiKS73TlvYiIBKVgERGRoBQsIiISlIJFRESCSvkCSRERSa+lb3X+GrdEbWGnr42bd17nr5fpDPVYREQkKAWLiIgEpWAREZGgFCwiIhKUgkVERIJSsIiISFAKFhERCUrBIiIiQSlYREQkKAWLiIgEpWAREZGgFCwiIhKUgkVERIJSsIiISFAKFhERCUrBIiIiQXUYLGb2iJnVmtnWJtNKzGyVmcWjxyHpLVNERPJFKj2WpcAVp027DVjj7jFgTTQuIiLScbC4+8vAodMmzwKWRcPLgKsC1yUiInmqq8dYytx9H0D0WBquJBERyWdFmVhJPB7PxGrSIh21J2oLgy8z2xK1iWyXkDd6e1tV7FibHK7vYN76Ix3NEV4+/n06W3O8oKFb64vFYu0+39VgSZhZubvvM7NyoLY7ReSqeDyeltrLTh4LvsxsStQmKCsty3YZeUFtBcXvFac0X/2ReooHpTZvSPn29+nKNhWLnZGmahp1dVfYM0BVNFwFPB2mHBERyXepnG78W+A/gfPM7F0zmw/8DJhhZnFgRjQuIiLS8a4wd7+2jaemB65FRER6AF15LyIiQWXkrDARkXw1+o3nmo1vG//VLFWSP9RjERGRoBQsIiISlHaFNbH0rebXlyRqC3vcNSci0r7Td31J56nHIiIiQSlYREQkKAWLiIgEpWAREZGgFCwiIhJUzp8VdvqZWiIi2dT0rDFdLNk69VhERCSonO+xiIjkKt3upXXqsYiISFAKFhERCUq7wkQkKN0SRdRjERGRoBQsIiISlIJFRESCUrCIiEhQChYREQmqW8FiZleY2VtmttPMbgtVlIiI5K8uB4uZFQIPAl8FLgCuNbMLQhUmIiL5ydy9ay80mwz82N3/WzR+O4C73wVQV1fXtQWLiEjeGDx4sJ0+rTu7woYB7zQZfzeaJiIivVh3gqVFSgHqpYiI9HLduaXLu8CIJuPDgfdPjbTWPRIRkZ6vOz2W14CYmZ1jZn2BOcAzYcoSEZF81eVgcfcTwCJgJbAdWO7u20IVlg1mNtvMtpnZSTOraGe+t81si5ltMrM/Z7LGXNCJdur1p6ObWYmZrTKzePQ4pI35GqLtaZOZ9Zp/0DraRsysn5k9ET1fY2YjM19lbkihreaZ2X812Y6+k406oRtnhfVEZnY+cBJ4GLjZ3VsNDTN7G6hw9wMZLC9npNJO0enofwVm0Ljb9DXgWnf/SyZrzTYz+zfgkLv/LPowGOLu1a3Md9TdB2a+wuxJZRsxs4XAGHdfYGZzgKvd/V+yUnAWpdhW82j8XFqUlSKb0JX3Tbj7dnd/K9t15LoU22kisNPdd7v7ceBxYFb6q8s5s4Bl0fAy4Kos1pJrUtlGmrbfCmC6mfXG47d59X5SsHSNA//PzF43s+9mu5gcpdPRG5W5+z6A6LG0jfn6m9mfzWy9mfWW8EllG0nOE+1+rwOGZqS63JLq++nrZvamma0wsxGtPJ8Rve6LvsxsNXB2K0/9b3d/OsXFXOLu75tZKbDKzHa4+8vhqsy+AO3Ua05Hb6+tOrGYz0Xb1CjgBTPb4u67wlSYs1LZRnrNdtSBVNrhD8Bv3f0TM1tAY09vWtora0WvCxZ3vzzAMt6PHmvN7Ckau6k9KlgCtFO7p6P3JO21lZklzKzc3feZWTlQ28YyTm1Tu81sLTAe6OnBkso2cmqed82sCBgMHMpMeTmlw7Zy94NNRn8BLM5AXa3SrrBOMrMzzGzQqWHgH4Gt2a0qJ+l09EbPAFXRcBXQordnZkPMrF80fBZwCdAbTnJIZRtp2n7XAC947zzjqMO2iv5xOWUmjWfrZoe76yf6Aa6m8T+DT4AEsDKa/lngj9HwKGBz9LONxl1DWa8919opGv8nGs9k2dUb2ylqg6HAGiAePZZE0yuA/4iGK4Et0Ta1BZif7boz2D4tthHgTmBmNNwf+B2wE9gAjMp2zTncVndFn0mbgReBL2SrVp1uLCIiQWlXmIiIBKVgERGRoBQsIiISlIJFRESCUrCIiEhQChYREQlKwSIiIkEpWEREJKj/DwEhud4ZpKbvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(ate_w_f, kde=False, label=\"PS W School\")\n",
    "sns.distplot(ate_wo_f, kde=False, label=\"PS W/O School\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the propensity score estimator that adds the feature school has a humongous variance, while the one without it is much more well behaved. Also, since school is not a confounder, the model without it is also not biased.\n",
    "\n",
    "As I've said, simply predicting the treatment is not what this is about. We actually need to construct the prediction in a way that controls for confounding, not in a way to predict the treatment. This leads to another problem often encountered in propensity score methods. In our mindset case, the data turned out to be very balanced. But this is not always the case. In some situations, the treated have a much higher probability of treatment than the untreated and the propensity score distribution doesn't overlap much."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEGCAYAAABGnrPVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZRV5Znv8e+PKhkMk2BXiUKrnVUmTm2itIKwiENHxbSQqzGLDC26sNNL1BiHBNKuXL3J9YrGXENnZTBqlKitkHQS1BsHghDwkqAGbON4S42tJaROBCkgAQHruX/sXdVFUdQ5VbVPnYHfZ61a5+y93733U2+ds59633cPigjMzMyyMqDUAZiZWXVxYjEzs0w5sZiZWaacWMzMLFNOLGZmliknFjMzy5QTi1UNSYdJCkmT85S7UNKuHm57uaQ7+hZh36S/2+eLvI+C6tCsO04s1m8k3Z0etELSLkn/KekHkkZntIu3gDHA6nR/Y9N9ndKp3ELgkB5u+1zgqrYJSb+SdHfvQ92dpKMl3SOpSdJ7ad38XNKpWe3DrL84sVh/W0ly8D8M+CJwHvDjLDYcEe9HxB8jYmeectsiormH294YEZv7FmHXJJ0JPAMcDPwTcBRwDvBb4LZi7NOsmJxYrL/tSA/+TRGxGPg2cJakIUpcI+l1STskvSbpSx1XljRd0lpJf5G0SdJTkj6aLuvcjfNW+rosnf9GWq69K0zS8HRbn+20nzGS3pd0Vjrd3hWWtlROB2Z2aIGdIunXkn7YaTtKf4/ru6oMSfuTJNblEXF6RDwSEa9FxHMRcRNwUqdVhqctmy2S3pL0lU7bq5V0vaQ/SNou6QVJ/9ypzFBJ307Xf0/SG5L+pav40vJzJL3bRcvPrEtOLFZq20g+h7XAbOAbwDzgaOCbwDxJswAkHQT8BLg/XT6RJDHtbbzk+PT1PJJW0t91LpC2QhYDMzst+hzQDCzpYrtXkLS8FqXbHQOsAn4AfEbS0A5lTyNpnf1oLzGeAdQBN3S1MCLe7TTrOmAF8BGS+rmpU3fZHSTddv8MHAl8PS3TVocCHgamAZenZS4A/tR535IGSPpO+vt+LCKW7+V3MNtNbakDsH2XpKOAS4HVEbFF0lzgOxHR9l9/o6QPAdcCd5IcwPcDFkXEG2mZl7rZRdvBcmNE/LGbcguAhyUdHBHr0nn/CNwXEe93LhwRLZJ2ANs6blfSvwP/CswgOcADXAw8GhFv7mXfR6SvL3YTX0cLI+L29P2/SppNkpyWSTqcJEkcFREvp2X+kNbh5SR1eBrwMeDvIuKZtMzrJMmqo0EkY1HHAhMj4j8LjM/MLRbrd6dI2ippG/A8yUHts5KGA2PZ8wD3a+CwtMvoOeAx4Pl0YPsKSeMyiGkJkCNppSDpOOBv6eHYT0TsAO4mGSchPSnhvwG3d7Oaehjrs52m3wbq0/fj0+09k9bxVklbgX8BGtIyJwDvdkgqe3MXSVKZ5KRiPeXEYv1tNUk3zpHAkIj4eES83mF559tttx9409bDVJL/up8m6eL6f5L+oS8Bpdu9j+S/fdLXtRHx+15s7jbg7yT9LUmrZyNJ19PevJK+HlXg9nd0mg7+63vc9noySR23/RxDkig7rpPP/wEOB84qMC6zdk4s1t+2RcSrEfFGRLzXNjMd62gi6abpaArwh4j4S1ouIuKpiPhfETGFpEVz0V721XYQrikgrgXAMZLGA59Jp7uzo6vtRsSrwBMkrZaLgbsiortrZh4naS1d29VCSQfkD73d79LXv07ruOPPax3KjEp/z+7cRzLu9CNJncefzLrlMRYrJzcC35LUCCwnaZlcQjIOg6STSc7GehxYT9K987ckYwddeQfYCpwh6QXgvS4GwwGIiOclrSXptvorkhMEuvMH4FRJHwRagJYOpznfBtxLMh50TncbiYi/SLoQ+LmkXwHfImnFfAA4kyRBfShPLG3belXSj4Db07PFfpNu5wTgr9KzzJ4gOfFgoaSrSLoXDwaOjIg7Om3vAUk7gfskDewwtmPWLbdYrJx8H/jvJGMCLwJzgLkR0ZY4WkjOBFsMNJKcaXUfyZlke4iIVpKk9GmSU4/X5tn/ApKuo0cjIpen7LdIEtd/kJwkMKnDsl+ksS6JiD/k2Q4R8QjJGWvNJEnyZZKuqJNJx2t64AvArSQtoBeBpSQtj9fTfQXwCeCXJGexvUKSBA/cS2z/TlJ/bScKmOUlP0HSLFuSRpEMqn8+PTCb7VPcFWaWEUn7kZyh9TVgHUnLxWyf48Rilp1JwDKS8ZcLuroGxmxf4K4wMzPLlAfvzcwsU0XrCmtpaXFTyMysyo0YMWKPu0e4xWJmZplyYjEzs0yVfWJpbGwsdQhlx3WyO9fHnlwnu3N97K7Y9VH2icXMzCqLE4uZmWXKF0ha2YgItm7dSmtra4/WGzx4MC0tLUWKqjL1tE4GDBjA0KFDSR4wadY3TixWNrZu3cqgQYMYOHBgj9YbNGgQgwcPLlJUlamndbJjxw62bt3KsGHDihiV7SvcFWZlo7W1tcdJxbIxcODAHrcUzfbGicXMzDLlxGLWwciRI7n22v96mON3vvMdbrzxxj5tc+PGjUyePJnJkydzxBFHcOSRR7ZP79jR+UnDhVu5ciWrV6/u8XrHHnssGzZs6PV+zfIp+zGW0WtWUNv0clG2vevUbh/uZyV29yt/Lqjczp072W+//DcSvvBDH8hbZtCgQTz00ENcddVVjB49uqD95zNq1CiefPJJAG688UaGDh3K5ZdfvluZXbt2UVvbs6/jk08+ydChQznppJMyidMsK2WfWMz6U21tLRdeeCHf+973+NrXvrbbsjfffJPLLruMd955hwMPPJDvfve7jBs3jksuuYRhw4bx7LPP0tzczNe//nWmT5+ed1+XXHIJBxxwAM899xzHHXccF198Mddccw3vvPMO+++/P/Pnz+eII47gkUce4ZZbbmHHjh2MGjWK22+/nW3btnHXXXdRU1PDwoULufnmmzniiCO48soraWpqorW1lZtuuokJEyawceNGZs2axTvvvMMJJ5yA72huxVZQV5ikKyW9IOl5SfdLGizpcEmrJTVKWijJo65WFS6++GIWLVq0x+m6X/7yl5kxYwarVq3i/PPPZ86cOe3LmpubefTRR1m4cCHXX399wft67bXXWLx4MTfccANXXHEFN998M7/+9a/5xje+wdVXXw3AxIkT+dWvfsXKlSs577zzmD9/PoceeigXXXQRs2fP5sknn+Tkk09m7ty5zJ49m2XLlnHnnXfyxS9+EYB58+YxYcIEVq5cydSpU2lqaup7JZl1I2+LRdIhwBeBoyJim6RFwAzgbODWiHhA0g+AWSTPLDeraMOHD2fGjBncdtttDBkypH3+008/zb333gvAjBkzuO6669qXfeITn2DAgAF8+MMf5k9/+lPB+5o+fTo1NTVs3bqVp556ipkzZ7Yvaxt/efvtt7noootobm5mx44dHHrooV1ua/ny5bz8ctJtHBFs2bKFLVu2sGrVqva4zzzzTEaOHFlwfJaNQrt18ymkO7ccFNoVVgsMkbQT2B9YD5wGfDZdvgC4HicWqxKzZ89mypQpfO5zn9trmY4XEw4aNKj9fU+6mj7wgeRA0drayogRI9rHYjr6yle+wqWXXsrZZ5/NypUrmTdvXpfbam1tZcmSJQwZMoTt27f72h4rmbyJJSLelnQL8CawDXgc+B2wKSJ2pcWagEP2to2+3PBsNJDL5Xq9fnc2FOlGbKPXrCjKdjccP6X9fTXeVG/w4MG7HaB37txZ8LqFlN2+fXveMhHB9u3bGTJkCOeccw733HMPM2bMYPv27YwfP54HHniA888/nwceeIATTzyR7du38/7777Njx4727bdtoyu7du1i586de6w3cOBAxo0bx6JFi5g2bRoRwYsvvsjRRx9NS0sLo0ePZvv27dx77720tra2J4533323fV8f+9jH+N73vsell14KwDPPPMMxxxzDSSedxP3338+VV17J0qVL2bRpE9u3b98jxs2bNxftu1YOSvmdac7VZLKdmzL685x7UN/qo6GhodvlhXSFHQBMBw4HNgE/AaZ2UXSv/6blC6I7G9esoK6urtfrd2dUH+LqTrHOYmuLt7GxsU91Wq5aWlp2+y+7kDO9oO2ssP3ylivkP3hJ7eW+9KUvcdddd1FbW8vgwYP55je/yWWXXcb3v//99sH7wYMHU1NTw8CBA9vX67iNzmpra9lvv/26XO+OO+7g6quvZv78+ezatYtzzz2XE044ga9+9at84Qtf4OCDD2b8+PE0NTUxePBgzjnnHC644AIef/xxbr75Zm655RauueYaTjvtNHbu3MnkyZO59dZbufbaa5k1axZnnHEGkyZNYuzYsQwePHiPGIcPH864cePy1lElKvV3pr41m66w7Kwran3kfea9pPOBsyJiVjp9ATAROB84KCJ2SZoIXB8RZ7atl9UTJDcuvLNoiaVYpxvXLnuoKNtti7fUX5JiaWlpYcSIET1ez90+e+pNnfS2/itBqb8zWY2xZGXSgOwSS2+fIPkmMEHS/ko6lU8HXgSWAZ9Ky8wEFmcSpZmZVbS8iSUiVgM/BdYAv0/X+SEwB7hK0qskQyF3FjFOMzOrEAWdFRYR1wHXdZr9OnBi5hGZmVlF873CzMwsU04sZmaWKd8rzMxsL8rtbK5K4cRiltq4cSPTpk0Dkotya2pq2u9w/MQTT/TqIWQrV65k4MCBPb4D8bHHHsvy5cszu8OyWX9yYrGyVej1QIN27qS2gAsk8123lO/29r61vVlhnFjMulHMW9tDksB8a3urNk4sZnm03dq+pqaGadOmceutt/LBD36QZ555hquvvpqHHnqo/db2kvjxj3/M/PnzueGGG7jooot2a/lcfPHFzJ49m4kTJ/LWW29x3nnn8dRTT7Xf2n7OnDk89thj3H333aX9pc36wInFLI9i3doe8K3trSo5sZjlUaxb25tVK1/HYlag4cOHc+ihh/KLX/wCSG6P//vf/x5Ibjl/8MEHA3D//fe3rzN06FC2bNnSPn3aaadx++23t08/99xzAJx88sn85Cc/AWDJkiVs2rSpuL+MWRE5sZj1wA9/+EPuueceJk2axIQJE/jlL38JwNy5c5k5cyZTp07d7RThqVOn8vDDDzN58mRWrVrFTTfdxNq1azn55JM56aSTuOuuu9rXX7VqFVOmTOGJJ55g7NixJfn9zLKQ97b5veXb5mfPt83vmm+bvyffNn93vf3OVOsFkuVw23wzM7OCObGYmVmmnFjMzCxTTixmZpapvIlF0ockPdvhZ7OkL0kaJWmJpMb09YD+CNiq14ABA9ovOLT+tWPHDgYM8P+Zlo28F0hGxCvARwAk1QBvAz8H5gJLI2KepLnp9Jwixpq5Yp29Zb0zdOhQtm7dyrZt23q03ubNmxk+fHiRoqpMPa2TAQMGMHTo0CJGZPuSnl55fzrwWkT8p6TpwCnp/AXAciossVh5kcSwYcN6vF4ul2PcuHFFiKhyuU6slHra9p0BtF1WXB8R6wHS1+JcbGJmZhWl4AskJQ0E1gFHR0SzpE0RMbLD8ncjon2cpeMFko2Njb0OcPSaFb1et9psOH5KqUMw26f87I81pQ6hKM496P0+rd/x4squLpDsSVfYVGBNRDSn082SxkTEekljgFwhQfTUxjUrinblfaUZldZjtV5531uujz25TnbX2/qob63OK+8huyvvu9KTrrDP8F/dYAAPAm33D58JLM4qKDMzq1wFJRZJ+wMfB37WYfY84OOSGtNlXd8n3MzM9ikFdYVFxF+A0Z3mbSA5S8zMzKydr4gyM7NMObGYmVmmnFjMzCxTTixmZpYpJxYzM8tUT+8VZmZW9jo/Urg5V1PFFzuWH7dYzMwsU04sZmaWKScWMzPLlBOLmZllyonFzMwy5cRiZmaZcmIxM7NMObGYmVmmnFjMzCxTTixmZpapQp8gOVLSTyW9LOklSRMljZK0RFJj+npAsYM1M7PyV2iLZT7waER8GDgOeAmYCyyNiAZgaTptZmb7uLyJRdJwYApwJ0BE7IiITcB0YEFabAHwyWIFaWZmlaOQFsvfAH8C7pK0VtIdkj4A1EfEeoD0ta6IcZqZWYVQRHRfQBoP/BaYFBGrJc0HNgOXR8TIDuXejYj2cZaWlpb2DTc2NvY6wNFrVvR63Wqz4fgppQ7BrCL87I81pQ6hrJ170Pt9Wr+hoaH9/YgRI9R5eSHPY2kCmiJidTr9U5LxlGZJYyJivaQxQK6QIHpq45oV1NW5MQQwKq3HxsbGPtVptXF97Glfr5POz15pzjVTX1dfomjK0bqifj7yJpaI+KOktyR9KCJeAU4HXkx/ZgLz0tfFRYvSAKhd9hAAo3M5aptezmy7u049J7NtmZkV+gTJy4H7JA0EXgcuIhmfWSRpFvAmcH5xQjQzs0pSUGKJiGeB8V0sOj3bcMzMrNL5ynszM8uUE4uZmWXKicXMzDLlxGJmZplyYjEzs0w5sZiZWaacWMzMLFNOLGZmliknFjMzy5QTi5mZZarQe4VZFWu7uWXWfHNLs32TWyxmZpYpJxYzM8uUE4uZmWXKicXMzDLlxGJmZpkq6KwwSW8AW4D3gV0RMV7SKGAhcBjwBvDpiHi3OGGamVml6EmL5dSI+EhEtD1Jci6wNCIagKXptJmZ7eP60hU2HViQvl8AfLLv4ZiZWaUrNLEE8Lik30n6QjqvPiLWA6SvdcUI0MzMKkuhV95Pioh1kuqAJZJe7slOGhsbex5ZajSQy+V6vX61qoQ62dCHv3tP9eUzVq325TppztV0Ma+5BJGUqYP69vloaGjodnlBiSUi1qWvOUk/B04EmiWNiYj1ksYAez3S5QuiOxvXrKCuzo2hjnK5XEXUyag+/N17orGxsU+fsWq0r9dJfeufd5tuzjVTX1dfomjK0bqifj7ydoVJ+oCkYW3vgTOA54EHgZlpsZnA4mIFaWZmlaOQFks98HNJbeX/LSIelfQ0sEjSLOBN4PzihWlmZpUib2KJiNeB47qYvwE4vRhBmZlZ5fKV92ZmliknFjMzy5QTi5mZZcqJxczMMuXEYmZmmXJiMTOzTDmxmJlZppxYzMwsU04sZmaWKScWMzPLlBOLmZllyonFzMwy5cRiZmaZcmIxM7NMObGYmVmmCn3mvZnZXt39yp/zF7J9hlssZmaWqYITi6QaSWslPZxOHy5ptaRGSQslDSxemGZmVil60mK5Anipw/RNwK0R0QC8C8zKMjAzM6tMBSUWSWOBTwB3pNMCTgN+mhZZAHyyGAGamVllKXTw/tvAV4Bh6fRoYFNE7Eqnm4BD9rZyY2NjrwMcDeRyuV6vX60qoU429OHv3p3Ra1bsPg1s7DSvtzYcPyWT7ZSDvnzveqo5V9Nv++qt5lxzqUMoHwf17fPR0NDQ7fK8iUXSPwC5iPidpFPaZndRNHobRHc2rllBXV1dr9evRrlcriLqZFQf/u7dqW16ebfpLOujWDH3t8bGxj5973qqvrW8zwprzjVTX1df6jDKyLqifj4KabFMAqZJOhsYDAwnacGMlFSbtlrGAuuKFqWZmVWMvGMsEfHViBgbEYcBM4AnIuJzwDLgU2mxmcDiokVpZmYVoy/XscwBrpL0Kkk3953ZhGRmZpWsR1feR8RyYHn6/nXgxOxDMjOzSuYr783MLFNOLGZmliknFjMzy5QTi5mZZcqJxczMMuXEYmZmmXJiMTOzTPkJkmb7MD/50YrBLRYzM8uUWyxWNLXLHip1CGZWAm6xmJlZppxYzMwsU04sZmaWKScWMzPLlBOLmZllyonFzMwylTexSBos6SlJ/yHpBUn/I51/uKTVkholLZQ0sPjhmplZuSukxfIecFpEHAd8BDhL0gTgJuDWiGgA3gVmFS9MMzOrFHkTSyS2ppP7pT8BnAb8NJ2/APhkUSI0M7OKUtAYi6QaSc8COWAJ8BqwKSJ2pUWagEOKE6KZmVWSgm7pEhHvAx+RNBL4OXBkV8X2tn5jY2PvogNGA7lcrtfrVyvXye6yqo8NffislptCvnfNuZp+iKQ8NOeaSx1C+Tiob8flhoaGbpf36F5hEbFJ0nJgAjBSUm3aahkLrOttEN3ZuGYFdXV1vV6/GuVyOddJB1nWx6g+fFbLSWNjY0Hfu/rWfePuxs25Zurr6ksdRhlZ16fjcj6FnBX2V2lLBUlDgL8HXgKWAZ9Ki80EFhcrSDMzqxyFtFjGAAsk1ZAkokUR8bCkF4EHJP1PYC1wZxHjNDOzCpE3sUTEc8BHu5j/OnBiMYIyM7PK5SvvzcwsU04sZmaWKScWMzPLlBOLmZllyonFzMwy1aMLJM2qXe2yh4qy3V2nnlOU7ZqVI7dYzMwsU04sZmaWKScWMzPLlBOLmZllyonFzMwy5cRiZmaZcmIxM7NMObGYmVmmnFjMzCxTTixmZpapQh5NPE7SMkkvSXpB0hXp/FGSlkhqTF8PKH64ZmZW7gq5V9gu4OqIWCNpGPA7SUuAC4GlETFP0lxgLjCneKGaGcDdr/w5b5nmXA31rfnLmRVD3hZLRKyPiDXp+y3AS8AhwHRgQVpsAfDJYgVpZmaVo0djLJIOAz4KrAbqI2I9JMkHqMs6ODMzqzwF3zZf0lDg34EvRcRmSQXvpLGxsRehJUYDuVyu1+tXK9fJ7sq9Pjb04TvQWXOupsByzZntsxq4Pjo4qG/H5YaGhm6XF5RYJO1HklTui4ifpbObJY2JiPWSxgB7/WbnC6I7G9esoK7OjaGOcrmc66SDSqiPUX34DnRWyNhJc66Z+rr6zPZZ6Vwfna3r03E5n7yJRUnT5E7gpYj43x0WPQjMBOalr4uLEqFZFcjyAWJHN7/X/v6Fj07NbLtmWSmkxTIJ+Efg95KeTef9C0lCWSRpFvAmcH5xQjQzs0qSN7FExJPA3gZUTs82HDMzq3S+8t7MzDLlxGJmZplyYjEzs0w5sZiZWaacWMzMLFNOLGZmliknFjMzy5QTi5mZZcqJxczMMuXEYmZmmXJiMTOzTDmxmJlZppxYzMwsU04sZmaWKScWMzPLlBOLmZllKm9ikfQjSTlJz3eYN0rSEkmN6esBxQ3TzMwqRSEtlruBszrNmwssjYgGYGk6bWZmVtCjiVdIOqzT7OnAKen7BcByYE6GcVkF+k3zeyXZ7+YtYnjsue+J9YNKEE3/OnrtI13OH7dlM8PfHt7r7b7w0am9Xtest2Ms9RGxHiB9rcsuJDMzq2R5WyxZaGxs7PW6o4FcLpddMFWiHOtk8xaVcN+b95j32JYSBFJGuqqTQjXnmjOMpDxU4+/Uawf17bjc0NDQ7fLeJpZmSWMiYr2kMUC3R7l8QXRn45oV1NW5QdRRLpcryzrpqjuqP2zespnhw3rf7VON+lon9XX1GUZTes255qr7nfpmXZ+Oy/n0tivsQWBm+n4msDibcMzMrNLlbbFIup9koP5ASU3AdcA8YJGkWcCbwPnFDNLM+tfeTgroK58UsG8o5Kywz+xl0ekZx2JmZlXAV96bmVmmnFjMzCxTTixmZpapfrmOxcwMindSAPjEgHLiFouZmWXKicXMzDLlxGJmZplyYjEzs0x58N7MqkJ3Jwb05TECPimg59xiMTOzTDmxmJlZppxYzMwsU04sZmaWKQ/em5l1w48Q6Dm3WMzMLFNusRi/aS7NI4XNrDq5xWJmZpnqU2KRdJakVyS9KmluVkGZmVnl6nVXmKQa4LvAx4Em4GlJD0bEi1kFV6362vW0eYsYHu6+MrPypIjo3YrSROD6iDgznf4qQETcCNDS0tK7DZuZWcUYMWKEOs/rS1fYIcBbHaab0nlmZrYP60ti2SNLAW6lmJnt4/pyunETMK7D9FhgXdtEV80jMzOrfn1psTwNNEg6XNJAYAbwYDZhmZlZpep1YomIXcBlwGPAS8CiiHiht9vLd+qypEGSFqbLV0s6rLf7qgQF1MdVkl6U9JykpZIOLUWc/anQ09slfUpSSBrfn/H1t0LqQ9Kn08/JC5L+rb9j7G8FfG/+WtIySWvT787ZpYizv0j6kaScpOf3slyS/jWtr+ckHZ/JjiOi5D9ADfAa8DfAQOA/gKM6lZkN/CB9PwNYWOq4S1wfpwL7p+8vqeb6KLRO0nLDgBXAb4HxpY67xJ+RBmAtcEA6XVfquMugTn4IXJK+Pwp4o9RxF7lOpgDHA8/vZfnZwCMkY+YTgNVZ7Ldcrrw/EXg1Il6PiB3AA8D0TmWmAwvS9z8FTpdUreM4eesjIpZFxF/Syd+SjHFVs0I+IwDfAG4GtvdncCVQSH38E/DdiHgXICJy/RxjfyukTgJoe5TkCDqMC1ejiFgBbOymyHTgx5H4LTBS0pi+7rdcEkshpy63l4mkG64FGN0v0fW/np7KPYvkv45qlrdOJH0UGBcRD/dnYCVSyGfkCOAISf9X0m8lndVv0ZVGIXVyPfB5SU3AL4HL+ye0slWUy0bK5SaUhZy6vC+d3lzw7yrp88B44GNFjaj0uq0TSQOAW4EL+yugEivkM1JL0h12CkmLdqWkYyJiU5FjK5VC6uQzwN0R8a30Iu970jppLX54Zakox9VyabF0e+py5zKSakmasd018SpZIfWBpL8HrgWmRVT9PV7y1ckw4BhguaQ3SPqLH6ziAfxCvzOLI2JnRPwBeIUk0VSrQupkFrAIICJ+AwwGDuyX6MpTQceaniqXxFLIqcsPAjPT958Cnoh09KkK5a2PtNvnNpKkUu1955CnTiKiJSIOjIjDIuIwknGnaRHxTGnCLbpCvjO/IDnJA0kHknSNvd6vUfavQurkTeB0AElHkiSWP/VrlOXlQeCC9OywCUBLRKzv60bLoissInZJajt1uQb4UUS8IOnrwDMR8SBwJ0mz9VWSlsqM0kVcXAXWxzeBocBP0nMY3oyIaSULusgKrJN9RoH18RhwhqQXgfeBL0fEhtJFXVwF1snVwO2SriTp8rmwiv9BRdL9JF2hB6bjStcB+wFExA9IxpnOBl4F/gJclMl+q7hOzcysBMqlK8zMzKqEE4uZmWXKicXMzDLlxGJmZplyYjEzs7ScRV8AAAAZSURBVEw5sZiZWaacWMzMLFNOLGZmlqn/D1HFy3tp6b1mAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(np.random.beta(4,1,500), kde=False, label=\"Non Treated\")\n",
    "sns.distplot(np.random.beta(1,3,500), kde=False, label=\"Treated\")\n",
    "plt.title(\"Positivity Check\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this happens, it means that positivity is not very strong. If a treated has a propensity score of, say, 0.9 and the maximum propensity score of the untreated is 0.7, we won't have any untreated to compare to the individual with the 0.9 propensity score. This lack of balancing can generate some bias, because we will have to extrapolate the treatment effect to unknown regions. Not only that, entities with very high or very low propensity scores have a very high weight, which increases variance. As a general rule of thumb, you are in trouble if any weight is higher than 20 (which happens with an untreated propensity score of 0.95 or a treated with a propensity score of 0.05). \n",
    "\n",
    "An alternative is clipping the weight to be at a maximum size of 20. This will decrease the variance, but it will actually generate more bias. To be honest, although this is a common practice to reduce variance, I don't really like it. You will never know if the bias you are inducing with clipping is too much. Also, if the distributions don't overlap, your data is probably not enough to make a causal conclusion anyway. To gain some further intuition about this, we can look at a technique that combines propensity score and matching\n",
    "\n",
    "## Propensity Score Matching\n",
    "\n",
    "As I've said before, you don't need to control for X when you have the propensity score. It suffices to control it. As such, you can think of the propensity score as performing a kind of dimensionality reduction on the feature space. It condenses all the features in X into a single treatment assignment dimension. For this reason, we can treat the propensity score as an input feature for other models. Take a regression, model for instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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.0759</td> <td>    0.065</td> <td>  -47.046</td> <td> 0.000</td> <td>   -3.204</td> <td>   -2.948</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>intervention</th>     <td>    0.3930</td> <td>    0.019</td> <td>   20.974</td> <td> 0.000</td> <td>    0.356</td> <td>    0.430</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>propensity_score</th> <td>    9.0519</td> <td>    0.200</td> <td>   45.298</td> <td> 0.000</td> <td>    8.660</td> <td>    9.444</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smf.ols(\"achievement_score ~ intervention + propensity_score\", data=data_ps).fit().summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we control for the propensity score, we now estimate a ATE of 0.39, which is lower than the 0.47 we did previously with a regression model without controlling for the propensity score. We can also use matching on the propensity score. This time, instead of trying to find matches that are similar in all the X features, we can find matches that just have the same propensity score.\n",
    "\n",
    "This is a huge improvement on top of the matching estimator, since it deals with the curse of dimensionality. Also, if a feature is unimportant for the treatment assignment, the propensity score model will learn that and give low importance to it when fitting the treatment mechanism. Matching on the features, on the other hand, would still try to find matches where individuals are similar on this unimportant feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Treatment Effect Estimates: Matching\n",
      "\n",
      "                     Est.       S.e.          z      P>|z|      [95% Conf. int.]\n",
      "--------------------------------------------------------------------------------\n",
      "           ATE      0.379      0.025     15.102      0.000      0.330      0.428\n",
      "           ATC      0.371      0.028     13.461      0.000      0.317      0.425\n",
      "           ATT      0.395      0.027     14.844      0.000      0.343      0.447\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from causalinference import CausalModel\n",
    "\n",
    "cm = CausalModel(\n",
    "    Y=data_ps[\"achievement_score\"].values, \n",
    "    D=data_ps[\"intervention\"].values, \n",
    "    X=data_ps[[\"propensity_score\"]].values\n",
    ")\n",
    "\n",
    "cm.est_via_matching(matches=1, bias_adj=True)\n",
    "\n",
    "print(cm.estimates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, we also get an ATE of 0.38, which is more in line with what we've seen before with propensity score weighting. Matching on the propensity score also give us some intuition about why is dangerous to have a small overlap in the propensity score between treated and untreated. If this happens, the matching on the propensity score discrepancy will be large, which will lead to bias. \n",
    "\n",
    "One final word of caution here is that the above standard errors are wrong, as they don't account for the uncertainty in the estimation of the propensity score. Unfortunately, [bootstrap doesn't work with matching](https://economics.mit.edu/files/11862). Also, the theory above is so recent that there are no Python implementations of propensity score methods with the correct standard errors. For this reason, we don't see a lot of propensity score matching in Python. \n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "Here, we've learned that the probability of getting the treatment is called the propensity score and that we can use this as a balancing score. What this means is that, if we have the propensity score, we don't need to control for the confounders directly. It is sufficient to control for the propensity score in order to identify the causal effect. We saw how the propensity scores acts as a dimensionality reduction on the confounder space.\n",
    "\n",
    "These proprieties allowed us to derive a weighting estimator for causal inference. Not only that, we saw how the propensity score can be used along other methods to control for confounding bias. \n",
    "\n",
    "Finally, we looked at some extrapolation problems that we might run into if we are unable to have a good overlap between the treated and untreated propensity score distribution. \n",
    "\n",
    "\n",
    "## References\n",
    "\n",
    "I like to think of this entire series is a tribute to Joshua Angrist, Alberto Abadie and Christopher Walters for their amazing Econometrics class. Most of the ideas here are taken from their classes at the American Economic Association. Watching them is what is keeping me sane during this tough year of 2020.\n",
    "* [Cross-Section Econometrics](https://www.aeaweb.org/conference/cont-ed/2017-webcasts)\n",
    "* [Mastering Mostly Harmless Econometrics](https://www.aeaweb.org/conference/cont-ed/2020-webcasts)\n",
    "\n",
    "I'll also like to reference the amazing books from Angrist. They have shown me that Econometrics, or 'Metrics as they call it, is not only extremely useful but also profoundly fun.\n",
    "\n",
    "* [Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com/)\n",
    "* [Mastering 'Metrics](https://www.masteringmetrics.com/)\n",
    "\n",
    "My final reference is Miguel Hernan and Jamie Robins' book. It has been my trustworthy companion in the most thorny causal questions I had to answer.\n",
    "\n",
    "* [Causal Inference Book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)\n",
    "\n",
    "The data that we used was taken from the article [Estimating Treatment Effects with Causal Forests: An Application](https://arxiv.org/pdf/1902.07409.pdf), by Susan Athey and Stefan Wager. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}