Application with 2019 ASEC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import lightgbm as lgb\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import confusion_matrix, roc_auc_score, f1_score, precision_recall_curve\n",
    "from DiscriminationMitigation import DiscriminationMitigator\n",
    "\n",
    "pd.set_option('display.max_columns', 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binary_metrics(y_true, y_pred):\n",
    "    '''\n",
    "    Calculates binary classification performance metrics for a given model.\n",
    "    :param y_true: array_like, truth values as int\n",
    "    :param y_pred: array_like, predicted values as int\n",
    "    :returns: dict, with keys for each metric: \n",
    "        accuracy - proportion of correct predictions out of total predictions\n",
    "        sensitivity - (aka recall), of all true positives reviews how many did we correctly predict as positive\n",
    "        specificity - (aka selectivity/TNR), of all true negatives how many did we correctly predict as negative\n",
    "        precision - of all predicted positive cases how many were actually positive\n",
    "        F-1 score - harmonic/weighted mean of precision and sensitivity scores\n",
    "        ROC-AUC - area under receiver operating characteristic curve\n",
    "        \n",
    "    '''\n",
    "    tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n",
    "    \n",
    "    metrics = {}\n",
    "    metrics['accuracy'] = round((tp + tn) / len(y_true), 4)\n",
    "    metrics['sensitivity/recall'] = round(tp / (fn + tp), 4) # aka recall\n",
    "    metrics['specificity'] = round(tn / (tn + fp), 4) # aka TNR\n",
    "    metrics['precision'] = round(tp / (tp + fp), 4)\n",
    "    metrics['f1'] = round(2 * (metrics['precision'] * metrics['sensitivity/recall']) \\\n",
    "                        / (metrics['precision'] + metrics['sensitivity/recall']), 4)\n",
    "    metrics['roc_auc'] = round(roc_auc_score(y_true, y_pred), 4)\n",
    "    \n",
    "    return metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def continuous_metrics(y_true, y_pred):\n",
    "    '''\n",
    "    Calculates performance metrics for a continuous outcome from a model.\n",
    "    :param y_true: array_like, truth values as float\n",
    "    :param y_pred: array_like, predicted values as float\n",
    "    :returns: dict, with keys for each metric:\n",
    "        mae - Mean Absolute Error\n",
    "        mse - Mean Squared Error\n",
    "        rmse - Root Mean Squared Error\n",
    "        r2 - Coefficient of Determination (r-squared)\n",
    "    '''\n",
    "    metrics = {}\n",
    "    metrics['mae'] = round(sum(abs(y_true - y_pred)) / len(y_true), 4)\n",
    "    metrics['mse'] = round(sum((y_true - y_pred)**2) / len(y_true), 4)\n",
    "    metrics['rmse'] = round(np.sqrt(metrics['mse']), 4)\n",
    "    metrics['r2'] = round(1 - (sum((y_true - y_pred)**2) / sum((y_true.mean() - y_true)**2)), 4)\n",
    "    \n",
    "    return metrics    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def combine_prediction(X_test, y_test, pred, outcome='50k'):\n",
    "    '''\n",
    "    Combines X_test dataframe with np.array of predictions\n",
    "    :param X_test: pd.DataFrame, X_test dataset\n",
    "    :param pred: np.array, predicted y values\n",
    "    :param outcome: str, column name to rename true y values\n",
    "    :returns:\n",
    "        pd.DataFrame\n",
    "    '''\n",
    "    X_test_combined = pd.concat([X_test, pd.DataFrame(y_test, columns=[outcome])], axis=1)\n",
    "    return pd.concat([X_test_combined.reset_index(drop=True), pd.DataFrame(pred, columns=['pred'])], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in [2019 ASEC (March CPS)](https://cps.ipums.org/cps/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('./data/asec_2019.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(180101, 18)\n"
     ]
    }
   ],
   "source": [
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   YEAR  MONTH  ASECFLAG   ASECWT  RELATE  AGE  SEX  RACE  MARST  POPSTAT  \\\n",
      "0  2019      3         1  2031.67     101   21    1   100      6        1   \n",
      "1  2019      3         1  1232.04     101   85    2   100      5        1   \n",
      "2  2019      3         1  1209.17     101   61    2   100      6        1   \n",
      "3  2019      3         1  1146.23     101   73    2   100      4        1   \n",
      "4  2019      3         1  1480.79     301   37    1   100      6        1   \n",
      "\n",
      "   HISPAN  LABFORCE  UHRSWORK1  EDUC  SCHLCOLL  WKSWORK1  UHRSWORKLY  INCWAGE  \n",
      "0       0         2       30.0    60       5.0        52          30  18000.0  \n",
      "1       0         1      999.0    73       0.0         0         999      0.0  \n",
      "2       0         2       44.0    73       0.0        52          44  12000.0  \n",
      "3       0         1      999.0    73       0.0         0         999      0.0  \n",
      "4       0         2       20.0    73       5.0        52          20  12000.0  \n"
     ]
    }
   ],
   "source": [
    "print(df.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Restrict dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop individuals with no employment or earnings last year\n",
    "df = df.loc[df['INCWAGE'] > 0] # No wage/salary last year\n",
    "df = df.loc[(df['WKSWORK1'] > 0) & (df['WKSWORK1'] <= 52)] # No weeks worked last year\n",
    "df = df.loc[(df['UHRSWORKLY'] > 0) & (df['UHRSWORKLY'] < 999)] # No usual hours/week worked last year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(85644, 18)\n"
     ]
    }
   ],
   "source": [
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to individuals aged 18-64\n",
    "df = df.loc[(df['AGE'] >=18) & (df['AGE'] <=64)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(78644, 18)\n"
     ]
    }
   ],
   "source": [
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to non-Hispanics\n",
    "df = df.loc[df['HISPAN'] == 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(63125, 18)\n"
     ]
    }
   ],
   "source": [
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to non-mixed-race Blacks and Whites\n",
    "df = df.loc[df.RACE.isin([100, 200])]\n",
    "df['blk'] = np.where(df['RACE'] == 200, 1, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(55508, 19)\n"
     ]
    }
   ],
   "source": [
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to adult civilians\n",
    "df = df.loc[df['POPSTAT'] == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(55069, 19)\n"
     ]
    }
   ],
   "source": [
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Engineer features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flag for part-time usual work\n",
    "df['pt'] = np.where(df['UHRSWORKLY'] < 35, 1, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Non-linearities for age\n",
    "df['AGE2'] = df['AGE'] ** 2\n",
    "df['AGE3'] = df['AGE'] ** 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure education coded correctly\n",
    "df = df.loc[df['EDUC'] <= 125] # 999 is missing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hourly wage\n",
    "df['hrwage'] = df['INCWAGE'] / (df['WKSWORK1'] * df['UHRSWORKLY'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    55069.000000\n",
       "mean        31.269112\n",
       "std        184.416973\n",
       "min          0.000962\n",
       "25%         14.375000\n",
       "50%         21.634615\n",
       "75%         34.188034\n",
       "max      25481.250000\n",
       "Name: hrwage, dtype: float64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['hrwage'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to people earning 1 < hr_wage < 100\n",
    "df = df.loc[(df['hrwage'] > 1) & (df['hrwage'] < 100)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "count    53671.000000\n",
      "mean        25.907302\n",
      "std         17.049376\n",
      "min          1.016667\n",
      "25%         14.202864\n",
      "50%         21.634615\n",
      "75%         33.333333\n",
      "max         99.759615\n",
      "Name: hrwage, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(df['hrwage'].describe())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Log hourly wage\n",
    "df['lnwage'] = np.log(df['hrwage'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count     53671.000000\n",
       "mean      52682.839876\n",
       "std       41242.872191\n",
       "min          25.000000\n",
       "25%       24617.500000\n",
       "50%       43680.000000\n",
       "75%       70000.000000\n",
       "max      420000.000000\n",
       "Name: INCWAGE, dtype: float64"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['INCWAGE'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flag for whether total earnings > 50,000 or not\n",
    "df['50k'] = np.where(df['INCWAGE'] > 50000, 1, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[['lnwage', 'hrwage', '50k', 'pt', 'INCWAGE', 'WKSWORK1', 'UHRSWORKLY', \n",
    "         'AGE', 'AGE2', 'AGE3', 'SEX', 'blk', 'MARST', 'SCHLCOLL', 'EDUC']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     lnwage     hrwage  50k  pt  INCWAGE  WKSWORK1  UHRSWORKLY  AGE  AGE2  \\\n",
      "0  2.445686  11.538462    0   1  18000.0        52          30   21   441   \n",
      "2  1.657229   5.244755    0   0  12000.0        52          44   61  3721   \n",
      "4  2.445686  11.538462    0   1  12000.0        52          20   37  1369   \n",
      "6  3.179655  24.038462    1   0  55000.0        52          44   53  2809   \n",
      "8  3.442059  31.251250    1   0  50002.0        40          40   62  3844   \n",
      "\n",
      "     AGE3  SEX  blk  MARST  SCHLCOLL  EDUC  \n",
      "0    9261    1    0      6       5.0    60  \n",
      "2  226981    2    0      6       0.0    73  \n",
      "4   50653    1    0      6       5.0    73  \n",
      "6  148877    2    0      4       5.0    73  \n",
      "8  238328    1    0      1       0.0   111  \n"
     ]
    }
   ],
   "source": [
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "categorical_features = ['SEX', 'MARST', 'SCHLCOLL', 'EDUC', 'pt', 'blk']\n",
    "for col in categorical_features:\n",
    "    df[col] = df[col].astype('category')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lnwage         float64\n",
       "hrwage         float64\n",
       "50k              int32\n",
       "pt            category\n",
       "INCWAGE        float64\n",
       "WKSWORK1         int64\n",
       "UHRSWORKLY       int64\n",
       "AGE              int64\n",
       "AGE2             int64\n",
       "AGE3             int64\n",
       "SEX           category\n",
       "blk           category\n",
       "MARST         category\n",
       "SCHLCOLL      category\n",
       "EDUC          category\n",
       "dtype: object"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Descriptive differences in earnings by race "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Association between race and likelihood of annual earnings surpassing 50K (i.e. \"high-income\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 1: Linear probability model regressing high/low earnings dummy on race dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>50k</td>       <th>  R-squared:         </th> <td>   0.012</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.012</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   738.0</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>2.04e-161</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:44</td>     <th>  Log-Likelihood:    </th> <td> -37646.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td> 53671</td>      <th>  AIC:               </th> <td>7.530e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 53669</td>      <th>  BIC:               </th> <td>7.531e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.4276</td> <td>    0.002</td> <td>  184.321</td> <td> 0.000</td> <td>    0.423</td> <td>    0.432</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.1486</td> <td>    0.005</td> <td>  -27.166</td> <td> 0.000</td> <td>   -0.159</td> <td>   -0.138</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>208719.027</td> <th>  Durbin-Watson:     </th> <td>   1.819</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>   <td> 0.000</td>   <th>  Jarque-Bera (JB):  </th> <td>8611.430</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>            <td> 0.375</td>   <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>        <td> 1.187</td>   <th>  Cond. No.          </th> <td>    2.85</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                    50k   R-squared:                       0.012\n",
       "Model:                            OLS   Adj. R-squared:                  0.012\n",
       "Method:                 Least Squares   F-statistic:                     738.0\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):          2.04e-161\n",
       "Time:                        10:04:44   Log-Likelihood:                -37646.\n",
       "No. Observations:               53671   AIC:                         7.530e+04\n",
       "Df Residuals:                   53669   BIC:                         7.531e+04\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.4276      0.002    184.321      0.000       0.423       0.432\n",
       "blk           -0.1486      0.005    -27.166      0.000      -0.159      -0.138\n",
       "==============================================================================\n",
       "Omnibus:                   208719.027   Durbin-Watson:                   1.819\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             8611.430\n",
       "Skew:                           0.375   Prob(JB):                         0.00\n",
       "Kurtosis:                       1.187   Cond. No.                         2.85\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = df['50k'].reset_index(drop=True)\n",
    "predictors = df['blk'].astype('int').reset_index(drop=True)\n",
    "predictors = sm.add_constant(predictors) # add constant\n",
    "sm.OLS(y, predictors, hasconst=True).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to Whites, Blacks have a 14.86% lower likelihood of earning >$50k per year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 2: Linear probability model regressing high/low earnings dummy on race dummy and control variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>50k</td>       <th>  R-squared:         </th> <td>   0.285</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.285</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   1259.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:45</td>     <th>  Log-Likelihood:    </th> <td> -28954.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td> 53671</td>      <th>  AIC:               </th> <td>5.797e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 53639</td>      <th>  BIC:               </th> <td>5.826e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    31</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "     <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>  <td>    0.3332</td> <td>    0.080</td> <td>    4.148</td> <td> 0.000</td> <td>    0.176</td> <td>    0.491</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>    <td>   -0.0733</td> <td>    0.005</td> <td>  -14.478</td> <td> 0.000</td> <td>   -0.083</td> <td>   -0.063</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AGE</th>    <td>   -0.0193</td> <td>    0.005</td> <td>   -3.989</td> <td> 0.000</td> <td>   -0.029</td> <td>   -0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AGE2</th>   <td>    0.0009</td> <td>    0.000</td> <td>    6.780</td> <td> 0.000</td> <td>    0.001</td> <td>    0.001</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AGE3</th>   <td>-8.663e-06</td> <td> 1.08e-06</td> <td>   -8.042</td> <td> 0.000</td> <td>-1.08e-05</td> <td>-6.55e-06</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>pt</th>     <td>   -0.2688</td> <td>    0.004</td> <td>  -66.484</td> <td> 0.000</td> <td>   -0.277</td> <td>   -0.261</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_SEX</th>  <td>   -0.1804</td> <td>    0.004</td> <td>  -49.260</td> <td> 0.000</td> <td>   -0.188</td> <td>   -0.173</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_MAR</th>  <td>   -0.0522</td> <td>    0.017</td> <td>   -2.995</td> <td> 0.003</td> <td>   -0.086</td> <td>   -0.018</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>1_MAR</th>  <td>   -0.0931</td> <td>    0.013</td> <td>   -6.923</td> <td> 0.000</td> <td>   -0.119</td> <td>   -0.067</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>2_MAR</th>  <td>   -0.0616</td> <td>    0.006</td> <td>   -9.560</td> <td> 0.000</td> <td>   -0.074</td> <td>   -0.049</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>3_MAR</th>  <td>   -0.0793</td> <td>    0.016</td> <td>   -4.995</td> <td> 0.000</td> <td>   -0.110</td> <td>   -0.048</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>4_MAR</th>  <td>   -0.0910</td> <td>    0.005</td> <td>  -17.442</td> <td> 0.000</td> <td>   -0.101</td> <td>   -0.081</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_SCH</th>  <td>    0.2098</td> <td>    0.018</td> <td>   11.639</td> <td> 0.000</td> <td>    0.174</td> <td>    0.245</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>1_SCH</th>  <td>    0.0819</td> <td>    0.045</td> <td>    1.810</td> <td> 0.070</td> <td>   -0.007</td> <td>    0.171</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>2_SCH</th>  <td>   -0.0258</td> <td>    0.013</td> <td>   -1.949</td> <td> 0.051</td> <td>   -0.052</td> <td>    0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>3_SCH</th>  <td>   -0.0316</td> <td>    0.016</td> <td>   -1.936</td> <td> 0.053</td> <td>   -0.064</td> <td>    0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>4_SCH</th>  <td>    0.0021</td> <td>    0.011</td> <td>    0.193</td> <td> 0.847</td> <td>   -0.019</td> <td>    0.023</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_EDU</th>  <td>    0.1145</td> <td>    0.102</td> <td>    1.122</td> <td> 0.262</td> <td>   -0.086</td> <td>    0.315</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>1_EDU</th>  <td>   -0.0084</td> <td>    0.082</td> <td>   -0.103</td> <td> 0.918</td> <td>   -0.170</td> <td>    0.153</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>2_EDU</th>  <td>   -0.0116</td> <td>    0.060</td> <td>   -0.195</td> <td> 0.846</td> <td>   -0.128</td> <td>    0.105</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>3_EDU</th>  <td>   -0.1063</td> <td>    0.055</td> <td>   -1.929</td> <td> 0.054</td> <td>   -0.214</td> <td>    0.002</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>4_EDU</th>  <td>   -0.0216</td> <td>    0.054</td> <td>   -0.400</td> <td> 0.689</td> <td>   -0.128</td> <td>    0.084</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>5_EDU</th>  <td>   -0.0028</td> <td>    0.053</td> <td>   -0.052</td> <td> 0.958</td> <td>   -0.107</td> <td>    0.101</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>6_EDU</th>  <td>   -0.0227</td> <td>    0.054</td> <td>   -0.421</td> <td> 0.674</td> <td>   -0.128</td> <td>    0.083</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>7_EDU</th>  <td>    0.0583</td> <td>    0.052</td> <td>    1.128</td> <td> 0.259</td> <td>   -0.043</td> <td>    0.160</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>8_EDU</th>  <td>    0.1411</td> <td>    0.052</td> <td>    2.726</td> <td> 0.006</td> <td>    0.040</td> <td>    0.243</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>9_EDU</th>  <td>    0.1733</td> <td>    0.052</td> <td>    3.316</td> <td> 0.001</td> <td>    0.071</td> <td>    0.276</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>10_EDU</th> <td>    0.1885</td> <td>    0.052</td> <td>    3.618</td> <td> 0.000</td> <td>    0.086</td> <td>    0.291</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>11_EDU</th> <td>    0.3491</td> <td>    0.052</td> <td>    6.744</td> <td> 0.000</td> <td>    0.248</td> <td>    0.451</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>12_EDU</th> <td>    0.4788</td> <td>    0.052</td> <td>    9.220</td> <td> 0.000</td> <td>    0.377</td> <td>    0.581</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>13_EDU</th> <td>    0.5770</td> <td>    0.053</td> <td>   10.830</td> <td> 0.000</td> <td>    0.473</td> <td>    0.681</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>14_EDU</th> <td>    0.5637</td> <td>    0.053</td> <td>   10.658</td> <td> 0.000</td> <td>    0.460</td> <td>    0.667</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>19550.706</td> <th>  Durbin-Watson:     </th> <td>   1.892</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>   <th>  Jarque-Bera (JB):  </th> <td>2717.034</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td> 0.112</td>   <th>  Prob(JB):          </th> <td>    0.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 1.921</td>   <th>  Cond. No.          </th> <td>1.88e+07</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)<br/>[2] The condition number is large, 1.88e+07. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                    50k   R-squared:                       0.285\n",
       "Model:                            OLS   Adj. R-squared:                  0.285\n",
       "Method:                 Least Squares   F-statistic:                     1259.\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):               0.00\n",
       "Time:                        10:04:45   Log-Likelihood:                -28954.\n",
       "No. Observations:               53671   AIC:                         5.797e+04\n",
       "Df Residuals:                   53639   BIC:                         5.826e+04\n",
       "Df Model:                          31                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.3332      0.080      4.148      0.000       0.176       0.491\n",
       "blk           -0.0733      0.005    -14.478      0.000      -0.083      -0.063\n",
       "AGE           -0.0193      0.005     -3.989      0.000      -0.029      -0.010\n",
       "AGE2           0.0009      0.000      6.780      0.000       0.001       0.001\n",
       "AGE3       -8.663e-06   1.08e-06     -8.042      0.000   -1.08e-05   -6.55e-06\n",
       "pt            -0.2688      0.004    -66.484      0.000      -0.277      -0.261\n",
       "0_SEX         -0.1804      0.004    -49.260      0.000      -0.188      -0.173\n",
       "0_MAR         -0.0522      0.017     -2.995      0.003      -0.086      -0.018\n",
       "1_MAR         -0.0931      0.013     -6.923      0.000      -0.119      -0.067\n",
       "2_MAR         -0.0616      0.006     -9.560      0.000      -0.074      -0.049\n",
       "3_MAR         -0.0793      0.016     -4.995      0.000      -0.110      -0.048\n",
       "4_MAR         -0.0910      0.005    -17.442      0.000      -0.101      -0.081\n",
       "0_SCH          0.2098      0.018     11.639      0.000       0.174       0.245\n",
       "1_SCH          0.0819      0.045      1.810      0.070      -0.007       0.171\n",
       "2_SCH         -0.0258      0.013     -1.949      0.051      -0.052       0.000\n",
       "3_SCH         -0.0316      0.016     -1.936      0.053      -0.064       0.000\n",
       "4_SCH          0.0021      0.011      0.193      0.847      -0.019       0.023\n",
       "0_EDU          0.1145      0.102      1.122      0.262      -0.086       0.315\n",
       "1_EDU         -0.0084      0.082     -0.103      0.918      -0.170       0.153\n",
       "2_EDU         -0.0116      0.060     -0.195      0.846      -0.128       0.105\n",
       "3_EDU         -0.1063      0.055     -1.929      0.054      -0.214       0.002\n",
       "4_EDU         -0.0216      0.054     -0.400      0.689      -0.128       0.084\n",
       "5_EDU         -0.0028      0.053     -0.052      0.958      -0.107       0.101\n",
       "6_EDU         -0.0227      0.054     -0.421      0.674      -0.128       0.083\n",
       "7_EDU          0.0583      0.052      1.128      0.259      -0.043       0.160\n",
       "8_EDU          0.1411      0.052      2.726      0.006       0.040       0.243\n",
       "9_EDU          0.1733      0.052      3.316      0.001       0.071       0.276\n",
       "10_EDU         0.1885      0.052      3.618      0.000       0.086       0.291\n",
       "11_EDU         0.3491      0.052      6.744      0.000       0.248       0.451\n",
       "12_EDU         0.4788      0.052      9.220      0.000       0.377       0.581\n",
       "13_EDU         0.5770      0.053     10.830      0.000       0.473       0.681\n",
       "14_EDU         0.5637      0.053     10.658      0.000       0.460       0.667\n",
       "==============================================================================\n",
       "Omnibus:                    19550.706   Durbin-Watson:                   1.892\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             2717.034\n",
       "Skew:                           0.112   Prob(JB):                         0.00\n",
       "Kurtosis:                       1.921   Cond. No.                     1.88e+07\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "[2] The condition number is large, 1.88e+07. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = df['50k'].reset_index(drop=True)\n",
    "predictors = df[['blk', 'AGE', 'AGE2', 'AGE3', 'pt']].reset_index(drop=True)\n",
    "predictors['blk'] = predictors['blk'].astype(int)\n",
    "predictors['pt'] = predictors['pt'].astype(int)\n",
    "for col in ['SEX', 'MARST', 'SCHLCOLL', 'EDUC']:\n",
    "    onehot_vector = pd.DataFrame(OneHotEncoder().fit_transform(df[[col]]).toarray()[:, 1:])\n",
    "    onehot_vector.columns = [str(i)+'_'+col[:3] for i in onehot_vector.columns]\n",
    "    predictors = pd.concat([predictors, onehot_vector], axis=1)\n",
    "predictors = sm.add_constant(predictors) # add constant\n",
    "\n",
    "sm.OLS(y, predictors, hasconst=True).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Holding other factors in the model constant, Blacks have about an 7.3% lower probability of earning >$50k per year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Association between race and log hourly wage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 3: OLS regression of log hourly wages including only race dummy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y = df['lnwage'].reset_index(drop=True)\n",
    "predictors = df['blk'].astype('int').reset_index(drop=True)\n",
    "predictors = sm.add_constant(predictors) # add constant\n",
    "sm.OLS(y, predictors, hasconst=True).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to Whites, Blacks earn on average 19.75% (0.1975 = 1-exp(-0.22)) per hour worked."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 4: OLS regression of log hourly wages using race dummy and control variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>lnwage</td>      <th>  R-squared:         </th> <td>   0.296</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.296</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   665.6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:45</td>     <th>  Log-Likelihood:    </th> <td> -45032.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td> 53671</td>      <th>  AIC:               </th> <td>9.013e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td> 53639</td>      <th>  BIC:               </th> <td>9.041e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>    31</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "     <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>  <td>    1.7981</td> <td>    0.142</td> <td>   12.662</td> <td> 0.000</td> <td>    1.520</td> <td>    2.076</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>    <td>   -0.1094</td> <td>    0.007</td> <td>  -15.589</td> <td> 0.000</td> <td>   -0.123</td> <td>   -0.096</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AGE</th>    <td>    0.0582</td> <td>    0.008</td> <td>    7.478</td> <td> 0.000</td> <td>    0.043</td> <td>    0.073</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AGE2</th>   <td>   -0.0008</td> <td>    0.000</td> <td>   -3.974</td> <td> 0.000</td> <td>   -0.001</td> <td>   -0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AGE3</th>   <td> 2.664e-06</td> <td> 1.62e-06</td> <td>    1.640</td> <td> 0.101</td> <td>-5.19e-07</td> <td> 5.85e-06</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>pt</th>     <td>   -0.1937</td> <td>    0.009</td> <td>  -22.032</td> <td> 0.000</td> <td>   -0.211</td> <td>   -0.176</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_SEX</th>  <td>   -0.2118</td> <td>    0.005</td> <td>  -42.615</td> <td> 0.000</td> <td>   -0.222</td> <td>   -0.202</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_MAR</th>  <td>   -0.1029</td> <td>    0.024</td> <td>   -4.304</td> <td> 0.000</td> <td>   -0.150</td> <td>   -0.056</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>1_MAR</th>  <td>   -0.1771</td> <td>    0.018</td> <td>   -9.582</td> <td> 0.000</td> <td>   -0.213</td> <td>   -0.141</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>2_MAR</th>  <td>   -0.0999</td> <td>    0.008</td> <td>  -12.037</td> <td> 0.000</td> <td>   -0.116</td> <td>   -0.084</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>3_MAR</th>  <td>   -0.1464</td> <td>    0.023</td> <td>   -6.490</td> <td> 0.000</td> <td>   -0.191</td> <td>   -0.102</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>4_MAR</th>  <td>   -0.1334</td> <td>    0.007</td> <td>  -19.561</td> <td> 0.000</td> <td>   -0.147</td> <td>   -0.120</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_SCH</th>  <td>    0.0819</td> <td>    0.040</td> <td>    2.027</td> <td> 0.043</td> <td>    0.003</td> <td>    0.161</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>1_SCH</th>  <td>    0.0736</td> <td>    0.135</td> <td>    0.546</td> <td> 0.585</td> <td>   -0.190</td> <td>    0.338</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>2_SCH</th>  <td>   -0.0823</td> <td>    0.021</td> <td>   -4.003</td> <td> 0.000</td> <td>   -0.123</td> <td>   -0.042</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>3_SCH</th>  <td>   -0.0635</td> <td>    0.023</td> <td>   -2.796</td> <td> 0.005</td> <td>   -0.108</td> <td>   -0.019</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>4_SCH</th>  <td>   -0.0218</td> <td>    0.014</td> <td>   -1.537</td> <td> 0.124</td> <td>   -0.050</td> <td>    0.006</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>0_EDU</th>  <td>   -0.0803</td> <td>    0.144</td> <td>   -0.556</td> <td> 0.578</td> <td>   -0.363</td> <td>    0.203</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>1_EDU</th>  <td>   -0.0342</td> <td>    0.133</td> <td>   -0.257</td> <td> 0.797</td> <td>   -0.295</td> <td>    0.226</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>2_EDU</th>  <td>   -0.1126</td> <td>    0.110</td> <td>   -1.023</td> <td> 0.306</td> <td>   -0.329</td> <td>    0.103</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>3_EDU</th>  <td>   -0.2498</td> <td>    0.104</td> <td>   -2.391</td> <td> 0.017</td> <td>   -0.454</td> <td>   -0.045</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>4_EDU</th>  <td>   -0.1348</td> <td>    0.101</td> <td>   -1.331</td> <td> 0.183</td> <td>   -0.333</td> <td>    0.064</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>5_EDU</th>  <td>   -0.0778</td> <td>    0.099</td> <td>   -0.783</td> <td> 0.433</td> <td>   -0.273</td> <td>    0.117</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>6_EDU</th>  <td>   -0.0978</td> <td>    0.101</td> <td>   -0.964</td> <td> 0.335</td> <td>   -0.297</td> <td>    0.101</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>7_EDU</th>  <td>    0.0561</td> <td>    0.097</td> <td>    0.579</td> <td> 0.563</td> <td>   -0.134</td> <td>    0.246</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>8_EDU</th>  <td>    0.1750</td> <td>    0.097</td> <td>    1.803</td> <td> 0.071</td> <td>   -0.015</td> <td>    0.365</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>9_EDU</th>  <td>    0.2232</td> <td>    0.097</td> <td>    2.291</td> <td> 0.022</td> <td>    0.032</td> <td>    0.414</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>10_EDU</th> <td>    0.2505</td> <td>    0.097</td> <td>    2.574</td> <td> 0.010</td> <td>    0.060</td> <td>    0.441</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>11_EDU</th> <td>    0.4950</td> <td>    0.097</td> <td>    5.102</td> <td> 0.000</td> <td>    0.305</td> <td>    0.685</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>12_EDU</th> <td>    0.6449</td> <td>    0.097</td> <td>    6.638</td> <td> 0.000</td> <td>    0.454</td> <td>    0.835</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>13_EDU</th> <td>    0.8200</td> <td>    0.099</td> <td>    8.261</td> <td> 0.000</td> <td>    0.625</td> <td>    1.015</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>14_EDU</th> <td>    0.8030</td> <td>    0.099</td> <td>    8.143</td> <td> 0.000</td> <td>    0.610</td> <td>    0.996</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>4391.268</td> <th>  Durbin-Watson:     </th> <td>   1.852</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>10292.104</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td>-0.509</td>  <th>  Prob(JB):          </th> <td>    0.00</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 4.889</td>  <th>  Cond. No.          </th> <td>1.88e+07</td> \n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)<br/>[2] The condition number is large, 1.88e+07. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 lnwage   R-squared:                       0.296\n",
       "Model:                            OLS   Adj. R-squared:                  0.296\n",
       "Method:                 Least Squares   F-statistic:                     665.6\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):               0.00\n",
       "Time:                        10:04:45   Log-Likelihood:                -45032.\n",
       "No. Observations:               53671   AIC:                         9.013e+04\n",
       "Df Residuals:                   53639   BIC:                         9.041e+04\n",
       "Df Model:                          31                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          1.7981      0.142     12.662      0.000       1.520       2.076\n",
       "blk           -0.1094      0.007    -15.589      0.000      -0.123      -0.096\n",
       "AGE            0.0582      0.008      7.478      0.000       0.043       0.073\n",
       "AGE2          -0.0008      0.000     -3.974      0.000      -0.001      -0.000\n",
       "AGE3        2.664e-06   1.62e-06      1.640      0.101   -5.19e-07    5.85e-06\n",
       "pt            -0.1937      0.009    -22.032      0.000      -0.211      -0.176\n",
       "0_SEX         -0.2118      0.005    -42.615      0.000      -0.222      -0.202\n",
       "0_MAR         -0.1029      0.024     -4.304      0.000      -0.150      -0.056\n",
       "1_MAR         -0.1771      0.018     -9.582      0.000      -0.213      -0.141\n",
       "2_MAR         -0.0999      0.008    -12.037      0.000      -0.116      -0.084\n",
       "3_MAR         -0.1464      0.023     -6.490      0.000      -0.191      -0.102\n",
       "4_MAR         -0.1334      0.007    -19.561      0.000      -0.147      -0.120\n",
       "0_SCH          0.0819      0.040      2.027      0.043       0.003       0.161\n",
       "1_SCH          0.0736      0.135      0.546      0.585      -0.190       0.338\n",
       "2_SCH         -0.0823      0.021     -4.003      0.000      -0.123      -0.042\n",
       "3_SCH         -0.0635      0.023     -2.796      0.005      -0.108      -0.019\n",
       "4_SCH         -0.0218      0.014     -1.537      0.124      -0.050       0.006\n",
       "0_EDU         -0.0803      0.144     -0.556      0.578      -0.363       0.203\n",
       "1_EDU         -0.0342      0.133     -0.257      0.797      -0.295       0.226\n",
       "2_EDU         -0.1126      0.110     -1.023      0.306      -0.329       0.103\n",
       "3_EDU         -0.2498      0.104     -2.391      0.017      -0.454      -0.045\n",
       "4_EDU         -0.1348      0.101     -1.331      0.183      -0.333       0.064\n",
       "5_EDU         -0.0778      0.099     -0.783      0.433      -0.273       0.117\n",
       "6_EDU         -0.0978      0.101     -0.964      0.335      -0.297       0.101\n",
       "7_EDU          0.0561      0.097      0.579      0.563      -0.134       0.246\n",
       "8_EDU          0.1750      0.097      1.803      0.071      -0.015       0.365\n",
       "9_EDU          0.2232      0.097      2.291      0.022       0.032       0.414\n",
       "10_EDU         0.2505      0.097      2.574      0.010       0.060       0.441\n",
       "11_EDU         0.4950      0.097      5.102      0.000       0.305       0.685\n",
       "12_EDU         0.6449      0.097      6.638      0.000       0.454       0.835\n",
       "13_EDU         0.8200      0.099      8.261      0.000       0.625       1.015\n",
       "14_EDU         0.8030      0.099      8.143      0.000       0.610       0.996\n",
       "==============================================================================\n",
       "Omnibus:                     4391.268   Durbin-Watson:                   1.852\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):            10292.104\n",
       "Skew:                          -0.509   Prob(JB):                         0.00\n",
       "Kurtosis:                       4.889   Cond. No.                     1.88e+07\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "[2] The condition number is large, 1.88e+07. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = df['lnwage'].reset_index(drop=True)\n",
    "predictors = df[['blk', 'AGE', 'AGE2', 'AGE3', 'pt']].reset_index(drop=True)\n",
    "predictors['blk'] = predictors['blk'].astype(int)\n",
    "predictors['pt'] = predictors['pt'].astype(int)\n",
    "for col in ['SEX', 'MARST', 'SCHLCOLL', 'EDUC']:\n",
    "    onehot_vector = pd.DataFrame(OneHotEncoder().fit_transform(df[[col]]).toarray()[:, 1:])\n",
    "    onehot_vector.columns = [str(i)+'_'+col[:3] for i in onehot_vector.columns]\n",
    "    predictors = pd.concat([predictors, onehot_vector], axis=1)\n",
    "predictors = sm.add_constant(predictors) # add constant\n",
    "sm.OLS(y, predictors, hasconst=True).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Holding other factors in the model constant, Blacks earn about 10.4% lower hourly wages (0.104 = 1-exp(-0.1094))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrimination in predicting high/low earners"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train-test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = df['50k']\n",
    "X = df[['blk', 'AGE', 'AGE2', 'AGE3', 'SEX', 'EDUC', 'SCHLCOLL', 'MARST', 'pt']]\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, random_state=999)\n",
    "X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.2, random_state=999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(38642, 9)\n",
      "(9661, 9)\n",
      "(5368, 9)\n"
     ]
    }
   ],
   "source": [
    "print(X_train.shape)\n",
    "print(X_val.shape)\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percent high-income Blacks in test set: 26.91%\n",
      "Percent high-income Whites in test set: 42.34%\n",
      "Difference: -15.43%\n"
     ]
    }
   ],
   "source": [
    "foo = pd.concat([X_test, y_test], axis=1)\n",
    "print(f\"Percent high-income Blacks in test set: {round(foo[foo.blk==1]['50k'].mean()*100, 2)}%\")\n",
    "print(f\"Percent high-income Whites in test set: {round(foo[foo.blk==0]['50k'].mean()*100, 2)}%\")\n",
    "print(f\"Difference: {round((foo[foo.blk==1]['50k'].mean() - foo[foo.blk==0]['50k'].mean())*100, 2)}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive model: Excluding race"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "model1 = lgb.LGBMClassifier(objective='binary',\n",
    "                           random_state=999,\n",
    "                           metric='logloss')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LGBMClassifier(metric='logloss', objective='binary', random_state=999)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train_mod = X_train.loc[:, X_train.columns != 'blk']\n",
    "X_val_mod = X_val.loc[:, X_val.columns != 'blk']\n",
    "X_test_mod = X_test.loc[:, X_test.columns != 'blk']\n",
    "model1.fit(X_train_mod, y_train, eval_set=[(X_val_mod, y_val)], eval_metric='logloss', early_stopping_rounds=10, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "naive_pred = model1.predict(X_test_mod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accuracy': 0.7493,\n",
       " 'sensitivity/recall': 0.6884,\n",
       " 'specificity': 0.7897,\n",
       " 'precision': 0.6852,\n",
       " 'f1': 0.6868,\n",
       " 'roc_auc': 0.7391}"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_metrics(y_test, naive_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "naive = combine_prediction(X_test, y_test, naive_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 5: Linear probability model regressing naive model predicted values on race dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>pred</td>       <th>  R-squared:         </th> <td>   0.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   27.24</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>1.87e-07</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:46</td>     <th>  Log-Likelihood:    </th> <td> -3776.6</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  5368</td>      <th>  AIC:               </th> <td>   7557.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  5366</td>      <th>  BIC:               </th> <td>   7570.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.4157</td> <td>    0.007</td> <td>   56.772</td> <td> 0.000</td> <td>    0.401</td> <td>    0.430</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.0927</td> <td>    0.018</td> <td>   -5.219</td> <td> 0.000</td> <td>   -0.128</td> <td>   -0.058</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>21127.127</td> <th>  Durbin-Watson:     </th> <td>   2.029</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>   <th>  Jarque-Bera (JB):  </th> <td> 884.818</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td> 0.399</td>   <th>  Prob(JB):          </th> <td>7.31e-193</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 1.178</td>   <th>  Cond. No.          </th> <td>    2.83</td> \n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   pred   R-squared:                       0.005\n",
       "Model:                            OLS   Adj. R-squared:                  0.005\n",
       "Method:                 Least Squares   F-statistic:                     27.24\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):           1.87e-07\n",
       "Time:                        10:04:46   Log-Likelihood:                -3776.6\n",
       "No. Observations:                5368   AIC:                             7557.\n",
       "Df Residuals:                    5366   BIC:                             7570.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.4157      0.007     56.772      0.000       0.401       0.430\n",
       "blk           -0.0927      0.018     -5.219      0.000      -0.128      -0.058\n",
       "==============================================================================\n",
       "Omnibus:                    21127.127   Durbin-Watson:                   2.029\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              884.818\n",
       "Skew:                           0.399   Prob(JB):                    7.31e-193\n",
       "Kurtosis:                       1.178   Cond. No.                         2.83\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "naive_reg = sm.add_constant(naive) # add constant\n",
    "naive_reg['blk'] = naive_reg['blk'].astype(int)\n",
    "sm.OLS(naive_reg['pred'], naive_reg[['const', 'blk']]).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to an observed differential of -0.1543 in the test set, the naive model predicts a smaller differential at -0.093"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discriminatory model: Including race"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "model2 = lgb.LGBMClassifier(objective='binary',\n",
    "                           random_state=999,\n",
    "                           metric='logloss')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LGBMClassifier(metric='logloss', objective='binary', random_state=999)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model2.fit(X_train, y_train, eval_set=[(X_val, y_val)], eval_metric='logloss', early_stopping_rounds=10, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "discrim_pred = model2.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accuracy': 0.7547,\n",
       " 'sensitivity/recall': 0.7206,\n",
       " 'specificity': 0.7773,\n",
       " 'precision': 0.6827,\n",
       " 'f1': 0.7011,\n",
       " 'roc_auc': 0.749}"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_metrics(y_test, discrim_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "discrim = combine_prediction(X_test, y_test, discrim_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 6: Linear probability model regressing discriminatory model predicted values on race dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>pred</td>       <th>  R-squared:         </th> <td>   0.027</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.027</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   184.8</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>2.07e-41</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:47</td>     <th>  Log-Likelihood:    </th> <td> -3756.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  5368</td>      <th>  AIC:               </th> <td>   7516.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  5366</td>      <th>  BIC:               </th> <td>   7530.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.4563</td> <td>    0.007</td> <td>   61.660</td> <td> 0.000</td> <td>    0.442</td> <td>    0.471</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.2231</td> <td>    0.016</td> <td>  -13.595</td> <td> 0.000</td> <td>   -0.255</td> <td>   -0.191</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>21106.787</td> <th>  Durbin-Watson:     </th> <td>   2.000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>   <th>  Jarque-Bera (JB):  </th> <td> 814.198</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td> 0.287</td>   <th>  Prob(JB):          </th> <td>1.58e-177</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 1.180</td>   <th>  Cond. No.          </th> <td>    2.83</td> \n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   pred   R-squared:                       0.027\n",
       "Model:                            OLS   Adj. R-squared:                  0.027\n",
       "Method:                 Least Squares   F-statistic:                     184.8\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):           2.07e-41\n",
       "Time:                        10:04:47   Log-Likelihood:                -3756.2\n",
       "No. Observations:                5368   AIC:                             7516.\n",
       "Df Residuals:                    5366   BIC:                             7530.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.4563      0.007     61.660      0.000       0.442       0.471\n",
       "blk           -0.2231      0.016    -13.595      0.000      -0.255      -0.191\n",
       "==============================================================================\n",
       "Omnibus:                    21106.787   Durbin-Watson:                   2.000\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              814.198\n",
       "Skew:                           0.287   Prob(JB):                    1.58e-177\n",
       "Kurtosis:                       1.180   Cond. No.                         2.83\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discrim_reg = sm.add_constant(discrim) # add constant\n",
    "discrim_reg['blk'] = discrim_reg['blk'].astype(int)\n",
    "sm.OLS(discrim_reg['pred'], discrim_reg[['const', 'blk']]).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to an observed differential of -0.1543 in the test set, the discriminatory model predicts a *larger* differential at -0.2231"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `DiscriminationMitigator`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "config = {'protected_class_features': ['blk']}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "mitigated = DiscriminationMitigator(df=X_test, model=model2, config=config).predictions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        unadj_pred  unif_wts\n",
      "94428            1       1.0\n",
      "145967           1       0.5\n",
      "18557            0       0.0\n",
      "175691           1       0.5\n",
      "132926           1       1.0\n"
     ]
    }
   ],
   "source": [
    "print(mitigated.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluating performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "mitigated['thresh_0.5'] = np.where(mitigated['unif_wts'] >= 0.5, 1, 0) # naive threshold of 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "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>unadj_pred</th>\n",
       "      <th>unif_wts</th>\n",
       "      <th>thresh_0.5</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>5368.000000</td>\n",
       "      <td>5368.000000</td>\n",
       "      <td>5368.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.421572</td>\n",
       "      <td>0.377981</td>\n",
       "      <td>0.440574</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.493857</td>\n",
       "      <td>0.451500</td>\n",
       "      <td>0.496502</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\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>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        unadj_pred     unif_wts   thresh_0.5\n",
       "count  5368.000000  5368.000000  5368.000000\n",
       "mean      0.421572     0.377981     0.440574\n",
       "std       0.493857     0.451500     0.496502\n",
       "min       0.000000     0.000000     0.000000\n",
       "25%       0.000000     0.000000     0.000000\n",
       "50%       0.000000     0.000000     0.000000\n",
       "75%       1.000000     1.000000     1.000000\n",
       "max       1.000000     1.000000     1.000000"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mitigated.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accuracy': 0.7513,\n",
       " 'sensitivity/recall': 0.7402,\n",
       " 'specificity': 0.7587,\n",
       " 'precision': 0.671,\n",
       " 'f1': 0.7039,\n",
       " 'roc_auc': 0.7494}"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binary_metrics(y_test, mitigated['thresh_0.5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 7: Linear probability model regressing discriminatory model predicted values on race dummy\n",
    "#### Adjusted using `DiscriminationMitigator`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>thresh_0.5</td>    <th>  R-squared:         </th> <td>   0.006</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   31.91</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>1.70e-08</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:48</td>     <th>  Log-Likelihood:    </th> <td> -3842.8</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  5368</td>      <th>  AIC:               </th> <td>   7690.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  5366</td>      <th>  BIC:               </th> <td>   7703.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.4565</td> <td>    0.007</td> <td>   61.687</td> <td> 0.000</td> <td>    0.442</td> <td>    0.471</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.1025</td> <td>    0.018</td> <td>   -5.649</td> <td> 0.000</td> <td>   -0.138</td> <td>   -0.067</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>19574.254</td> <th>  Durbin-Watson:     </th> <td>   2.006</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th>  <td> 0.000</td>   <th>  Jarque-Bera (JB):  </th> <td> 876.238</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>           <td> 0.236</td>   <th>  Prob(JB):          </th> <td>5.34e-191</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>       <td> 1.078</td>   <th>  Cond. No.          </th> <td>    2.83</td> \n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:             thresh_0.5   R-squared:                       0.006\n",
       "Model:                            OLS   Adj. R-squared:                  0.005\n",
       "Method:                 Least Squares   F-statistic:                     31.91\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):           1.70e-08\n",
       "Time:                        10:04:48   Log-Likelihood:                -3842.8\n",
       "No. Observations:                5368   AIC:                             7690.\n",
       "Df Residuals:                    5366   BIC:                             7703.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.4565      0.007     61.687      0.000       0.442       0.471\n",
       "blk           -0.1025      0.018     -5.649      0.000      -0.138      -0.067\n",
       "==============================================================================\n",
       "Omnibus:                    19574.254   Durbin-Watson:                   2.006\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              876.238\n",
       "Skew:                           0.236   Prob(JB):                    5.34e-191\n",
       "Kurtosis:                       1.078   Cond. No.                         2.83\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mitigated_reg = sm.add_constant(mitigated) # add constant\n",
    "mitigated_reg['blk'] = X_test['blk']\n",
    "mitigated_reg['blk'] = mitigated_reg['blk'].astype(int)\n",
    "sm.OLS(mitigated_reg['thresh_0.5'], mitigated_reg[['const', 'blk']]).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrimination in predicting high-wage workers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train-test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = df['lnwage']\n",
    "X = df[['blk', 'AGE', 'AGE2', 'AGE3', 'SEX', 'EDUC', 'SCHLCOLL', 'MARST', 'pt']]\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, random_state=999)\n",
    "X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.2, random_state=999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(38642, 9)\n",
      "(9661, 9)\n",
      "(5368, 9)\n"
     ]
    }
   ],
   "source": [
    "print(X_train.shape)\n",
    "print(X_val.shape)\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean hourly wage of Blacks in test set: $21.26\n",
      "Mean hourly wage of Whites in test set: $26.56\n",
      "Difference: $-5.3\n",
      "Log difference: -0.22\n"
     ]
    }
   ],
   "source": [
    "foo = pd.concat([X_test, y_test], axis=1)\n",
    "print(f\"Mean hourly wage of Blacks in test set: ${round(np.exp(foo[foo.blk==1]['lnwage']).mean(), 2)}\")\n",
    "print(f\"Mean hourly wage of Whites in test set: ${round(np.exp(foo[foo.blk==0]['lnwage']).mean(), 2)}\")\n",
    "print(f\"Difference: ${round(np.exp(foo[foo.blk==1]['lnwage']).mean() - np.exp(foo[foo.blk==0]['lnwage']).mean(), 2)}\")\n",
    "print(f\"Log difference: {round(foo[foo.blk==1]['lnwage'].mean() - foo[foo.blk==0]['lnwage'].mean(), 2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive model: Excluding race"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "model3 = lgb.LGBMRegressor(random_state=999,\n",
    "                           metric='mean_squared_error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LGBMRegressor(metric='mean_squared_error', random_state=999)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train_mod = X_train.loc[:, X_train.columns != 'blk']\n",
    "X_val_mod = X_val.loc[:, X_val.columns != 'blk']\n",
    "X_test_mod = X_test.loc[:, X_test.columns != 'blk']\n",
    "model3.fit(X_train_mod, y_train, eval_set=[(X_val_mod, y_val)], eval_metric='mean_squared_error', early_stopping_rounds=10, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "naive_pred_cont = model3.predict(X_test_mod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mae 0.4183\n",
      "mse 0.306\n",
      "rmse 0.5532\n",
      "r2 0.2979\n"
     ]
    }
   ],
   "source": [
    "for key, val in continuous_metrics(y_test, naive_pred_cont).items():\n",
    "    print(key, val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "naive_cont = combine_prediction(X_test, y_test, naive_pred_cont, outcome='lnwage')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 8: OLS regression of naive model predicted values on race dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>pred</td>       <th>  R-squared:         </th> <td>   0.012</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.011</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   64.78</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>1.02e-15</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:48</td>     <th>  Log-Likelihood:    </th> <td> -2089.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  5368</td>      <th>  AIC:               </th> <td>   4183.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  5366</td>      <th>  BIC:               </th> <td>   4197.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    3.0644</td> <td>    0.005</td> <td>  575.472</td> <td> 0.000</td> <td>    3.054</td> <td>    3.075</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.1066</td> <td>    0.013</td> <td>   -8.048</td> <td> 0.000</td> <td>   -0.133</td> <td>   -0.081</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>155.295</td> <th>  Durbin-Watson:     </th> <td>   2.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>  90.652</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.162</td>  <th>  Prob(JB):          </th> <td>2.07e-20</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.452</td>  <th>  Cond. No.          </th> <td>    2.83</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   pred   R-squared:                       0.012\n",
       "Model:                            OLS   Adj. R-squared:                  0.011\n",
       "Method:                 Least Squares   F-statistic:                     64.78\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):           1.02e-15\n",
       "Time:                        10:04:48   Log-Likelihood:                -2089.7\n",
       "No. Observations:                5368   AIC:                             4183.\n",
       "Df Residuals:                    5366   BIC:                             4197.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          3.0644      0.005    575.472      0.000       3.054       3.075\n",
       "blk           -0.1066      0.013     -8.048      0.000      -0.133      -0.081\n",
       "==============================================================================\n",
       "Omnibus:                      155.295   Durbin-Watson:                   2.005\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               90.652\n",
       "Skew:                          -0.162   Prob(JB):                     2.07e-20\n",
       "Kurtosis:                       2.452   Cond. No.                         2.83\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "naive_reg_cont = sm.add_constant(naive_cont) # add constant\n",
    "naive_reg_cont['blk'] = naive_reg_cont['blk'].astype(int)\n",
    "sm.OLS(naive_reg_cont['pred'], naive_reg_cont[['const', 'blk']]).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to a observed difference of -0.22 log percentage points in between Blacks and Whites in the test set, a naive model predicts -0.1066 log percentage point difference. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discriminatory model: Including race"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "model4 = lgb.LGBMRegressor(random_state=999,\n",
    "                           metric='mean_squared_error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LGBMRegressor(metric='mean_squared_error', random_state=999)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model4.fit(X_train, y_train, eval_set=[(X_val, y_val)], eval_metric='MSE', early_stopping_rounds=10, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "discrim_pred_cont = model4.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mae 0.4166\n",
      "mse 0.3032\n",
      "rmse 0.5506\n",
      "r2 0.3043\n"
     ]
    }
   ],
   "source": [
    "for key, val in continuous_metrics(y_test, discrim_pred_cont).items():\n",
    "    print(key, val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "discrim_cont = combine_prediction(X_test, y_test, discrim_pred_cont, outcome='lnwage')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 9: OLS regression of discriminatory model predicted values on race dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>pred</td>       <th>  R-squared:         </th> <td>   0.040</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.040</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   264.8</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>3.70e-58</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:49</td>     <th>  Log-Likelihood:    </th> <td> -1992.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  5368</td>      <th>  AIC:               </th> <td>   3988.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  5366</td>      <th>  BIC:               </th> <td>   4002.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    3.0775</td> <td>    0.005</td> <td>  580.198</td> <td> 0.000</td> <td>    3.067</td> <td>    3.088</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.1971</td> <td>    0.012</td> <td>  -16.273</td> <td> 0.000</td> <td>   -0.221</td> <td>   -0.173</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>149.378</td> <th>  Durbin-Watson:     </th> <td>   2.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>  95.898</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.198</td>  <th>  Prob(JB):          </th> <td>1.50e-21</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.478</td>  <th>  Cond. No.          </th> <td>    2.83</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                   pred   R-squared:                       0.040\n",
       "Model:                            OLS   Adj. R-squared:                  0.040\n",
       "Method:                 Least Squares   F-statistic:                     264.8\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):           3.70e-58\n",
       "Time:                        10:04:49   Log-Likelihood:                -1992.2\n",
       "No. Observations:                5368   AIC:                             3988.\n",
       "Df Residuals:                    5366   BIC:                             4002.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          3.0775      0.005    580.198      0.000       3.067       3.088\n",
       "blk           -0.1971      0.012    -16.273      0.000      -0.221      -0.173\n",
       "==============================================================================\n",
       "Omnibus:                      149.378   Durbin-Watson:                   2.005\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               95.898\n",
       "Skew:                          -0.198   Prob(JB):                     1.50e-21\n",
       "Kurtosis:                       2.478   Cond. No.                         2.83\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "discrim_reg_cont = sm.add_constant(discrim_cont) # add constant\n",
    "discrim_reg_cont['blk'] = discrim_reg_cont['blk'].astype(int)\n",
    "sm.OLS(discrim_reg_cont['pred'], discrim_reg_cont[['const', 'blk']]).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to a observed difference of -0.22 log percentage points in between Blacks and Whites in the test set, a discriminatory model predicts a difference of -0.1971, which is considerably larger than the naive model prediction at -0.1066. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `DiscriminationMitigator`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "mitigated_cont = DiscriminationMitigator(df=X_test, model=model4, config=config).predictions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        unadj_pred  unif_wts\n",
      "94428     3.301243  3.289145\n",
      "145967    3.146659  3.047996\n",
      "18557     2.801680  2.762560\n",
      "175691    3.143289  3.040616\n",
      "132926    3.291120  3.220549 \n",
      "\n",
      "        unadj_pred     unif_wts\n",
      "count  5368.000000  5368.000000\n",
      "mean      3.046788     3.009338\n",
      "std       0.357940     0.338295\n",
      "min       2.047689     2.093242\n",
      "25%       2.808305     2.794492\n",
      "50%       3.046140     3.018227\n",
      "75%       3.301243     3.270785\n",
      "max       3.823983     3.704577\n"
     ]
    }
   ],
   "source": [
    "print(mitigated_cont.head(), '\\n')\n",
    "print(mitigated_cont.describe())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mae 0.4213\n",
      "mse 0.308\n",
      "rmse 0.555\n",
      "r2 0.2933\n"
     ]
    }
   ],
   "source": [
    "for key, val in continuous_metrics(y_test, mitigated_cont['unif_wts']).items():\n",
    "    print(key, val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table 10: OLS regression of discriminatory model predicted values on race dummy\n",
    "#### Adjusted using `DiscriminationMitigator`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>        <td>unif_wts</td>     <th>  R-squared:         </th> <td>   0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   58.00</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 14 Apr 2021</td> <th>  Prob (F-statistic):</th> <td>3.08e-14</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>10:04:50</td>     <th>  Log-Likelihood:    </th> <td> -1770.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  5368</td>      <th>  AIC:               </th> <td>   3545.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  5366</td>      <th>  BIC:               </th> <td>   3558.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC3</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    3.0241</td> <td>    0.005</td> <td>  602.816</td> <td> 0.000</td> <td>    3.014</td> <td>    3.034</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>blk</th>   <td>   -0.0951</td> <td>    0.012</td> <td>   -7.616</td> <td> 0.000</td> <td>   -0.120</td> <td>   -0.071</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>151.197</td> <th>  Durbin-Watson:     </th> <td>   2.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td> 105.460</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.234</td>  <th>  Prob(JB):          </th> <td>1.26e-23</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.498</td>  <th>  Cond. No.          </th> <td>    2.83</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors are heteroscedasticity robust (HC3)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:               unif_wts   R-squared:                       0.010\n",
       "Model:                            OLS   Adj. R-squared:                  0.010\n",
       "Method:                 Least Squares   F-statistic:                     58.00\n",
       "Date:                Wed, 14 Apr 2021   Prob (F-statistic):           3.08e-14\n",
       "Time:                        10:04:50   Log-Likelihood:                -1770.3\n",
       "No. Observations:                5368   AIC:                             3545.\n",
       "Df Residuals:                    5366   BIC:                             3558.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:                  HC3                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          3.0241      0.005    602.816      0.000       3.014       3.034\n",
       "blk           -0.0951      0.012     -7.616      0.000      -0.120      -0.071\n",
       "==============================================================================\n",
       "Omnibus:                      151.197   Durbin-Watson:                   2.005\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              105.460\n",
       "Skew:                          -0.234   Prob(JB):                     1.26e-23\n",
       "Kurtosis:                       2.498   Cond. No.                         2.83\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC3)\n",
       "\"\"\""
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mitigated_reg_cont = sm.add_constant(mitigated_cont) # add constant\n",
    "mitigated_reg_cont['blk'] = X_test['blk']\n",
    "mitigated_reg_cont['blk'] = mitigated_reg_cont['blk'].astype(int)\n",
    "sm.OLS(mitigated_reg_cont['unif_wts'], mitigated_reg_cont[['const', 'blk']]).fit(cov_type='HC3', use_t=True).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to a predicted of differential at -0.1066 according to the naive model, the mitigated model yields a smaller predicted differential at -0.0951."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}