5B.Machine_Learning_Prediction.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, Benjamin Feder, Nathan Barrett </center>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6412954\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6412954.svg\" alt=\"DOI\"></a>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Machine Learning Part 2: Prediction Models and Model Evaluation**\n",
    "\n",
    "Now that the features and label have been created for supervised machine learning, we can deploy and evaluate various algorithms in trying to predict future employment for a cohort of graduates. Recall that we will build our model on 2017 bachelor's degree recipients in Texas and then test the model on the 2018 cohort.\n",
    "\n",
    "**If you have not done so, please review the [features creation notebook](./5A.ML_Feature_Creation.ipynb) prior to running any code in this notebook.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **1. Read in the Data**\n",
    "\n",
    "As always, let's load R packages and establish the database connection first. Note that in this notebook, we include additional packages `caret`, `rpart`, `rpart.plot`, and `randomForest` to get functions for calculating evalution metrics and fitting models."
   ]
  },
  {
   "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",
    "# For faster date conversions\n",
    "library(lubridate)\n",
    "\n",
    "# Classification and regression training package. For streamlining the process for creating predictive models\n",
    "library(caret)\n",
    "\n",
    "# For Decision Tree model and to plot the tree\n",
    "library(rpart)\n",
    "library(rpart.plot)\n",
    "\n",
    "# For Random Forest model\n",
    "library(randomForest)\n",
    "\n",
    "library(scales)"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the training set that we created in the first ML notebook\n",
    "# Adding in major name instead of CIP code\n",
    "query <- \"\n",
    "SELECT nts.gradid, nts.gradtypi, SUBSTRING(cl.CIPTitle2020, 1, 10) gradmaj, nts.instregion, nts.transfer_ind, nts.total_sems, nts.total_hours, nts.employed_count, nts.employed_prop, nts.wage_ind label\n",
    "FROM tr_tx_2021.dbo.nb_training_set nts \n",
    "LEFT JOIN ds_public_1.dbo.cip_lookup cl\n",
    "ON nts.gradmaj = cl.CIPCode2020\n",
    "\" \n",
    "\n",
    "training_set <- dbGetQuery(con, query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(training_set)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the testing set that we created in the first ML notebook\n",
    "query <- \"\n",
    "SELECT nts.gradid, nts.gradtypi, SUBSTRING(cl.CIPTitle2020, 1, 10) gradmaj, nts.instregion, nts.transfer_ind, nts.total_sems, nts.total_hours, nts.employed_count, nts.employed_prop, nts.wage_ind label\n",
    "FROM tr_tx_2021.dbo.nb_testing_set nts \n",
    "LEFT JOIN ds_public_1.dbo.cip_lookup cl\n",
    "ON nts.gradmaj = cl.CIPCode2020\n",
    "\" \n",
    "\n",
    "testing_set <- dbGetQuery(con, query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(testing_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **2. Prepare the Data**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to clustering, a dataset must have certain properties for running a prediction model. We will clean our data frames by following 5 steps:\n",
    "\n",
    "1. Remove overlap\n",
    "1. Remove non-explanatory features\n",
    "1. Set variable types\n",
    "1. Analyze missingness\n",
    "1. Examine scales across variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Remove Overlap**\n",
    "\n",
    "When running and assessing prediction algorithms, there should not be any units (individuals, in this case) present in both the training and testing sets, as we want to provide the algorithm with unseen test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see amount of individuals before removal\n",
    "testing_set %>%\n",
    "    summarize(\n",
    "        n_distinct(gradid)\n",
    "    )\n",
    "\n",
    "# remove all gradid values from testing_set that are in training_set\n",
    "testing_set <- testing_set %>% \n",
    "    anti_join(training_set, by = 'gradid')\n",
    "\n",
    "# see amount of individuals after removal\n",
    "testing_set %>%\n",
    "    summarize(\n",
    "        n_distinct(gradid)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Remove non-explanatory features**\n",
    "\n",
    "Like we did with clustering, we will want to remove any non-explanatory features before running a supervised machine learning model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove gradid\n",
    "training_set <- training_set %>% select(-gradid)\n",
    "\n",
    "testing_set <- testing_set %>% select(-gradid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Set Variable Types**\n",
    "\n",
    "Take a look at the type of each column. Note that the categorical features are still in character type, such as `instregion` and `gradmaj`, or integer type, such as `gradtypi` and `transfer_ind`. We need to convert them to factors so that they can be used as dummy variables in ML models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# look at structure of data frame\n",
    "str(training_set)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change all character variables to factors\n",
    "training_set <- training_set %>% mutate_if(is.character, as.factor)\n",
    "testing_set <- testing_set %>% mutate_if(is.character, as.factor) \n",
    "\n",
    "# convert some numeric variables into factor type. We will leave the label as numeric\n",
    "training_set$gradtypi <- factor(training_set$gradtypi)\n",
    "training_set$transfer_ind <- factor(training_set$transfer_ind)\n",
    "\n",
    "testing_set$gradtypi <- factor(testing_set$gradtypi)\n",
    "testing_set$transfer_ind <- factor(testing_set$transfer_ind)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save data frames as df_training and df_testing b/c will need original data frames later\n",
    "df_training <- training_set\n",
    "df_testing <- testing_set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Analyze Missingness**\n",
    "\n",
    "Before we run a machine learning model, we want to make sure that there are no missing values in our data so that all rows are preserved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check if columns have missing values\n",
    "sapply(df_training, function(x) sum(is.na(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Examine Scales Across Variables**\n",
    "\n",
    "Just like we did prior to clustering, we will want to confirm that all of the numerical features are on similar scales to prevent inappropriate overweighting of certain features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get descriptions of each variable using \"summary\" function\n",
    "summary(df_training)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scale the numeric variables for both sets\n",
    "df_training <- df_training %>%\n",
    "    mutate(\n",
    "        total_sems = scale(total_sems)[,1],\n",
    "        total_hours = scale(total_hours)[,1],\n",
    "        employed_count = scale(employed_count)[,1],\n",
    "        employed_prop = scale(employed_prop)[,1]\n",
    "    )\n",
    "\n",
    "df_testing <- df_testing %>%\n",
    "    mutate(\n",
    "        total_sems = scale(total_sems)[,1],\n",
    "        total_hours = scale(total_hours)[,1],\n",
    "        employed_count = scale(employed_count)[,1],\n",
    "        employed_prop = scale(employed_prop)[,1]\n",
    "    )\n",
    "\n",
    "# see evidence of scaling\n",
    "head(df_testing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **3. Create Functions**\n",
    "\n",
    "Some code will be used many times in this notebook, such as the code to calculate precision and recall, the code to get precision at K, and the code to compare a model's results with baselines. We can put these code in functions so that we can easily repeat these calculation processes by using one line of code.\n",
    "\n",
    "The description of these functions is available below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. `precision_recall(test_data, label, pscore)`: calculates precision and recall \n",
    "\n",
    "    After we run an ML model and use its results to predict outcomes for the testing data, we can use the `precision_recall()` function to calculate precision and recall. **This function returns a data frame which contains precision and recall at various k (0<k<1).** It has three arguments:\n",
    "    - `test_data`: the data frame of the testing data\n",
    "    - `label`: the name of the label column. It should be a **string**.\n",
    "    - `pscore`: the name of the predicted score column. It should be a **string**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. `precision_at_k(k, test_data, label, pscore)`: calculate precision at K. \n",
    "\n",
    "    We can use the `precision_at_k()` function to get the precision at the k% we choose. **This function returns a number.** It has four arguments:\n",
    "    - `k`: the percent of population we have enough resources to intervene (e.g., help graduates get employed)\n",
    "    - `test_data`: the data frame of the testing data\n",
    "    - `label`: the name of the label column. It should be a **string**.\n",
    "    - `pscore`: the name of the predicted score column. It should be a **string**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. `compare_w_baseline(k, model, test_data, label, pscore)`: compare a model's precision at k with baselines\n",
    "\n",
    "    This function compares our ML model with two baseline models. In the first baseline model, we randomly assign labels (0 or 1) to each person. In the second baseline model, we assign 1 to every person. **This function returns a data frame which includes precision at k (we choose the k) of our ML model and the two baseline models.** It has five arguments: \n",
    "    - `k`: the percent of population we have enough resources to intervene (e.g., help graduates get employed)\n",
    "    - `model`: the name of the ML model. It should be a **string**\n",
    "    - `test_data`: the data frame of the testing data\n",
    "    - `label`: the name of the label column. It should be a **string**.\n",
    "    - `pscore`: the name of the predicted score column. It should be a **string**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we just need to run the next three blocks of code so that they will be ready for us to use in the model evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a function, which returns a DataFrame containing precision and recall at K% of the population\n",
    "precision_recall <- function(test_data, label, pscore) {\n",
    "    \n",
    "    # Get the actual label and predicted score in the testing data\n",
    "    df_temp <- df_testing[, c(label, pscore)] \n",
    "    \n",
    "    #Calculate Precision and Recall at K\n",
    "    df_temp <- df_temp %>%\n",
    "        arrange(desc(!!sym(pscore))) %>% # Sort the rows descendingly based on the predicted score\n",
    "        mutate(rank = row_number()) %>% # Add rank to each row\n",
    "        mutate(recall = cumsum(label == 1)/sum(label == 1),  # Calculate Recall at K\n",
    "               precision = cumsum(label == 1)/rank, # Calculate precision at k\n",
    "               k = rank/(nrow(test_data))) # Percent of population\n",
    "    \n",
    "    return(df_temp)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a function, which returns precision at K\n",
    "precision_at_k <- function(k=0.1, test_data, label, pscore) {\n",
    "    \n",
    "    # Get the Precision-Recall at K DataFrame\n",
    "    df <- precision_recall(test_data, label, pscore)\n",
    "    \n",
    "    # Assign a few parameters\n",
    "    pct_pop <- k  # Percent of population the resource can cover\n",
    "    test_pop <- nrow(df_testing) # Total number of people in the testing data\n",
    "    pop_at_k <- as.integer(pct_pop * test_pop) # At K percent of the population, how many people the recourse can cover\n",
    "    \n",
    "    # Get precision at K% from the Precision-Recall DataFrame\n",
    "    prec_at_k <- df$precision[pop_at_k]\n",
    "    \n",
    "    return(prec_at_k)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a function, which returns a DataFrame containing measures for baseline models \n",
    "# Baseline model 1: Randomly assign the label (0 or 1)\n",
    "# Baseline model 2: Guess everyone is employed\n",
    "\n",
    "compare_w_baseline <- function(k=0.1, model, test_data, label, pscore) {\n",
    "    # Set a seed so we get consistent results\n",
    "    set.seed(42)\n",
    "    \n",
    "    # Assign a few parameters\n",
    "    pct_pop <- k  # Percent of population the resource can cover\n",
    "    test_pop <- nrow(df_testing) # Total number of people in the testing data\n",
    "    pop_at_k <- as.integer(pct_pop * test_pop) # At K percent of the population, how many people the recourse can cover\n",
    "    \n",
    "    # Get the Precision-Recall at K DataFrame\n",
    "    df <- precision_recall(test_data, label, pscore)\n",
    "    \n",
    "    # Generate Precision-Recall at K for baseline model 1\n",
    "    df_random <- df %>%\n",
    "        mutate(random_score = runif(nrow(df))) %>% # Generate a row of random scores\n",
    "        arrange(desc(random_score)) %>% # Sort the data by the random scores\n",
    "        mutate(random_rank = row_number()) %>% # Add rank to each row\n",
    "        mutate(random_recall = cumsum(label==1)/sum(label==1), # Calculate Recall at K\n",
    "           random_precision = cumsum(label==1)/random_rank, # Calculate Precision at K\n",
    "           random_k = random_rank/(nrow(df)))\n",
    "    \n",
    "    # Precision at K of the model\n",
    "    model_precision_at_k <- precision_at_k(k, test_data, label, pscore)\n",
    "    \n",
    "    # Precision at K of baseline model 1\n",
    "    random_precision_at_k <- df_random$random_precision[pop_at_k]\n",
    "    \n",
    "    # Precision at K of baseline model 2\n",
    "    allemp_precision_at_k <- sum(test_data$label)/test_pop\n",
    "    \n",
    "    # Create a DataFrame which shows all measures\n",
    "    df_compare_prec <- data.frame(\"model\" = c(model, \"Random\", \"All Employed\"),\n",
    "                              \"precision\" = c(model_precision_at_k, random_precision_at_k, allemp_precision_at_k))\n",
    "    \n",
    "    return(df_compare_prec)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **4. Logistic Regression**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the logit regression with the training dataset\n",
    "# label is outcome variable and we will predict on all other features\n",
    "# family = binomial (link = 'logit') specifies functional form and error distribution\n",
    "lr_model <- glm(label ~ ., family = binomial(link = 'logit'), data = df_training)\n",
    "\n",
    "# Show feature importance\n",
    "summary(lr_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above results, the column **Estimate** shows the importance of each feature in predicting whether a graduate will be employed in the 4th quarter after graduation. The stars at the end of each row indicates whether a feature coefficient is statistically significant. Note that our model automatically leaves out one category for each group of dummy variables. For example, of all the `gradmaj` dummies, agriculture is the omitted category. \n",
    "\n",
    "We can see that based on our logistic regression model, most of our features are important predictors. Additionally, we can interpret the estimate for factors, as a comparison to the feature left out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Model Evaluation**\n",
    "\n",
    "Next, let's evaluate the logistic regression's performace on the testing set. There are a variety of evaluation metrics we can use, and we will introduce the following:\n",
    "- Confusion matrix\n",
    "- Precision at k\n",
    "- Comparison to baseline models\n",
    "- Precision-recall curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Confusion Matrix**\n",
    "\n",
    "First, let's use the function `predict()` to predict outcomes of graduates in our testing data based on the results of our logistic regression. Recall that the predicted scores don't tell us the predicted outcome of a graduate. We can say that a graduate with score 0.9 is predicted to be less likely to be employed in the 4th quarter after graduation than a graduate with score 0.7. But whether 0.9 implies the graduate will be employed or not depends on our choice of the threshold. For example, if we choose a threshold of 0.5, any graduate with predicted scores greater than 0.5 will be defined as employed and any graduate with predicted scores less than 0.5 will be defined as unemployed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict on testing set\n",
    "# type = \"response\" ensures we get predicted probabilities\n",
    "df_testing$predict_score <- predict(lr_model, df_testing, type = \"response\")\n",
    "\n",
    "# Set a threshold for the predicted score\n",
    "# Assume people who get more than 0.5 predicted score will be employed 4th quarter after graduation\n",
    "threshold <- .5\n",
    "\n",
    "# Add the predicted outcome to a new column in df_testing\n",
    "# If predicted score is greater than the threshold, then predict label = 1, otherwise, predict label =0.\n",
    "df_testing$predict_label <- ifelse(df_testing$predict_score > threshold, 1, 0) \n",
    "\n",
    "# Confusion matrix\n",
    "confusionMatrix(factor(df_testing$predict_label), factor(df_testing$label))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the results of `confusionMatrix()`, you can view the number of true positives, false positives, true negatives, and false negatives predicted with the model.  Additionally, you can see the specificity, sensitivity, and accuracy of the models. We will discuss some of these metrics in the following sections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Precision at K**\n",
    "\n",
    "Most of the time, we may not care about how accurate our model is in predicting both positive (employed) and negative (not employed) outcomes. Instead, we may want to check how accurate our model is in predicting positive outcomes, which is captured by **precision**.\n",
    "\n",
    "$$Precision = \\frac{True Positive}{True Positive + False Positive}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of arbitrarily controlling the threshold like we did before, we may want to think about that given our resources, such as funding, time, and staff, what percentage of the population we can help. In the context of this project, we need to decide what percentage of graduates we can provide assistance to so that they have a better opportunity to be employed after graduation. Suppose that we have enough resources to cover 10% of the graduates. We can use the function we created in section 2, `precision_at_k()`, to get the precision at 10% of the population. \n",
    "\n",
    "> Note: In the case of intervention, it may be more useful to assign `label = 1` as the negative scenario, thus isolating the most likely to not find employment, in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check precision at K%\n",
    "k_pct <- .1 # Here, we are checking precision at 10% of the population, change the value to check precision at different K\n",
    "\n",
    "lr_prec_at_10 <- precision_at_k(k_pct, df_testing, \"label\", \"predict_score\")\n",
    "\n",
    "print(paste0(\"In the Logistic Regression Model, precision at \", label_percent()(k_pct), \n",
    "             \" of the population is: \", (lr_prec_at_10)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result implies that at 10% of the population, among all the graduates the logistic regression model predicts to be employed, REDACTED were actually employed in the 4th quarter after graduation. <font color=red> You can change the value assigned to `k_pct` in the first line of the code to explore precision at other percent of the population.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Compare our model with baselines**\n",
    "\n",
    "How good is the logistic regression's precision at 10% of the population? Is it accurate enough or not? Recall the function `compare_w_baseline()`, which allows us to compare our model with two baselines. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare precision at K% of the population with baseline models\n",
    "df_compare <- compare_w_baseline(k=0.1, \"Logistic Regression\", df_testing, \"label\", \"predict_score\")\n",
    "\n",
    "df_compare"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the logistic regression model's precision is higher than the two baseline models' precision. This implies that at 10% of the population, logistic regression model is better at predicting who will be employed compared to randomly guessing who will be employed or by guessing that everyone will be employed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a bar plot to show the comparision of the Logistic Regression model and the baselines (random guess or guess everyone stay)\n",
    "\n",
    "# For easier reading, increase base font size\n",
    "theme_set(theme_gray(base_size = 16))\n",
    "# Adjust repr.plot.width and repr.plot.height to change the size of graphs\n",
    "options(repr.plot.width = 10, repr.plot.height = 5)\n",
    "\n",
    "# Specify source dataset and x and y variables\n",
    "lr_baseline_plot <- ggplot(df_compare, aes(x = model, y = precision)) + \n",
    "    geom_col() + # Plots bars on the graph\n",
    "    geom_text(size = 5, aes(label = format(precision, digit = 3), vjust = -0.5)) + # Show values on top of the bar\n",
    "    scale_y_continuous(breaks = seq(0, 1, 0.2), limits = c(0, 1)) + # Adjust the y scale to set the interval for tick marks\n",
    "    labs(title = \"Precision at 10% against the baseline, Logistic Regression\", # Add graph title\n",
    "         x = \" \", y = 'Precision at 10%') + \n",
    "    theme(axis.title.x = element_text(face=\"bold\"), # Adjust the style of X-axis label\n",
    "          axis.title.y = element_text(face=\"bold\"), # Adjust the styles of the two Y-axes labels\n",
    "          axis.text.x = element_text(face=\"bold\", size = 16),\n",
    "          plot.title = element_text(hjust = 0.5))  # Center the graph title\n",
    "\n",
    "print(lr_baseline_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Precision-Recall Curve**\n",
    "\n",
    "Another measure we often use to evaluate the performance of a ML model is recall. It shows us what percentage of people with actual positive outcomes our model can capture. In the context of this project, recall tells us what percentage of those who found employment our ML model can accurately predict. \n",
    "\n",
    "Again, we do not need to set a specific threshold, as we can plot their values on a line chart so that it is easier for us to see how precision and recall change with our choice of k."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Precision = \\frac{True Positive}{True Positive + False Positive}$$\n",
    "\n",
    "$$Recall = \\frac{True Positive}{True Positive + False Negative}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get Precision-Recall at K data frame\n",
    "df_measure_at_k <- precision_recall(df_testing, \"label\", \"predict_score\")\n",
    "\n",
    "# See the top records of the DataFrame\n",
    "head(df_measure_at_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For easier reading, increase base font size\n",
    "theme_set(theme_gray(base_size = 16))\n",
    "# Adjust repr.plot.width and repr.plot.height to change the size of graphs\n",
    "options(repr.plot.width = 10, repr.plot.height = 5)\n",
    "\n",
    "# Create the precision-recall curve\n",
    "lr_pr_curve <- ggplot(df_measure_at_k, aes(x=k)) + # Plot percent of population (k) on the x-axis\n",
    "    geom_line(aes(y=precision), color = 'blue') + # Add the precision curve\n",
    "    geom_line(aes(y=recall), color = 'red') + # Add the recall curve\n",
    "    scale_y_continuous(       # We need to create a dual-axis graph, so we need to define two y axes\n",
    "        name = \"Precision\",     # Label of the first axis\n",
    "        sec.axis = sec_axis(~.*1,name=\"Recall\"), # Add a second axis and specify its label\n",
    "        breaks = seq(0, 1, 0.1)) +  # Adjust the tick mark on Y-axis\n",
    "    scale_x_continuous(breaks = seq(0, 1, 0.2)) + # Adjust the tick mark on X-axis\n",
    "    labs(title = \"Precision-Recall Curve, Logistic Regression\", # Add graph title\n",
    "         x = \"Percent of Population\") + # Add X-axis label\n",
    "    theme(axis.title.x = element_text(face=\"bold\"), # Adjust the style of X-axis label\n",
    "          axis.title.y.left = element_text(face=\"bold\", color=\"blue\"), # Adjust the styles of the two Y-axes labels\n",
    "          axis.title.y.right = element_text(face='bold', color = 'red'),\n",
    "          plot.title = element_text(hjust = 0.5))  # Center the graph title\n",
    "\n",
    "# Display the graph that we just created\n",
    "print(lr_pr_curve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above graph, the blue line represents precision. We can see that it stays relatively constant when we increase the choice of k. The red line shows recall. Its values increases as we increase the choice of k. At 10% of the population, precision is high and recall is low.\n",
    "\n",
    "We would expect to see this because the precision curve starts at a high point, because when k is low, only graduates with the highest predicted scores will be defined as employed. When k increases (meaning we are selecting a higher percent of the population), we relax the threshold and predict graduates with relatively low predicted scores to be not employed. Therefore, the precision decreases, which is what we see. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **5. Decision Tree**\n",
    "\n",
    "There are other supervised machine learning algorithms besides logistic regressions. Here, we will try a decision tree model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the decision tree model\n",
    "# method = 'class' is to designate the classification tree method (label is binary)\n",
    "dt_model <- rpart(label ~ ., method = 'class', data = df_training)\n",
    "\n",
    "# Print results\n",
    "printcp(dt_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above result, we can see that the only variable the decision tree model actually used in prediction was `total_sems`. In the logistic regression model, we saw that majority of the variables were statistically significant. \n",
    "\n",
    "We can also view the tree graphically. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the tree\n",
    "prp(dt_model, # Your decision tree model\n",
    "    type = 0, # type of trees\n",
    "    extra = 100, # what information to show in each node\n",
    "    main = \"Decision Tree Model\") # Add a title to your decision tree graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can evaluate the decision tree model in a similar fashion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the Decision Tree model predicted score and save it in a column in the testing DataFrame\n",
    "df_testing$dt_predict_score <- predict(dt_model, df_testing, type = 'prob')[,2]\n",
    "\n",
    "# Get the Decision Tree model Precision-Recall at K% DataFrame\n",
    "df_dt_measure_at_k <- precision_recall(df_testing, \"label\", \"dt_predict_score\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For easier reading, increase base font size\n",
    "theme_set(theme_gray(base_size = 16))\n",
    "# Adjust repr.plot.width and repr.plot.height to change the size of graphs\n",
    "options(repr.plot.width = 10, repr.plot.height = 5)\n",
    "\n",
    "# Create the precision-recall curve\n",
    "dt_pr_curve <- ggplot(df_dt_measure_at_k, aes(x=k)) + \n",
    "    geom_line(aes(y=precision), color = 'blue') + # Add the precision curve\n",
    "    geom_line(aes(y=recall), color = 'red') + # Add the recall curve\n",
    "    scale_y_continuous(       # We need to create a dual-axis graph, so we need to define two y axes\n",
    "        name = \"Precision\",     # Label of the first axis\n",
    "        sec.axis = sec_axis(~.*1,name=\"Recall\"), # Add a second axis and specify its label\n",
    "        breaks = seq(0, 1, 0.1)) +  # Adjust the tick mark on Y-axis\n",
    "    scale_x_continuous(breaks = seq(0, 1, 0.2)) + # Adjust the tick mark on X-axis\n",
    "    labs(title = \"Precision-Recall Curve, Decision Tree\", # Add graph title\n",
    "         x = \"Percent of Population\") + # Add X-axis label\n",
    "    theme(axis.title.x = element_text(face=\"bold\"), # Adjust the style of X-axis label\n",
    "          axis.title.y.left = element_text(face=\"bold\", color=\"blue\"), # Adjust the styles of the two Y-axes labels\n",
    "          axis.title.y.right = element_text(face = 'bold', color = 'red'),\n",
    "          plot.title = element_text(hjust = 0.5))  # Center the graph title\n",
    "\n",
    "# Display the graph that we just created\n",
    "print(dt_pr_curve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the precision-recall curve for the decision tree model looks relatively similar to the curve for the logistic regression model. However, the precision at low values of *k* is lower than that of the logistic regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **6. Compare Multiple Models**\n",
    "\n",
    "In this section, we provide you code to compare evaluation metrics for several ML models all at once. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \"refresh\" the training data and the testing data to remove predicted scores\n",
    "\n",
    "df_training <- training_set\n",
    "df_testing <- testing_set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to a logistic regression and decision tree, we will also try inputting our data into a random forest. This list can be expanded to include other ML models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a list to include all the models we want to compare\n",
    "# LR: Logistic Regresion\n",
    "# DT: Decision Tree\n",
    "# RF: Random Forest\n",
    "model_list <- c(\"LR\", \"DT\", \"RF\")\n",
    "\n",
    "# Define percent of population the resource can cover\n",
    "k <- 0.1\n",
    "pct_pop <- k  # Percent of population the resource can cover\n",
    "test_pop <- nrow(df_testing) # Total number of people in the testing data\n",
    "pop_at_k <- as.integer(pct_pop * test_pop) # At K percent of the population, how many people the resource can cover\n",
    "\n",
    "# make data frame number of rows to be number of models\n",
    "n <- length(model_list)\n",
    "df_compare_models <- data.frame(\n",
    "    Model = model_list, \n",
    "    accuracy = double(n),\n",
    "    precision_at_k = double(n),\n",
    "    recall_at_k = double(n)\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then loop through the models to run them all in the same series of commands.\n",
    "\n",
    "> Note: The code cell below will take approximately 8 minutes to execute. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (model in model_list) {\n",
    "    \n",
    "    # Logististic Regression Model\n",
    "    if (model==\"LR\") {\n",
    "        fit <- glm(label ~ ., family = binomial(link = 'logit'), data = df_training) # Fit the model\n",
    "        df_testing$predict_score <- predict(fit, df_testing, type = 'response') # Predict scores\n",
    "    }\n",
    "    \n",
    "    # Decision Tree Model\n",
    "    if (model==\"DT\") {\n",
    "        fit <- rpart(label ~ ., method = 'class', data = df_training) # Fit the model\n",
    "        df_testing$predict_score <- predict(fit, df_testing, type = 'prob')[,2] # Predict scores\n",
    "    }\n",
    "    \n",
    "    # Random Forest Model\n",
    "    if (model == \"RF\"){\n",
    "        df_training$label <- factor(df_training$label) # for the algorithm to run we need to convert the label to a factor.\n",
    "        df_testing$label <- factor(df_testing$label)\n",
    "        # ntree determines number of trees\n",
    "        # mtry is number of variables randomly sampled as candidates at each split\n",
    "        # importance = TRUE to include feature importances\n",
    "        fit <- randomForest(label ~ ., data = df_training, type = 'class', ntree = 500, mtry = 6, importance = TRUE) # Fit the model\n",
    "        df_testing$predict_score <- predict(fit, df_testing) # Predict scores\n",
    "        df_testing$predict_score <- as.numeric(as.character(df_testing$predict_score)) # for our previously defined functions to work, we need to convert the predict_score and label to numeric.\n",
    "        df_testing$label <- as.numeric(as.character(df_testing$label))\n",
    "    }\n",
    "    \n",
    "    \n",
    "    # Get Precision-Recall DataFrame\n",
    "    df_prec_rec <- precision_recall(df_testing, \"label\", \"predict_score\")\n",
    "    \n",
    "    # Calculate accuracy\n",
    "    threshold <- df_prec_rec$predict_score[pop_at_k] # Get the predicted score at K%\n",
    "    df_testing$predict_label <- ifelse(df_testing$predict_score >= threshold, 1, 0) # Predict the label, if > threshold, then 1; if < threshold, then 0\n",
    "    df_testing <- df_testing %>% mutate(accurate = 1*(label == predict_label)) # Generate an indicate of whether the prediction is correct\n",
    "    acc = (sum(df_testing$accurate)/nrow(df_testing)) # Calculate accuracy\n",
    "    df_compare_models$accuracy <- ifelse(df_compare_models$Model == model,acc,df_compare_models$accuracy) # Save accuracy to the DataFrame\n",
    "    \n",
    "    # Calculate precision and save it in the df_compare_models DataFrame\n",
    "    prec_at_k <- precision_at_k(k, df_testing, \"label\", \"predict_score\")\n",
    "    df_compare_models$precision_at_k <- ifelse(df_compare_models$Model == model, prec_at_k, df_compare_models$precision_at_k)\n",
    "    \n",
    "    # Calculate Recall and save it in the df_compare_models DataFrame\n",
    "    rec_at_k <- df_prec_rec$recall[pop_at_k]\n",
    "    df_compare_models$recall_at_k <- ifelse(df_compare_models$Model == model, rec_at_k, df_compare_models$recall_at_k)\n",
    "    \n",
    "}\n",
    "\n",
    "# check results\n",
    "df_compare_models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that at the 10% of the population, the decision tree model is the most accurate. However, the logistic regression model is better at predicting if a graduate will be employed and can capture more graduates that will be employed than the other two models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. References\n",
    "\n",
    "Lou, Tian. (2022, March 18). Machine Learning Model Deployment and Evaluation Using Illinois Unemployment Insurance Data. Zenodo. https://doi.org/10.5281/zenodo.6369160"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Footnotes:**\n",
    "<span id=\"fn1\"> 1. For more information, see <a href='http://www.milbo.org/rpart-plot/prp.pdf'>plotting rpart trees with the rpart.plot package</a>. </span>  \n",
    "\n",
    "[[Go back]](#13)\n",
    "\n",
    "<span id=\"fn2\"> 2. For more information, see <a href='https://cran.r-project.org/web/packages/randomForest/randomForest.pdf'>Breiman and Cutler's Random Forests for Classification and\n",
    "Regression</a>. </span>  \n",
    "\n",
    "[[Go back]](#14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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
}