4.Clustering.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> Nathan Barrett, Benjamin Feder, Joshua Edelmann </center>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6412741\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6412741.svg\" alt=\"DOI\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unsupervised Machine Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes, we do not have a clear outcome variable but nevertheless want to explore and discover any inherent groupings or configurations in the data. Unsupervised machine learning methods can help tackle these problems. Clustering is the most common unsupervised machine learning technique, but you might also be aware of principal components analysis (PCA) or neural networks implementations such as self-organizing maps (SOM). This notebook will provide an introduction to unsupervised machine learning through a clustering example.\n",
    "\n",
    "In this example, we will try to identify patterns within Texas' labor market - namely patterns in the types of employers. Can we isolate specific employers based on derived features addressing their employees' experiences as measured by stability, opportunity, quality, and firm characteristics? We hope to do so using an unsupervised machine learning algorithm. Once we cluster Texas' employers, we can see how Texas' bachelor's degree recipients in the 2015 calendar year fit within the scope of the Texas labor market."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we will focus on **K-Means clustering** (*k* defines the number of clusters), which is considered to be the most commonly-used clustering method. The algorithm works as follows:\n",
    "1. Select *k* (the number of clusters you want to generate).\n",
    "2. Initialize by selecting k points as centroids of the *k* clusters. This is typically done by selecting k points uniformly at random.\n",
    "3. Assign each point a cluster according to the nearest centroid.\n",
    "4. Recalculate cluster centroids based on the assignment in **(3)** as the mean of all data points belonging to that cluster.\n",
    "5. Repeat **(3)** and **(4)** until convergence. Convergence will be further discussed in the sections below.\n",
    "\n",
    "Please reference Chapter 7 of the Big Data and Social Science Textbook and the accompanying class videos for more information on unsupervised machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates how *k*-means clustering can be used to better understand Texas' labor market in 2017. We've already developed a handful of employer measures in a supplemental notebook. We will experiment with a few different values of *k* to see how we can best understand the labor market by looking for differentiation between each of the clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Packages and Set Up\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main R package that we will use for clustering is called `cluster`. We also import all our usual packages for database connection and data manipulation/visualization."
   ]
  },
  {
   "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",
    "\n",
    "# clustering\n",
    "library(cluster, warn.conflicts=F, quietly=T)"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code adjusting overall graph attributes\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 = 12, repr.plot.height = 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Read in the Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will read in a table `employer_yearly_agg`, which contains characteristics of Texas' labor market from 2015-2018, and limit it to information solely from 2017. This nearly mimics the years included in the employment outcomes analyses in the data exploration and visualization notebooks. For more information as to how this table was created, please refer to the Supplemental notebook [\"Supplemental_Employer_Measures.ipynb.\"](Supplemental/Supplemental_Employer_Measures.ipynb)\n",
    "\n",
    "> Note: It is possible to cluster on all employers in all years at the same time as well—just keep in mind that there are a different set of assumptions that are prevalent with that decision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`employer_yearly_agg` contains variables tracking the following characteristics on a quarterly basis, with the average (`avg_`) quarterly values reported within each year. The measures can be broken up within four separate categories:\n",
    "\n",
    "Stability:\n",
    "\n",
    "    - Number of new hires who become full quarter employees (hired in t-1 whom we see in t+1)\n",
    "    - Ratio of full quarter employees to all employees \n",
    "\n",
    "Opportunity:\n",
    "\n",
    "    - Employment growth ((Employment t – employment t-1)/abs( employment at t +employment t-1)/2) \n",
    "    - Number of new hires\n",
    "\n",
    "Job Quality:\n",
    "\n",
    "    - Average earnings per employee \n",
    "    - Average full quarter earnings per employee\n",
    "    - Earnings per employee at 25th percentile\n",
    "    - Earnings per employee at 75th percentile\n",
    "    \n",
    "Firm characteristics:\n",
    "\n",
    "    - Total Payroll\n",
    "    - Total full quarter employment\n",
    "    - Total employment\n",
    "    \n",
    "In addition, the variable `num_quarters` tracks the number of quarters in which the associated `empr_no` existed in Texas' Unemployment Insurance wage records with at least five employees in the given year.\n",
    "\n",
    "> Note: For example, the variable `avg_num_employed` contains the average number of quarterly employees for a given year/employer combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read into R\n",
    "qry <- \"\n",
    "SELECT *\n",
    "FROM tr_tx_2021.dbo.employer_yearly_agg\n",
    "WHERE year = 2017\n",
    "\"\n",
    "\n",
    "emp <- dbGetQuery(con, qry)\n",
    "\n",
    "# see employers\n",
    "head(emp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see all available columns\n",
    "names(emp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Clean the Data\n",
    "\n",
    "A dataset must contain specific attributes in order to run a k-means clustering algorithm. We will examine `emp` through three steps:\n",
    "\n",
    "1. Remove non-explanatory and non-continuous features\n",
    "2. Examine scales across variables\n",
    "3. Analyze missingness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove Non-Explanatory and Non-Continuous Features\n",
    "\n",
    "For any algorithm, non-explanatory variables such as unique ids should be removed from the data set. Additionally, k-means algorithms only work properly with continuous features. This is because k-means calculates its distance measure using euclidean distance, which is the distance between each data point and the centroid of a cluster. It is hard to assign positions for categorical variables in the euclidean space.\n",
    "\n",
    "> There are more sophisticated clustering algorithms that do not use Euclidean distances and thus allow categorical variables in the model. If you are interested in them, you can take a look at the functions `kmodes` and `gower.dist` in R - you will need to download their respective libraries first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove features without explanatory power\n",
    "emp_ml <- emp %>%\n",
    "    select(-c(empr_no, year, num_quarters))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check data type of all variables - make sure all of them are numeric\n",
    "str(emp_ml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine Scales Across Variables\n",
    "\n",
    "If the different explanatory metrics are on a variety of numerical scales, the k-means algorithm will overweigh variables on larger scales due to its distance metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get descriptions of each variable using \"summary\" function\n",
    "summary(emp_ml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert all numeric variables to numeric type otherwise integer64 won't scale\n",
    "emp_ml_num <- emp_ml %>%\n",
    "    sapply(as.numeric)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scale the features since variables like avg_emp_rate are much smaller than avg_total_earnings\n",
    "emp_ml_scale <- scale(emp_ml_num)\n",
    "\n",
    "# View first rows after scaling\n",
    "emp_ml_scale %>% \n",
    "   head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyze Missingness\n",
    "\n",
    "If an employer has missing information in any of the columns, the row will be dropped in the clustering method.\n",
    "\n",
    "> Note that you should **never remove data** if possible - in a real world setting you would likely want to fill any missing data with an imputation method or a baseline assumption. We will discuss missing data during the Inference lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check number of rows (where each row is a unique employer/year combination)\n",
    "nrow(emp_ml_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# na.omit will remove any rows with any NA values\n",
    "emp_ml_scale <- na.omit(emp_ml_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check number of rows after dropping rows with any NA values\n",
    "# see that there is no missing data\n",
    "nrow(emp_ml_scale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Choose the Number of Clusters, *K*\n",
    "\n",
    "Running a *k*-means model is simple: we just need to use `kmeans()` and choose the number of clusters (called `centers`). What number should we choose? Here, we have a bunch of features, so it is hard to visualize the data and decide the proper number.\n",
    "\n",
    "Luckily, there are a few procedures we can use to help us find an appropriate *k* value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elbow Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One technique that can be used to help determine *k* in the elbow method. In the elbow method, we try different k values and calculate the sum of squared errors (`SSE`) after the model converges. Then we plot all the `SSE` by K in a line-chart. The line-chart should resemble an arm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set seed to ensure work is reproducible because k-means has random starting points\n",
    "set.seed(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to compute total within-cluster sum of squares\n",
    "# note that we are using small \"nstart\" value to minimize time of code to run\n",
    "# will discuss nstart later\n",
    "wss <- function(k) {\n",
    "    kmeans(emp_ml_scale, centers=k, nstart = 1)$tot.withinss\n",
    "}\n",
    "\n",
    "# compute and plot wss for k = 1 to k = 15\n",
    "k.values <- 1:15\n",
    "\n",
    "# extract wss values for each k\n",
    "wss_values <- map_dbl(k.values, wss)\n",
    "\n",
    "# plot the resulting SSE for each value of k\n",
    "plot(\n",
    "    k.values, \n",
    "    wss_values, \n",
    "    type = \"b\", pch=19, frame=FALSE,\n",
    "    xlab = \"Number of clusters K\", \n",
    "    ylab = \"Total within-clusters sum of squares\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We try to choose the number around the inflection point, where the change in SSE becomes negligible, indicating that there is little room to improve the model by increasing k (the bend in the elbow). On our graph, the elbow curve begins to flatten around k = 4.\n",
    "\n",
    "> Note: The jaggedness in the curve is due to the default `nstart` value being `1`. `nstart` specifies a number of initial configurations and reports on the best one. However, due to the size of the data, if you increase `nstart` in the function above, it may not converge based on the number of clusters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Try Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model and run on emp_ml_scale\n",
    "set.seed(1)\n",
    "k4 <- kmeans(emp_ml_scale, centers = 4, nstart = 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> An optimal number is usually somewhere between 20 and 50. (See more information in the Resources section - Professor Steorts, Duke University)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see structure\n",
    "str(k4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the `kmeans` function returns the following components, with the most useful for us now as:\n",
    "- `cluster` - an integer indicating a cluster to which each point is allocated\n",
    "- `centers` - a matrix of cluster centers\n",
    "- `totss` - the total sum of squares\n",
    "- `withinss` - vector of within-cluster sum of squares, one component per cluster.\n",
    "- `tot.withinss` - total within-cluster sum of squares, i.e. `sum(withinss)`\n",
    "- `betweenss` - the between-cluster sum of squares, i.e. `totss-tot.withinss`\n",
    "- `size` - the number of points in each cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see size of cluster\n",
    "k4$size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that most of the employers are concentrated in cluster 2. In the perfect world, we would want them to be distributed more evenly across clusters, but in some cases, it may make sense that they wouldn't. Most importantly, we are looking for **high intra-cluster similarity and low inter-cluster similarity**. We will continue to run our model through a battery of tests to further inform our clustering decision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Describe Features across Clusters\n",
    "\n",
    "To start to tease out the potential for high intra-cluster and low inter-cluster similarity, we can take a look at basic descriptives of the employers in these clusters. This will allow us to start to get a sense of some of the important differences across the clusters in the hopes of categorizing each cluster as a separate \"type\" of employer relative to those in other clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove missing values (none here)\n",
    "emp <- na.omit(emp) \n",
    "\n",
    "# add cluster number to the original dataframe\n",
    "frame_4 <- emp %>%                     \n",
    "    mutate(k4.cluster = k4$cluster)  \n",
    "\n",
    "# remove empr_nbr, year, and num_quarters columns\n",
    "frame_4 <- frame_4 %>%\n",
    "    select(-c(empr_no, year, num_quarters))\n",
    "\n",
    "# summarize and add in sizes of each cluster\n",
    "# add suffix \"by_employer\" to each summarize variable\n",
    "frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(across(everything(), list(by_employer=mean))) %>%\n",
    "    mutate(\n",
    "        size = k4$size\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, we can see that our biggest cluster, cluster 2, contains relatively small employers with medium stability but low opportunity and job quality. Cluster 3 also has smaller employers, and on average, they also have low job quality, low stability, and high opportunity. In contrast, cluster 1 contains some of the larger employers in the state that tend to have better job quality while cluster 4 has midsize employers with the best job quality and stability.\n",
    "\n",
    "> When clustering, be cognizant that small numbers of employers per cluster may result in additional redaction due to disclosure regulations concerning the number of employers contributing to each cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare Key Variables of Interest\n",
    "\n",
    "We can also visualize the differences between clusters in more detail by including the standard deviation in conjunction with the average for some of the differentiating features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# focus on avg_avg_earnings, avg_bottom_25_pctile, avg_top_75_pctile (all measures of job quality)\n",
    "# need to convert summarize data frame into a long format with each variable/cluster combination corresponding to a single row\n",
    "frame_4_mean <- frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    select(k4.cluster, avg_avg_earnings, avg_bottom_25_pctile, avg_top_75_pctile) %>%\n",
    "    summarise(across(everything(), mean)) %>%\n",
    "    pivot_longer(-k4.cluster, names_to = \"variable\", values_to = \"mean\")\n",
    "\n",
    "# Save results with standard deviation to a dataframe\n",
    "frame_4_sd <- frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    select(k4.cluster, avg_avg_earnings, avg_bottom_25_pctile, avg_top_75_pctile) %>%\n",
    "    summarise(across(everything(), sd)) %>%\n",
    "    pivot_longer(-k4.cluster, names_to = \"variable\", values_to = \"sd\") %>%\n",
    "    select(-c(k4.cluster, variable))\n",
    "\n",
    "# Bind two dataframes together\n",
    "df <- cbind(frame_4_mean,frame_4_sd)\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare mean + sd across all four clusters for variables of interest\n",
    "ggplot(df, aes(x=k4.cluster, y=mean, fill=k4.cluster)) +\n",
    "    geom_bar(stat=\"identity\", position = position_dodge()) +   # plot bars for the mean values\n",
    "    geom_errorbar(aes(ymax = mean + sd, ymin = mean),          # add standard deviation bars\n",
    "                  width=.2,\n",
    "                  position = position_dodge(.9)) +\n",
    "    facet_grid(. ~ variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can really see the differences in job quality across the clusters as well as where there may be potential overlap. There are additional steps you can take to further justify a clustering choice, but in clustering there may not be a single right answer - every time we run a different number of clusters, interesting patterns about our data can be exposed. Sometimes it may be useful to try running the algorithm on different numbers of clusters and comparing the outputs between models.\n",
    "\n",
    "We really want to know *whether the clusters that we find represent true subgroups in our data*. This could be a crucial input toward choosing the right number of clusters. (See more information on additional methods for selecting `k` in the Resources section - Professor Steorts, Duke University)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare industries\n",
    "\n",
    "In theory, you can also compare clusters by looking at the most common industries within each cluster. When examing industries by cluster, be careful that some employers may have multiple associated NAICS codes within a calendar year."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Relating Back to the Cohort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Armed with results of describing the different types of employers existing within Texas' labor force, we can now relate the results of our unsupervised machine learning model back to our cohort of interest.  In this section we will take a look at how our original cohort of 2015 bachelor's degree recipients fit within these clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read earnings of cohort into R subset to 2017 jobs\n",
    "qry <- \"\n",
    "SELECT *\n",
    "FROM tr_tx_2021.dbo.nb_cohort_wages_link\n",
    "WHERE year(job_date) = '2017'\n",
    "\"\n",
    "df_wages <- dbGetQuery(con, qry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add cluster number to the original dataframe\n",
    "frame_4 <- emp %>%                     \n",
    "    mutate(k4.cluster = k4$cluster)  \n",
    "\n",
    "# Join wages table with frame_4 clustering results\n",
    "df_wages_clus <- df_wages %>%\n",
    "    inner_join(frame_4, by='empr_no')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Employers who Employed a Member of the Cohort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see number of employers by cluster that employed someone in the cohort\n",
    "df_wages_clus %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_cohort = n_distinct(empr_no))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare within cohort employers to all employers in original clusters\n",
    "# Get number of unique employers per cluster in the full dataframe (all employers)\n",
    "frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_all = n_distinct(empr_no))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare with percentages\n",
    "\n",
    "cohort_emp <- df_wages_clus %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_cohort = n_distinct(empr_no))\n",
    "\n",
    "emp_all <- frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_all = n_distinct(empr_no))\n",
    "\n",
    "# Join cohort employers with all employers, andd find percentage\n",
    "cohort_emp %>%\n",
    "    inner_join(emp_all, by = 'k4.cluster') %>%\n",
    "    mutate(percentage = (emp_cohort / emp_all) * 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a reminder, one limitation of our original employers file is that it doesn't include quarterly information for employers with less than 5 employees in a specific quarter. We can see that outside of the employers in Cluster 1 (our smallest cluster by amount of employers), the majority of employers in other clusters did not employ anyone in our cohort in 2017."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Average Earnings and Number of Individuals for Cohort by Cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate average quarterly earnings and number individuals employers in at least one quarter for cohort by cluster\n",
    "df_wages_clus %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(\n",
    "        mean_earnings_cohort = mean(wage),\n",
    "        num_individuals = n_distinct(gradid)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate average earnings for all employees by cluster\n",
    "frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(mean_earnings_all = mean(avg_avg_earnings))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that among our 2015 graduating cohort, those who worked for employers in clusters 1 and 4 earned REDACTED, on average. We can see similar trends in average quarterly wages amongst all employees and those in our cohort acrosss the clusters, except that those in the cohort made less on average, relative to the average employee, especially in the highest-paying cluster. Additionally, we can see that out of those in the cohort who were employed in Texas in 2017, the majority were employed by employers in Clusters 2 and 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Breakdown by Most Common Majors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also see how recipients of the most common degrees within the graduating cohort placed into these clusters.  Recall we have the most common majors for the overall cohort saved as the csv file `common_major.csv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert ADRF username Firstname.Lastname.UserID\n",
    "username <- \"____\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# most common majors\n",
    "common_major <- read_csv(sprintf(\"U:\\\\%s\\\\TX Training\\\\Results\\\\common_major.csv\", username))\n",
    "\n",
    "common_major"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# limit df_wages_clus to only common majors\n",
    "\n",
    "#Create a 2 digit CIP program code from the full CIP code in `gradmaj` and change to numeric\n",
    "df_wages_clus <- df_wages_clus %>%\n",
    "    mutate(\n",
    "        CIP_Program = as.numeric(substring(gradmaj, 1, 2))\n",
    "    )\n",
    "\n",
    "# limit to common majors\n",
    "df_wages_clus_major <- df_wages_clus %>%\n",
    "    inner_join(common_major, by = 'CIP_Program')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find proportions employment by cluster within the most common majors\n",
    "prop_by_cip_cluster <- df_wages_clus_major %>%\n",
    "    group_by(CIPTitle2010, k4.cluster) %>%\n",
    "    summarize(\n",
    "        n_inds = n_distinct(gradid), \n",
    "        n_employers = n_distinct(empr_no)\n",
    "    ) %>%\n",
    "    ungroup() %>%\n",
    "    group_by(CIPTitle2010) %>%\n",
    "    mutate(\n",
    "        prop_by_cluster = n_inds / sum(n_inds)\n",
    "    ) %>%\n",
    "    ungroup()\n",
    "\n",
    "head(prop_by_cip_cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show in bar plot\n",
    "prop_by_cip_cluster %>%\n",
    "    ggplot(aes(x=word(CIPTitle2010, 1), y = prop_by_cluster, fill=as.character(k4.cluster))) +\n",
    "    geom_col(position=position_dodge()) +\n",
    "    labs(\n",
    "        x = \"Major\",\n",
    "        y = \"Proportion by Cluster\",\n",
    "        fill = \"Cluster\"\n",
    "    ) +\n",
    "    ylim(0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As somewhat expected given the sorting of the overall cohort into primarily employers in clusters 2 and 4, all degree recipients within the most common major groupings were most often employed by employers in clusters 2 and 4. However, you can note that those receiving degrees in Engineering, Business, and Health were more likely to be employed by employers in Cluster 4, the cluster with the highest average quarterly pay, whereas those majoring in Biological Sciences and Multi/Interdisciplinary Studies were more likely to be employed by employers in Cluster 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red><h3> Checkpoint 1: By Institution Type </h3></font> \n",
    "\n",
    "Repeat this breakdown by institution type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the file with hints and solutions\n",
    "source(\"nb4_hints_and_solutions.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove the `#` symbol in the checkpoint hints and solutions cells to get hints and solutions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#checkpoint_1.hint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#checkpoint_1.solution()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find proportions employment by cluster within institution types\n",
    "# replace ___\n",
    "prop_by_cip_cluster_inst <- df_wages_clus %>%\n",
    "    group_by(___, k4.cluster) %>%\n",
    "    summarize(\n",
    "        n_inds = n_distinct(gradid), \n",
    "        n_employers = n_distinct(empr_no)\n",
    "    ) %>%\n",
    "    ungroup() %>%\n",
    "    group_by(___) %>%\n",
    "    mutate(\n",
    "        prop_by_cluster = n_inds / sum(n_inds)\n",
    "    ) %>%\n",
    "    ungroup()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# graph using bar plot\n",
    "# replace ___\n",
    "___ %>%\n",
    "    ggplot(aes(x=as.character(___), y = prop_by_cluster, fill=as.character(___))) +\n",
    "    geom_col(position=position_dodge()) +\n",
    "    labs(\n",
    "        x = \"___\",\n",
    "        y = \"Proportion by Cluster\",\n",
    "        fill = \"Cluster\"\n",
    "    ) +\n",
    "    ylim(0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Chappell, Joseph, Feder, Benjamin, & Barrett, Nathan. (2022, April 1). Cluster Analysis Using Tennessee Employment Data. Zenodo. https://doi.org/10.5281/zenodo.6407271"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources\n",
    "- UC Business Analytics R Programming Guide: https://uc-r.github.io/kmeans_clustering\n",
    "- Rebecca Steorts, Assistant Professor, Duke University, Department of Statistical Science, Data Mining and Machine Learning course: https://github.com/resteorts/data-mine/tree/master/lectures_2018/10-unsupervise/10-kmeans.pdf"
   ]
  }
 ],
 "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
}