Instrumental Variables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Instrumental Variables for Treatment Effects in Quasi-Experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an example from Chapter 11 of the book [Methods Matter: Improving Causal Inference in Educational and Social Science](https://www.amazon.com/Methods-Matter-Improving-Inference-Educational-ebook/dp/B00HNSNBO4) by Richard Murnane and John Willett. This chapter explains how we can use Instrumental Variables Estimation to get a treatment effect using observational data.\n",
    "\n",
    "In this example, we use IV to recover the treatment effect in cases where random assignment is applied to an offer to participate, where not every participates (compliers AND never takers allowed!), and where other people participate through some other means. Specifically, we would like to evaluate the effectiveness of a financial aid offer on the likelihood of a student finishing 8th grade in Bogotá, Colombia.\n",
    "\n",
    "The data is sampled from [Vouchers for Private Schooling in Colombia: Evidence from a Randomized Natural Experiment](http://www.nber.org/papers/w8343) (2002) by Joshua Angrist and others.\n",
    "\n",
    "Python implementation credit goes to J. Nathan Matias, who provides implementations of these methods in a [github repository](https://github.com/natematias/research_in_python)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import math\n",
    "import statsmodels.formula.api as smf  # for doing statistical regression\n",
    "import statsmodels.api as sm      # access to the wider statsmodels library, including R datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "def coefficient_to_odds_ratio(coefficient):\n",
    "    \"\"\"\n",
    "    This function interprets a coefficient from a Logistic Regression model in as an odds ratio\n",
    "    \n",
    "    Given:\n",
    "    coefficient -- coefficient from loistic regression output\n",
    "    \n",
    "    Returns:\n",
    "    odds ratio (i.e. p/(1-p))\n",
    "    \"\"\"\n",
    "    odds = math.exp(coefficient)\n",
    "    print(\"Those using financial aid are {:2f} times more likely to finish the 8th grade!\".format(odds))\n",
    "    \n",
    "def ate_odds_ratio(prob0, prob1):\n",
    "    \"\"\"\n",
    "    This function calculates the average treatment effects as an odds ratio\n",
    "    when the outcome is binary.\n",
    "    \n",
    "    Given:\n",
    "    prob0 -- probability given no treatment\n",
    "    prob1 -- probability given treatment\n",
    "    \n",
    "    Returns:\n",
    "    odds ratio (i.e. p/(1-p))\n",
    "    \"\"\"\n",
    "    odds = (prob1/(1-prob1)) - (prob0/(1-prob0))\n",
    "    print(\"Odds ratio: {}\".format(odds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the Dataset\n",
    "\n",
    "The dataset includes the following variables:\n",
    "- `finish8th`: did the student finish 8th grade or not (outcome variable)\n",
    "- `won_lottry`: won the lottery to receive offer of financial aid\n",
    "- `use_fin_aid`: did the student use financial aid of any kind (not exclusive to the lottery) or not\n",
    "- `base_age`: student age\n",
    "- `male`: is the student male or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>won_lottry</th>\n",
       "      <th>male</th>\n",
       "      <th>base_age</th>\n",
       "      <th>finish8th</th>\n",
       "      <th>use_fin_aid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>10.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     id  won_lottry  male  base_age  finish8th  use_fin_aid\n",
       "0   3.0         1.0   0.0      11.0        1.0          1.0\n",
       "1   4.0         0.0   1.0      11.0        1.0          1.0\n",
       "2   5.0         0.0   1.0      11.0        1.0          0.0\n",
       "3   6.0         0.0   0.0       9.0        0.0          0.0\n",
       "4  10.0         1.0   1.0      11.0        1.0          1.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Import data set\n",
    "voucher_df = pd.read_sas('colvoucher.sas7bdat') # reading in sas file\n",
    "voucher_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================================================================\n",
      "                              OVERALL SUMMARY\n",
      "==============================================================================\n",
      "                id   won_lottry         male     base_age    finish8th  \\\n",
      "count  1171.000000  1171.000000  1171.000000  1171.000000  1171.000000   \n",
      "mean   1357.010248     0.505551     0.504697    12.004270     0.681469   \n",
      "std     890.711584     0.500183     0.500192     1.347038     0.466106   \n",
      "min       3.000000     0.000000     0.000000     7.000000     0.000000   \n",
      "25%     616.000000     0.000000     0.000000    11.000000     0.000000   \n",
      "50%    1280.000000     1.000000     1.000000    12.000000     1.000000   \n",
      "75%    1982.500000     1.000000     1.000000    13.000000     1.000000   \n",
      "max    4030.000000     1.000000     1.000000    17.000000     1.000000   \n",
      "\n",
      "       use_fin_aid  \n",
      "count  1171.000000  \n",
      "mean      0.581554  \n",
      "std       0.493515  \n",
      "min       0.000000  \n",
      "25%       0.000000  \n",
      "50%       1.000000  \n",
      "75%       1.000000  \n",
      "max       1.000000  \n",
      "==============================================================================\n",
      "                         USE FINANCIAL AID = 0\n",
      "==============================================================================\n",
      "                id  won_lottry        male    base_age   finish8th  \\\n",
      "count   490.000000  490.000000  490.000000  490.000000  490.000000   \n",
      "mean   1472.571429    0.102041    0.518367   12.097959    0.604082   \n",
      "std     931.468109    0.303011    0.500173    1.388897    0.489547   \n",
      "min       5.000000    0.000000    0.000000    7.000000    0.000000   \n",
      "25%     707.250000    0.000000    0.000000   11.000000    0.000000   \n",
      "50%    1412.500000    0.000000    1.000000   12.000000    1.000000   \n",
      "75%    2103.250000    0.000000    1.000000   13.000000    1.000000   \n",
      "max    4030.000000    1.000000    1.000000   17.000000    1.000000   \n",
      "\n",
      "       use_fin_aid  \n",
      "count        490.0  \n",
      "mean           0.0  \n",
      "std            0.0  \n",
      "min            0.0  \n",
      "25%            0.0  \n",
      "50%            0.0  \n",
      "75%            0.0  \n",
      "max            0.0  \n",
      "==============================================================================\n",
      "                         USE FINANCIAL AID = 1\n",
      "==============================================================================\n",
      "                id  won_lottry        male    base_age   finish8th  \\\n",
      "count   681.000000  681.000000  681.000000  681.000000  681.000000   \n",
      "mean   1273.860499    0.795888    0.494860   11.936858    0.737151   \n",
      "std     851.223236    0.403347    0.500341    1.312992    0.440505   \n",
      "min       3.000000    0.000000    0.000000    9.000000    0.000000   \n",
      "25%     564.000000    1.000000    0.000000   11.000000    0.000000   \n",
      "50%    1168.000000    1.000000    0.000000   12.000000    1.000000   \n",
      "75%    1751.000000    1.000000    1.000000   13.000000    1.000000   \n",
      "max    4023.000000    1.000000    1.000000   16.000000    1.000000   \n",
      "\n",
      "       use_fin_aid  \n",
      "count        681.0  \n",
      "mean           1.0  \n",
      "std            0.0  \n",
      "min            1.0  \n",
      "25%            1.0  \n",
      "50%            1.0  \n",
      "75%            1.0  \n",
      "max            1.0  \n"
     ]
    }
   ],
   "source": [
    "print(\"==============================================================================\")\n",
    "print(\"                              OVERALL SUMMARY\"                                 )\n",
    "print(\"==============================================================================\")\n",
    "\n",
    "print(voucher_df.describe())\n",
    "\n",
    "for i in range(2):\n",
    "    print(\"==============================================================================\")\n",
    "    print(\"                         USE FINANCIAL AID = %(i)d\" % {\"i\":i}                  )\n",
    "    print(\"==============================================================================\")\n",
    "    print(voucher_df[voucher_df['use_fin_aid']==i].describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average Treatment Effect\n",
    "\n",
    "We can calculate the ATE empirically simply using the summary stats above, however this is an unconditional average treatment effect.\n",
    "\n",
    "Use the fact that \n",
    "\n",
    "$ ATE = E[Y=1|T=1] - E[Y=1| T=0] $.\n",
    "\n",
    "In case of binary outcome, we use\n",
    "\n",
    "$ ATE = Pr[Y=1|T=1] - Pr[Y=1|T=0] $\n",
    "\n",
    "and subsitute the conditional expectations/probabilities with their empirical estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.133069"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate unconditional average treatment effect\n",
    "ate = 0.737151 - 0.604082\n",
    "ate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Odds ratio: 1.2786901424379509\n"
     ]
    }
   ],
   "source": [
    "# Calculate this in terms of odds ratio\n",
    "ate_odds_ratio(0.604082, 0.737151)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Naive Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================================================================\n",
      "                            LOGISTIC REGRESSION\n",
      "==============================================================================\n",
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:              finish8th   No. Observations:                 1171\n",
      "Model:                            GLM   Df Residuals:                     1167\n",
      "Model Family:                Binomial   Df Model:                            3\n",
      "Link Function:                  logit   Scale:                             1.0\n",
      "Method:                          IRLS   Log-Likelihood:                -694.99\n",
      "Date:                Wed, 27 Jun 2018   Deviance:                       1390.0\n",
      "Time:                        15:40:56   Pearson chi2:                 1.17e+03\n",
      "No. Iterations:                     4                                         \n",
      "===============================================================================\n",
      "                  coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------\n",
      "Intercept       4.2496      0.596      7.128      0.000       3.081       5.418\n",
      "use_fin_aid     0.5773      0.130      4.441      0.000       0.323       0.832\n",
      "male           -0.4192      0.130     -3.216      0.001      -0.675      -0.164\n",
      "base_age       -0.2967      0.048     -6.156      0.000      -0.391      -0.202\n",
      "===============================================================================\n"
     ]
    }
   ],
   "source": [
    "print(\"==============================================================================\")\n",
    "print(\"                            LOGISTIC REGRESSION\"                               )\n",
    "print(\"==============================================================================\")\n",
    "result = smf.glm(formula = \" finish8th ~ use_fin_aid + male + base_age\", \n",
    "                 data=voucher_df,\n",
    "                 family=sm.families.Binomial()).fit()\n",
    "print(result.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Those using financial aid are 1.781220 times more likely to finish the 8th grade!\n"
     ]
    }
   ],
   "source": [
    "# Log odd units to odds ratio units\n",
    "coefficient_to_odds_ratio(result.params[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two Stage Least Squares Estimation\n",
    "\n",
    "Recall that using 2SLS estimation for the instrumental variables model requires estimating two equations:\n",
    "\n",
    "__Stage 1:__ \n",
    "\n",
    "$ T = Z\\alpha + X\\psi + \\nu $\n",
    "\n",
    "__Stage 2:__\n",
    "\n",
    "$ Y = {\\hat{T}} \\gamma + X\\beta + \\epsilon $\n",
    "\n",
    "The instrument we choose must satisfy the two conditions:\n",
    " 1. __Instrument relevance:__ $ Cov(Z, T) \\neq 0 $\n",
    " 2. __Instrument exogeneity:__ $ Cov(Z, \\epsilon) = 0 $\n",
    " \n",
    "Why is there endogeneity and what instrument is available in our dataset to combat this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==============================================================================\n",
      "                                  FIRST STAGE\n",
      "==============================================================================\n",
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:            use_fin_aid   No. Observations:                 1171\n",
      "Model:                            GLM   Df Residuals:                     1167\n",
      "Model Family:                Binomial   Df Model:                            3\n",
      "Link Function:                  logit   Scale:                             1.0\n",
      "Method:                          IRLS   Log-Likelihood:                -488.00\n",
      "Date:                Wed, 27 Jun 2018   Deviance:                       975.99\n",
      "Time:                        15:40:56   Pearson chi2:                 1.16e+03\n",
      "No. Iterations:                     5                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept      0.3455      0.731      0.472      0.637      -1.088       1.779\n",
      "won_lottry     3.5514      0.178     19.934      0.000       3.202       3.901\n",
      "male          -0.1622      0.164     -0.992      0.321      -0.483       0.158\n",
      "base_age      -0.1184      0.061     -1.946      0.052      -0.238       0.001\n",
      "==============================================================================\n",
      "\n",
      "\n",
      "==============================================================================\n",
      "                                  SECOND STAGE\n",
      "==============================================================================\n",
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:              finish8th   No. Observations:                 1171\n",
      "Model:                            GLM   Df Residuals:                     1167\n",
      "Model Family:                Binomial   Df Model:                            3\n",
      "Link Function:                  logit   Scale:                             1.0\n",
      "Method:                          IRLS   Log-Likelihood:                -696.65\n",
      "Date:                Wed, 27 Jun 2018   Deviance:                       1393.3\n",
      "Time:                        15:40:56   Pearson chi2:                 1.17e+03\n",
      "No. Iterations:                     4                                         \n",
      "======================================================================================\n",
      "                         coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------------\n",
      "Intercept              4.0756      0.604      6.753      0.000       2.893       5.258\n",
      "use_fin_aid_fitted     0.7743      0.192      4.036      0.000       0.398       1.150\n",
      "male                  -0.4175      0.130     -3.208      0.001      -0.673      -0.162\n",
      "base_age              -0.2919      0.048     -6.077      0.000      -0.386      -0.198\n",
      "======================================================================================\n"
     ]
    }
   ],
   "source": [
    "print(\"==============================================================================\")\n",
    "print(\"                                  FIRST STAGE\"                                 )\n",
    "print(\"==============================================================================\")\n",
    "result = smf.glm(formula = \"use_fin_aid ~ won_lottry + male + base_age\", \n",
    "                 data=voucher_df,\n",
    "                 family=sm.families.Binomial()).fit()\n",
    "voucher_df['use_fin_aid_fitted']= result.predict()\n",
    "print(result.summary())\n",
    "\n",
    "print()\n",
    "print()\n",
    "print(\"==============================================================================\")\n",
    "print(\"                                  SECOND STAGE\"                                )\n",
    "print(\"==============================================================================\")\n",
    "result = smf.glm(formula = \" finish8th ~ use_fin_aid_fitted + male + base_age\", \n",
    "                 data=voucher_df,\n",
    "                 family=sm.families.Binomial()).fit()\n",
    "print(result.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Those using financial aid are 2.168986 times more likely to finish the 8th grade!\n"
     ]
    }
   ],
   "source": [
    "# Transforming the log-odds units into odds ratio\n",
    "coefficient_to_odds_ratio(result.params[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting the Local Average Treatment Effect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we use IV to get a causal effect, what are our results actually telling us? According to Murnane and Willett, \"an estimate of a treatment effect obtained by IV methods should be regarded as an __estimated local average treatment effect (LATE)__. As mentioned earlier in the lecture slides, the LATE estimate can depend on your choice of instruments.\n",
    "\n",
    "As stated in the book:\n",
    "\n",
    "- Compliers \"are willing to have their behavior determined by the outcomes of the lottery, regardless of the particular experimental conditions to which they were assigned\" (278).\n",
    "- Always-Takers \"are families who will find and make use of financial aid to pay private-school fees\" regardless of the lottery. They may find aid outside the lottery\n",
    "- Never-takers are the mirror image of always-takers: \"they will not make use of financial aid to pay childrens' fees at a private secondary school under any circumstances\" (278)\n",
    "- (there are other possible groups, like \"defiers\" (Gennetian et all, 2005) who always do the opposite of what investigators ask them to do, but we make the assumption of \"no defiers\" in this dataset)\n",
    "\n",
    "In this context, IV estimates of the __local average treatment effect (LATE) for this quasi-experiment only applies to \"compliers\"__--and not to never-takers or always-takers.\n",
    "\n",
    "__Tests for instrument validity/strength of instruments is an established line of research that we can point you to if wanted!__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}