intermediate_notebooks/E2E/census/census_education2income_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Census Notebook\n",
    "**Authorship**<br />\n",
    "Original Author: Taurean Dyer<br />\n",
    "Last Edit: Taurean Dyer, 9/26/2019<br />\n",
    "\n",
    "**Test System Specs**<br />\n",
    "Test System Hardware: GV100<br />\n",
    "Test System Software: Ubuntu 18.04<br />\n",
    "RAPIDS Version: 0.10.0a - Docker Install<br />\n",
    "Driver: 410.79<br />\n",
    "CUDA: 10.0<br />\n",
    "\n",
    "\n",
    "**Known Working Systems**<br />\n",
    "RAPIDS Versions:0.8, 0.9, 0.10\n",
    "\n",
    "# Intro\n",
    "Held every 10 years, the US census gives a detailed snapshot in time about the makeup of the country.  The last census in 2010 surveyed nearly 309 million people.  IPUMS.org provides researchers an open source data set with 1% to 10% of the census data set.  In this notebook, we want to see how education affects total income earned in the US based on data from each census from the 1970 to 2010 and see if we can predict some results if the census was held today, according to the national average.  We will go through the ETL, training the model, and then testing the prediction.  We'll make every effort to get as balanced of a dataset as we can.  We'll also pull some extra variables to allow for further self-exploration of gender based education and income breakdowns.  On a single Titan RTX, you can run the whole notebook workflow on the 4GB dataset of 14 million rows by 44 columns in less than 3 minutes.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Let's begin!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import cuml\n",
    "import cudf\n",
    "import dask_cudf\n",
    "import sys\n",
    "import os\n",
    "from pprint import pprint\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get your data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ipums dataset is in our S3 bucket and zipped.  \n",
    "1. We'll need to create a folder for our data in the `/data` folder\n",
    "1. Download the zipped data into that folder from S3\n",
    "1. Load the zipped data quickly into cudf using it's read_csv() parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "\n",
    "data_dir = '../../../data/census/'\n",
    "if not os.path.exists(data_dir):\n",
    "    print('creating census data directory')\n",
    "    os.system('mkdir ../../../data/census')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download the IPUMS dataset\n",
    "base_url = 'https://rapidsai-data.s3.us-east-2.amazonaws.com/datasets/'\n",
    "fn = 'ipums_education2income_1970-2010.csv.gz'\n",
    "if not os.path.isfile(data_dir+fn):\n",
    "    print(f'Downloading {base_url+fn} to {data_dir+fn}')\n",
    "    urllib.request.urlretrieve(base_url+fn, data_dir+fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data(cached = data_dir+fn):\n",
    "    if os.path.exists(cached):\n",
    "        print('use ipums data')\n",
    "        X = cudf.read_csv(cached, compression='infer')\n",
    "    else:\n",
    "        print(\"No data found!  Please check your that your data directory is ../../../data/census/ and that you downloaded the data.  If you did, please delete the `../../../data/census/` directory and try the above 2 cells again\")\n",
    "        X = null\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = load_data(data_dir+fn)\n",
    "print('data',df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.head(5).to_pandas())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_counts = df.YEAR.value_counts()\n",
    "print(original_counts) ### Remember these numbers!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ETL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning Income data\n",
    "First, let's focus on cleaning out the bad values for Total Income `INCTOT`. First, let's see if there are an `N/A` values, as when we did `head()`, we saw some in other columns, like CBSERIAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['INCTOT_NA'] = df['INCTOT'].isna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.INCTOT_NA.value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, great, there are no `N/A`s...or are there?  Let's drop `INCTOT_NA` and see what our value counts look like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=df.drop('INCTOT_NA')\n",
    "print(df.INCTOT.value_counts().to_pandas())  ### Wow, look how many people in America make $10,000,000!  Wait a minutes... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not that many people make $10M a year. Checking https://usa.ipums.org/usa-action/variables/INCTOT#codes_section, `9999999`is INCTOT's code for `N/A`.  That was why when we ran `isna`, RAPIDS won't find any.  Let's first create a new dataframe that is only NA values, then let's pull those encoded `N/A`s out of our working dataframe!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('data',df.shape)\n",
    "tdf = df.query('INCTOT == 9999999')\n",
    "df = df.query('INCTOT != 9999999')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('working data',df.shape)\n",
    "print('junk count data',tdf.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're down by nearly 1/4 of our original dataset size.  For the curious, now we should be able to get accurate Total Income data, by year, not taking into account inflation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.groupby('YEAR')['INCTOT'].mean()) # without that cleanup, the average would have bene in the millions...."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Normalize Income for inflation\n",
    "Now that we have reduced our dataframe to a baseline clean data to answer our question, we should normalize the amounts for inflation.  `CPI99`is the value that IPUMS uses to contian the inflation factor.  All we have to do is multipy by year.  Let's see how that changes the Total Income values from just above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.groupby('YEAR')['CPI99'].mean()) ## it just returns the CPI99\n",
    "df['INCTOT'] = df['INCTOT'] * df['CPI99']\n",
    "print(df.groupby('YEAR')['INCTOT'].mean()) ## let's see what we got!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning Education Data\n",
    "Okay, great!  Now we have income cleaned up, it should also have cleaned much of our next sets of values of interes, namely Education and Education Detailed.  However, there are still some `N/A`s in key variables to worry about, which can cause problmes later.  Let's create a list of them..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suspect = ['CBSERIAL','EDUC', 'EDUCD', 'EDUC_HEAD', 'EDUC_POP', 'EDUC_MOM','EDUCD_MOM2','EDUCD_POP2', 'INCTOT_MOM','INCTOT_POP','INCTOT_MOM2','INCTOT_POP2', 'INCTOT_HEAD']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(0, len(suspect)):\n",
    "    df[suspect[i]] = df[suspect[i]].fillna(-1)\n",
    "    print(suspect[i], df[suspect[i]].value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get drop any rows of any `-1`s in Education and Education Detailed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "totincome = ['EDUC','EDUCD']\n",
    "for i in range(0, len(totincome)):\n",
    "    query = totincome[i] + ' != -1'\n",
    "    df = df.query(query)\n",
    "    print(totincome[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.shape)\n",
    "df.head().to_pandas().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, the good news is that we lost no further rows, start to normalize the data so when we do our OLS, one year doesn't unfairly dominate the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalize the Data\n",
    "The in the last step, need to keep our data at about the same ratio as we when started (1% of the population), with the exception of 1980, which was a 5% and needs to be reduced.  This is why we kept the temp dataframe `tdf` - to get the counts per year.   we will find out just how many have to realize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Working data: \\n', df.YEAR.value_counts())\n",
    "print('junk count data: \\n', tdf.YEAR.value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now, so that we can do MSE, let's make all the dtypes the same. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "keep_cols = ['YEAR', 'DATANUM', 'SERIAL', 'CBSERIAL', 'HHWT', 'GQ', 'PERNUM', 'SEX', 'AGE', 'INCTOT', 'EDUC', 'EDUCD', 'EDUC_HEAD', 'EDUC_POP', 'EDUC_MOM','EDUCD_MOM2','EDUCD_POP2', 'INCTOT_MOM','INCTOT_POP','INCTOT_MOM2','INCTOT_POP2', 'INCTOT_HEAD', 'SEX_HEAD']\n",
    "df = df.loc[:, keep_cols]\n",
    "#df = df.drop(col for col in df.columns if col not in keep_cols)\n",
    "for i in range(0, len(keep_cols)):\n",
    "    df[keep_cols[i]] = df[keep_cols[i]].fillna(-1)\n",
    "    print(keep_cols[i], df[keep_cols[i]].value_counts())\n",
    "    df[keep_cols[i]]= df[keep_cols[i]].astype('float64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## I WANTED TO REDUCE THE 1980 SAMPLE HERE, BUT .SAMPLE() IS NEEDED AND NOT WORKING, UNLESS THERE IS A WORK AROUND..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the important data now clean and normalized, let's start doing the regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ridge Regression\n",
    "We have 44 variables.  The other variables may provide important predictive information.  The Ridge Regression technique with cross validation to identify the best hyperparamters may be the best way to get the most accurate model.  We'll have to \n",
    "\n",
    "* define our performance metrics\n",
    "* split our data into train and test sets\n",
    "* train and test our model\n",
    "\n",
    "Let's begin and see what we get!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As our performance metrics we'll use a basic mean squared error and coefficient of determination implementation\n",
    "def mse(y_test, y_pred):\n",
    "    return ((y_test.reset_index(drop=True) - y_pred.reset_index(drop=True)) ** 2).mean()\n",
    "\n",
    "def cod(y_test, y_pred):\n",
    "    y_bar = y_test.mean()\n",
    "    total = ((y_test - y_bar) ** 2).sum()\n",
    "    residuals = ((y_test.reset_index(drop=True) - y_pred.reset_index(drop=True)) ** 2).sum()\n",
    "    return 1 - (residuals / total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cuml.preprocessing.model_selection import train_test_split\n",
    "trainsize = .9\n",
    "yCol = \"EDUC\"\n",
    "from cuml.preprocessing.model_selection import train_test_split\n",
    "from cuml.linear_model.ridge import Ridge\n",
    "\n",
    "def train_and_score(data, clf, train_frac=0.8, n_runs=20):\n",
    "    mse_scores, cod_scores = [], []\n",
    "    for _ in range(n_runs):\n",
    "        X_train, X_test, y_train, y_test = cuml.preprocessing.model_selection.train_test_split(df, yCol, train_size=.9)\n",
    "        y_pred = clf.fit(X_train, y_train).predict(X_test)\n",
    "        mse_scores.append(mse(y_test, y_pred))\n",
    "        cod_scores.append(cod(y_test, y_pred))\n",
    "    return mse_scores, cod_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## Results\n",
    " **Moment of truth!  Let's see how our regression training does!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "n_runs = 20\n",
    "clf = Ridge()\n",
    "mse_scores, cod_scores = train_and_score(df, clf, n_runs=n_runs)\n",
    "print(f\"median MSE ({n_runs} runs): {np.median(mse_scores)}\")\n",
    "print(f\"median COD ({n_runs} runs): {np.median(cod_scores)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Fun fact:** if you made INCTOT the y axis, your prediction results would not be so pretty!  It just shows that your education level can be an indicator for your income, but your income is NOT a great predictor for your education level.  You have better odds flipping a coin!\n",
    "\n",
    "* median MSE (50 runs): 518189521.07548225\n",
    "* median COD (50 runs): 0.425769113846303"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next Steps/Self Study\n",
    "* You can pickle the model and use it in another workflow\n",
    "* You can redo the workflow with based on head of household using `EDUC`, `SEX`, and `INCTOT` for X in `X`_HEAD\n",
    "* You can see the growing role of education with women in their changing role in the workforce and income with \"EDUC_MOM\" and \"EDUC_POP"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}