{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-03T03:21:00.563796Z",
     "start_time": "2019-09-03T03:20:54.260623Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pathlib\n",
    "from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score, mean_squared_error\n",
    "%matplotlib inline\n",
    "from plotnine import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.189356Z",
     "start_time": "2019-09-02T22:22:17.177815Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def predict_g(X, y):\n",
    "    clf=RandomForestClassifier(random_state=42, oob_score=True, n_estimators=10, max_depth=3, class_weight='balanced')\n",
    "    clf.fit(X,y)\n",
    "    print('oob_score: '+ str(clf.oob_score_))\n",
    "    return clf.predict_proba(X)[:,1]\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.210529Z",
     "start_time": "2019-09-02T22:22:17.192456Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def predict_Q(X,y):\n",
    "    reg=RandomForestRegressor(random_state=42, n_estimators=10, max_depth=3)\n",
    "    reg.fit(X,y)\n",
    "    return reg.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.230529Z",
     "start_time": "2019-09-02T22:22:17.210529Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def calc_delta(alpha, g, bias):\n",
    "    alpha_term=(1/alpha)-1\n",
    "    #check broadcasting\n",
    "    bias_term=sp.special.digamma(g*alpha_term+1)\\\n",
    "        -sp.special.digamma((1-g)*alpha_term)\\\n",
    "        -sp.special.digamma(g*alpha_term)\\\n",
    "        +sp.special.digamma((1-g)*alpha_term+1)\n",
    "    delta=bias/np.mean(bias_term)\n",
    "    return delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.258529Z",
     "start_time": "2019-09-02T22:22:17.234535Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def calc_beta_shapes(g, alpha):\n",
    "    alpha_shape=g*((1/alpha)-1)\n",
    "    beta_shape=(1-g)*((1/alpha)-1)\n",
    "    return alpha_shape, beta_shape\n",
    "    \n",
    "\n",
    "def calc_Rsq(alpha, delta, calculated_params, treatment, response):\n",
    "    alpha_shape, beta_shape=calc_beta_shapes(calculated_params['g'], alpha)\n",
    "    #TODO vv: does broadcasting work as expected for second argument of sp.special.polygamma ?\n",
    "    Rsq_num=delta**2*np.mean(sp.special.polygamma(1, alpha_shape+treatment)\\\n",
    "                             +sp.special.polygamma(1, beta_shape+(1-treatment)))\n",
    "    Rsq_den=mean_squared_error(response, calculated_params['Q'])\n",
    "    return Rsq_num/Rsq_den\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.274529Z",
     "start_time": "2019-09-02T22:22:17.258529Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def calc_Rsq_hat(variable_col, response_col, df, Q):\n",
    "    y=df[response_col].values\n",
    "#     if Q=='calc':\n",
    "#         Q=predict_Q(df.drop([response_col], axis=1).values, y)\n",
    "    Qhat=predict_Q(df.drop([response_col, variable_col], axis=1).values, y)\n",
    "    Rsqhat=(mean_squared_error(y, Qhat)-mean_squared_error(y, Q))/(mean_squared_error(y, Qhat))\n",
    "    return Rsqhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.294530Z",
     "start_time": "2019-09-02T22:22:17.278530Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def calc_alpha_hat(variable_col, response_col, treatment_col, df, g):\n",
    "    #call y t\n",
    "    y=df[treatment_col].values\n",
    "#     if g=='calc':\n",
    "#         g=predict_g(df.drop([response_col, treatment_col], axis=1), y)\n",
    "    ghat=predict_g(df.drop([response_col, treatment_col, variable_col], axis=1), y)\n",
    "    ahat_numerator=np.var(g)-np.var(ghat)\n",
    "    ahat_denominator=np.mean(ghat*(1-ghat))\n",
    "    ahat=ahat_numerator/ahat_denominator\n",
    "    return ahat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.330530Z",
     "start_time": "2019-09-02T22:22:17.302532Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def plot_sensitivity_graph(df, treatment_col, response_col, variable_cols, bias):\n",
    "    #Calculate g and Q\n",
    "    # TODO: vv what's going on here? why the drop?\n",
    "    g=predict_g(df.drop([response_col, treatment_col], axis=1).values, df[treatment_col].values)\n",
    "    Q=predict_Q(df.drop([response_col], axis=1).values, df[response_col].values)\n",
    "    calculated_params=pd.DataFrame({'g': g, 'Q': Q})\n",
    "    \n",
    "    #Calculate alpha, delta, Rsq\n",
    "    alpha_series=pd.Series(np.arange(0.01,1,0.01))\n",
    "    delta_series=alpha_series.apply(calc_delta, g=calculated_params['g'], bias=bias)\n",
    "    sensitivity_df=pd.DataFrame({'alphas':alpha_series, 'deltas':delta_series})\n",
    "    Rsq=[]\n",
    "    for n in range(len(alpha_series)):\n",
    "        Rsq.append(calc_Rsq(alpha_series[n], delta_series[n], calculated_params, df[treatment_col], df[response_col]))\n",
    "    sensitivity_df['Rsq']=Rsq\n",
    "    \n",
    "    #Plot observed co-variates\n",
    "    alpha_hat=[]\n",
    "    Rsq_hat=[]\n",
    "    \n",
    "    for var in variable_cols:\n",
    "        alpha_hat.append(calc_alpha_hat(var, response_col, treatment_col, df, calculated_params['g']))\n",
    "        Rsq_hat.append(calc_Rsq_hat(var, response_col, df, calculated_params['Q']))\n",
    "    \n",
    "    #call variable importances\n",
    "    observed_covariates=pd.DataFrame({'alpha_hat': alpha_hat, 'Rsq_hat':Rsq_hat})\n",
    "    #print(observed_covariates)\n",
    "    return (ggplot(data=sensitivity_df, mapping=aes(x='alphas', y='Rsq'))\n",
    "            +geom_line()\n",
    "            +theme_light()\n",
    "            +geom_point(data=observed_covariates, mapping=aes(x='alpha_hat', y='Rsq_hat'), color='red')\n",
    "            +labs(x='alpha')\n",
    "           )\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-03T03:21:00.623796Z",
     "start_time": "2019-09-03T03:21:00.615798Z"
    }
   },
   "outputs": [],
   "source": [
    "col_names=['treatment', 'age', 'education', 'black', 'hispanic', 'married', 'nodegree', 'RE74', 'RE75', 'RE78']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-03T03:21:20.113741Z",
     "start_time": "2019-09-03T03:21:20.030383Z"
    }
   },
   "outputs": [],
   "source": [
    "control=pd.read_csv(pathlib.Path.cwd().parents[0].joinpath('data', 'nswre74_control.txt'), header=None, sep=r\"\\s\\s\", names=col_names, engine='python')\n",
    "treatment=pd.read_csv(pathlib.Path.cwd().parents[0].joinpath('data', 'nswre74_treated.txt'), header=None, sep=r\"\\s\\s\", names=col_names, engine='python')\n",
    "psid=pd.read_csv(pathlib.Path.cwd().parents[0].joinpath('data', 'psid_controls.txt'), header=None, sep=r\"\\s\\s\", names=col_names, engine='python')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.482448Z",
     "start_time": "2019-09-02T22:22:17.470529Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(260, 10) (185, 10) (2490, 10)\n"
     ]
    }
   ],
   "source": [
    "print(control.shape, treatment.shape, psid.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.542441Z",
     "start_time": "2019-09-02T22:22:17.486442Z"
    }
   },
   "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>treatment</th>\n",
       "      <th>age</th>\n",
       "      <th>education</th>\n",
       "      <th>black</th>\n",
       "      <th>hispanic</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",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>0.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12383.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>10740.08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>0.0</td>\n",
       "      <td>45.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11796.47</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   treatment   age  education  black  hispanic  married  nodegree  RE74  RE75  \\\n",
       "0        0.0  23.0       10.0    1.0       0.0      0.0       1.0   0.0   0.0   \n",
       "1        0.0  26.0       12.0    0.0       0.0      0.0       0.0   0.0   0.0   \n",
       "2        0.0  22.0        9.0    1.0       0.0      0.0       1.0   0.0   0.0   \n",
       "3        0.0  18.0        9.0    1.0       0.0      0.0       1.0   0.0   0.0   \n",
       "4        0.0  45.0       11.0    1.0       0.0      0.0       1.0   0.0   0.0   \n",
       "\n",
       "       RE78  \n",
       "0      0.00  \n",
       "1  12383.68  \n",
       "2      0.00  \n",
       "3  10740.08  \n",
       "4  11796.47  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imbens_1=pd.concat([control,treatment]).reset_index(drop=True)\n",
    "imbens_1.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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>treatment</th>\n",
       "      <th>age</th>\n",
       "      <th>education</th>\n",
       "      <th>black</th>\n",
       "      <th>hispanic</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",
       "      <td>count</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.00000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "      <td>445.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>mean</td>\n",
       "      <td>0.415730</td>\n",
       "      <td>25.370787</td>\n",
       "      <td>10.195506</td>\n",
       "      <td>0.833708</td>\n",
       "      <td>0.08764</td>\n",
       "      <td>0.168539</td>\n",
       "      <td>0.782022</td>\n",
       "      <td>2102.265311</td>\n",
       "      <td>1377.138368</td>\n",
       "      <td>5300.763699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>std</td>\n",
       "      <td>0.493402</td>\n",
       "      <td>7.100282</td>\n",
       "      <td>1.792119</td>\n",
       "      <td>0.372762</td>\n",
       "      <td>0.28309</td>\n",
       "      <td>0.374766</td>\n",
       "      <td>0.413337</td>\n",
       "      <td>5363.582400</td>\n",
       "      <td>3150.960771</td>\n",
       "      <td>6631.491695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>min</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>17.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.00000</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",
       "      <td>25%</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>20.000000</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>50%</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>24.000000</td>\n",
       "      <td>10.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>3701.812000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>75%</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>28.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>824.388900</td>\n",
       "      <td>1220.836000</td>\n",
       "      <td>8124.715000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>max</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>55.000000</td>\n",
       "      <td>16.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>39570.680000</td>\n",
       "      <td>25142.240000</td>\n",
       "      <td>60307.930000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        treatment         age   education       black   hispanic     married  \\\n",
       "count  445.000000  445.000000  445.000000  445.000000  445.00000  445.000000   \n",
       "mean     0.415730   25.370787   10.195506    0.833708    0.08764    0.168539   \n",
       "std      0.493402    7.100282    1.792119    0.372762    0.28309    0.374766   \n",
       "min      0.000000   17.000000    3.000000    0.000000    0.00000    0.000000   \n",
       "25%      0.000000   20.000000    9.000000    1.000000    0.00000    0.000000   \n",
       "50%      0.000000   24.000000   10.000000    1.000000    0.00000    0.000000   \n",
       "75%      1.000000   28.000000   11.000000    1.000000    0.00000    0.000000   \n",
       "max      1.000000   55.000000   16.000000    1.000000    1.00000    1.000000   \n",
       "\n",
       "         nodegree          RE74          RE75          RE78  \n",
       "count  445.000000    445.000000    445.000000    445.000000  \n",
       "mean     0.782022   2102.265311   1377.138368   5300.763699  \n",
       "std      0.413337   5363.582400   3150.960771   6631.491695  \n",
       "min      0.000000      0.000000      0.000000      0.000000  \n",
       "25%      1.000000      0.000000      0.000000      0.000000  \n",
       "50%      1.000000      0.000000      0.000000   3701.812000  \n",
       "75%      1.000000    824.388900   1220.836000   8124.715000  \n",
       "max      1.000000  39570.680000  25142.240000  60307.930000  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imbens_1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.558441Z",
     "start_time": "2019-09-02T22:22:17.546444Z"
    }
   },
   "outputs": [],
   "source": [
    "imbens_2=pd.concat([treatment, psid]).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:17.622447Z",
     "start_time": "2019-09-02T22:22:17.562442Z"
    }
   },
   "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>treatment</th>\n",
       "      <th>age</th>\n",
       "      <th>education</th>\n",
       "      <th>black</th>\n",
       "      <th>hispanic</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",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>37.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9930.0460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3595.8940</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>1.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>24909.4500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7506.1460</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>1.0</td>\n",
       "      <td>33.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</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": [
       "   treatment   age  education  black  hispanic  married  nodegree  RE74  RE75  \\\n",
       "0        1.0  37.0       11.0    1.0       0.0      1.0       1.0   0.0   0.0   \n",
       "1        1.0  22.0        9.0    0.0       1.0      0.0       1.0   0.0   0.0   \n",
       "2        1.0  30.0       12.0    1.0       0.0      0.0       0.0   0.0   0.0   \n",
       "3        1.0  27.0       11.0    1.0       0.0      0.0       1.0   0.0   0.0   \n",
       "4        1.0  33.0        8.0    1.0       0.0      0.0       1.0   0.0   0.0   \n",
       "\n",
       "         RE78  \n",
       "0   9930.0460  \n",
       "1   3595.8940  \n",
       "2  24909.4500  \n",
       "3   7506.1460  \n",
       "4    289.7899  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imbens_2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:19.941377Z",
     "start_time": "2019-09-02T22:22:17.626442Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "oob_score: 0.5483146067415731\n",
      "oob_score: 0.5752808988764045\n",
      "oob_score: 0.5191011235955056\n",
      "oob_score: 0.5280898876404494\n",
      "oob_score: 0.5123595505617977\n",
      "oob_score: 0.503370786516854\n",
      "oob_score: 0.5191011235955056\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 17044 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ggplot: (-9223363244695129852)>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_sensitivity_graph(imbens_1.drop(['RE74', 'RE75'], axis=1), 'treatment', 'RE78', ['age', 'education', 'black', 'hispanic', 'married', 'nodegree'], 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "oob_score: 0.5438202247191011\n",
      "oob_score: 0.5595505617977528\n",
      "oob_score: 0.46741573033707867\n",
      "oob_score: 0.5303370786516854\n",
      "oob_score: 0.5370786516853933\n",
      "oob_score: 0.5573033707865168\n",
      "oob_score: 0.5730337078651685\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 16812 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ggplot: (-9223363244695186951)>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_sensitivity_graph(imbens_1, 'treatment', 'RE78', ['age', 'education', 'black', 'hispanic', 'married', 'nodegree'], 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-02T22:22:21.449561Z",
     "start_time": "2019-09-02T22:22:19.941377Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "oob_score: 0.8762616822429906\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:458: UserWarning: Some inputs do not have OOB scores. This probably means too few trees were used to compute any reliable oob estimates.\n",
      "  warn(\"Some inputs do not have OOB scores. \"\n",
      "/home/victor/anaconda3/envs/sns/lib/python3.7/site-packages/sklearn/ensemble/forest.py:463: RuntimeWarning: invalid value encountered in true_divide\n",
      "  predictions[k].sum(axis=1)[:, np.newaxis])\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "oob_score: 0.8871028037383177\n",
      "oob_score: 0.8904672897196262\n",
      "oob_score: 0.8557009345794393\n",
      "oob_score: 0.8736448598130842\n",
      "oob_score: 0.7958878504672897\n",
      "oob_score: 0.8485981308411215\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 17032 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ggplot: (-9223363244695263060)>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_sensitivity_graph(imbens_2.drop(['RE74', 'RE75'], axis=1), 'treatment', 'RE78', ['age', 'education', 'black', 'hispanic', 'married', 'nodegree'], 1000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
