5.Inference_Imputation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img style=\"float: center;\" src=\"images/CI_horizontal.png\" width=\"400\">\n",
    "</center>\n",
    "<center>\n",
    "    <span style=\"font-size: 1.5em;\">\n",
    "        <a href='https://www.coleridgeinitiative.org'>Website</a>\n",
    "    </span>\n",
    "</center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> Joshua Edelmann, Rukhshan Mian, Benjamin Feder </center>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6407277\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6407277.svg\" alt=\"DOI\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imputation, Outcome Measurement, and Inference \n",
    "\n",
    "## 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 2015-16 Tennessee 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: Wages notebook, you initially examined the earnings distribution for all members of this cohort who had positive earnings in this time period in Tennessee. To evaluate the earnings outcomes of all 2015-16 Tennessee graduates, you need to decide what to do when you cannot find their earnings in the Tennessee Unemployment Insurance (UI) wage records. A person may not appear in Tennessee'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 Tennessee.\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 Tennessee 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\n",
    "- Adding in Ohio, Kentucky, and Missouri UI wage records for the cohort in question"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R Setup and Server 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)\n",
    "                                 "
   ]
  },
  {
   "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 Data Exploration: Wages notebook, we will create a new column, `quarter_number`, to track the quarter after graduation 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_tn_2021.dbo.nb_cohort_wages_link;\n",
    "\"\n",
    "\n",
    "df_wages <- dbGetQuery(con, qry)"
   ]
  },
  {
   "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 limited to 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we will want to estimate the total earnings for each `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(SSN) %>%\n",
    "    summarize(tot_wages = sum(wge_amt)) %>%\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(SSN)\n",
    "    )\n",
    "\n",
    "cat('The total graduates with positive earnings during their first quarter after graduation:', q1_num$n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compare the total graduate with positive earnings with the full cohort."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in cohort\n",
    "qry <- \"\n",
    "SELECT * \n",
    "FROM tr_tn_2021.dbo.grads1516_dated;\n",
    "\"\n",
    "\n",
    "df <- dbGetQuery(con, qry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat('The total number of graduates with positive earnings during their first quarter after graduation 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",
    "# filter by quarter number 4\n",
    "q4_wages <- df_wages %>%\n",
    "    filter(quarter_number == __)\n",
    "\n",
    "# replace __ in the code below with the appropriate function\n",
    "# grouping by SSN and summing wages\n",
    "q4_wages <- q4_wages %>%\n",
    "    ___(SSN) %>%\n",
    "    summarize(tot_wages = sum(wge_amt)) %>%\n",
    "    ungroup()\n",
    "\n",
    "# replace __ in the code below with the appropriate variable\n",
    "# count total number of distinct SSN\n",
    "q4_num <- __ %>%\n",
    "    summarize(n=n_distinct(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 Tennessee. Let's add in members of our cohort who did not appear in Tennessee'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 conduct different imputation methods on the cohort's earnings in their first quarter after graduation as we progress throughout this notebook.\n",
    "\n",
    "We can add in those with missing wages by using a `left_join()` of the original cohort, `df`, to `q1_wages`, as this will add in one row for each `SSN` 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(\"SSN\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a quick check, we can see if the number of individuals in `q1_all_wages` that have wages is equal to the total number of individuals in the cohort that were in `q1_wages`. We can do this by adding in an indicator variable, `wage_ind`, 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(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",
    "        Gender = 'U',\n",
    "        BirthYear ='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",
    "# add grads without positive earnings\n",
    "\n",
    "q4_all_wages <- df %>%\n",
    "   ___(q4_wages, c(\"SSN\"))\n",
    "\n",
    "# see number of na values by column\n",
    "colSums(is.na(____))\n",
    "\n",
    "\n",
    "# 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",
    "    )\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 five 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 Tennessee UI data\n",
    "3. Filling in missing values with the average Tennessee UI earnings of people who are in the same degree fields and have the same gender\n",
    "4. Regression\n",
    "5. Adding in Ohio, Kentucky, and Missouri UI wage records for the cohort in question"
   ]
  },
  {
   "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))\n"
   ]
  },
  {
   "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",
    "# drop missing values\n",
    "q4_no_missing <- ___ %>% \n",
    "    filter(!is.na(tot_wages))\n",
    "\n",
    "# see earnings distribution\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",
    "# fill all null tot_wages with 0\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",
    "# Take a look at the distribution. How does it vary from the distribution you get in method 1?\n",
    "cat('Average earnings if missing wages are dropped is $', round(mean(q4_no_missing$tot_wages), 2), sep = '', '.')\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 `Gender` and `CIP_Family`."
   ]
  },
  {
   "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",
    " - Create an additional column `mean_wages`\n",
    "- Recode so that missing values are populated with mean earnings\n",
    " - Store data 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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean earnings by gender/CIP_Family\n",
    "q1_all_wages %>%\n",
    "    group_by(Gender, CIP_Family) %>%\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/CIP_Family\n",
    "q1_major_gend <- q1_all_wages %>%\n",
    "    group_by(Gender, CIP_Family) %>%\n",
    "    summarize(mean_wages = mean(tot_wages, na.rm=T)) %>%\n",
    "    ungroup()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now merge the two data frames, `q1_major_gend` and `q1_all_wages` using `inner_join`. Here, we are adding in the imputed wages from `q1_major_gend` to the data frame `q1_all_wages`\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('Gender', 'CIP_Family')) %>%\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('Gender', 'CIP_Family'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 Tennessee 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": [
    "It appears we do not have available earnings for every combination of gender and primary degree."
   ]
  },
  {
   "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=TRUE), 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 (`CIP_Family`) 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",
    "\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",
    "    _____(q4_major_gend, by=c('Gender', '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(____(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",
    "    ____(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, na.rm=TRUE), 2), sep = '', '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Regression imputation\n",
    "\n",
    "We can also use regression to try to get more accurate imputed earnings. We will build a regression equation from the obervations 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 information of graduates such as birth year, gender, major, and time of graduation as input variables to build out the regression.\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(SSN, tot_wages, BirthYear, Gender, CIP_Family, grad_date)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ease of interpreting the linear regression results, 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(CIP_Family) %>% \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 other.\n",
    "\n",
    "q1_reg <- q1_reg  %>%\n",
    "    mutate(\n",
    "        major_group = ifelse(q1_reg$CIP_Family %in% top_5_majors$CIP_Family, CIP_Family, \"Other\")\n",
    "    ) %>%\n",
    "    select(-CIP_Family)\n",
    "    "
   ]
  },
  {
   "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",
    "    )\n",
    "\n",
    "nrow(q1_reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: NAs introduced by coercion means when we change `BirthYear` from a character vector to a numeric variable, R forced all `BirthYear` values that equal 'U' to equal NA. (Previously, we recoded all missing values to equal 'U'). This is because R did not know how turn the 'U' value into a numeric variable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1_reg %>% summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to filter out all `BirthYear` values that are NA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1_reg <- q1_reg %>% \n",
    "    filter(!is.na(BirthYear))\n",
    "nrow(q1_reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: In this example, all instances where `BirthYear` is unknown `Gender` is unknown as well. Because of this you will not see any `U` in `Gender` when looking at coefficients. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we are creating two data frames to predict the missing wages. `q1_wages_na` is the data frame that we need to predict wages for and `q1_wages_pred` is the data frame that is used to create the linear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1_wages_na <- q1_reg %>%\n",
    "    filter(is.na(tot_wages)) %>%\n",
    "# don't need tot_wages because they are null \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 + Gender + major_group, data = q1_wages_pred)"
   ]
  },
  {
   "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 suggets that younger graduates tend to earn less, males tend to earn more, and engineering majors and health professionals earn more than the comparison group (business management) 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 data frame `q1_wages_pred`, which already has `tot_wages`, with our data frame `q1_wages_na_w_earnings` 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": [
    "### 5. Add in Ohio, Kentucky, and Missouri UI data\n",
    "\n",
    "Finally, let's see how the earnings distribution changes when we add in some bordering states' UI wage records. We will leverage wage records from Ohio, Missouri, and Kentucky to find earnings in these states for members of our cohort. From there, we can combine the tables and analyze the resulting changes to the overall earnings distribution. \n",
    "\n",
    "Recall that in the `2.Data_Exploration_Wages` notebook, we already found wage records for those in our cohort in Kentucky. The resulting data frame has been saved as the permanent table `ky_cohort_link`, which was created from the following code:\n",
    "\n",
    "    SELECT cohort.*, w.naics, w.wages, w.coleridge_id, w.fein, w.job_date, 'KY' AS state\n",
    "    INTO tr_tn_2021.dbo.ky_cohort_link\n",
    "    FROM tr_tn_2021.dbo.grads1516_dated cohort\n",
    "    JOIN tr_tn_2021.dbo.ky_ui_wages_dated w\n",
    "    on cohort.SSN = w.ssn \n",
    "    WHERE w.job_date >= cohort.grad_date and DATEADD(quarter, 13, cohort.grad_date) >= w.job_date and w.wages > 0;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also performed a similar process to link to Ohio's and Missouri's wage records, and saved the results to the permanent tables `oh_wages` and `mo_wages`, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    SELECT cohort.*, oh_wages.ssn_hash, oh_wages.wages, \n",
    "    convert(datetime, concat(oh_wages.quarter*3-2, '/', '01', '/', oh_wages.year)) AS job_date, 'OH' AS state\n",
    "    INTO tr_tn_2021.dbo.oh_wages \n",
    "    FROM tr_tn_2021.dbo.grads1516_dated cohort\n",
    "    JOIN ds_oh_djfs.dbo.ui_wage_by_employer_olda oh_wages\n",
    "    ON cohort.SSN = oh_wages.ssn_hash;\n",
    "    \n",
    "\n",
    "    SELECT tn_cohort.*, mo_wages.wage,\n",
    "    convert(datetime, concat(mo_wages.quarter*3-2, '/', '01', '/', mo_wages.year)) AS job_date, 'MO' AS state\n",
    "    INTO tr_tn_2021.dbo.mo_wages \n",
    "    FROM tr_tn_2021.dbo.grads1516_dated tn_cohort\n",
    "    JOIN ds_mo_des.dbo.ui_wages mo_wages\n",
    "    ON tn_cohort.SSN = mo_wages.ssn;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With wage information subset to only records from our cohort, we can combine the three tables in SQL using `UNION ALL`. Once combined, we can explore the resulting earnings distribution for our cohort."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine mo, oh, and ky wage records\n",
    "qry = \"\n",
    "SELECT SSN, CIP_Family, MajorGroups, Gender, BirthYear, grad_date, wage as wages, job_date, state\n",
    "FROM tr_tn_2021.dbo.mo_wages\n",
    "WHERE job_date >= grad_date and DATEADD(quarter, 13, grad_date) >= job_date and wage > 0\n",
    "UNION ALL\n",
    "SELECT SSN, CIP_Family, MajorGroups, Gender, BirthYear, grad_date, wages, job_date, state\n",
    "FROM  tr_tn_2021.dbo.oh_wages\n",
    "WHERE job_date >= grad_date and DATEADD(quarter, 13, grad_date) >= job_date and wages > 0\n",
    "UNION ALL\n",
    "SELECT SSN, CIP_Family, MajorGroups, Gender, BirthYear, grad_date, wages, job_date, state \n",
    "FROM tr_tn_2021.dbo.ky_cohort_link;\n",
    "\"\n",
    "combined_wages <- dbGetQuery(con, qry)\n",
    "head(combined_wages)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how many distinct `SSN` values we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_wages %>%\n",
    "    summarize(n=n_distinct(SSN))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have combined wage records for contiguous states together, we can extend this to include Tennessee's wage records as well. Before doing so, we need to rename the `wge_amt` variable to `wages` because we need to have common variable names when we combine the data across the 4 states. Otherwise, `rbind` will not work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename the wge_amt column to wages\n",
    "df_wages <- df_wages %>%\n",
    "    rename(wages = wge_amt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also add in the `state` variable to the Tennessee wages (`df_wages`), to designate that the wage records are from Tennessee. We will also select for key variables that are in the `combined_wages` data frame that we created above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add in state for tn wages\n",
    "df_wages$state = \"TN\"\n",
    "\n",
    "# select key variables\n",
    "df_wages <- df_wages %>%\n",
    "select(SSN, CIP_Family, MajorGroups, Gender, BirthYear, grad_date, wages, job_date, state, quarter_number)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `combined_wages` does not have `quarter_number`, which tracks the quarter relative to graduation of employment. We can add it in using the same code we applied to `df_wages`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add in quarter_number variable into combined_wages\n",
    "combined_wages <- combined_wages %>%\n",
    "    mutate(\n",
    "        quarter_number = round(as.double(difftime(as.Date(job_date), as.Date(grad_date), units = \"weeks\")/13), 0)\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine our two data frames \n",
    "combined_wages <- rbind(combined_wages, df_wages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter for first quarter\n",
    "q1_combined_wages <- combined_wages %>%\n",
    "    filter(quarter_number == 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before looking into the earnings distribution, let's see the amount of individuals for which we have wage records in their first quarter after graduation by state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of individuals by state first quarter after graduation\n",
    "q1_combined_wages %>%\n",
    "    group_by(state) %>%\n",
    "    summarize(n=n_distinct(SSN)) %>%\n",
    "    arrange(desc(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we do not have many records across states, the additional information will provide more known earnings to our cohort. We can find the aggregate earnings for every individual in the cohort across these four states in their first quarter after graduation by aggregating `q1_combined_wages` by each `SSN`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's see the earnings distribution after we add UI records from other states\n",
    "q1_combined_agg_wages <- q1_combined_wages %>%\n",
    "    group_by(SSN) %>%\n",
    "    summarize(tot_wages = sum(wages)) %>%\n",
    "    ungroup()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how many distinct `SSN` values are in `q1_combined_agg_wages`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see number of distinct coleridge_id values \n",
    "n_distinct(q1_combined_agg_wages$SSN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now look at how the earnings distribution has changed with the addition of wage records from contiguous states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(q1_combined_agg_wages$tot_wages)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you may recall, the number of unique `SSN` values in `q1_combined_wages` is not the same number as in the original cohort. Although we combined all earnings observations for these four states, we are still missing a portion of the original cohort that did not appear in any of these states' UI wage records in their first quarter after graduation.\n",
    "\n",
    "To allow for reasonable comparison, let's add in these individuals using another `left_join()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add in those present in df but not q1_combined_wages\n",
    "q1_all_combined_wages <- df %>%\n",
    "    left_join(q1_combined_agg_wages, 'SSN')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see number of unique coleridge_ids\n",
    "n_distinct(q1_all_combined_wages$SSN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat('By adding in UI wage records from a handful of bordering states, we have managed to find wage records for', \n",
    "   n_distinct(q1_combined_agg_wages$SSN) - n_distinct(q1_wages$SSN),\n",
    "   'more people, as well as augmented earnings for some others.')"
   ]
  },
  {
   "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",
    "- `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(SSN, tot_wages) %>% head()\n",
    "\n",
    "q1_no_missing <- q1_no_missing %>%\n",
    "    select(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(SSN, tot_wages) %>% head()\n",
    "\n",
    "q1_reg_earnings <- q1_reg_earnings %>%\n",
    "    select(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(SSN, tot_wages) %>% head()\n",
    "\n",
    "q1_wages_zero <- q1_wages_zero %>%\n",
    "    select(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(SSN, imputed_wages) %>% rename(tot_wages = imputed_wages) %>% head()\n",
    "\n",
    "q1_major_gend_impute <- q1_major_gend_impute %>%\n",
    "    select(SSN, tot_wages) %>%\n",
    "    mutate(method = 'mean')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adapt to other states\n",
    "\n",
    "q1_combined_agg_wages <- q1_combined_agg_wages %>%\n",
    "    select(SSN, tot_wages) %>%\n",
    "    mutate(method = \"combined states\")\n"
   ]
  },
  {
   "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, q1_combined_agg_wages)\n",
    "max(all_methods$tot_wages, na.rm=TRUE)"
   ]
  },
  {
   "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 outliers that have a tot_wages values greater then REDACTED\n",
    "\n",
    "all_methods <- all_methods %>% \n",
    "    filter(tot_wages < REDACTED)\n",
    "\n",
    "# check the max value\n",
    "max(all_methods$tot_wages, na.rm=TRUE)"
   ]
  },
  {
   "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 Q1 Earnings Distribution's Quartiles are moderately affected across \\n imputation methods\",\n",
    "        x='Imputation Method',\n",
    "        y='Quarter Earnings',\n",
    "        caption = 'Source: TN/KY/OH/MO 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 significant change on the overall earnings distribution',\n",
    "        y='Number of Workers',\n",
    "        x='Q1 After Graduation Wages',\n",
    "        caption = 'Source: TN/KY/OH/MO 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": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Feder, Benjamin, & Barrett, Nathan. (2022, March 30). Imputation Using New Jersey Unemployment Insurance Wage Records. Zenodo. https://doi.org/10.5281/zenodo.6399344\n",
    "\n",
    "Peugh, J. L., & Enders, C. K. (2004). Missing Data in Educational Research: A Review of Reporting Practices and Suggestions for Improvement. _Review of Educational Research_, 74(4), 525-556. doi: 10.3102/00346543074004525\n",
    "\n",
    "Rubin, D. B. (1976). Inference and Missing Data. _Biometrika_, 63(3), 581-592. doi:10.2307/2335739"
   ]
  },
  {
   "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
}