05_Inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: center;\" src=\"images/CI_horizontal.png\" width=\"600\">\n",
    "<center>\n",
    "    <span style=\"font-size: 1.5em;\">\n",
    "        <a href='https://www.coleridgeinitiative.org'>Website</a>\n",
    "    </span>\n",
    "</center>\n",
    "\n",
    "<center>Brian Kim, Allison Nunez</center>\n",
    "<br>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6426187\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6426187.svg\" alt=\"DOI\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inference\n",
    "\n",
    "## Introduction \n",
    "\n",
    "Surveys are often used designed to glean estimates of an entire population based on the responses from a subset of individuals. If the survey population is just a subsample of the overall population, the survey results will be accompanied by weights to mirror the population in question. This notebook will cover how to properly incorporate survey weights into an analysis and the importance of considering the data-generating process in doing so.\n",
    "\n",
    "Additionally, how should missing values be treated in a dataset? Unfortunately, there is usually no *right* answer. However, it is possible to impute these missing values, providing a best guess for each missing point's true value. Here, the second part of the notebook will describe common imputation methods to use when approaching missing values in the future.\n",
    "\n",
    "### Learning Objectives\n",
    "\n",
    "* Leverage survey weights to get accurate statistics from the sample.\n",
    "\n",
    "* Explore options for imputing missing values.\n",
    "\n",
    "The first section of this notebook, [Survey Weights](06_01_Inference.ipynb/##Survey-Weights), will focus on the SDR data for the 2015 cohort of graduates, particularly on using the survey weights to estimate the overall post-graduation earnings distribution. The second part of the notebook, [Missing Data](06_01_Inference.ipynb/##Missing-Data), will shift to the SED data, as it will walk through imputing missing ages at the time of dissertation completion for the 2015 cohort.\n",
    "\n",
    "The imputation methods covered in this notebook include:\n",
    "   - Dropping all \"missing\" values\n",
    "   - Substituting missing values with the average age at dissertation completion for those in the same PhD field\n",
    "   - Regression imputation\n",
    "   - Mode imputation (for categorical variables)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R Setup\n",
    "\n",
    "As always, start by importing the required libraries, as well as creating the connection to the database.\n",
    "\n",
    "> Recall that the `survey` package can be used to work with survey weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#database interaction imports\n",
    "library(odbc)\n",
    "\n",
    "# for data manipulation/visualization\n",
    "library(tidyverse)\n",
    "\n",
    "# scaling data, calculating percentages, overriding default graphing\n",
    "library(scales)\n",
    "\n",
    "# adding weights \n",
    "library(survey)\n",
    "\n",
    "# to better view images\n",
    "# For easier viewing of graphs\n",
    "# Adjust repr.plot.width and repr.plot.height to change the size of graphs\n",
    "theme_set(theme_gray(base_size = 24))\n",
    "options(repr.plot.width = 20, repr.plot.height = 12)\n",
    "options(scipen=999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connect to the database\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": [
    "## Survey Weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned previously, the SDR data includes survey weights, which allows for the production of estimates for the total population of SED PhD recipients. In general, survey weights are used because the sample is not necessarily taken evenly from the population. Sometimes, researchers decide to intentionally oversample from certain subpopulations in order to make sure they have enough people from that group. Re-weighting is also done afterwards to adjust for non-response or other factors that may reduce the sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following section covers computing statistics with and without sampling weights to show how the results differ and demonstrate why using weights is so important."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Survey of Doctorate Recipients (SDR)**\n",
    "\n",
    "As the SDR data contains the sub-samples of the SED population, survey weights must always be used in the calculations.\n",
    "\n",
    "The example in this section will focus on estimating the distribution of earnings for the 2015 SED cohort. In the SDR data, the variable `sdryr` (the year of first award of a U.S. PhD degree) can be used to subset to those in the 2015 cohort. The `salary` and `wtsurvy` variables from the dataset will be needed for this computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the relevant variables from the SDR data to find the earnings \n",
    "# among the 2015 cohort\n",
    "\n",
    "query <- \"\n",
    "SELECT salary, wtsurvy\n",
    "FROM ds_nsf_ncses.dbo.nsf_sdr_2017\n",
    "WHERE sdryr = '2015' \n",
    "\"\n",
    "\n",
    "earnings_2015 <- dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the head of the table\n",
    "\n",
    "head(earnings_2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consult the data dictionary for the values that can be expected in the `salary` variable.  The Survey of Doctorate Recipients (SDR) data dictionary states that `9999998` is a value reserved for a 'Logical Skip'. Therefore, these values can be removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the rows with the logical skip value\n",
    "\n",
    "earnings_2015 <- earnings_2015 %>%\n",
    "    filter(salary != '9999998')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, as done in the Data Exploration notebook, the survey weights can be added using `svydesign` function from `survey` package to calculate the weighted earnings distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add weights using \"svydesign\" function from \"survey\" library\n",
    "\n",
    "weighted_earnings_2015 <- svydesign(ids=~1, data=earnings_2015, weights=earnings_2015$wtsurvy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earnings percentiles can be found using the `svyquantile` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find percentiles on the weighted data \n",
    "\n",
    "svyquantile(~salary, weighted_earnings_2015, c(0, .25, .50, .75, 1), na.rm=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution can then be compared with that of the non-weighted earnings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "as.data.frame(quantile(earnings_2015$salary))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red><h3> Checkpoint 1: Find the distribution of earnings for a different cohort</h3></font>\n",
    "\n",
    "Using the `svydesign` function above, find the distribution of earnings for the 2012 cohort."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the relevant variables from the SDR data to find the earnings \n",
    "# among the 2012 cohort\n",
    "\n",
    "query <- \"\n",
    "SELECT salary, wtsurvy\n",
    "FROM ds_nsf_ncses.dbo.nsf_sdr_2017\n",
    "WHERE sdryr = '2012' \n",
    "\"\n",
    "\n",
    "earnings_2012 <- dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the first rows of the table\n",
    "head(earnings_2012)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove the rows with the logical skip value\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the weighted estimates\n",
    "# Change the names of the variables\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the percentiles\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, it is important to keep in mind that when using survey data (due to the specificities with which every particular survey is designed), the weights always need to be applied in order to be able to draw accurate conclusions about the general population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Missing Data\n",
    "\n",
    "Datasets sometimes contain variables with missing (or unknown) data. There are many options for approaching these missing values, such as ignoring them or leveraging various methods to fill those in, in order to be able to use the data. This example focuses on the `age_at_dissertation` variable in the `nsf_sed` table, where there are a lot of missing ages.\n",
    "\n",
    "Keep in mind that these imputed values will be **approximations**, and must be treated as such. If choosing to impute missing values in the project or future work, acknowledge the process and clearly state for which variables the values are imputed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, see how many missing ages there are\n",
    "# Get the count of non-null age at dissertation and count of all rows\n",
    "\n",
    "query <- \"\n",
    "select count(age_at_dissertation) as age_at_diss_count, count(*) as total_count \n",
    "from ds_nsf_ncses.dbo.nsf_sed;\n",
    "\"\n",
    "dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that there is a sizable amount of observations with missing ages, imputation might be most suitable. The following subsections will walk through different methods of imputation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 1. Dropping missing values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oftentimes, missing values are treated by ignoring them. However, ignoring potentially non-missing values assumes that they represent the same distribution as the present one. This method should **never, ever, ever** be used in practice.\n",
    "\n",
    "> Deleting missing values is often called listwise deletion and essentially assumes that missing values are missing completely at random (MCAR)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's get the SED table\n",
    "query <- \"\n",
    "SELECT drf_id, age_at_dissertation, phdfield_name, srceprim\n",
    "FROM ds_nsf_ncses.dbo.nsf_sed\n",
    "\"\n",
    "sed <- dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the overall mean \n",
    "overall_mean = mean(sed$age_at_dissertation, na.rm=TRUE)\n",
    "\n",
    "overall_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 2. Mean Imputation\n",
    "\n",
    "One of the simplest ways of imputing values is by taking the mean and filling it in for the missing values. It's possible to do this by using the overall mean, as well as by certain subgroups. Given that the overall mean when ignoring missing values was found above and stored in `overall_mean`, this imputation technique can be completed in a single code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set missing means to the overall mean\n",
    "complete_sed <- sed %>%\n",
    "    mutate(age_at_dissertation = ifelse(is.na(age_at_dissertation), overall_mean, age_at_dissertation))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(complete_sed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead, as mentioned at the beginning of the subsection, mean imputation can also be used on more granular subgroups, such as within each `phdfield_name`. This assumes that the missing `age_at_dissertation` values are conditional on the `phdfield_name`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate non-missing ages within each phdfield_name\n",
    "mean_by_field <- sed %>%\n",
    "                group_by(phdfield_name) %>%\n",
    "                summarise(mean_by_field = mean(age_at_dissertation, na.rm=TRUE))\n",
    "mean_by_field"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After calculating the mean non-missing ages within each PhD field, these means can be joined to the original `sed` data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# join mean_by_field with sed by phdfield_name\n",
    "sed_phd <- sed %>%\n",
    "    inner_join(mean_by_field, by = 'phdfield_name')\n",
    "\n",
    "head(sed_phd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The imputed ages can now be calculated where missing `age_at_dissertation` values will take on the appropriate `mean_by_field` calculation and non-missing ones will remain the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create imputed age column\n",
    "sed_phd <- sed_phd %>%\n",
    "    mutate(\n",
    "        imp_age = ifelse(is.na(age_at_dissertation), mean_by_field, age_at_dissertation)\n",
    "    )\n",
    "\n",
    "head(sed_phd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see means after imputation\n",
    "mean(sed_phd$imp_age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method 3. Regression Imputation\n",
    "\n",
    "Regression can also be used to try to compute more accurate values. A regression equation will be built with the observations for which the age is known, and then it will be applied to essentially predict the missing values. This is, in effect, an extension of the mean imputation by subgroup. Here, the primary source of support as well as the PhD field will be used as predictors to generate the regression.\n",
    "\n",
    "> Note: Here the assumptions associated with linear regressions are not checked, as this example is aimed at merely displaying how to use a linear regression for imputation. If planning on using regression imputation, please check all assumptions before employing a predictive model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the needed variables\n",
    "query <- \"\n",
    "SELECT drf_id, age_at_dissertation, phdfield_name, srceprim\n",
    "FROM ds_nsf_ncses.dbo.nsf_sed\n",
    "\"\n",
    "sed <- dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the model will be built based on information only from the members of the cohort with non-missing age values, the `sed` dataset can be split into two datasets - one for training (`sed_nona`, without any missing data for the training of the model) and one for testing (`sed_na`, with only missing data to use for prediction)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataset with no missing data, for the training of the model\n",
    "sed_nona <- sed %>%\n",
    "            drop_na()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataset with only missing data, to use for prediction\n",
    "# Remove a column with age as it is always null\n",
    "sed_na <- sed %>%\n",
    "    filter(is.na(age_at_dissertation)) %>%\n",
    "    select(-c(age_at_dissertation))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model creation process for a linear regression can be done using the `lm()` function. The variable to predict is on the left-hand side of `lm()` before the `~`, and the predictors (`phdfield_name` and `srceprim` in this case) are on the right-hand side of the `~`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model and fit coefficients\n",
    "model <- lm(age_at_dissertation ~ phdfield_name + srceprim, data = sed_nona)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the model has been fit, the missing ages can be predicted using the coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict age for test set\n",
    "pred_age <- data.frame(age_at_dissertation = predict(model, newdata=sed_na))\n",
    "\n",
    "head(pred_age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the output for `predict()` retains the same order of rows from `sed`, the `age_at_dissertation` variable from `pred_age` can be easily added into the existing `sed` data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see updated data frame with predicted age\n",
    "cbind(sed_na, pred_age) %>% \n",
    "    head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save updated data frame\n",
    "sed_na_w_age <- cbind(sed_na, pred_age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, before seeing the effects of the imputation method, the training set should be combined to get the full set of `age_at_dissertation` values, both imputed and previously known ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save combined training and testing sets\n",
    "sed_all <- rbind(sed_na_w_age, sed_nona)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entire age distribution can now be viewed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see age distribution for full cohort\n",
    "summary(sed_all$age_at_dissertation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see age distribution for imputed portion of cohort\n",
    "summary(sed_na_w_age$age_at_dissertation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red><h3> Checkpoint: Use another variable for a regression</h3></font> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think of another variable that can be used in the regression to predict the age at dissertation. Include it in the model, and try running it again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imputation for Categorical Missing Values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Impute using mode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For categorical variables, it is impossible to apply something like mean imputation, because there is not a mean to calculate. A method for imputing missing categorical values is to find the most frequent value (mode), and impute using the mode.\n",
    "\n",
    "This example will cover using mode imputation for missing `marital` status based on each `age_at_dissertation` value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's get the marital status and age\n",
    "query <- \"\n",
    "SELECT marital, age_at_dissertation, drf_id\n",
    "FROM ds_nsf_ncses.dbo.nsf_sed\n",
    "\"\n",
    "marital_status <- dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, as exhibited in previous methods, the missing values must be removed. Once they are removed, the most frequent marital status per age can be found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only non-missing marital statuses\n",
    "no_miss_marital <- marital_status %>% \n",
    "    filter(marital != 'NA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `group_by` function, in combination with `top_n`, can be used together to isolate the most common `marital` value within each `age_at_dissertation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get only the most frequent marital status values per age\n",
    "marital_mode <- no_miss_marital %>%\n",
    "                 group_by(age_at_dissertation) %>%\n",
    "                count(marital) %>%\n",
    "                top_n(1) %>%\n",
    "                select(-c(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(marital_mode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Marital_mode` can now function as a lookup table with the most frequent value of a marital status per age. This lookup table can now be joined with the original `marital_status` table. This will create two columns: `marital.x` (original marital status value) and `marital.y` (marital status value from the lookup table)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge original and lookup tables\n",
    "merged <- inner_join(marital_status, marital_mode, by=c('age_at_dissertation'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(merged)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the NA values in the `marital.x` column can be replaced with the known values in the column `marital.y`, and then these two columns (`marital.y` and `marital.x`) can be removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replace the NA values in the `marital.x` column with the known values in the column `marital.y`\n",
    "merged <- merged %>%\n",
    "    mutate(\n",
    "        marital_imp = ifelse(marital.x == 'NA', marital.y, marital.x)\n",
    "    ) %>%\n",
    "    select(-c(marital.x, marital.y))\n",
    "\n",
    "head(merged)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Optional) Advanced: Using machine learning to impute values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To impute values, the machine learning algorithm can also be used called the `K-nearest Neighbors`. The principle behind it 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 there were cases where the data on a marital status was completely missing per age, their marital status could be approximated by referring to other characteristics which could be shared by that age 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 age group)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Optional) Replicate Weights\n",
    "\n",
    "The SDR also comes with 104 replicate weights. The following code cells show how replicate weights can be used to estimate the median salary for those in 2015 SED using the results from the 2015 SDR data.\n",
    "\n",
    "First, the main SDR table must be joined to the replicate weights table based on matching `refid` values.\n",
    "\n",
    "> There are non-optional sections below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the relevant variables from the SDR data to find \n",
    "# the earnings among the 2015 cohort\n",
    "\n",
    "query <- \"\n",
    "SELECT salary, wtsurvy, rw.*\n",
    "FROM ds_nsf_ncses.dbo.nsf_sdr_2017 sdr \n",
    "JOIN ds_nsf_ncses.dbo.nsf_sdr_rw_2017 rw\n",
    "ON sdr.refid = rw.refid\n",
    "WHERE sdryr = '2015'\n",
    "\"\n",
    "earnings_rw_2015 <- dbGetQuery(con,query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the head of the table\n",
    "\n",
    "head(earnings_rw_2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The replicate weight values are all columns that start with \"rw\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separate out replicate weights\n",
    "rw_2015 <- earnings_rw_2015 %>%\n",
    "    select(starts_with(\"rw\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(rw_2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the sample weights have been separated into their own data frame, the median salary can be calculated using the same method as above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate salary estimates\n",
    "rw_earnings_2015 <- svydesign(ids=~1, data=earnings_rw_2015, weights=earnings_rw_2015$wtsurvy)\n",
    "\n",
    "# extract median salary estimate\n",
    "earnings_median <- as.data.frame(svyquantile(~salary, rw_earnings_2015, .50, na.rm=TRUE))$'0.5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earnings_median"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data documentation, as detailed in the SDR 2017 Replicate Weight User Guide, describes how to find the variance of this estimate. This formula is given by\n",
    "\n",
    "$$ v_{rep}(\\hat{\\theta}) = \\sum^R_{r=1} 0.038461 * (\\hat{\\theta}_r - \\hat{\\theta})^2. $$\n",
    "\n",
    "> The constant term (0.038461) is taken from the Replicate Weight User Guide.\n",
    "\n",
    "After initializing a vector to store $\\hat{\\theta}_r$ values, these values can be found using each set of replicate weights in order to find the estimate 104 times. Do this with a `for` loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize a vector with NaN first\n",
    "thetas <- rep(NA, 104)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thetas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change weights to numeric\n",
    "\n",
    "rw_2015 <- mutate_all(rw_2015, function(x) as.numeric(as.character(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loop through and calculate\n",
    "for(i in 1:104){\n",
    "    weighted <- svydesign(ids=~1, data=earnings_rw_2015, weights=rw_2015[,i])\n",
    "    thetas[i] <- as.data.frame(svyquantile(~salary, weighted, .50))$'0.5'\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see thetas\n",
    "thetas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, use array operations to calculate the variance. That is, for example, if subtracting a scalar from the array, it will do the operation on each element of the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate variance\n",
    "med_var = sum(0.038461 * (thetas - earnings_median) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see variance\n",
    "med_var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard error\n",
    "sqrt(med_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Close the database connection\n",
    "dbDisconnect(con)"
   ]
  }
 ],
 "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"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}