5.Inference_Imputation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> <img style=\"float: center;\" src=\"images/CI_horizontal.png\" width=\"450\">\n",
    "<center>\n",
    "    <span style=\"font-size: 1.5em;\">\n",
    "        <a href='https://www.coleridgeinitiative.org'>Website</a>\n",
    "    </span> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Imputation, Outcome Measurement, and Inference</center>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6399344\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6399344.svg\" alt=\"DOI\"></a>\n",
    "\n",
    "<center> Benjamin Feder, Nathan Barrett </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "What should you do when you encounter missing values in your data? Unfortunately, there is usually no *right* answer. However, different decisions on how to address missing values can have different implications on the inferences you draw from any analysis. One option is to omit observations that have missing values, but you may also choose to retain all observations and instead impute these missing values using various techniques. Imputation provides your best guess for each missing point's true value. Here, you will learn how to implement common imputation methods for missing data and how those methods affect outcome measures and inferences.\n",
    "\n",
    "### Learning Objectives\n",
    "\n",
    "* Explore options for imputing missing values\n",
    "\n",
    "* Visualize estimate changes following imputation\n",
    "\n",
    "\n",
    "In this notebook, you will focus on 2012-13 New Jersey graduates' earnings during their first year after graduation. The notebook provides demonstrations using the first quarter after graduation and you will replicate using the fourth quarter. Recall that in the [Data Exploration](2.Dataset_Exploration.ipynb) notebook, you initially examined the earnings distribution for all members of this cohort who had positive earnings in this time period in New Jersey. To evaluate the earnings outcomes of all 2012-13 New Jersey graduates, you need to decide what to do when you cannot find their earnings in the New Jersey Unemployment Insurance (UI) wage records. A person may not appear in New Jersey's UI wage records for several reasons:\n",
    "- The person is unemployed. \n",
    "- The person is out of labor force, e.g., schooling, childcare, etc...\n",
    "- The person was employed outside of New Jersey.\n",
    "- The person's job is not covered in UI wage records, e.g.,self-employed, independent contractors, federal government works, etc. <a href='https://www.nap.edu/read/10206/chapter/11#294'>(Hotz and Scholz, 2002)</a>\n",
    "\n",
    "You will explore the resulting earnings outcomes after applying different earnings imputation methods. The methods covered in this notebook include:\n",
    "- Dropping all \"missing\" values\n",
    "- Filling in zero for people who do not have records in New Jersey UI wage records data \n",
    "- Substituting missing values with the average earnings of people who are in the same degree fields and have the same gender\n",
    "- Regression imputation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R Setup and Database Connection\n",
    "\n",
    "Before you begin, you need to run the code cells below to import the libraries and connect to the proper server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Database interaction imports\n",
    "library(odbc, warn.conflicts=F, quietly=T)\n",
    "\n",
    "# For data manipulation/visualization\n",
    "library(tidyverse, warn.conflicts=F, quietly=T)\n",
    "\n",
    "# For faster date conversions\n",
    "library(lubridate, warn.conflicts=F, quietly=T)\n",
    "\n",
    "# Use percent() function\n",
    "library(scales, warn.conflicts=F, quietly=T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Establish a Connection to the Server**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we are ready to connect to the server. We will create the connection using the `DBI`  and `ODBC` libraries. \n",
    "\n",
    "> **Loading R libraries** and **establishing connection** should always be the first step in your Jupyter Notebooks. Make sure you copy these code chunks when you create a new notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connect to the server\n",
    "con <- DBI::dbConnect(odbc::odbc(),\n",
    "                     Driver = \"SQL Server\",\n",
    "                     Server = \"msssql01.c7bdq4o2yhxo.us-gov-west-1.rds.amazonaws.com\",\n",
    "                     Trusted_Connection = \"True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brief Manipulation: Isolating Earnings during first quarter after graduation\n",
    "\n",
    "Before we start performing imputation, we need to do some quick data manipulation to isolate earnings from the first quarter after each individual's graduation. To do so, using the same approach as we did in the last [section](2.Data_Exploration.ipynb) of the Data Exploration notebook, we will create a new column, `quarter_number`, by taking the difference between `job_date` and `grad_date` then dividing by 13 and rounding to the nearest whole number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in earnings table\n",
    "qry <- \"\n",
    "SELECT * \n",
    "FROM tr_nj_2021.dbo.nb_cohort_wages_link;\n",
    "\"\n",
    "\n",
    "df_wages <- dbGetQuery(con, qry)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add in quarter number\n",
    "df_wages <- df_wages %>%\n",
    "    mutate(\n",
    "        quarter_number = round(as.double(difftime(as.Date(job_date), as.Date(grad_date), units = \"weeks\")/13), 0)\n",
    "    )\n",
    "\n",
    "# see evidence\n",
    "df_wages %>%\n",
    "    select(grad_date, job_date, quarter_number) %>%\n",
    "    head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can simply create a data frame with just first quarter post-graduation wages using `filter()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter quarter 1 after graduation\n",
    "q1_wages <- df_wages %>%\n",
    "    filter(quarter_number == 1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we will want to estimate the total earnings for each `hashed_ssn` in this quarter, not necessarily their wages per employer, let's aggregate `q1_wages` to find the total earnings for each member of this cohort in the entire quarter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get total earnings in quarter\n",
    "q1_wages <- q1_wages %>%\n",
    "    group_by(hashed_ssn) %>%\n",
    "    summarize(tot_wages = sum(wage)) %>%\n",
    "    ungroup()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of graduates with positive earnings\n",
    "q1_num <- q1_wages %>%\n",
    "    summarize(\n",
    "        n=n_distinct(hashed_ssn)\n",
    "    )\n",
    "\n",
    "cat('The total graduates with positive earnings during their first quarter after graduation:', q1_num$n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in cohort\n",
    "qry <- \"\n",
    "SELECT * \n",
    "FROM tr_nj_2021.dbo.nb_cohort_dated;\n",
    "\"\n",
    "\n",
    "df <- dbGetQuery(con, qry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in major information\n",
    "qry <- \"\n",
    "SELECT code_2010 as major, title_2010 as major_title\n",
    "FROM ds_nj_oshe.dbo.supplements_cipcode;\n",
    "\"\n",
    "\n",
    "df_cip <- dbGetQuery(con,qry)\n",
    "\n",
    "df <- df %>% \n",
    "    left_join(df_cip, by.y = \"major\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat('That is', percent(q1_num$n/nrow(df), .01), 'of the study cohort.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checkpoint 1: Identifying Earnings in the Fourth Quarter after Graduation\n",
    "\n",
    "Given the code above, create a data subset `q4_wages` that contains all earnings for the cohort in their fourth quarter after graduation. How many members of our cohort had positive earnings in this quarter? Do you expect this number to be higher or lower than the number in the first quarter?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace __ in the code below with the appropriate quarter\n",
    "\n",
    "q4_wages <- df_wages %>%\n",
    "    filter(quarter_number == __)\n",
    "\n",
    "# replace __ in the code below with the appropriate function\n",
    "q4_wages <- q4_wages %>%\n",
    "    ___(hashed_ssn) %>%\n",
    "    summarize(tot_wages = sum(wage)) %>%\n",
    "    ungroup()\n",
    "\n",
    "# replace __ in the code below with the appropriate variable\n",
    "q4_num <- __ %>%\n",
    "    summarize(n=n_distinct(hashed_ssn))\n",
    "\n",
    "cat('The total graduates with positive earnings during their fourth quarter after graduation:', q4_num$n)\n",
    "cat('\\nThat is', percent(q4_num$n/nrow(df), .01), 'of the study cohort.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add graduates without positive earnings for Q1\n",
    "\n",
    "Our current data frame, `q1_wages`, only contains individuals with positive earnings in their first quarter after graduation in New Jersey. Let's add in members of our cohort who did not appear in New Jersey's wage records during this time period, as well the additional variables from the original cohort table to better describe the individuals. This will let us easily analyze different earnings distributions in the cohort's first quarter after graduation as we progress throughout this notebook.\n",
    "\n",
    "We can do so by using a `left_join()` of the original cohort, `df`, to `q1_wages`, as this will add in one row for each `coleridge_id` in the original cohort that was not included in `q1_wages`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add grads without positive earnings\n",
    "q1_all_wages <- df %>%\n",
    "    left_join(q1_wages, c(\"hashed_ssn\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a quick check, we can see if the number of individuals in `q1_all_wages` that either have or do not have null wages makes sense given the total number of individuals in the cohort that were in `q1_wages`. We can do so by adding in an indicator variable if the `wages` column was null for each potential wage record in `q1_all_wages`, and then counting the number of distinct individuals based on this new variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# employment outcomes for all of those in our original cohort\n",
    "q1_all_wages %>%\n",
    "    mutate(wage_ind = ifelse(is.na(tot_wages), 'no', 'yes')) %>%\n",
    "    group_by(wage_ind) %>%\n",
    "    summarize(n=n_distinct(hashed_ssn))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1_num$n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that these numbers make sense. If they did not add up, chances are there was an issue with the details of your join."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to confirm, we can check to see if the number of rows in `q1_all_wages` is equal to the number of rows in `df`, the original cohort, as each individual in the original cohort should correspond to a single row regardless of employment status."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrow(df) == nrow(q1_all_wages)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also check to see if we have any missing values for our demographic variables. If so, let's fill these in as `unknown` so they won't be dropped in future analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see number of na values by column\n",
    "colSums(is.na(q1_all_wages))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1_all_wages<-q1_all_wages %>%\n",
    "    replace_na(list(\n",
    "        hashed_njsmartid = 'U',\n",
    "        hashed_instid ='U',\n",
    "        birthyear = 'U',\n",
    "        admstatus = 'U',\n",
    "        yr_matriculation = 'U',\n",
    "        sem_matriculation = 'U',\n",
    "        accumulatedcredit = 'U',\n",
    "        accumlatedgpa = 'U',\n",
    "        hisp = 'U',\n",
    "        race_ai_alask = 'U',\n",
    "        race_as = 'U',\n",
    "        race_aa = 'U',\n",
    "        race_pi = 'U',\n",
    "        race_wh = 'U',\n",
    "        race_single = 'U'\n",
    "    )\n",
    "              )\n",
    "\n",
    "# see na distribution now\n",
    "colSums(is.na(q1_all_wages))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Theoretically, you could apply these imputation methods to these missing demographic values. However, for the purposes of this notebook, we will focus our imputation techniques on missing earnings values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checkpoint 2: Replicate for Q4 ###\n",
    "\n",
    "Create a data frame `q4_all_wages` that mirrors `q1_all_wages` except for Q4. Feel free to add in as many code cells as you deem necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace __ in the code below with the appropriate word\n",
    "\n",
    "q4_all_wages <- df %>%\n",
    "   ___(q4_wages, c(\"hashed_ssn\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see number of na values by column by replacing __ with the appropriate variable\n",
    "colSums(is.na(___))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for each variable with missing data fill in with unknown, as shown above. \n",
    "\n",
    "q4_all_wages<-q4_all_wages %>%\n",
    "    replace_na(list(\n",
    "        ___ = 'U',\n",
    "        ___ ='U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U',\n",
    "        ___ = 'U'\n",
    "    )\n",
    "              )\n",
    "\n",
    "# see na distribution now\n",
    "colSums(is.na(q4_all_wages))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Impute Wage Values\n",
    "\n",
    "Now that we have confirmed that our `q1_all_wages` dataframe is ready to use for testing our imputation methods, we can get started. To recall, here are the four methods we will be trying out in this notebook:\n",
    "1. Dropping all people with \"missing\" values on the variable of interest (Q1 wages)\n",
    "2. Filling in zero for people who do not have records in New Jersey UI data\n",
    "3. Filling in missing values with the average New Jersey UI earnings of people who are in the same degree fields and have the same gender\n",
    "4. Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Drop All Missing Values\n",
    "\n",
    "First, let's look at the earnings outcomes during first quarter after graduation when we drop all missing earnings values. Here, by ignoring potentially non-missing values, we are hoping that they mirror the same distribution as the present one. Although this is fairly common, you should **never, ever, ever** use this method in practice. \n",
    "\n",
    "> Deleting missing values is often called listwise deletion and essentially assumes that missing values are missing completely at random (MCAR). For a scholarly treatment of this issue, see (amongst others): \n",
    "> - Rubens (1976) \"Inference and Missing Data\" for the initial presentation, or\n",
    "> - Peugh and Enders (2004) \"Missing Data in Educational Research: A Review of Reporting Practices and Suggestions for Improvement\" for a more recent discussion.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop missing values\n",
    "q1_no_missing <- q1_all_wages %>% \n",
    "    filter(!is.na(tot_wages))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see earnings distribution\n",
    "summary(q1_no_missing$tot_wages)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checkpoint 3: Replicate for Q4\n",
    "\n",
    "What does the earnings distribution look like for Q4 when you drop missing values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace __ with the appropriate variable\n",
    "\n",
    "q4_no_missing <- ___ %>% \n",
    "    filter(!is.na(tot_wages))\n",
    "\n",
    "summary(q4_no_missing$tot_wages)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Fill in Missing Values with Zero\n",
    "\n",
    "Next, let's see how the earnings distribution shifts when we encode all missing earnings outcomes as 0. Here, we are assuming that all missing earnings are due to unemployment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill all null tot_wages with 0\n",
    "q1_wages_zero <- q1_all_wages %>%\n",
    "    mutate(tot_wages = ifelse(is.na(tot_wages), 0, tot_wages)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take a look at the distribution. How does it vary from the distribution you get in method 1?\n",
    "summary(q1_wages_zero$tot_wages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat('Average earnings if missing wages are dropped is $', round(mean(q1_no_missing$tot_wages), 2), sep = '', '.')\n",
    "\n",
    "cat('\\nAverage earnings if missing wages are imputed as 0 is $', round(mean(q1_wages_zero$tot_wages), 2), sep = '', '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checkpoint 4: Replicate for Q4 ####\n",
    "\n",
    "What does the earnings distribution look like for Q4 when you fill missing values with zero?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace __ with the appropriate function\n",
    "\n",
    "q4_wages_zero <- q4_all_wages %>%\n",
    "    mutate(tot_wages = ifelse(___(tot_wages), 0, tot_wages)) \n",
    "summary(q4_wages_zero$tot_wages)\n",
    "\n",
    "cat('Average earnings if missing wages are dropped is $', round(mean(q4_no_missing$tot_wages), 2), sep = '', '.')\n",
    "\n",
    "cat('\\nAverage earnings if missing wages are imputed as 0 is $', round(mean(q4_wages_zero$tot_wages), 2), sep = '', '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Fill in Missing Values with Major/Gender Mean Earnings\n",
    "\n",
    "Now, instead of either ignoring missing values or assuming the earnings are 0, we will try imputing missing earnings for each individual as the average quarterly earnings of the other individuals in our cohort of the same `sex` and `major_title`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, our strategy is as follows:\n",
    "- Using populated wages, find mean earnings for each major by gender\n",
    "- Merge the mean earnings, based on major and gender, to the overall cohort\n",
    " - creates an additional column `mean_wages`\n",
    "- Recode so that missing values are populated with mean earnings\n",
    " - data stored in a new column `imputed_wages`\n",
    "\n",
    "\n",
    ">Note: This process is frequently termed mean imputation. Implementing this method will compress the variance and covariance of the imputed variable, resulting in biased parameter estimates for all parameters except the mean (Peugh & Enders, 2004, p.529). In this example, we are assuming that the missing values in wages are conditional on both gender and major. We also assume that the missingness is not truly indicative of lack of wages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the sex column, since group REDACTED has so few values, we will drop group REDACTED to create a binary variable and recode the remaining groups so Male = 1 and Female = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter out sex = 0 and recode\n",
    "q1_all_wages <- q1_all_wages %>%\n",
    "    filter(sex != 0) %>%\n",
    "    mutate(\n",
    "        sex = ifelse(sex == '1', '1', '0')\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean earnings by gender/kpeds_major1 grouping\n",
    "q1_all_wages %>%\n",
    "    group_by(sex, major_title) %>%\n",
    "    summarize(mean_wages = mean(tot_wages, na.rm=T)) %>%\n",
    "    head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mean earnings by gender/kpeds_major1 grouping saved\n",
    "q1_major_gend <- q1_all_wages %>%\n",
    "    group_by(sex, major_title) %>%\n",
    "    summarize(mean_wages = mean(tot_wages, na.rm=T)) %>%\n",
    "    ungroup()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will merge the two DataFrames, `q1_major_gend` and `q1_all_wages` using `inner_join`.\n",
    "> Note: `left_join()` would also work in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see if join works\n",
    "q1_all_wages %>%\n",
    "    inner_join(q1_major_gend, by=c('sex', 'major_title')) %>%\n",
    "    head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save join results to q1_joined_major_gend\n",
    "q1_joined_major_gend <- q1_all_wages %>%\n",
    "    inner_join(q1_major_gend, by=c('sex', 'major_title'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can add a new column to `q1_joined_major_gend` to include the mean wage, based on gender and major, *if* the individual did not appear in the New Jersey UI wage records data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see if mutation works as designed\n",
    "q1_joined_major_gend %>%\n",
    "    mutate(imputed_wages = ifelse(is.na(tot_wages), mean_wages, tot_wages)) %>%\n",
    "    select(tot_wages, mean_wages, imputed_wages) %>%\n",
    "    head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save mutation to q1_major_gend_impute\n",
    "q1_major_gend_impute <- q1_joined_major_gend %>%\n",
    "    mutate(imputed_wages = ifelse(is.na(tot_wages), mean_wages, tot_wages))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In using this method, there is a chance we cannot impute missing values for all individuals in the cohort. If `imputed_wages` is still `NA`, we can assume there were no individuals in the cohort with non-missing earnings with the same major/gender combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see if any still don't have imputed earnings\n",
    "q1_major_gend_impute %>%\n",
    "    filter(is.na(imputed_wages)) %>%\n",
    "    summarize(n=n())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, it seems as though we do not have available earnings for every combination of gender and primary degree. For the sake of the exercise, we will ignore the earnings of those whose we could not impute using this method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see distribution\n",
    "summary(q1_major_gend_impute$imputed_wages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat('Average earnings if missing wages are dropped is $', round(mean(q1_no_missing$tot_wages), 2), sep = '', '.')\n",
    "\n",
    "cat('\\nAverage earnings if missing wages are imputed as 0 is $', round(mean(q1_wages_zero$tot_wages), 2), sep = '', '.')\n",
    "\n",
    "cat('\\nAverage earnings if missing wages are imputed using major/gender means earnings is $', round(mean(q1_major_gend_impute$imputed_wages, na.rm=T), 2), sep = '', '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Checkpoint 5: Replicate for Q4 ####\n",
    "\n",
    "Impute missing earnings values as the mean earnings of individuals in the cohort with the same gender (`gender`) and degree designation (`kpeds_major1`) in quarter 4. What does the earnings distribution look like? For how many individuals could you not impute values using this method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill in __ with appropriate variables \n",
    "\n",
    "# recode sex\n",
    "q4_all_wages <- q4_all_wages %>%\n",
    "    filter(sex != 0)\n",
    "\n",
    "q4_all_wages <- q4_all_wages %>% \n",
    "    mutate(\n",
    "        sex=ifelse(sex == 1, '1', '0')\n",
    "    )\n",
    "\n",
    "#mean earnings by gender/kpeds_major1 grouping saved\n",
    "q4_major_gend <- q4_all_wages %>%\n",
    "    group_by(___, ___) %>%\n",
    "    summarize(mean_wages = mean(tot_wages, na.rm=T)) %>%\n",
    "    ungroup()\n",
    "\n",
    "\n",
    "# save join results to q1_joined_major_gend\n",
    "q4_joined_major_gend <- q4_all_wages %>%\n",
    "    inner_join(q4_major_gend, by=c('sex', 'major_title'))\n",
    "\n",
    "# save mutation to q1_major_gend_impute\n",
    "q4_major_gend_impute <- q4_joined_major_gend %>%\n",
    "    mutate(imputed_wages = ifelse(is.na(tot_wages), mean_wages, tot_wages))\n",
    "\n",
    "\n",
    "\n",
    "# see if any still don't have imputed earnings\n",
    "q4_major_gend_impute %>%\n",
    "    filter(is.na(imputed_wages)) %>%\n",
    "    summarize(n=n())\n",
    "q4_major_gend_impute <- q4_major_gend_impute %>% drop_na()\n",
    "\n",
    "# see earnings distribution\n",
    "summary(q4_major_gend_impute$imputed_wages)\n",
    "\n",
    "cat('Average earnings if missing wages are dropped is $', round(mean(q4_no_missing$tot_wages), 2), sep = '', '.')\n",
    "\n",
    "cat('\\nAverage earnings if missing wages are imputed as 0 is $', round(mean(q4_wages_zero$tot_wages), 2), sep = '', '.')\n",
    "\n",
    "cat('\\nAverage earnings if missing wages are imputed using major/gender means earnings is $', round(mean(q4_major_gend_impute$imputed_wages), 2), sep = '', '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Regression imputation\n",
    "\n",
    "We can also use regression to try to get more accurate earnings values. We will build a regression equation from the observations for which we know the earnings, then use the equation to predict the missing earnings values. This is, in effect, an extension of the mean imputation by subgroup. Here, we will use demographic information of graduates such as birth year, sex, major, and time of graduation.\n",
    "\n",
    "> Note: We will not be checking the assumptions associated with linear regressions, as this example is aimed at merely displaying how to use a linear regression for imputation. If you plan on using regression imputation, please check all assumptions before employing a predictive model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset to variables included in regression analysis\n",
    "q1_reg <- q1_all_wages %>%\n",
    "    select(hashed_ssn, tot_wages, birthyear, sex, major_title, grad_date)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make sense of the linear regression, we will filter for the top 5 majors and then group all other majors together. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# finding top 5 majors\n",
    "top_5_majors <- q1_reg %>%\n",
    "    count(major_title) %>% \n",
    "    arrange(desc(n)) %>%\n",
    "    head(5)\n",
    "\n",
    "top_5_majors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating an additional column with the top 5 majors and all others assigned to 6.\n",
    "\n",
    "q1_reg <- q1_reg  %>%\n",
    "    mutate(\n",
    "        major_group = ifelse(q1_reg$major_title %in% top_5_majors$major_title, major_title, \"Other\")\n",
    "    )\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see variable types and outputs\n",
    "q1_reg %>%\n",
    "    glimpse()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, it may make sense for `birthyear` to be a numeric variable rather than a character vector, as there may be some predictive power in numerically analyzing the ages of the graduates. Let's change `birthyear` to a `numeric` variable.\n",
    "\n",
    "> As you saw above, there are a few individuals in the cohort with an unknown birthdate, so they will not be included in the predictive portion of this analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert birthyear to numeric\n",
    "q1_reg <- q1_reg %>%\n",
    "    mutate(\n",
    "        birthyear = as.numeric(birthyear)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# don't need tot_wages because they are null \n",
    "q1_wages_na <- q1_reg %>%\n",
    "    filter(is.na(tot_wages)) %>%\n",
    "    select(-c(tot_wages))\n",
    "\n",
    "# removing NAs from tot_wages\n",
    "q1_wages_pred <- q1_reg %>%\n",
    "    filter(!is.na(tot_wages))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model creation process for a linear regression can be done using the `lm()` function. The variable we are trying to predict is on the left-hand side of `lm()` before the `~`, and the predictors are all of the variables on the right-hand side of the `~`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model and fit coefficients\n",
    "q1_wages_model <- lm(tot_wages ~ birthyear + sex + major_group, data = q1_reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see model summary\n",
    "summary(q1_wages_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Part of regression-based imputation is to evaluate your model for any unusual relationships. Examining the above result suggest that younger graduates tend to earn less, males tend to earn more, and business majors and registered nurses earn more than the comparison group (accounting majors) while the other major grouping earns less. While there is certainly more we could add to inform this model the sign of these coefficients make theoretical sense. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have fit coefficients for each of the predictors in the model, we can predict the `tot_wages` variable for the missing data using `predict()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict earnings\n",
    "pred_earnings <- data.frame(tot_wages = predict(q1_wages_model, newdata=q1_wages_na))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(pred_earnings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the output for `predict()` retains the same order of rows from `q1_wages_na`, we can add the `tot_wages` variable from `pred_earnings` into the existing `q1_wages_na` data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see updated data frame with predicted earnings\n",
    "cbind(q1_wages_na, pred_earnings) %>% \n",
    "    head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save updated data frame with predicted earnings\n",
    "q1_wages_na_w_earnings <- cbind(q1_wages_na, pred_earnings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, before we can see the effects of the imputation method, we need to combine our training set, which already has `tot_wages`, with our testing set and its predicted `tot_wages`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine the q1_wages_na_w_earning and q1_wages_pred\n",
    "rbind(q1_wages_na_w_earnings, q1_wages_pred) %>% \n",
    "    head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the combined dataframes\n",
    "q1_reg_earnings <- rbind(q1_wages_na_w_earnings, q1_wages_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see the entire earnings distribution for the cohort after applying regression imputation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see earnings distribution for full cohort\n",
    "summary(q1_reg_earnings$tot_wages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see earnings distribution for imputed portion of cohort\n",
    "summary(q1_wages_na_w_earnings$tot_wages)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Earnings Distributions\n",
    "\n",
    "We can quickly determine if these different imputation methods significantly altered the pre-imputation wage distribution by visualizing the overall earnings distribution. Plotting side-by-side boxplots can be an effective choice. To do so, we need to bind the earnings from all of these methods by rows, meaning they must have the same columns. For the sake of simplicity, we will have three columns in this data frame:\n",
    "\n",
    "- `hashed_ssn`, the person identifier\n",
    "- `tot_wages`, cumulative earnings in first quarter post-graduation\n",
    "- `method`, type of imputation method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adapt q1_no_missing\n",
    "q1_no_missing %>%\n",
    "    select(hashed_ssn, tot_wages) %>% head()\n",
    "\n",
    "q1_no_missing <- q1_no_missing %>%\n",
    "    select(hashed_ssn, tot_wages) %>%\n",
    "    mutate(method = 'remove missing')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adapt q1_reg_earnings\n",
    "q1_reg_earnings%>%\n",
    "    select(hashed_ssn, tot_wages) %>% head()\n",
    "\n",
    "q1_reg_earnings <- q1_reg_earnings %>%\n",
    "    select(hashed_ssn, tot_wages) %>%\n",
    "    mutate(method = 'regression')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adapt q1_wages_zero\n",
    "q1_wages_zero %>%\n",
    "    select(hashed_ssn, tot_wages) %>% head()\n",
    "\n",
    "q1_wages_zero <- q1_wages_zero %>%\n",
    "    select(hashed_ssn, tot_wages) %>%\n",
    "    mutate(method = 'zero')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adapt q1_major_gend_impute\n",
    "q1_major_gend_impute %>% select(hashed_ssn, imputed_wages) %>% rename(tot_wages = imputed_wages) %>% head()\n",
    "\n",
    "q1_major_gend_impute <- q1_major_gend_impute %>%\n",
    "    select(hashed_ssn, tot_wages) %>%\n",
    "    mutate(method = 'mean')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that these methods all have the same column names, we can feed them into `rbind()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine earnings from all methods\n",
    "all_methods <- rbind(q1_major_gend_impute, q1_reg_earnings, q1_no_missing, q1_wages_zero)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before visualizing these distributions, we will filter out extreme outliers that would affect the rest of the visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# removing rows that have a tot_wages values greater then 50000\n",
    "\n",
    "all_methods <- all_methods %>% \n",
    "    filter(tot_wages < 50000)\n",
    "\n",
    "# check the max value\n",
    "max(all_methods$tot_wages, na.rm=T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# boxplot of all methods\n",
    "all_methods %>%\n",
    "    ggplot(aes(x=method, y = tot_wages)) +\n",
    "    geom_boxplot() + \n",
    "    labs(\n",
    "        title = \"The REDACTED Earnings Distribution's Quartiles are REDACTED REDACTED across \\n imputation methods\",\n",
    "        x='Imputation Method',\n",
    "        y='Quarter Earnings',\n",
    "        caption = 'Source: NJ UI wage records data'\n",
    "    ) +\n",
    "    theme_minimal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple histograms\n",
    "\n",
    "We can also look at the differences in the earnings distribution by looking at side-by-side histograms. Instead of using the `geom_` layer `geom_boxplot()`, we will use `geom_histogram()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_methods %>%\n",
    "    ggplot(aes(x=tot_wages)) +\n",
    "    geom_histogram() + \n",
    "    facet_grid(method ~ .) +\n",
    "    labs(\n",
    "        title = 'Zero Imputation has a REDACTED REDACTED on the overall earnings distribution',\n",
    "        y = 'Number of Workers',\n",
    "        x='Quarterly Wages',\n",
    "        caption = 'Source: NJ UI wage records data'\n",
    "    ) +\n",
    "    theme_minimal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Optional) Advanced: Using machine learning to impute values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To impute values, we can also use machine learning algorithms such as `K-nearest Neighbors` and `Decision Trees`. The principle behind `K-nearest Neighbors` is quite simple: the missing values can be imputed by values of \"closest neighbors\" - as approximated by other, known, features. \n",
    "\n",
    "For example, if we had cases where the data on earnings of some graduates was completely missing, we could approximate their earnings by referring to other characteristics which could be shared by major group (their 'closest neighbors' in terms of characteristics).\n",
    "\n",
    "The algorithm calculates the distance between the input values (the missing values) and helps to identify the nearest possible value based on other features (such as known characteristics of the closest major group). Imputing missing data using machine learning has become a research hotbed, and there are plenty of papers covering the various algorithms if you are curious."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}