Lalonde dataset - Causal Inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Causal Inference in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is an exploration of causal inference in python using the famous Lalonde dataset.\n",
    "* Causal inference is a technique to estimate the effect of one variable onto another, given the presence of other influencing variables (confonding factors) that we try to keep 'controlled'.\n",
    "\n",
    "The study looked at the effectiveness of a job training program (the treatment) on the real earnings of an individual, a couple years after completion of the program.\n",
    "\n",
    "The data consists of a number of demographic variables (age, race, academic background, and previous real earnings), as well as a treatment indicator, and the real earnings in the year 1978 (the response).\n",
    "\n",
    "*Robert Lalonde, \"Evaluating the Econometric Evaluations of Training Programs\", American Economic Review, Vol. 76, pp. 604-620* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we'll be using the tools provided by Laurence Wong in the Package *CausalInference*. Comments on what each function does come from the very good package documentation: http://laurence-wong.com/software/\n",
    "\n",
    "This package relies heavily on Rubin causal model, and so will this analysis https://en.wikipedia.org/wiki/Rubin_causal_model\n",
    "\n",
    "The reason why several models exist is that it is impossible to observe the causal effect on a single unit, and so assumptions must be made to estimate the missing counterfactuals. We'll explain what all that means in this post."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://pypi.org/project/CausalInference/\n",
    "from causalinference import CausalModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>treat</th>\n",
       "      <th>age</th>\n",
       "      <th>educ</th>\n",
       "      <th>black</th>\n",
       "      <th>hispan</th>\n",
       "      <th>married</th>\n",
       "      <th>nodegree</th>\n",
       "      <th>re74</th>\n",
       "      <th>re75</th>\n",
       "      <th>re78</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>NSW1</th>\n",
       "      <td>1</td>\n",
       "      <td>37</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9930.0460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSW2</th>\n",
       "      <td>1</td>\n",
       "      <td>22</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3595.8940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSW3</th>\n",
       "      <td>1</td>\n",
       "      <td>30</td>\n",
       "      <td>12</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>24909.4500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSW4</th>\n",
       "      <td>1</td>\n",
       "      <td>27</td>\n",
       "      <td>11</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7506.1460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSW5</th>\n",
       "      <td>1</td>\n",
       "      <td>33</td>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>289.7899</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      treat  age  educ  black  hispan  married  nodegree  re74  re75  \\\n",
       "NSW1      1   37    11      1       0        1         1   0.0   0.0   \n",
       "NSW2      1   22     9      0       1        0         1   0.0   0.0   \n",
       "NSW3      1   30    12      1       0        0         0   0.0   0.0   \n",
       "NSW4      1   27    11      1       0        0         1   0.0   0.0   \n",
       "NSW5      1   33     8      1       0        0         1   0.0   0.0   \n",
       "\n",
       "            re78  \n",
       "NSW1   9930.0460  \n",
       "NSW2   3595.8940  \n",
       "NSW3  24909.4500  \n",
       "NSW4   7506.1460  \n",
       "NSW5    289.7899  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lalonde = pd.read_csv('lalonde.csv', index_col=0)\n",
    "lalonde.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>treat</th>\n",
       "      <th>age</th>\n",
       "      <th>educ</th>\n",
       "      <th>black</th>\n",
       "      <th>hispan</th>\n",
       "      <th>married</th>\n",
       "      <th>nodegree</th>\n",
       "      <th>re74</th>\n",
       "      <th>re75</th>\n",
       "      <th>re78</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "      <td>614.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.301303</td>\n",
       "      <td>27.363192</td>\n",
       "      <td>10.268730</td>\n",
       "      <td>0.395765</td>\n",
       "      <td>0.117264</td>\n",
       "      <td>0.415309</td>\n",
       "      <td>0.630293</td>\n",
       "      <td>4557.546569</td>\n",
       "      <td>2184.938207</td>\n",
       "      <td>6792.834483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.459198</td>\n",
       "      <td>9.881187</td>\n",
       "      <td>2.628325</td>\n",
       "      <td>0.489413</td>\n",
       "      <td>0.321997</td>\n",
       "      <td>0.493177</td>\n",
       "      <td>0.483119</td>\n",
       "      <td>6477.964479</td>\n",
       "      <td>3295.679043</td>\n",
       "      <td>7470.730792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>16.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>20.000000</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>238.283425</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>25.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1042.330000</td>\n",
       "      <td>601.548400</td>\n",
       "      <td>4759.018500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>32.000000</td>\n",
       "      <td>12.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>7888.498250</td>\n",
       "      <td>3248.987500</td>\n",
       "      <td>10893.592500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>55.000000</td>\n",
       "      <td>18.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>35040.070000</td>\n",
       "      <td>25142.240000</td>\n",
       "      <td>60307.930000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            treat         age        educ       black      hispan     married  \\\n",
       "count  614.000000  614.000000  614.000000  614.000000  614.000000  614.000000   \n",
       "mean     0.301303   27.363192   10.268730    0.395765    0.117264    0.415309   \n",
       "std      0.459198    9.881187    2.628325    0.489413    0.321997    0.493177   \n",
       "min      0.000000   16.000000    0.000000    0.000000    0.000000    0.000000   \n",
       "25%      0.000000   20.000000    9.000000    0.000000    0.000000    0.000000   \n",
       "50%      0.000000   25.000000   11.000000    0.000000    0.000000    0.000000   \n",
       "75%      1.000000   32.000000   12.000000    1.000000    0.000000    1.000000   \n",
       "max      1.000000   55.000000   18.000000    1.000000    1.000000    1.000000   \n",
       "\n",
       "         nodegree          re74          re75          re78  \n",
       "count  614.000000    614.000000    614.000000    614.000000  \n",
       "mean     0.630293   4557.546569   2184.938207   6792.834483  \n",
       "std      0.483119   6477.964479   3295.679043   7470.730792  \n",
       "min      0.000000      0.000000      0.000000      0.000000  \n",
       "25%      0.000000      0.000000      0.000000    238.283425  \n",
       "50%      1.000000   1042.330000    601.548400   4759.018500  \n",
       "75%      1.000000   7888.498250   3248.987500  10893.592500  \n",
       "max      1.000000  35040.070000  25142.240000  60307.930000  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's have an overview of the data\n",
    "lalonde.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Here is the *raw* difference in earning between the control group and the treated group:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>median</th>\n",
       "      <th>mean</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>treat</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4975.505</td>\n",
       "      <td>6984.169742</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4232.309</td>\n",
       "      <td>6349.143530</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         median         mean\n",
       "treat                       \n",
       "0      4975.505  6984.169742\n",
       "1      4232.309  6349.143530"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lalonde.groupby('treat')['re78'].agg(['median','mean'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The control group has higher earning that the treatment group - does this mean the treatment had a negative impact?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "treat\n",
       "0    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "1    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "Name: re78, dtype: object"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD8CAYAAABthzNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFp1JREFUeJzt3X20XWV94PHvzyQQEWwIJDbmkiapWVZwaUkvDA4OOFIHSDHYmTqTLFcNbw0tGRuHmaWgXUXXGtfS2gG0TDWxUsFaEPEFpIhGKlhXKzHhRYKASYGBS6gJsUgVEYi/+eM8gZPLvveee+/Z5+Xe72ets87ez95nP7/Hu82P53n2S2QmkiQN95JuByBJ6k0mCElSJROEJKmSCUKSVMkEIUmqZIKQJFUyQUiSKpkgJEmVTBCSpEozux3AZBx++OG5ePHibochSX1l69atj2fmvLH26+sEsXjxYrZs2dLtMCSpr0TE/2tlP4eYJEmVTBCSpEomCElSpb6eg5Ckbnn22WcZGhri6aef7nYoI5o9ezYDAwPMmjVrQr+vLUFExOXAacCuzHztsG3/C/goMC8zH4+IAD4GrACeAs7IzNvrik2SJmtoaIhDDjmExYsX0/gnrLdkJnv27GFoaIglS5ZM6Bh1DjF9BjhleGFEHAG8BXi4qfhUYFn5rAU+UWNckjRpTz/9NIcddlhPJgeAiOCwww6bVA+ntgSRmd8Gflyx6RLgPUDzq+xOB67Mhu8CcyJiQV2xSVI79Gpy2Gey8XV0kjoiVgKPZuZdwzYtBB5pWh8qZVXHWBsRWyJiy+7du2uKVJLUsUnqiDgIeD/wn6o2V5RVviw7MzcCGwEGBwd9obaknvDWv/hOW4/31Xe9ccx9brrpJtavX8/evXs555xzuOCCC9oaQyevYvp1YAlwV+n2DAC3R8SxNHoMRzTtOwDsrDOY0f6YrfxhJKmb9u7dy7p169i0aRMDAwMcc8wxrFy5kiOPPLJtdXRsiCkz787M+Zm5ODMX00gKyzPzX4DrgXdGw3HATzLzsU7FJkn9ZvPmzbzqVa9i6dKlHHDAAaxatYrrrruurXXUliAi4irgn4BXR8RQRJw9yu43Ag8AO4BPAefVFZckTQWPPvooRxzxwsDLwMAAjz76aFvrqG2IKTNXj7F9cdNyAuvqikWSpprGP5v7a/dVVT5qQ5L60MDAAI888sLFn0NDQ7zyla9sax0mCEnqQ8cccwzbt2/nwQcf5JlnnuHqq69m5cqVba3DZzFJUht0+urHmTNnctlll3HyySezd+9ezjrrLI466qj21tHWo0mSOmbFihWsWLGituM7xCRJqmSCkCRVMkFIkiqZICRJlUwQkqRKJghJUiUvc5WkdthwYnuPd+6tY+5y1llnccMNNzB//ny2bdvW3vqxByFJfeuMM87gpptuqu34JghJ6lMnnHACc+fOre34JghJUiUThCSpkglCklTJBCFJquRlrpLUDi1cltpuq1ev5pZbbuHxxx9nYGCAD37wg5x99mhvdx4fE4Qk9amrrrqq1uM7xCRJqlRbgoiIyyNiV0Rsayr7aETcFxHfj4gvR8Scpm0XRsSOiLg/Ik6uKy5JUmvq7EF8BjhlWNkm4LWZ+Trgh8CFABFxJLAKOKr85i8jYkaNsUnSpGVmt0MY1WTjqy1BZOa3gR8PK/tGZj5XVr8LDJTl04GrM/MXmfkgsAM4tq7YJGmyZs+ezZ49e3o2SWQme/bsYfbs2RM+Rjcnqc8CPl+WF9JIGPsMlTJJ6kkDAwMMDQ2xe/fubocyotmzZzMwMDD2jiPoSoKIiPcDzwGf21dUsVtlWo6ItcBagEWLFtUSnySNZdasWSxZsqTbYdSq41cxRcQa4DTgHflC32wIOKJptwFgZ9XvM3NjZg5m5uC8efPqDVaSprGOJoiIOAV4L7AyM59q2nQ9sCoiDoyIJcAyYHMnY5Mk7a+2IaaIuAp4E3B4RAwBF9G4aulAYFNEAHw3M/8wM++JiGuAH9AYelqXmXvrik2SNLbaEkRmrq4o/vQo+38I+FBd8UiSxsc7qSVJlUwQkqRKJghJUiUThCSpkglCklTJBCFJqmSCkCRVMkFIkiqZICRJlUwQkqRKJghJUiUThCSpkglCklTJBCFJqmSCkCRVMkFIkiqZICRJlUwQkqRKJghJUiUThCSpkglCklSptgQREZdHxK6I2NZUNjciNkXE9vJ9aCmPiPh4ROyIiO9HxPK64pIktabOHsRngFOGlV0A3JyZy4CbyzrAqcCy8lkLfKLGuCRJLagtQWTmt4EfDys+HbiiLF8BvK2p/Mps+C4wJyIW1BWbJGlsnZ6DeEVmPgZQvueX8oXAI037DZUySVKX9MokdVSUZeWOEWsjYktEbNm9e3fNYUnS9NXpBPGjfUNH5XtXKR8CjmjabwDYWXWAzNyYmYOZOThv3rxag5Wk6azTCeJ6YE1ZXgNc11T+znI103HAT/YNRUmSumNmXQeOiKuANwGHR8QQcBHwYeCaiDgbeBh4e9n9RmAFsAN4CjizrrgkSa2pLUFk5uoRNp1UsW8C6+qKRZI0fr0ySS1J6jEmCElSJROEJKmSCUKSVMkEIUmqZIKQJFUyQUiSKpkgJEmVTBCSpEomCElSJROEJKlSSwkiIl5bdyCSpN7Sag/ikxGxOSLOi4g5tUYkSeoJLSWIzHwj8A4aL/XZEhF/GxFvqTUySVJXtTwHkZnbgT8B3gucCHw8Iu6LiP9cV3CSpO5pdQ7idRFxCXAv8GbgrZn5mrJ8SY3xSZK6pNUXBl0GfAp4X2b+fF9hZu6MiD+pJTJJUle1miBWAD/PzL0AEfESYHZmPpWZn60tOklS17Q6B/FN4KVN6weVMknSFNVqgpidmT/dt1KWD6onJElSL2g1QfwsIpbvW4mI3wJ+Psr+kqQ+1+ocxLuBL0TEzrK+APhvE600Iv4HcA6QwN3AmeWYVwNzgduB38/MZyZahyRpclq9Ue57wG8AfwScB7wmM7dOpMKIWAj8MTCYma8FZgCrgI8Al2TmMuBfgbMncnxJUnuM52F9xwCvA44GVkfEOydR70zgpRExk8ZcxmM07qm4tmy/AnjbJI4vSZqkloaYIuKzwK8DdwJ7S3ECV463wsx8NCL+HHiYxjzGN4CtwBOZ+VzZbQhYON5jS5Lap9U5iEHgyMzMyVYYEYcCpwNLgCeALwCnVuxaWVdErAXWAixatGiy4UiSRtDqENM24FfbVOdvAw9m5u7MfBb4EvDvgTllyAlgANhZ9ePM3JiZg5k5OG/evDaFJEkartUexOHADyJiM/CLfYWZuXICdT4MHBcRB9EYYjoJ2AJ8C/g9GlcyrQGum8CxJUlt0mqC+EC7KszM2yLiWhqXsj4H3AFsBP4OuDoi/ncp+3S76qxy8RPrR9644eCxD3Dure0LRpJ6UEsJIjNvjYhfA5Zl5jfLf/3PmGilmXkRcNGw4geAYyd6TElSe7X6uO8/oHEJ6oZStBD4Sl1BSZK6r9VJ6nXA8cCT8PzLg+bXFZQkqftaTRC/aH7sRbnaaNKXvEqSelerCeLWiHgfjbuf30Lj3oWv1heWJKnbWk0QFwC7aTxY71zgRhrvp5YkTVGtXsX0SxqvHP1UveFIknpFq89iepCKOYfMXNr2iCRJPWE8z2LaZzbwdhrvbZAkTVGtvg9iT9Pn0cy8lMbjuSVJU1SrQ0zLm1ZfQqNHcUgtEUmSekKrQ0z/p2n5OeAh4L+2PRpJUs9o9Sqm/1h3IJKk3tLqENP5o23PzIvbE44kqVeM5yqmY4Dry/pbgW8Dj9QRlCSp+8bzwqDlmflvABHxAeALmXlOXYFJkrqr1UdtLAKeaVp/Bljc9mgkST2j1R7EZ4HNEfFlGndU/y5wZW1RSZK6rtWrmD4UEV8D/kMpOjMz76gvLElSt7U6xARwEPBkZn4MGIqIJTXFJEnqAa2+cvQi4L3AhaVoFvA3dQUlSeq+VnsQvwusBH4GkJk78VEbkjSltZognsnMpDzyOyJeNplKI2JORFwbEfdFxL0R8YaImBsRmyJie/k+dDJ1SJImp9UEcU1EbADmRMQfAN9kci8P+hhwU2b+BvB64F4ab627OTOXATeXdUlSl7R6FdOfl3dRPwm8GvjTzNw0kQoj4uXACcAZ5djPAM9ExOnAm8puVwC30Jj3kCR1wZgJIiJmAF/PzN8GJpQUhllK4/3Wfx0Rrwe2AuuBV2TmYwCZ+VhEzB8hnrXAWoBFixa1IRxJUpUxh5gycy/wVET8SpvqnAksBz6RmUfTmPhueTgpMzdm5mBmDs6bN69NIUmShmv1TuqngbsjYhPlSiaAzPzjCdQ5BAxl5m1l/VoaCeJHEbGg9B4WALsmcGxJUpu0miD+rnwmLTP/JSIeiYhXZ+b9wEnAD8pnDfDh8n1dO+qTJE3MqAkiIhZl5sOZeUWb630X8LmIOAB4ADiTxnDXNRFxNvAw8PY21ylJGoexehBfoTFfQER8MTP/Szsqzcw7abxjYriT2nF8SdLkjTVJHU3LS+sMRJLUW8ZKEDnCsiRpihtriOn1EfEkjZ7ES8syZT0z8+W1RidJ6ppRE0RmzuhUIJKk3jKe90FIkqYRE4QkqZIJQpJUyQQhSapkgpAkVTJBSJIqmSAkSZVMEJKkSiYISVIlE4QkqZIJQpJUyQQhSapkgpAkVTJBSJIqmSAkSZVMEJKkSiYISVKlriWIiJgREXdExA1lfUlE3BYR2yPi8xFxQLdikyR1twexHri3af0jwCWZuQz4V+DsrkQlSQK6lCAiYgD4HeCvynoAbwauLbtcAbytG7FJkhq61YO4FHgP8MuyfhjwRGY+V9aHgIVVP4yItRGxJSK27N69u/5IJWma6niCiIjTgF2ZubW5uGLXrPp9Zm7MzMHMHJw3b14tMUqSYGYX6jweWBkRK4DZwMtp9CjmRMTM0osYAHZ2ITZJUtHxHkRmXpiZA5m5GFgF/H1mvgP4FvB7Zbc1wHWdjk2S9IJeug/ivcD5EbGDxpzEp7scjyRNa90YYnpeZt4C3FKWHwCO7WY8kqQX9FIPQpLUQ0wQkqRKJghJUiUThCSpkglCklTJBCFJqtTVy1x71fZdP60sXzb/4A5HIkndYw9CklTJHsREbThx4r8999b2xSFJNbEHIUmqZIKQJFUyQUiSKjkH0UVv/YvvVJZ/9V1v7HAkkvRi9iAkSZXsQdSs6p6K80foOUhSL7EHIUmqZIKQJFVyiGkcRnoEB4zvMRwXP7F+9B02jHIsb7KT1CH2ICRJlUwQkqRKHU8QEXFERHwrIu6NiHsiYn0pnxsRmyJie/k+tNOxSZJe0I0exHPA/8zM1wDHAesi4kjgAuDmzFwG3FzWJUld0vEEkZmPZebtZfnfgHuBhcDpwBVltyuAt3U6NknSC7p6FVNELAaOBm4DXpGZj0EjiUTE/C6G1rt8zLikDulagoiIg4EvAu/OzCcjotXfrQXWAixatKi+AMdptEtgJakfdSVBRMQsGsnhc5n5pVL8o4hYUHoPC4BdVb/NzI3ARoDBwcHsSMAd1q77LSRpMrpxFVMAnwbuzcyLmzZdD6wpy2uA6zodmyTpBd3oQRwP/D5wd0TcWcreB3wYuCYizgYeBt7ehdimnObeyPCHBPpYcUmj6XiCyMzvACNNOJzUyVj60UjDTw49jcz3bkgT453UkqRKJghJUiUThCSpkglCklTJ90HoRUaa1B2NE77S1GOCUO28ikjqTw4xSZIq2YOYRl70qtMRXm168RPV91qcP+dj7Q5JUg8zQUwRPixQUrs5xCRJqmSCkCRVcohJbTGRS2NH+41XOEndZw9CklTJHsQ05sS2pNHYg5AkVTJBSJIqOcSklr3oRrtx8CY7qf+YINRXvPJJ6hyHmCRJlexBqCPGPTzV/Jyoc29tbzCSWmKCUE9qvgT3/AnchPe8DSeO+PDBkR5W+DwTk6a5nhtiiohTIuL+iNgRERd0Ox5Jmq56qgcRETOA/wu8BRgCvhcR12fmD7obmbqp5eGpsXoEksalpxIEcCywIzMfAIiIq4HTAROExjSl7gzfcGJ36nVYrbMm83fuwN+q14aYFgKPNK0PlTJJUof1Wg8iKspyvx0i1gJry+pPI+L+CdZ1OPD4BH/ba6ZKW3qrHX9YdTq2rLfa0qoXt7k/21FtqrSl0Y7JnZ+/1spOvZYghoAjmtYHgJ3NO2TmRmDjZCuKiC2ZOTjZ4/SCqdKWqdIOmDptmSrtgKnTlk62o9eGmL4HLIuIJRFxALAKuL7LMUnStNRTPYjMfC4i/jvwdWAGcHlm3tPlsCRpWuqpBAGQmTcCN3agqkkPU/WQqdKWqdIOmDptmSrtgKnTlo61IzJz7L0kSdNOr81BSJJ6xLRMEL34OI+IuDwidkXEtqayuRGxKSK2l+9DS3lExMdL/N+PiOVNv1lT9t8eEWuayn8rIu4uv/l4REzqGrlR2nFERHwrIu6NiHsiYn0ft2V2RGyOiLtKWz5YypdExG0lrs+XCyqIiAPL+o6yfXHTsS4s5fdHxMlN5R07FyNiRkTcERE39Hk7Hip//zsjYksp68fza05EXBsR95X/v7yh59qRmdPqQ2Py+5+BpcABwF3AkT0Q1wnAcmBbU9mfAReU5QuAj5TlFcDXaNw3chxwWymfCzxQvg8ty4eWbZuBN5TffA04taZ2LACWl+VDgB8CR/ZpWwI4uCzPAm4rMV4DrCrlnwT+qCyfB3yyLK8CPl+Wjyzn2YHAknL+zej0uQicD/wtcENZ79d2PAQcPqysH8+vK4BzyvIBwJxea0ctf8Be/pT/wb7etH4hcGG34yqxLGb/BHE/sKAsLwDuL8sbgNXD9wNWAxuayjeUsgXAfU3l++1Xc5uuo/Fsrb5uC3AQcDvw72jcbDVz+PlE4+q7N5TlmWW/GH6O7duvk+cijXuKbgbeDNxQ4uq7dpTjP8SLE0RfnV/Ay4EHKfPAvdqO6TjE1E+P83hFZj4GUL7nl/KR2jBa+VBFea3K0MTRNP7Luy/bUoZl7gR2AZto/JfyE5n5XEX9z8dctv8EOIzxt7EOlwLvAX5Z1g+jP9sBjacrfCMitkbjyQrQf+fXUmA38Ndl2O+vIuJlvdaO6ZggxnycRx8YqQ3jLa9NRBwMfBF4d2Y+OdquFWU905bM3JuZv0njv8CPBV4zSv092ZaIOA3YlZlbm4tHqbsn29Hk+MxcDpwKrIuIE0bZt1fbMpPGkPInMvNo4Gc0hpRG0pV2TMcEMebjPHrIjyJiAUD53lXKR2rDaOUDFeW1iIhZNJLD5zLzS6W4L9uyT2Y+AdxCY/x3TkTsu4eouf7nYy7bfwX4MeNvY7sdD6yMiIeAq2kMM13ah+0AIDN3lu9dwJdpJO5+O7+GgKHMvK2sX0sjYfRWO+oaJ+zVD43M/QCNSbZ9E2pHdTuuEtti9p+D+Cj7T1j9WVn+HfafsNpcyufSGNc8tHweBOaWbd8r++6bsFpRUxsCuBK4dFh5P7ZlHjCnLL8U+AfgNOAL7D+5e15ZXsf+k7vXlOWj2H9y9wEaE7sdPxeBN/HCJHXftQN4GXBI0/I/Aqf06fn1D8Cry/IHSht6qh21nYi9/KFxRcAPaYwnv7/b8ZSYrgIeA56lkf3PpjHuezOwvXzv+8MHjRcr/TNwNzDYdJyzgB3lc2ZT+SCwrfzmMoZNjrWxHW+k0ZX9PnBn+azo07a8DrijtGUb8KelfCmNK0R20PhH9sBSPrus7yjblzYd6/0l3vtpupqk0+ci+yeIvmtHifmu8rlnX119en79JrClnF9fofEPfE+1wzupJUmVpuMchCSpBSYISVIlE4QkqZIJQpJUyQQhSapkgpAkVTJBSJIqmSAkSZX+P5zGrEgG2Vh5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1af724aec50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lalonde.groupby('treat')['re78'].plot(kind='hist', bins=20, alpha=0.8, legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This dataset is not a balanced trial.** Indeed people in the control group are very different from people in the test (treatment) group. Below is a plot of the different income distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "treat\n",
       "0    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "1    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "Name: re74, dtype: object"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD8CAYAAABthzNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFsxJREFUeJzt3X2QZXV95/H312FgREh4mMaMNJ2ZUWIClhmxh3VLRaPJgrMRQmpjZiqVoKADK2RxzVYJaCluyirXiKghCwyRBYwZQAmCLqIju0CsjY6DjjCKOMNDpIcpBgYVDc/Nd/+4v4Y7zenu29333HN7+v2qutXn/O55+Nwzt/s75/zOQ2QmkiSN96KmA0iS+pMFQpJUyQIhSapkgZAkVbJASJIqWSAkSZUsEJKkShYISVIlC4QkqdJeTQeYjcWLF+fSpUubjiFJc8ptt932cGYOTDXdnC4QS5cuZdOmTU3HkKQ5JSL+tZPpPMQkSapkgZAkVbJASJIqzek+CElqytNPP83IyAhPPPFE01EmtGjRIgYHB1m4cOGM5rdASNIMjIyMsP/++7N06VIiouk4L5CZ7Nq1i5GREZYtWzajZXiISZJm4IknnuDggw/uy+IAEBEcfPDBs9rDsUBI0gz1a3EYM9t8FghJUiX7ICSpC97+t9/q6vK+8pdvmHKaG2+8kTPPPJPR0VHe/e53c9ZZZ3U1w7wuENP5B+3kH0uSemV0dJTTTz+dDRs2MDg4yMqVKzn++OM54ogjurYODzFJ0hy0ceNGXvGKV7B8+XL23ntvVq9ezXXXXdfVdVggJGkO2r59O4cddthz44ODg2zfvr2r66itQETEpRGxMyK2tLVdFRGby+u+iNhc2pdGxONt711UVy5J2hNk5gvaun1WVZ19EJcBFwBXjDVk5p+ODUfEecAv2qa/OzNX1JhHkvYYg4OD3H///c+Nj4yM8LKXvayr66htDyIzbwUeqXovWmXuHcD6utYvSXuylStXsnXrVu69916eeuoprrzySo4//viurqOps5jeCDyYmVvb2pZFxPeBR4EPZeY/NxNNkqav12c67rXXXlxwwQUce+yxjI6OcvLJJ3PkkUd2dx1dXVrn1rD73sMOYCgzd0XEa4EvR8SRmfno+BkjYi2wFmBoaKgnYSWpH61atYpVq1bVtvyen8UUEXsBfwxcNdaWmU9m5q4yfBtwN/BbVfNn5rrMHM7M4YGBKZ+YJ0maoSZOc/194MeZOTLWEBEDEbGgDC8HDgfuaSCbJKmo8zTX9cC/AK+MiJGIOKW8tZoXdk4fA9weET8AvgSclpmVHdySpN6orQ8iM9dM0P7OirZrgGvqyiJJmj6vpJYkVbJASJIqzeu7uUpS11z8pu4u79Rbppzk5JNP5qtf/SqHHHIIW7ZsmXL66XIPQpLmqHe+853ceOONtS3fAiFJc9QxxxzDQQcdVNvyLRCSpEoWCElSJQuEJKmSBUKSVMnTXCWpGzo4LbXb1qxZw80338zDDz/M4OAgH/3oRznllFOmnrFDFghJmqPWr6/3mWseYpIkVbJASJIqWSAkaYYys+kIk5ptPguEJM3AokWL2LVrV98Wicxk165dLFq0aMbLsJNakmZgcHCQkZERHnrooaajTGjRokUMDg7OeH4LhCTNwMKFC1m2bFnTMWrlISZJUiULhCSpkgVCklSptgIREZdGxM6I2NLWdm5EbI+IzeW1qu29syNiW0TcFRHH1pVLktSZOvcgLgOOq2g/PzNXlNcNABFxBLAaOLLM8z8jYkGN2SRJU6itQGTmrcAjHU5+AnBlZj6ZmfcC24Cj68omSZpaE30QZ0TE7eUQ1IGl7VDg/rZpRkrbC0TE2ojYFBGb+vn8Y0ma63pdIC4EXg6sAHYA55X2qJi28vLEzFyXmcOZOTwwMFBPSklSbwtEZj6YmaOZ+SxwCc8fRhoBDmubdBB4oJfZJEm762mBiIglbaMnAmNnOF0PrI6IfSJiGXA4sLGX2SRJu6vtVhsRsR54M7A4IkaAjwBvjogVtA4f3QecCpCZP4yIq4EfAc8Ap2fmaF3ZJElTq61AZOaaiubPTTL9x4CP1ZVHkjQ9XkktSapkgZAkVbJASJIqWSAkSZXm9QODPvXzMzuf+OL9dh8/9ZbuhpGkPuMehCSpkgVCklTJAiFJqmSBkCRVskBIkipZICRJlSwQkqRKFghJUiULhCSpkgVCklTJAiFJqmSBkCRVskBIkipZICRJlWorEBFxaUTsjIgtbW1/ExE/jojbI+LaiDigtC+NiMcjYnN5XVRXLklSZ+rcg7gMOG5c2wbgVZn5auAnwNlt792dmSvK67Qac0mSOlBbgcjMW4FHxrV9IzOfKaPfBgbrWr8kaXaa7IM4Gfha2/iyiPh+RNwSEW9sKpQkqaWRR45GxAeBZ4AvlKYdwFBm7oqI1wJfjogjM/PRinnXAmsBhoaGehVZkuadnu9BRMRJwB8Cf5aZCZCZT2bmrjJ8G3A38FtV82fmuswczszhgYGBXsWWpHmnpwUiIo4DPgAcn5mPtbUPRMSCMrwcOBy4p5fZJEm7q+0QU0SsB94MLI6IEeAjtM5a2gfYEBEA3y5nLB0D/PeIeAYYBU7LzEcqFyxJ6onaCkRmrqlo/twE014DXFNXFknS9HkltSSpkgVCklTJAiFJqmSBkCRVskBIkip1VCAi4lV1B5Ek9ZdO9yAuioiNEfHesVt0S5L2bB0ViMx8A/BnwGHApoj4x4j4g1qTSZIa1XEfRGZuBT5E61YZbwI+Wx7+88d1hZMkNafTPohXR8T5wJ3AW4C3Z+bvlOHza8wnSWpIp7fauAC4BDgnMx8fa8zMByLiQ7UkkyQ1qtMCsQp4PDNHASLiRcCizHwsMz9fWzpJUmM67YP4JvDitvF9S5skaQ/VaYFYlJm/Ghspw/vWE0mS1A86LRD/FhFHjY2Ux4I+Psn0kqQ5rtM+iPcBX4yIB8r4EuBP64kkSeoHHRWIzPxuRPw28EoggB9n5tO1JpMkNWo6T5RbCSwt87wmIsjMK2pJJUlqXEcFIiI+D7wc2EzrmdEACVggJGkP1ekexDBwRGZmnWEkSf2j07OYtgC/Md2FR8SlEbEzIra0tR0UERsiYmv5eWBpj4j4bERsi4jb28+akiT1XqcFYjHwo4j4ekRcP/bqYL7LgOPGtZ0F3JSZhwM3lXGAtwGHl9da4MIOs0mSatDpIaZzZ7LwzLw1IpaOaz4BeHMZvhy4mdYdYk8AriiHsb4dEQdExJLM3DGTdUuSZqfT01xviYjfBA7PzG9GxL7Aghmu86Vjf/Qzc0dEHFLaDwXub5tupLTtViAiYi2tPQyGhoZmGEGSNJVOb/f9HuBLwMWl6VDgy13OEhVtL+gUz8x1mTmcmcMDAwNdjiBJGtNpH8TpwOuBR+G5hwcdMukcE3swIpYAlJ87S/sIrSfWjRkEHkCS1IhOC8STmfnU2EhE7EXF/+47dD1wUhk+Cbiurf0vytlMrwN+Yf+DJDWn007qWyLiHODF5VnU7wW+MtVMEbGeVof04ogYAT4CfBy4OiJOAX4K/EmZ/AZaz53YBjwGvGsan0OS1GWdFoizgFOAO4BTaf0x//upZsrMNRO89daKaZPWoSxJUh/o9CymZ2k9cvSSeuNIkvpFp/diupfqM4qWdz2RJKkvTOdeTGMW0eo3OKj7cSRJ/aKjs5gyc1fba3tmfhp4S83ZJEkN6vQQU/uN815Ea49i/1oSSZL6QqeHmM5rG34GuA94R9fTSJL6RqdnMf1e3UEkSf2l00NM75/s/cz8VHfiSJL6xXTOYlpJ63YYAG8HbmX3u69KkvYgnRaIxcBRmflLgIg4F/hiZr67rmCSpGZ1erO+IeCptvGngKVdTyNJ6hud7kF8HtgYEdfSuqL6ROCK2lJJkhrX6VlMH4uIrwFvLE3vyszv1xdLktS0Tg8xAewLPJqZnwFGImJZTZkkSX2g00eOfgT4AHB2aVoI/ENdoSRJzet0D+JE4Hjg3wAy8wG81YYk7dE6LRBPlQf6JEBEvKS+SJKkftBpgbg6Ii4GDoiI9wDfxIcHSdIerdOzmD5ZnkX9KPBK4MOZuaHWZJKkRk1ZICJiAfD1zPx9YNZFISJeCVzV1rQc+DBwAPAe4KHSfk5m3jDb9UmSZmbKApGZoxHxWET8emb+YrYrzMy7gBXwXPHZDlwLvAs4PzM/Odt1SJJmr9MrqZ8A7oiIDZQzmQAy87/Mcv1vBe7OzH+NiFkuSpLUTZ0WiP9dXt22GljfNn5GRPwFsAn4q8z8WQ3rlCR1YNICERFDmfnTzLy82yuOiL1pXVsxdvHdhcBf0zqV9q9pPcXu5Ir51gJrAYaGhrodS5JUTHWa65fHBiLimi6v+23A9zLzQYDMfDAzRzPzWVqn0B5dNVNmrsvM4cwcHhgY6HIkSdKYqQpEe8fA8i6vew1th5ciYknbeycCW7q8PknSNEzVB5ETDM9KROwL/AFwalvzJyJiRVnPfePekyT12FQF4ncj4lFaexIvLsOU8czMX5vJSjPzMeDgcW1/PpNlSZLqMWmByMwFvQoiSeov03kehCRpHrFASJIqWSAkSZUsEJKkShYISVIlC4QkqZIFQpJUyQIhSapkgZAkVbJASJIqWSAkSZUsEJKkShYISVIlC4QkqZIFQpJUyQIhSapkgZAkVbJASJIqWSAkSZUmfSZ1nSLiPuCXwCjwTGYOR8RBwFXAUuA+4B2Z+bOmMkrSfNb0HsTvZeaKzBwu42cBN2Xm4cBNZVyS1ICmC8R4JwCXl+HLgT9qMIskzWtNFogEvhERt0XE2tL20szcAVB+HtJYOkma5xrrgwBen5kPRMQhwIaI+HEnM5VishZgaGioznySNK81tgeRmQ+UnzuBa4GjgQcjYglA+bmzYr51mTmcmcMDAwO9jCxJ80ojBSIiXhIR+48NA/8B2AJcD5xUJjsJuK6JfJKk5g4xvRS4NiLGMvxjZt4YEd8Fro6IU4CfAn/SUD5JmvcaKRCZeQ/wuxXtu4C39j6RJGm8fjvNVZLUJywQkqRKFghJUiULhCSpkgVCklTJAiFJqmSBkCRVskBIkipZICRJlSwQkqRKFghJUiULhCSpUpMPDJpTtu781W7j7//bb0047Vf+8g11x5Gk2rkHIUmqZIGQJFWyQEiSKlkgJEmVLBCSpEoWCElSJQuEJKlSz6+DiIjDgCuA3wCeBdZl5mci4lzgPcBDZdJzMvOGXufrhrdPco3EeP1yzcRczCypXk1cKPcM8FeZ+b2I2B+4LSI2lPfOz8xPNpBJkjROzwtEZu4AdpThX0bEncChvc4hSZpco30QEbEUeA3wndJ0RkTcHhGXRsSBE8yzNiI2RcSmhx56qGoSSVIXNFYgImI/4BrgfZn5KHAh8HJgBa09jPOq5svMdZk5nJnDAwMDPcsrSfNNIwUiIhbSKg5fyMx/AsjMBzNzNDOfBS4Bjm4imySppecFIiIC+BxwZ2Z+qq19SdtkJwJbep1NkvS8Js5iej3w58AdEbG5tJ0DrImIFUAC9wGnNpBNklQ0cRbTt4CoeGtOXvMgSXsqr6SWJFXyiXKq1XSu0Aav0pb6iQWiYdP9Azod/rGVNBseYpIkVbJASJIqWSAkSZUsEJKkShYISVIlz2LStNV55pWk/mGB0LzgE/Ok6fMQkySpknsQe7C5eChoT8/s3onmEvcgJEmV3IOYoU/9/MwZz/v+Az7TxSSSVA8LxBxjYVI/83DbnsUC0YDZ/JFX/eZiP4hUBwuE1EP+D1tziQVC2gNYeFQHC4Q6NhcPjc3lfpe6DnX5ECd1qu8KREQcB3wGWAD8fWZ+vOFI0rxmn8z81VfXQUTEAuDvgLcBRwBrIuKIZlNJ0vzUb3sQRwPbMvMegIi4EjgB+FGjqfYQc/EQ0WzN9jPP5UNU0mz1W4E4FLi/bXwE+HcNZZEau+5krl7vMq3cF++3+/ipt3Q3zFxw8ZtmPm8Ptle/FYioaMvdJohYC6wto7+KiLtmsb7FwMOzmL+XzFqfmvK+ccZzfnXitzrIOvP1zlZF7s637WlVv/49Nbe+t6fFbPL+ZicT9VuBGAEOaxsfBB5onyAz1wHrurGyiNiUmcPdWFbdzFqfuZR3LmWFuZV3LmWF3uTtq05q4LvA4RGxLCL2BlYD1zecSZLmpb7ag8jMZyLiDODrtE5zvTQzf9hwLEmal/qqQABk5g3ADT1aXVcOVfWIWeszl/LOpawwt/LOpazQg7yRmVNPJUmad/qtD0KS1CfmZYGIiOMi4q6I2BYRZzWY476IuCMiNkfEptJ2UERsiIit5eeBpT0i4rMl8+0RcVTbck4q02+NiJO6mO/SiNgZEVva2rqWLyJeWz7/tjLvjM9znCDruRGxvWzfzRGxqu29s8t674qIY9vaK78b5cSJ75TPcFU5iWKmWQ+LiP8bEXdGxA8j4szS3q/bdqK8fbd9I2JRRGyMiB+UrB+dbPkRsU8Z31beXzrTz9DlvJdFxL1t23ZFae/tdyEz59WLVuf33cByYG/gB8ARDWW5D1g8ru0TwFll+Czgf5ThVcDXaF0r8jrgO6X9IOCe8vPAMnxgl/IdAxwFbKkjH7AR+Pdlnq8Bb+ty1nOB/1Yx7RHl330fYFn5PiyY7LsBXA2sLsMXAf95FlmXAEeV4f2Bn5RM/bptJ8rbd9u3fN79yvBC4Dtlm1UuH3gvcFEZXg1cNdPP0OW8lwH/qWL6nn4X5uMexHO388jMp4Cx23n0ixOAy8vw5cAftbVfkS3fBg6IiCXAscCGzHwkM38GbACO60aQzLwVeKSOfOW9X8vMf8nWt/iKtmV1K+tETgCuzMwnM/NeYBut70Xld6P8j+stwJcqPvdMsu7IzO+V4V8Cd9K6i0C/btuJ8k6kse1bttGvyujC8spJlt++zb8EvLXkmdZnmEnWKfJOpKffhflYIKpu5zHZl71OCXwjIm6L1hXiAC/NzB3Q+sUEDintE+Xu9efpVr5Dy/D49m47o+yKXzp2yGYGWQ8Gfp6Zz3Q7azmk8Rpa/3Ps+207Li/04faNiAURsRnYSesP5d2TLP+5TOX9X5Q8Pft9G583M8e27cfKtj0/IvYZn7fDXLP6LszHAjHl7Tx66PWZeRStu9eeHhHHTDLtRLn75fNMN18vcl8IvBxYAewAzivtfZE1IvYDrgHel5mPTjbpNHP1Km9fbt/MHM3MFbTuxHA08DuTLL/xbTs+b0S8Cjgb+G1gJa3DRh9oIu98LBBT3s6jVzLzgfJzJ3AtrS/zg2W3kPJzZ5l8oty9/jzdyjdShse3d01mPlh++Z4FLqG1fWeS9WFau/J7jWufsYhYSOuP7Rcy859Kc99u26q8/bx9S76fAzfTOlY/0fKfy1Te/3Vahyp7/vvWlve4clgvM/NJ4H8x8207u+9Cp50Ve8qL1sWB99DqeBrrZDqygRwvAfZvG/5/tPoO/obdOyo/UYb/I7t3Tm3M5zun7qXVMXVgGT6oizmXsnvHb9fy0bq1yut4vvNsVZezLmkb/q+0jikDHMnuHZD30Op8nPC7AXyR3Ts53zuLnEHrWPCnx7X35badJG/fbV9gADigDL8Y+GfgDydaPnA6u3dSXz3Tz9DlvEvatv2ngY838V3o6R/FfnnROhPgJ7SOTX6woQzLy5frB8APx3LQOv55E7C1/Bz7Rw5aD1O6G7gDGG5b1sm0OtG2Ae/qYsb1tA4dPE3rfyKndDMfMAxsKfNcQLlws4tZP1+y3E7rnl7tf9A+WNZ7F21ndUz03Sj/XhvLZ/gisM8ssr6B1m7+7cDm8lrVx9t2orx9t32BVwPfL5m2AB+ebPnAojK+rby/fKafoct5/0/ZtluAf+D5M516+l3wSmpJUqX52AchSeqABUKSVMkCIUmqZIGQJFWyQEiSKlkgJEmVLBCSpEoWCElSpf8PPf8JzyAO7OMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1af725fdb70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lalonde.groupby('treat')['re74'].plot(kind='hist', bins=20, alpha=0.8, legend=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "treat\n",
       "0    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "1    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "Name: age, dtype: object"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFClJREFUeJzt3X2QXXV9x/H3FxJcUNqYkGBkiRsko0SnPHTD0AFRwSqNErAjFGqdKMFoixgHpxIYp8iMzuBMBVJx1ChqRCRElAZSG40RdPgHSIQKEjoopLBJSkKEAR8gJHz7xz3rbMM+3H0499y75/2a2bnnnHvOPZ+cnbvf/H6/8xCZiSSpvg6oOoAkqVoWAkmqOQuBJNWchUCSas5CIEk1ZyGQpJqzEEhSzVkIJKnmLASSVHNTqg7QjMMOOyx7enqqjiFJHWXz5s1PZebMkdbriELQ09PDpk2bqo4hSR0lIv6nmfXsGpKkmrMQSFLNWQgkqeY6YoxAkqrw4osv0tfXx/PPP191lGF1dXXR3d3N1KlTx7S9hUCShtDX18ehhx5KT08PEVF1nEFlJrt376avr4+5c+eO6TNK7RqKiGkRcUtEPBwRWyLiryJiekRsiIhHitdXl5lBksbq+eefZ8aMGW1bBAAighkzZoyr1VL2GMEKYH1mvhE4FtgCLAc2ZuY8YGMxL0ltqZ2LQL/xZiytEETEnwGnAtcDZOaezHwGOAtYVay2Cji7rAySpJGVOUZwFLAL+GZEHAtsBpYBh2fmDoDM3BERs0rMIEkT5swv3jWhn3f7xac0td769etZtmwZ+/bt48ILL2T58ontSCmzEEwBTgAuzsy7I2IFo+gGioilwFKAOXPmlJNwBOP5pTf7C5ak4ezbt4+LLrqIDRs20N3dzYIFC1i0aBHz58+fsH2UOUbQB/Rl5t3F/C00CsOTETEboHjdOdjGmbkyM3szs3fmzBFvlSFJk9I999zD0UcfzVFHHcVBBx3Eeeedx9q1ayd0H6UVgsz8X+CJiHhDseh04CHgNmBxsWwxMLH/IkmaRLZt28aRRx75p/nu7m62bds2ofso+zqCi4EbI+Ig4FHgQzSKz5qIWAI8DpxTcgZJ6liZ+bJlE30mU6mFIDPvB3oHeev0MvcrSZNFd3c3TzzxxJ/m+/r6eO1rXzuh+/BeQ5LUxhYsWMAjjzzCY489xp49e1i9ejWLFi2a0H14iwlJalIVZwNOmTKF6667jne9613s27ePCy64gDe96U0Tu48J/TRJ0oRbuHAhCxcuLO3z7RqSpJqzEEhSzVkIJKnmLASSVHMWAkmqOQuBJNWcp49KUrO++taJ/byP/GzEVS644ALWrVvHrFmzePDBByd2/wVbBJLUxj74wQ+yfv36UvdhIZCkNnbqqacyffr0Uvcx6buGJvqJQpI02dgikKSasxBIUs1ZCCSp5ib9GIEkTZgmTvecaOeffz533nknTz31FN3d3Vx55ZUsWbJkQvdhIZCkNnbTTTeVvg+7hiSp5iwEklRzFgJJGkZmVh1hROPNaCGQpCF0dXWxe/futi4Gmcnu3bvp6uoa82c4WCxJQ+ju7qavr49du3ZVHWVYXV1ddHd3j3l7C4EkDWHq1KnMnTu36hils2tIkmrOQiBJNVdq11BEbAWeA/YBezOzNyKmAzcDPcBW4NzMfLrMHJKkobWiRfD2zDwuM3uL+eXAxsycB2ws5iVJFamia+gsYFUxvQo4u4IMkqRC2WcNJfDjiEjgq5m5Ejg8M3cAZOaOiJg12IYRsRRYCjBnzpySY0688T4Q5/aLT5mgJJI0vLILwcmZub34Y78hIh5udsOiaKwE6O3tbd+rOSSpw5XaNZSZ24vXncCtwInAkxExG6B43VlmBknS8EorBBHxyog4tH8aeCfwIHAbsLhYbTGwtqwMkqSRldk1dDhwa0T07+e7mbk+Iu4F1kTEEuBx4JwSM0iSRlBaIcjMR4FjB1m+Gzi9rP1KkkbHK4slqeYsBJJUcxYCSao5C4Ek1ZyFQJJqzkIgSTVnIZCkmrMQSFLNWQgkqeYsBJJUcxYCSao5C4Ek1ZyFQJJqzkIgSTVnIZCkmrMQSFLNWQgkqeYsBJJUcxYCSao5C4Ek1ZyFQJJqzkIgSTVnIZCkmrMQSFLNWQgkqeZKLwQRcWBE3BcR64r5uRFxd0Q8EhE3R8RBZWeQJA2tFS2CZcCWAfOfB67JzHnA08CSFmSQJA2h1EIQEd3Au4GvF/MBnAbcUqyyCji7zAySpOGV3SK4FvgU8FIxPwN4JjP3FvN9wBElZ5AkDaO0QhAR7wF2ZubmgYsHWTWH2H5pRGyKiE27du0qJaMkqdwWwcnAoojYCqym0SV0LTAtIqYU63QD2wfbODNXZmZvZvbOnDmzxJiSVG+lFYLMvCwzuzOzBzgP+Glmvh+4A3hfsdpiYG1ZGSRJI6viOoJLgUsi4tc0xgyuryCDJKkwZeRVxi8z7wTuLKYfBU5sxX4lSSNrqkUQEW8uO4gkqRrNdg19JSLuiYh/iohppSaSJLVUU4UgM08B3g8cCWyKiO9GxF+XmkyS1BJNDxZn5iPAp2kM9r4V+LeIeDgi/rascJKk8jU7RvAXEXENjXsGnQacmZnHFNPXlJhPklSyZs8aug74GnB5Zv6xf2Fmbo+IT5eSTJLUEs0WgoXAHzNzH0BEHAB0ZeYfMvOG0tJJkkrX7BjBT4CDB8wfUiyTJHW4ZgtBV2b+rn+mmD6knEiSpFZqthD8PiJO6J+JiL8E/jjM+pKkDtHsGMEngO9FRP+dQmcDf1dOJElSKzVVCDLz3oh4I/AGGs8UeDgzXyw1mSSpJUZz07kFQE+xzfERQWZ+u5RUkqSWaaoQRMQNwOuB+4F9xeIELASS1OGabRH0AvMzc9DHSkqSOlezZw09CLymzCCSpGo02yI4DHgoIu4BXuhfmJmLSkklSWqZZgvBZ8oMIUmqTrOnj/4sIl4HzMvMn0TEIcCB5UaTJLVCs2cNfRhYCkyncfbQEcBXgNPLi1ZvZ37xrkr2e/vFp1SyX0nVaXaw+CLgZOBZ+NNDamaVFUqS1DrNFoIXMnNP/0xETKFxHYEkqcM1Wwh+FhGXAwcXzyr+HnB7ebEkSa3SbCFYDuwCHgA+AvyQxvOLJUkdrtmzhl6i8ajKr5UbR5LUas2eNfQYg4wJZOZRE55IktRSo7nXUL8u4Bwap5IOKSK6gJ8Dryj2c0tmXhERc4HVxfa/AD4wcCBaktRaTY0RZObuAT/bMvNa4LQRNnsBOC0zjwWOA86IiJOAzwPXZOY84GlgyTjyS5LGqdmuoRMGzB5Ao4Vw6HDbFHcq7X/O8dTiJ2kUkL8vlq+icfuKLzedWJI0oZrtGvrCgOm9wFbg3JE2iogDgc3A0cCXgN8Az2Tm3mKVPhpXKQ+27VIaVzMzZ86cJmNKkkar2bOG3j6WD8/MfcBxETENuBU4ZrDVhth2JbASoLe314vXJKkkzXYNXTLc+5l59QjvPxMRdwInAdMiYkrRKugGtjeZVZJUgmYvKOsF/pFGN84RwEeB+TTGCQYdK4iImUVLgIg4GHgHsAW4A3hfsdpiYO1Yw0uSxm80D6Y5ITOfA4iIzwDfy8wLh9lmNrCqGCc4AFiTmesi4iFgdUR8FrgPuH7M6SVJ49ZsIZgDDDzXfw/QM9wGmflL4PhBlj8KnNjkfiVJJWu2ENwA3BMRt9IY3H0v8O3SUkmSWqbZs4Y+FxH/CbylWPShzLyvvFiqyngeiONDbaTO1OxgMcAhwLOZuQLoK24VIUnqcE0Vgoi4ArgUuKxYNBX4TlmhJEmt02yL4L3AIuD3AJm5nRFuMSFJ6gzNDhbvycyMiASIiFeWmEkluvqZZaNa/5JpK0pKIqldNNsiWBMRX6VxVfCHgZ/gQ2okaVJo9qyhfy2eVfws8AbgXzJzQ6nJJEktMWIhKK4M/lFmvgOYdH/87SqRVHcjdg0VdxD9Q0T8eQvySJJarNnB4ueBByJiA8WZQwCZ+fFSUkmSWqbZQvAfxY8kaZIZthBExJzMfDwzV7UqkCSptUYaI/j3/omI+H7JWSRJFRipEMSA6aPKDCJJqsZIhSCHmJYkTRIjDRYfGxHP0mgZHFxMU8xnZv5ZqekkSaUbthBk5oGtCiJJqkazp4+qRbzSWVKrjebBNJKkSchCIEk1ZyGQpJpzjEATZjwPvge4/eJTJiiJpNGwRSBJNWchkKSaK61rKCKOBL4NvAZ4CViZmSsiYjpwM9ADbAXOzcyny8ox0Ty9U9JkU2aLYC/wycw8BjgJuCgi5gPLgY2ZOQ/YWMxLkipSWiHIzB2Z+Yti+jlgC3AEcBbQf1vrVcDZZWWQJI2sJWMEEdEDHA/cDRyemTugUSyAWa3IIEkaXOmFICJeBXwf+ERmPjvS+gO2WxoRmyJi065du8oLKEk1V2ohiIipNIrAjZn5g2LxkxExu3h/NrBzsG0zc2Vm9mZm78yZM8uMKUm1VlohiIgArge2ZObVA966DVhcTC8G1paVQZI0sjKvLD4Z+ADwQETcXyy7HLgKWBMRS4DHgXNKzCBJGkFphSAz7+L/P+pyoNPL2q8kaXS8sliSas5CIEk1ZyGQpJqzEEhSzVkIJKnmLASSVHMWAkmqOQuBJNWchUCSas6H12tYo30iG/hUNqnT2CKQpJqzEEhSzdk1pAk32u6k/q6kM79415j3efvFp4x52/EYT+bxqurfrMnHFoEk1ZyFQJJqzq4h1V6V3Tvj0YldaWpPtggkqeYsBJJUcxYCSao5C4Ek1ZyFQJJqzkIgSTXn6aOaFDr1FNC6Ge/vydNey2GLQJJqzkIgSTVnIZCkmittjCAivgG8B9iZmW8ulk0HbgZ6gK3AuZn5dFkZJA2ujrenqOO/uVlltgi+BZyx37LlwMbMnAdsLOYlSRUqrRBk5s+B3+63+CxgVTG9Cji7rP1LkprT6jGCwzNzB0DxOmuoFSNiaURsiohNu3btallASaqbth0szsyVmdmbmb0zZ86sOo4kTVqtLgRPRsRsgOJ1Z4v3L0naT6sLwW3A4mJ6MbC2xfuXJO2nzNNHbwLeBhwWEX3AFcBVwJqIWAI8DpxT1v7bxWgf5N5uny9p8iutEGTm+UO8dXpZ+5QkjV7bDhZLklrDu4+q44y2O+ySaStKSqJW68S7zHbCFc22CCSp5iwEklRzFgJJqjnHCCRpBJ04NjEatggkqeYsBJJUc3YNSRqVyd5NUke2CCSp5iwEklRzFgJJqjnHCDTpeUuK4Xl8ZItAkmrOQiBJNWfXkKTS2f3U3mwRSFLNWQgkqebsGlLl2u25y2V3Y9hNonZji0CSas5CIEk1ZyGQpJpzjEAap7LHONptDEWTjy0CSao5C4Ek1VwlXUMRcQawAjgQ+HpmXlVFDkmj14quKk+xba2Wtwgi4kDgS8DfAPOB8yNifqtzSJIaqugaOhH4dWY+mpl7gNXAWRXkkCRRTSE4AnhiwHxfsUySVIEqxghikGX5spUilgJLi9nfRcR/D/F5hwFPjTXMurFu2JxxZSuZ2cbGbGNTcra3jGfjtj1u8fFxZ3tdMytVUQj6gCMHzHcD2/dfKTNXAitH+rCI2JSZvRMXb+KYbWzMNjZmGxuzVdM1dC8wLyLmRsRBwHnAbRXkkCRRQYsgM/dGxMeAH9E4ffQbmfmrVueQJDVUch1BZv4Q+OEEfdyI3UcVMtvYmG1szDY2tc8WmS8bp5Uk1Yi3mJCkmuuoQhAR34iInRHx4IBln4mIbRFxf/GzsKJsR0bEHRGxJSJ+FRHLiuXTI2JDRDxSvL66jbJVfuwioisi7omI/yqyXVksnxsRdxfH7ebixIJ2yfatiHhswHE7rtXZBmQ8MCLui4h1xXzlx22IXO10zLZGxANFjk3Fssq/p8NkK/172lGFAPgWcMYgy6/JzOOKn4kaexitvcAnM/MY4CTgouLWGcuBjZk5D9hYzLdLNqj+2L0AnJaZxwLHAWdExEnA54ts84CngSVtlA3gnwcct/sryNZvGbBlwHw7HDd4eS5on2MG8PYiR/+pme3wPe23fzYo+XvaUYUgM38O/LbqHIPJzB2Z+Yti+jkaX4IjaNw+Y1Wx2irg7DbKVrls+F0xO7X4SeA04JZieVXHbahsbSEiuoF3A18v5oM2OG775+oQlX9Pq9RRhWAYH4uIXxZdR5U06QaKiB7geOBu4PDM3AGNP8jArOqSvSwbtMGxK7oR7gd2AhuA3wDPZObeYpXKbkOyf7bM7D9unyuO2zUR8YoqsgHXAp8CXirmZ9Aex23/XP3a4ZhBo5j/OCI2F3cwgPb5ng6WDUr+nk6GQvBl4PU0mu47gC9UGSYiXgV8H/hEZj5bZZb9DZKtLY5dZu7LzONoXGV+InDMYKu1NlWx0/2yRcSbgcuANwILgOnApa3OFRHvAXZm5uaBiwdZtaXHbYhc0AbHbICTM/MEGndAvigiTq0wy/4Gy1b697TjC0FmPll8WV8CvkbjD0klImIqjT+0N2bmD4rFT0bE7OL92TT+Z9kW2drp2BV5ngHupDGOMS0i+q9zGfQ2JK00INsZRVdbZuYLwDep5ridDCyKiK007uB7Go3/iVd93F6WKyK+0ybHDIDM3F687gRuLbK0xfd0sGyt+J52fCHo/+UV3gs8ONS6JecI4HpgS2ZePeCt24DFxfRiYG27ZGuHYxcRMyNiWjF9MPAOGmMYdwDvK1ar6rgNlu3hAX8wgkZfcsuPW2ZelpndmdlD4zYtP83M91PxcRsi1z+0wzEr9v/KiDi0fxp4Z5GlHb6ng2Zrxfe0ox5eHxE3AW8DDouIPuAK4G3FqWgJbAU+UlG8k4EPAA8UfcoAlwNXAWsiYgnwOHBOG2U7vw2O3WxgVTQeWHQAsCYz10XEQ8DqiPgscB+NQtYu2X4aETNpdMXcD3y0gmxDuZTqj9tgbmyTY3Y4cGujHjEF+G5mro+Ie6n+ezpUthvK/p56ZbEk1VzHdw1JksbHQiBJNWchkKSasxBIUs1ZCCSp5iwEklRzFgJJqjkLgSTV3P8BBUlnXAY0jSsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1af724a95f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lalonde.groupby('treat')['age'].plot(kind='hist', bins=20, alpha=0.8, legend=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "treat\n",
       "0    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "1    AxesSubplot(0.125,0.125;0.775x0.755)\n",
       "Name: educ, dtype: object"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD8CAYAAABthzNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFX9JREFUeJzt3X2QHHWdx/H3VxJcUTTkAQ0smKA5NFjnyS0Upxxy4gnmJKClXihLgwlGS0Q87uoIaolelVVyKopypUThDB7yID4EOEQjipZVRzAgSiBgIiBsiGSJPKiIIfF7f0wvGZfe7GR3enp29/2qmtp++M30t5phPvl1/7o7MhNJkoZ6Rt0FSJK6kwEhSSplQEiSShkQkqRSBoQkqZQBIUkqZUBIkkoZEJKkUgaEJKnUlLoLGIuZM2fmnDlz6i5DksaVm2+++aHMnDVSu3EdEHPmzGHt2rV1lyFJ40pE/LqVdh5ikiSVMiAkSaUqC4iIuCgitkTEuqZln4yIOyPiFxHxrYiY1rTurIjYGBF3RcSxVdUlSWpNlecgvgKcD1zctGw1cFZmbo+Ic4CzgDMjYj6wCDgE2A/4fkT8VWbuqLA+SRq1J598kv7+fp544om6SxlWT08Pvb29TJ06dVTvrywgMvPHETFnyLLvNc3eCLy5mD4BuCwz/wTcExEbgcOB/6uqPkkai/7+fvbee2/mzJlDRNRdztNkJlu3bqW/v5+5c+eO6jPqPAexBPhOMb0/cH/Tuv5i2dNExLKIWBsRawcGBiouUZLKPfHEE8yYMaMrwwEgIpgxY8aYeji1BEREfAjYDlwyuKikWemj7jJzRWb2ZWbfrFkjDuOVpMp0azgMGmt9Hb8OIiIWA28AjsmdzzvtBw5oatYLPNDp2iRJO3U0ICLiOOBM4NWZ+XjTqquAr0XEuTROUs8DbupkbZI0Fsd//idt/byrTztyxDbXXXcdp59+Ojt27OCUU05h+fLlba2hsoCIiEuBo4GZEdEPnE1j1NIzgdVF1+fGzHxPZt4eEVcAd9A49HSqI5ikzhvLj1wrP2hqnx07dnDqqaeyevVqent7Oeyww1i4cCHz589v2zaqHMV0UsniC3fR/uPAx6uqR5ImkptuuokXv/jFHHTQQQAsWrSIVatWtTUgvJJaksahTZs2ccABO0/d9vb2smnTprZuw4CQpHFo5xifndo9qsqAkKRxqLe3l/vv33n5WH9/P/vtt19bt2FASNI4dNhhh7Fhwwbuuecetm3bxmWXXcbChQvbuo1x/TwISeoWnR7FNWXKFM4//3yOPfZYduzYwZIlSzjkkEPau422fpokqWMWLFjAggULKvt8DzFJkkoZEJKkUgaEJKmUASFJKmVASJJKGRCSpFIOc5Wkdrjg1e39vHf/aMQmS5Ys4ZprrmHfffdl3bp17d0+9iAkadw6+eSTue666yr7fANCksapo446iunTp1f2+QaEJKmUASFJKmVASJJKGRCSpFIOc5WkdmhhWGq7nXTSSdxwww089NBD9Pb28rGPfYylS5e27fMNCEkapy699NJKP99DTJKkUgaEJKmUASFJo5SZdZewS2Otz4CQpFHo6elh69atXRsSmcnWrVvp6ekZ9Wd4klqSRqG3t5f+/n4GBgbqLmVYPT099Pb2jvr9lQVERFwEvAHYkpkvK5ZNBy4H5gD3Am/NzIcjIoDzgAXA48DJmXlLVbVJ0lhNnTqVuXPn1l1Gpao8xPQV4Lghy5YD12fmPOD6Yh7g9cC84rUM+EKFdUmSWlBZQGTmj4HfDll8ArCymF4JnNi0/OJsuBGYFhGzq6pNkjSyTp+kfn5mbgYo/u5bLN8fuL+pXX+x7GkiYllErI2Itd187E+SxrtuGcUUJctKhwZk5orM7MvMvlmzZlVcliRNXp0OiAcHDx0Vf7cUy/uBA5ra9QIPdLg2SVKTTgfEVcDiYnoxsKpp+Tui4Qjg0cFDUZKkelQ5zPVS4GhgZkT0A2cDnwCuiIilwH3AW4rm19IY4rqRxjDXd1ZVlySpNZUFRGaeNMyqY0raJnBqVbVIknZft5ykliR1GQNCklTKgJAklTIgJEmlDAhJUikDQpJUyoCQJJUyICRJpQwISVIpA0KSVMqAkCSVMiAkSaUMCElSKQNCklTKgJAklTIgJEmlDAhJUikDQpJUyoCQJJUyICRJpQwISVIpA0KSVMqAkCSVMiAkSaUMCElSKQNCklSqloCIiH+JiNsjYl1EXBoRPRExNyLWRMSGiLg8IvasozZJUkPHAyIi9gfeD/Rl5suAPYBFwDnAZzJzHvAwsLTTtUmSdqrrENMU4FkRMQXYC9gMvAa4sli/EjixptokSdQQEJm5CfgUcB+NYHgUuBl4JDO3F836gf3L3h8RyyJibUSsHRgY6ETJkjQp1XGIaR/gBGAusB/wbOD1JU2z7P2ZuSIz+zKzb9asWdUVKkmTXB2HmF4L3JOZA5n5JPBN4JXAtOKQE0Av8EANtUmSCnUExH3AERGxV0QEcAxwB/BD4M1Fm8XAqhpqkyQV6jgHsYbGyehbgNuKGlYAZwJnRMRGYAZwYadrkyTtNGXkJu2XmWcDZw9ZfDdweA3lSJJKeCW1JKmUASFJKtVSQETEy6ouRJLUXVrtQXwxIm6KiPdGxLRKK5IkdYWWAiIzjwTeBhwArI2Ir0XEP1ZamSSpVi2fg8jMDcCHaQxHfTXwuYi4MyLeVFVxkqT6tHoO4q8j4jPAeho31Ts+M19aTH+mwvokSTVp9TqI84EvAR/MzD8OLszMByLiw5VUJkmqVasBsQD4Y2buAIiIZwA9mfl4Zn61suokSbVp9RzE94FnNc3vVSyTJE1QrQZET2b+fnCmmN6rmpIkSd2g1YD4Q0QcOjgTEX8L/HEX7SVJ41yr5yA+AHw9Igaf0TAb+OdqSpIkdYOWAiIzfxoRLwEOBgK4s3jYjyRpgtqd230fBswp3vOKiCAzL66kKklS7VoKiIj4KvAi4FZgR7E4AQNCkiaoVnsQfcD8zMwqi5EkdY9WRzGtA15QZSGSpO7Sag9iJnBHRNwE/GlwYWYurKQqSVLtWg2Ij1ZZhCSp+7Q6zPVHEfFCYF5mfj8i9gL2qLY0SVKdWr3d97uAK4ELikX7A9+uqihJUv1aPUl9KvAq4DF46uFB+1ZVlCSpfq0GxJ8yc9vgTERMoXEdhCRpgmo1IH4UER8EnlU8i/rrwNXVlSVJqlurAbEcGABuA94NXEvj+dSSpAmq1VFMf6bxyNEvtWOjETEN+DLwMhqHqpYAdwGX07jf073AWzPz4XZsT5K0+1odxXRPRNw99DWG7Z4HXJeZLwFeDqyn0Uu5PjPnAdcX85KkmuzOvZgG9QBvAaaPZoMR8VzgKOBkgOLk97aIOAE4umi2ErgBOHM025AkjV1LPYjM3Nr02pSZnwVeM8ptHkTjfMZ/R8TPIuLLEfFs4PmZubnY3mYcRitJtWr1dt+HNs0+g0aPYu8xbPNQ4LTMXBMR57Ebh5MiYhmwDODAAw8cZQmSpJG0eojp003T2ylOIo9ym/1Af2auKeavpBEQD0bE7MzcHBGzgS1lb87MFcAKgL6+Pq/FkKSKtDqK6R/atcHM/E1E3B8RB2fmXcAxwB3FazHwieLvqnZtU5K0+1o9xHTGrtZn5rm7ud3TgEsiYk/gbuCdNA5dXRERS4H7aJwIlyTVZHdGMR0GXFXMHw/8GLh/NBvNzFv5y5FRg44ZzedJktpvdx4YdGhm/g4gIj4KfD0zT6mqMElSvVq91caBwLam+W00rniWJE1QrfYgvgrcFBHfonFrjDcCF1dWlTTJHf/5n4z6vVefdmQbK9Fk1uoopo9HxHeAvy8WvTMzf1ZdWZKkurV6iAlgL+CxzDwP6I+IuRXVJEnqAq3erO9sGvdFOqtYNBX4n6qKkiTVr9UexBuBhcAfADLzAUZ/qw1J0jjQakBsy8ykeMxocXM9SdIE1mpAXBERFwDTIuJdwPdp08ODJEndqdVRTJ8qnkX9GHAw8JHMXF1pZZKkWo0YEBGxB/DdzHwtYChIXW4s11BIzUY8xJSZO4DHI+J5HahHktQlWr2S+gngtohYTTGSCSAz319JVZKk2rUaEP9bvCRJk8QuAyIiDszM+zJzZacKkiR1h5HOQXx7cCIivlFxLZKkLjJSQETT9EFVFiJJ6i4jBUQOMy1JmuBGOkn98oh4jEZP4lnFNMV8ZuZzK61OklSbXQZEZu7RqUIkSd1ld54HIUmaRAwISVIpA0KSVKrVK6mlcWksN667+rQj21hJtc595PSW2p0x7byKK9FEYkBIFfGuqhrvPMQkSSpVW0BExB4R8bOIuKaYnxsRayJiQ0RcHhF71lWbJKneHsTpwPqm+XOAz2TmPOBhYGktVUmSgJoCIiJ6gX8CvlzMB/Aa4MqiyUrgxDpqkyQ11NWD+Czw78Cfi/kZwCOZub2Y7wf2r6MwSVJDxwMiIt4AbMnMm5sXlzQtvTlgRCyLiLURsXZgYKCSGiVJ9QxzfRWwMCIWAD3Ac2n0KKZFxJSiF9ELPFD25sxcAawA6Ovr8w6zUgW8rkJQQw8iM8/KzN7MnAMsAn6QmW8Dfgi8uWi2GFjV6dokSTt103UQZwJnRMRGGuckLqy5Hkma1Gq9kjozbwBuKKbvBg6vsx5J0k7d1IOQJHURA0KSVMqAkCSVMiAkSaUMCElSKQNCklTKgJAklTIgJEmlDAhJUikDQpJUqtZbbUjatVbvqipVwR6EJKmUASFJKmVASJJKGRCSpFIGhCSplKOYpHa54NV/MXvuI78vbeZznDVe2IOQJJUyICRJpQwISVIpA0KSVMqT1Op6x3/+J3WXIE1K9iAkSaUMCElSKQNCklTKgJAklep4QETEARHxw4hYHxG3R8TpxfLpEbE6IjYUf/fpdG2SpJ3q6EFsB/41M18KHAGcGhHzgeXA9Zk5D7i+mJck1aTjAZGZmzPzlmL6d8B6YH/gBGBl0WwlcGKna5Mk7VTrOYiImAO8AlgDPD8zN0MjRIB966tMklRbQETEc4BvAB/IzMd2433LImJtRKwdGBiorkBJmuRquZI6IqbSCIdLMvObxeIHI2J2Zm6OiNnAlrL3ZuYKYAVAX19fdqRgSSMa6xXvV592ZJsqUbvUMYopgAuB9Zl5btOqq4DFxfRiYFWna5Mk7VRHD+JVwNuB2yLi1mLZB4FPAFdExFLgPuAtNdSmSeTcR07fdYMLntP4++4fVV+M1IU6HhCZ+RMghll9TCdrkSQNzyupJUmlDAhJUimfByGpK4xlFJQjoKphD0KSVMqAkCSVMiAkSaUMCElSKQNCklTKgJAklTIgJEmlvA5CHTHWO31K6jx7EJKkUgaEJKmUASFJKuU5CLXM8wjqVt7HqRr2ICRJpexBSKrciE/vK5wx7byKK9HusAchSSplD0LSpOb5i+HZg5AklTIgJEmlPMQ0zjjUVBNZqyezwRPanWAPQpJUyh6EAIchSno6exCSpFKTtgcx1mP5E31423DsaUiTR9cFREQcB5wH7AF8OTM/UXNJpTxZ3L1250SnNBZ1/g504h+pXXWIKSL2AP4LeD0wHzgpIubXW5UkTU7d1oM4HNiYmXcDRMRlwAnAHbVWpco4rFFVa3ePcjJ9D7uqBwHsD9zfNN9fLJMkdVhkZt01PCUi3gIcm5mnFPNvBw7PzNOa2iwDlhWzBwN3jXJzM4GHxlBuJ42XWq2z/cZLrdbZXlXX+cLMnDVSo247xNQPHNA03ws80NwgM1cAK8a6oYhYm5l9Y/2cThgvtVpn+42XWq2zvbqlzm47xPRTYF5EzI2IPYFFwFU11yRJk1JX9SAyc3tEvA/4Lo1hrhdl5u01lyVJk1JXBQRAZl4LXNuBTY35MFUHjZdarbP9xkut1tleXVFnV52kliR1j247ByFJ6hITPiAi4riIuCsiNkbE8pL1z4yIy4v1ayJiTg01HhARP4yI9RFxe0Q87cqeiDg6Ih6NiFuL10c6XWdTLfdGxG1FHWtL1kdEfK7Yp7+IiENrqPHgpn11a0Q8FhEfGNKmtn0aERdFxJaIWNe0bHpErI6IDcXffYZ57+KizYaIWFxDnZ+MiDuL/7bfiohpw7x3l9+TDtT50YjY1PTfd8Ew793lb0QH6ry8qcZ7I+LWYd7bsf35lMycsC8aJ7p/BRwE7An8HJg/pM17gS8W04uAy2uoczZwaDG9N/DLkjqPBq6pe58WtdwLzNzF+gXAd4AAjgDWdMH34Dc0xn53xT4FjgIOBdY1LftPYHkxvRw4p+R904G7i7/7FNP7dLjO1wFTiulzyups5XvSgTo/CvxbC9+NXf5GVF3nkPWfBj5S9/4cfE30HsRTt+7IzG3A4K07mp0ArCymrwSOiYjoYI1k5ubMvKWY/h2wnvF9BfkJwMXZcCMwLSJm11jPMcCvMvPXNdbwFzLzx8Bvhyxu/i6uBE4seeuxwOrM/G1mPgysBo7rZJ2Z+b3M3F7M3kjjeqVaDbM/W9HKb0Tb7KrO4nfnrcClVW1/d030gGjl1h1PtSm+9I8CMzpSXYniENcrgDUlq/8uIn4eEd+JiEM6WthfSuB7EXFzcWX7UN12y5RFDP8/XbfsU4DnZ+ZmaPyjAdi3pE237dslNHqLZUb6nnTC+4pDYRcNc8ium/bn3wMPZuaGYdZ3fH9O9IAo6wkMHbbVSpuOiIjnAN8APpCZjw1ZfQuNQyQvBz4PfLvT9TV5VWYeSuOuu6dGxFFD1nfTPt0TWAh8vWR1N+3TVnXTvv0QsB24ZJgmI31PqvYF4EXA3wCbaRy+Gapr9idwErvuPXR8f070gBjx1h3NbSJiCvA8RtdVHZOImEojHC7JzG8OXZ+Zj2Xm74vpa4GpETGzw2UO1vJA8XcL8C0a3fRmrez3Tnk9cEtmPjh0RTft08KDg4fiir9bStp0xb4tTo6/AXhbFgfIh2rhe1KpzHwwM3dk5p+BLw2z/W7Zn1OANwGXD9emjv050QOilVt3XAUMjgR5M/CD4b7wVSmOPV4IrM/Mc4dp84LBcyMRcTiN/3ZbO1flU3U8OyL2HpymccJy3ZBmVwHvKEYzHQE8OnjopAbD/qusW/Zpk+bv4mJgVUmb7wKvi4h9ikMmryuWdUw0Hup1JrAwMx8fpk0r35NKDTnv9cZhtt8tt/d5LXBnZvaXraxtf3byjHgdLxojan5JY6TCh4pl/0Hjyw3QQ+Pww0bgJuCgGmo8kka39hfArcVrAfAe4D1Fm/cBt9MYZXEj8Mqa9udBRQ0/L+oZ3KfNtQaNBz/9CrgN6Kup1r1o/OA/r2lZV+xTGqG1GXiSxr9il9I493U9sKH4O71o20fj6YqD711SfF83Au+soc6NNI7bD35XB0cB7gdcu6vvSYfr/Grx/fsFjR/92UPrLOaf9hvRyTqL5V8Z/F42ta1tfw6+vJJaklRqoh9ikiSNkgEhSSplQEiSShkQkqRSBoQkqZQBIUkqZUBIkkoZEJKkUv8PKewwTbu7idoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1af7507f400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lalonde.groupby('treat')['educ'].plot(kind='hist', bins=20, alpha=0.8, legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# notation, aims and assumptions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notations.**\n",
    "* Y represents the *response*, here is is 1978 earnings ('re78')\n",
    "* D represents the *treatment*: the job training program ('treat')\n",
    "* X represents the *confounding variables*, here it likely is age, education, race and marital status. X is also called a *covariate* or the *counter factual*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Aims.**\n",
    "What we want to know here is the Average Treatment Effect (**ATE**):\n",
    "$$\\Delta = E[Y_{1} - Y_{0}]$$\n",
    "\n",
    "However, as we saw, if we try to estimate this quantity from the row observational distribution, we get:\n",
    "$$\\Delta_{raw} = E[Y|D=1] - E[Y|D=0] = E[Y_{1}|D=1] - E[Y_{0}|D=0] \\\\\n",
    "\\neq \\Delta = E[Y_{1} - Y_{0}]$$\n",
    "\n",
    "because: \n",
    "$$E[Y_{i}|D=i] \\neq E[Y_{i}]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**General problem.** If we believe that age, education, race, and marital status all have a likely influence on earnings Y, we need a way to disentangle the effect of D on Y from the perturbative effect of X on Y."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assumptions.**\n",
    "The Causalinference package is based on a typical assumption called *unconfoundedness* or *ignorability*:\n",
    "$$(Y(0), Y(1)) \\perp D \\; | \\; X$$\n",
    "\n",
    "Indeed we saw that the treatment assignment is probably not independent of each subject's potential outcomes, e.g. poorer people are more represented in the treatment group than in the control group.\n",
    "\n",
    "However the treatment is assumed to be *unconfounded* in the sense that the dependence between the treatment assignment and the outcomes is *only through* something we observe, namely the covariates X.\n",
    "\n",
    "What this means is that if we control for X, i.e. look across people with similar levels of X, then the difference between treated and control should be attributable to the treatment itself, just as a randomized experiment would be.\n",
    "\n",
    "This is the assumption, and if it doesn't hold our results could be completely wrong."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest type of model we can use is a linear model:\n",
    "\n",
    "$$Y_0 = \\alpha + \\beta X + \\epsilon$$\n",
    "$$Y_1 = Y_0 + \\gamma D $$\n",
    "\n",
    "If this is accurate, fitting the following model to the data using linear regression will give us an estimate of the Average Treatment Effect (ATE):\n",
    "$$Y = \\alpha + \\beta X + \\gamma D$$\n",
    "\n",
    "$\\epsilon$ is called a residual and represents the noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "covariates = ['age', 'educ', 'black', 'hispan', 'married', 'nodegree', 're74', 're75']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Treatment Effect Estimates: OLS\n",
      "\n",
      "                     Est.       S.e.          z      P>|z|      [95% Conf. int.]\n",
      "--------------------------------------------------------------------------------\n",
      "           ATE   1548.244    734.521      2.108      0.035    108.584   2987.904\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\romai\\Anaconda3\\lib\\site-packages\\causalinference\\estimators\\ols.py:21: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  olscoef = np.linalg.lstsq(Z, Y)[0]\n"
     ]
    }
   ],
   "source": [
    "# we use the CausalModel method from the causalinference package\n",
    "\n",
    "causal = CausalModel(\n",
    "    Y=lalonde['re78'].values, \n",
    "    D=lalonde['treat'].values, \n",
    "    X=lalonde[covariates].values)\n",
    "\n",
    "causal.est_via_ols(adj=1)\n",
    "# adj=1 corresponds to the simplicity of the model we entered\n",
    "# This is called a \"constant treatment effect\"\n",
    "\n",
    "print(causal.estimates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This model predicts that the Average Treatment Effect (ATE, the job training) is $1548 extra annual earnings. This is very different from our previous raw results predicting that the job training had negative effects on earnings!**\n",
    "\n",
    "Assuming that our model accurately describes the counterfactual X, CausalModel provides the 95% confidence interval. What this means is that, if we were to repeat this treatment experiment, in 95% of the cases the Average Treatment Effect would be within that interval. That doesn't mean that the true value is within that interval.\n",
    "\n",
    "Based on the assumption that the residuals are normally distributed, the 95% confidence interval is calculated as:\n",
    "$$ AVG \\pm 1.96 * STD  / \\sqrt{n}$$\n",
    "\n",
    "In practice, as the confidence interval is very large, my interpretation is that the experiment should have had more people if a better estimate of the extra earnings was desired. Ways to control the standard deviation could also be explored. \n",
    "\n",
    "Overall, assuming that we controlled for all the effects and did it well, it seems that the job training had a positive effect on earnings. Indeed, although the standard deviation is very large, the p value of 0.035 rejects the null hypothesis (no effect) with a confidence level of 97.5%. However, the truth is that we don't know if we modelled the counterfactual well, and this could change everything... As we will see later, estimators such as the Ordinary Least Square (OLS) estimator can behave poorly when there is not enough covariate overlap, and that's because the estimator needs to extrapolate too much from one group to another.\n",
    "\n",
    "A more structured approach as we will see below can allow us to increase our confidence that the covariants are well controlled for. We will see many steps, but one simple idea is the technique of matching: the idea is to find for each sample which received the treatment a *similar* sample in the control group, and to directly compare these values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Structure for a more complete approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Pre-processing phase:**\n",
    "1. assess covariate balance\n",
    "2. estimate propensity score\n",
    "3. trim sample\n",
    "4. stratify sample\n",
    "\n",
    "**Estimation phase:**\n",
    "1. blocking estimator\n",
    "or/and\n",
    "2. matching estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-processing phase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the pre-processing phase, the data is inspected and manipulated to allow credible analysis to be conducted on it.\n",
    "\n",
    "As we discussed in the previous section, one key method for disantangling the treatment effect from the covariant effects is the *matching technique*. In this technique we compare subjects that have similar covariate values (i.e. same age, rage, income etc). However, our ability to compare such pairs depends heavily on the degree of overlap of the covariates between the treatment and control group. This is called *covariate balance*. \n",
    "\n",
    "Said otherwise, to control the effect of education, one way is to look at people in the tested group and in the non-tested group that all have the same level of education, say a bachelor degree. However, if nobody in the test group has a bachelor degree while many do in the non-test group, this procedure is impossible.\n",
    "\n",
    "**(1) assess covariate balance** to assess whether how easily people can be matched. If there is too much unbalance, direct matching will rarely be possible, and we may need to use more complex techniques, if at all possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['treat', 'age', 'educ', 'black', 'hispan', 'married', 'nodegree',\n",
       "       're74', 're75', 're78'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lalonde.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Summary Statistics\n",
      "\n",
      "                       Controls (N_c=429)         Treated (N_t=185)             \n",
      "       Variable         Mean         S.d.         Mean         S.d.     Raw-diff\n",
      "--------------------------------------------------------------------------------\n",
      "              Y     6984.170     7294.162     6349.144     7867.402     -635.026\n",
      "\n",
      "                       Controls (N_c=429)         Treated (N_t=185)             \n",
      "       Variable         Mean         S.d.         Mean         S.d.     Nor-diff\n",
      "--------------------------------------------------------------------------------\n",
      "             X0       28.030       10.787       25.816        7.155       -0.242\n",
      "             X1       10.235        2.855       10.346        2.011        0.045\n",
      "             X2        0.203        0.403        0.843        0.365        1.668\n",
      "             X3        0.142        0.350        0.059        0.237       -0.277\n",
      "             X4        0.513        0.500        0.189        0.393       -0.719\n",
      "             X5        0.597        0.491        0.708        0.456        0.235\n",
      "             X6     5619.237     6788.751     2095.574     4886.620       -0.596\n",
      "             X7     2466.484     3291.996     1532.055     3219.251       -0.287\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(causal.summary_stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Raw-diff** is the raw difference between the means of the control and treatment groups.\n",
    "\n",
    "As we saw previously, the treated group (trained) is earning $635 less than the control group, which is surprising."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Nor-diff** in this package is Imbens and Rubin's normalized differences (2015) in average covariates, defined as:\n",
    "$$\\frac{\\bar{X}_{k,t}-\\bar{X}_{k,c}}{\\sqrt{(s_{k,t}^2 + s_{k,c}^2) / 2}}$$\n",
    "\n",
    "Here $\\bar{X}_{k,t}$ and $s_{k,t}$ are the sample mean and sample standard deviation of the kth covariate of the treatment group, and $\\bar{X}_{k,c}$ and $s_{k,c}$ are the analogous statistics for the control group.\n",
    "\n",
    "The aim here is to assess the overlap between the control and treatment groups. It can be seen that X2, X4, and X6 (black, married, revenue in 1974) have a large normalized difference, beyond 0.5. This can be interpreted as an imbalance. Concretely, there are way more black people, less married people and lower income in 1974 in the treatment group than in the control group.\n",
    "\n",
    "The impact of imbalance is to make the matching technique harder to apply.\n",
    "We'll see later how we can try to correct for it (however, ideally the study would be more balanced!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(2) Propensity Score** - the probability of receiving the treatment, conditional on the covariates.\n",
    "\n",
    "Propensity is useful for assessing and improving covariate balance. Indeed a theorem by Rosenbaum and Rubin in 1983, proves that, for subjects that share the same propensity score (even if their covariate vectors are different), the difference between the treated and the control units actually identifies a conditional average treatment effect.\n",
    "\n",
    "**Thus, instead of matching on the covariate vectors X themselves, we can also match on the single-dimensional propensity score p(X), aggregate across subjects, and still arrive at a valid estimate of the overall average treatment effect.**\n",
    "\n",
    "$$ E[Y(1) - Y(0) | p(X)] \\approx E[Y(1) - Y(0)]$$\n",
    "\n",
    "This is if $p(X) = P(D=1 | X)$, which the CausalInference package estimates for us using a sequence of likelihood ratio tests.\n",
    "\n",
    "reference: http://laurence-wong.com/software/propensity-score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Estimated Parameters of Propensity Score\n",
      "\n",
      "                    Coef.       S.e.          z      P>|z|      [95% Conf. int.]\n",
      "--------------------------------------------------------------------------------\n",
      "     Intercept    -21.096      2.687     -7.851      0.000    -26.363    -15.829\n",
      "            X2      2.635      0.367      7.179      0.000      1.915      3.354\n",
      "            X4     -3.026      0.717     -4.222      0.000     -4.431     -1.621\n",
      "            X6      0.000      0.000      0.847      0.397     -0.000      0.000\n",
      "            X3      5.137      1.845      2.785      0.005      1.521      8.753\n",
      "            X1      1.175      0.316      3.713      0.000      0.555      1.796\n",
      "            X5      0.376      0.450      0.836      0.403     -0.505      1.258\n",
      "            X7      0.000      0.000      1.496      0.135     -0.000      0.000\n",
      "            X0      0.988      0.142      6.983      0.000      0.711      1.266\n",
      "         X0*X0     -0.015      0.002     -6.524      0.000     -0.019     -0.010\n",
      "         X1*X1     -0.064      0.018     -3.539      0.000     -0.099     -0.028\n",
      "         X6*X0     -0.000      0.000     -2.876      0.004     -0.000     -0.000\n",
      "         X6*X6      0.000      0.000      2.420      0.016      0.000      0.000\n",
      "         X3*X0     -0.223      0.081     -2.752      0.006     -0.382     -0.064\n",
      "         X4*X3      2.845      1.071      2.656      0.008      0.746      4.945\n",
      "         X2*X4      1.525      0.781      1.952      0.051     -0.006      3.055\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#this function estimates the propensity score, so that propensity methods can be employed\n",
    "causal.est_propensity_s()\n",
    "print(causal.propensity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(3) Trim sample.** This excludes subjects with extreme propensity scores. Indeed it will be very hard to match those extreme subjects, so the usual strategy is to focus attention on the remaining units that exhibit a higher degree of covariate balance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extreme propensity is a very high probability to be either in the control group or the treatment group\n",
    "# that makes matching difficult\n",
    "\n",
    "#by default, causal.cutoff is set to 1\n",
    "# the trim function will drop units whose estimated propensity score <= 0.1 or >= 0.9\n",
    "#causal.cutoff = 0.1\n",
    "#causal.trim()\n",
    "\n",
    "#however, there is a procedure that tried to select an optimal cutoff value\n",
    "causal.trim_s()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Summary Statistics\n",
      "\n",
      "                       Controls (N_c=157)         Treated (N_t=140)             \n",
      "       Variable         Mean         S.d.         Mean         S.d.     Raw-diff\n",
      "--------------------------------------------------------------------------------\n",
      "              Y     5476.633     6020.122     6351.987     6397.833      875.353\n",
      "\n",
      "                       Controls (N_c=157)         Treated (N_t=140)             \n",
      "       Variable         Mean         S.d.         Mean         S.d.     Nor-diff\n",
      "--------------------------------------------------------------------------------\n",
      "             X0       23.777        7.418       24.986        7.510        0.162\n",
      "             X1       10.210        2.405       10.329        2.177        0.052\n",
      "             X2        0.471        0.501        0.836        0.372        0.826\n",
      "             X3        0.248        0.433        0.071        0.258       -0.496\n",
      "             X4        0.261        0.441        0.221        0.417       -0.093\n",
      "             X5        0.650        0.479        0.664        0.474        0.031\n",
      "             X6     2673.838     4479.406     2110.413     4215.220       -0.130\n",
      "             X7     1906.682     3082.253     1635.661     3414.688       -0.083\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(causal.summary_stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this new subset, the normal difference for most variables is rather balanced.\n",
    "Only X2 (number of black people) is still unbalanced.\n",
    "\n",
    "It is worth noting that the initial sample of 614 people (429 controls, 185 treated) has been drastically trimmed to 297 people (157 controls, 140 treated).\n",
    "\n",
    "In this more balanced sub-sample, without using any model, the average earnings in 1978 is more like what we would expect: populations that received training (treated) earn in average $875 more than the control group."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(4) Stratify sample** - group similar subjects together. People are grouped in layers of similar propensity scores. These bins should have an improved covariate balance, and we should be able to compare and match samples within those bins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the default is to have 5 bins with equal number of samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Stratification Summary\n",
      "\n",
      "              Propensity Score         Sample Size     Ave. Propensity   Outcome\n",
      "   Stratum      Min.      Max.  Controls   Treated  Controls   Treated  Raw-diff\n",
      "--------------------------------------------------------------------------------\n",
      "         1     0.089     0.201        66         9     0.136     0.171  1398.600\n",
      "         2     0.205     0.463        48        26     0.323     0.384  1614.014\n",
      "         3     0.465     0.674        27        47     0.555     0.572  -505.292\n",
      "         4     0.676     0.909        16        58     0.779     0.812  2210.672\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\romai\\Anaconda3\\lib\\site-packages\\causalinference\\core\\summary.py:110: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return (mean_t-mean_c) / np.sqrt((sd_c**2+sd_t**2)/2)\n"
     ]
    }
   ],
   "source": [
    "# however, it is possible to split the sample in a more data-driven way.\n",
    "# The larger the sample, the more bins we can afford, so that samples can be increasingly similar within smaller bins\n",
    "# the limit in dividing too much is that there are too few datapoints in each bin for the bins to be statistically different (t-test)\n",
    "\n",
    "causal.stratify_s()\n",
    "print(causal.strata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within bins, the raw difference in outcome should be a good representation of the real treatment effect. For example:\n",
    "\n",
    "- People in group 1 are unlikely to be in the treatment group (well off?). For them, the training improved their earnings by $1399 in average.\n",
    "\n",
    "- People in group 4 are likely to be in the treatment group (poor?). For them, the training improved their earnings even more, with a mean of $2211 for that year 1978.\n",
    "\n",
    "Something that looks quite bad is that outcomes for the group 3 are totally different from that of the other groups. The trend seems to be that the higher the propensity score, the higher the raw difference in outcome for each stratum. but this one shows opposite results... This may be a sign that we haven't controlled for enough factors (or that the propensity calculation is wrong?). Or it might also be a true representation or reality: some people may benefit from the job training, while other may not. It might also be random and the reflection that we are working with a relatively small sample (74 elements in bin 3).\n",
    "\n",
    "Let's see in the analysis phase if regressions within each stratum will be able to control for confounding variables better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation phase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the estimation phase, treatment effects of the training can be estimated in several ways."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(1) The blocking estimator** - although each layer of the stratum is pretty balanced and gives reasonable raw results, this estimator goes further and controls for the confounding factors within each layer of the stratum. More precisely, this estimator uses a least square estimate within each propensity bin, and from this produces an overall average treatment effect estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#causal.est_via_blocking()\n",
    "#print(causal.estimates)\n",
    "\n",
    "# for some reason I'm having a singular matrix when calculating this blocking estimator\n",
    "# on one of the stratum\n",
    "# I've tried changing the stratum structure and the set of variables,\n",
    "# however, the singularity persists when calculating the covariance matrix\n",
    "# this would need a closer look at the dataset, which I haven't taken the time to do yet\n",
    "\n",
    "# this is one of the issue of this causalinference package:\n",
    "# it needs to invert large matrixes, which can fail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(2) The matching estimator** - although each layer of the stratum is pretty balanced and gives reasonable raw results, this matching estimator controls for the confounding factors by matching even more thinely samples within each layer of the stratum. More precisely, this pairing is done via nearest-neighborhood matching. If the matching is imperfect, biias correction is recommended."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If issues arrive with least square, such as excessive extrapolation, this matching estimator pushes until the end the unconfoundedness assumption and nonparametrically matches subjects with similar covariate values together. In other words, if the confounding factors are equal for both element of a pair, the difference between the two will be the real treatment effect. In the causalinference package, samples are weighted by the inverse of the standard deviation of the sample covariate, so as to normalize. \n",
    "\n",
    "Where matching discrepancy exist, least square will be used, but very locally, so large extrapolations should be less of a problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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    383.508   1206.472      0.318      0.751  -1981.177   2748.193\n",
      "           ATC    594.156   1504.175      0.395      0.693  -2354.027   3542.338\n",
      "           ATT    147.281   1368.636      0.108      0.914  -2535.245   2829.807\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\romai\\Anaconda3\\lib\\site-packages\\causalinference\\estimators\\matching.py:100: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  return np.linalg.lstsq(X, Y)[0][1:]  # don't need intercept coef\n"
     ]
    }
   ],
   "source": [
    "causal.est_via_matching(bias_adj=True)\n",
    "print(causal.estimates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model provides estimates of three quantities: ATE, ATT and ATC:\n",
    "* ATE is the Average Treatment Effect, and this is what we are most interested in. $ATE = E[Y_1-Y_0] \\approx E[Y_1-Y_0 | X]$\n",
    " * Here is seems that the average effect of the treatment (job training) was to increase earnings by \\$384.\n",
    " * However, this effect may just be a random variation, and the treatment may well not have any impact (the null hypothesis). The probability to reject the null hypothesis is 25%. The most common interpretation of this number is that the treatment of job trainings did not have a statistically significant impact on earnings, given the models and data processing we did\n",
    "* ATT is the Average Treatment effect of the Treated $ATT = E[Y_1-Y_0 | D=1]$\n",
    "* ATC is the Average Treatment effect of the Control $ATC = E[Y_1-Y_0 | D=0]$"
   ]
  },
  {
   "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   1027.087    883.785      1.162      0.245   -705.131   2759.305\n",
      "           ATC    754.217   1006.164      0.750      0.453  -1217.865   2726.298\n",
      "           ATT   1333.092    937.670      1.422      0.155   -504.741   3170.925\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\romai\\Anaconda3\\lib\\site-packages\\causalinference\\estimators\\matching.py:100: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  return np.linalg.lstsq(X, Y)[0][1:]  # don't need intercept coef\n"
     ]
    }
   ],
   "source": [
    "# allowing several matches\n",
    "causal.est_via_matching(bias_adj=True, matches=4)\n",
    "print(causal.estimates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allowing several matches attributes \\$1027 of revenue increase to the treatment, with 75% probability to be significant. A common interpretation would be still to reject this as proof of statistical significance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The effect of training is hard to establish firmly. Although it seems the sample from Lalonde had positive effects, it is actually quite likely to be without any effect.\n",
    "\n",
    "This isn't so far from what Lalonde concluded: http://people.hbs.edu/nashraf/LaLonde_1986.pdf\n",
    "By glancing at it, Lalonde seemed to know the gender of participants, which does not seem to be in this dataset, or may be hidden in the NSW vs AFDC.\n",
    "\n",
    "\n",
    "More work could be done to better estimate the counterfactual. For instance we could introduce polynomial variables to capture non-linear effects and/or introduce categorical variables to bin numerical variables such aseducation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This was an example of how the CausalInference package could be used, and our conclusions are attached to those models. This package relies heavily on propensity matching and its ignorability / confoundedness assumption.\n",
    "\n",
    "Other models exist, e.g. Bayesian models. This will be for another time for us. Meanwhile, the curious can have a look at this other post: https://engl.is/causal-analysis-introduction-examples-in-python-and-pymc.html"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}