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> Joseph Chappell, Benjamin Feder, Nathan Barrett </center>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6407271\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6407271.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 patterns 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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clustering is used to group data points together that are similar to each other. Optimally, a given clustering method will produce groupings with high intra-cluster (within) similarity and low inter-cluster (between) similarity. Clustering algorithms typically require a distance or similarity metric to generate clusters. They take a dataset and a distance metric (and sometimes additional parameters), and they generate clusters based on that distance metric. The most common distance metric used is Euclidean distance, but other commonly-used metrics are Manhattan, Minkowski, Chebyshev, Cosine, Hamming, Pearson, and Mahalanobis.\n",
    "\n",
    "Most clustering algorithms also require the user to specify the number of clusters (or some other parameter that indirectly determines the number of clusters) in advance as a parameter. This is often difficult to do a priori and typically makes clustering an iterative and interactive task. Another aspect of clustering that makes it interactive is often the difficulty in automatically evaluating the quality of the clusters. While various analytical clustering metrics have been developed, the best clustering is task-dependent and thus must be evaluated by the user. There may be different clusterings that can be generated with the same data. You can imagine clustering similar news stories based on the topic content, on the writing style or sentiment. The right set of clusters depends on the user and the task at hand. Clustering is therefore typically used for exploring the data, generating clusters, exploring the clusters, and then rerunning the clustering method with different parameters or modifying the clusters (by splitting or merging the previous set of clusters). Interpreting a cluster can be nontrivial: you can look at the centroid of a cluster, look at frequency distributions of different features (and compare them to the prior distribution of each feature), or other aspects.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm stops or \"converges\" when the assignments do not change from one iteration to the next. The final set of clusters, however, depends on the starting points. If initialized differently, it is possible that different clusters are obtained. One common practical trick is to run *k*-means several times, each with different (random) starting points. The *k*-means algorithm is fast, simple, and easy to use, and is often a good first clustering algorithm to try and see if it fits your needs. When the mean of the data points cannot be computed, a related method called *K-medoids* can be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates using *k*-means clustering to better understand Tennessee's labor market in 2017. We've already developed a handful of employer measures in a supplemental notebook. We will try 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 Tennessee's labor market from 2015-2019. This 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`employer_yearly_agg` contains variables tracking the following characteristics on a quarterly basis, with the average (`avg_`) quuarterly values reported within each year. The measures can be broken up within three separate categories:\n",
    "\n",
    "Stability:\n",
    "\n",
    "    - Number of full quarter employees\n",
    "    - Separation growth rate\n",
    "\n",
    "Opportunity:\n",
    "\n",
    "    - Hiring growth rate\n",
    "    - Number of employees\n",
    "    - Employment growth rate\n",
    "\n",
    "Job Quality:\n",
    "\n",
    "    - Total payroll\n",
    "    - Earnings per employee at the 25th percentile\n",
    "    - Earnings per employee at the 75th percentile\n",
    "    - Total payroll for full quarter employees\n",
    "    - Average earnings per employee\n",
    "    - Average earnings per full quarter employee\n",
    "\n",
    "    \n",
    "In addition, the variable `num_quarters` tracks the number of quarters in which the associated `empr_nbr` existed in Tennessee's 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_tn_2021.dbo.employer_yearly_agg\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": [
    "To get a glimpse at the labor market in a specific calendar year, we will subset our data frame to include information only for one year: 2017.\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subset a dataframe by rows from 2017\n",
    "\n",
    "emp <- emp %>%\n",
    "    filter(year == 2017)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that we only have 2017\n",
    "emp %>%\n",
    "    distinct(year)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clean the Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to remove the `empr_nbr`, `year`, and `num_quarters` variables from our data frame since these features do not provide any explanatory power for our k-means algorithm. 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 empr_nbr, year, and num_quarters\n",
    "emp_ml <- emp %>%\n",
    "    select(-c(empr_nbr, 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": [
    "**It is important that we consider scaling these features** before we compute *k*-means clustering, especially if the metrics are on a variety of numerical scales. Let's see if they are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get descriptions of each variable using \"summary\" function\n",
    "summary(emp_ml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that we have variables on different numerical scales - we can scale them using `scale()` function on our data frame `emp_ml`. Before we do so, though, let's convert all variables to `numeric` types, as the `scale()` function will return `Inf` and `-Inf` for the variables that are of `integer64` type."
   ]
  },
  {
   "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\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": [
    "Before running a clustering algorithm, we need to make sure that there are no missing values. Here we will use `na.omit()` to remove all rows with any NA values. If an employer has missing information in any of the columns, a row will be dropped.\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": [
    "## 2. Choose the Number of Clusters, *K*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 11 features, so it is hard to visualize the data and decide the proper number by using our eyes. Let's start with a small number, such as 3, and see how the results look.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because *k*-means clustering will generate different results (due to different starting points), we will set a seed so that the work in this notebook can be reproducible using the `set.seed()`. To get the same results, you must use the same seed before running the clustering algorithm every time. Luckily, if you set the same seed as your collaborators and are running the same *k*-means algorithm, you will see the same results, even if you are working in different environments, such as jupyter notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### k = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model and run on emp_ml\n",
    "set.seed(1)\n",
    "k3 <- kmeans(emp_ml_scale, centers=3, nstart=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> `nstart` specifies a number of initial configurations and reports on the best one - 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(k3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`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": "markdown",
   "metadata": {},
   "source": [
    "Let's check the size of each cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k3$size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that most of the employers are concentrated in cluster 1. 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**.\n",
    "\n",
    "Are there major differences in the characteristics of employers in each cluster?\n",
    "\n",
    "We can take a look at basic descriptives of the employers in these clusters by adding our clustering results to the original dataframe, `emp`, and call this dataframe `frame_3`. 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\n",
    "emp <- na.omit(emp) \n",
    "\n",
    "# add cluster number to the original dataframe\n",
    "frame_3 <- emp %>%                     \n",
    "    mutate(k3.cluster = k3$cluster)  \n",
    "\n",
    "# remove empr_nbr, year, and num_quarters columns\n",
    "frame_3 <- frame_3 %>%\n",
    "    select(-c(empr_nbr, year, num_quarters))\n",
    "\n",
    "# summarize and add in sizes of each cluster\n",
    "frame_3 %>%\n",
    "    group_by(k3.cluster) %>%\n",
    "    summarize_all(\"mean\") %>%\n",
    "    mutate(\n",
    "        size = k3$size\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, we can see that our biggest cluster, cluster 1, contains relatively small employers that pay their employees lower wages and have less turnover. Cluster 3 also has relatively small employers, and on average, they pay their employees less and have higher turnover rates than those in the other two clusters. In contrast, cluster 2 contains some of the larger employers in the state that tend to pay better and have higher job stability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate clusters\n",
    "\n",
    "One simple way to characterize the resulting clusters is to compare the summary stats between key variables of interest. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualize the differences between the clusters in more detail by finding the mean and standard deviation for the following variables: `avg_avg_earnings`, `avg_bottom_25_pctile`, and `avg_top_75_pctile`. First, we will need to convert our summarized data frame of the mean and standard deviations for the desired variables into a long format, with each variable/cluster combination corresponding to a distinct row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save results with mean to a dataframe\n",
    "frame_3_mean <- frame_3 %>%\n",
    "    group_by(k3.cluster) %>%\n",
    "    select(k3.cluster, avg_avg_earnings, avg_bottom_25_pctile, avg_top_75_pctile) %>%\n",
    "    summarize_all(\"mean\") %>%\n",
    "    mutate(\n",
    "        avg_avg_earnings = as.numeric(avg_avg_earnings)\n",
    "        ) %>%\n",
    "    pivot_longer(-k3.cluster, names_to = \"variable\", values_to = \"mean\")\n",
    "\n",
    "# Save results with standard deviation to a dataframe\n",
    "frame_3_sd <- frame_3 %>%\n",
    "    group_by(k3.cluster) %>%\n",
    "    select(k3.cluster, avg_avg_earnings, avg_bottom_25_pctile, avg_top_75_pctile) %>%\n",
    "    summarize_all(\"sd\") %>%\n",
    "    mutate(\n",
    "        avg_avg_earnings = as.numeric(avg_avg_earnings)\n",
    "    ) %>%\n",
    "    pivot_longer(-k3.cluster, names_to = \"variable\", values_to = \"sd\") %>%\n",
    "    select(-c(k3.cluster, variable))\n",
    "\n",
    "# Bind two dataframes together\n",
    "df <- cbind(frame_3_mean,frame_3_sd)\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use this data frame to visualize the means and standard deviations in our 3 clusters for these 3 variables: `avg_avg_earnings`, `avg_bottom_25_pctile`, and `avg_top_75_pctile` using a bar plot, separating the mean and standard deviations values for each variable into different grids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(df, aes(x=k3.cluster, y=mean, fill=k3.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) +                                  # plot by 3 variables of interest\n",
    "    ggtitle(\"Overlap in mean + SD between clusters 1 and 3\") +  # add title\n",
    "    xlab(\"Clusters\") +                                          # add label for x-axis\n",
    "    ylab(\"Mean\") +                                              # add label for y-axis\n",
    "    theme(text = element_text(size=16),                         # increase text font\n",
    "          axis.text.x = element_text(size=18, face=\"bold\"),     # increase text font on x-axis and make it bold\n",
    "          legend.position = \"none\")                             # remove legend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization functions\n",
    "\n",
    "We can also create a function to facilitate visualizing different columns in a similar way. The function takes the mean, standard deviation, and title of the visualization as its arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save it to a dataframe\n",
    "frame_3_mean_sd <- frame_3 %>%\n",
    "    group_by(k3.cluster) %>%\n",
    "    select(k3.cluster, avg_avg_earnings, avg_bottom_25_pctile, avg_top_75_pctile) %>%\n",
    "    summarise_all(list(mean = \"mean\", sd = \"sd\"))\n",
    "\n",
    "# Visualize average earnings by cluster\n",
    "viz <- function(mean, sd, title) {\n",
    "    ggplot(frame_3_mean_sd, aes(x=k3.cluster, y=mean, fill=k3.cluster)) +\n",
    "    geom_bar(position = position_dodge(), stat=\"identity\", fill=\"gray\") +\n",
    "    geom_errorbar(aes(ymax= mean + sd, ymin = mean),\n",
    "                  width=.2,\n",
    "                  position = position_dodge(.9)) +\n",
    "    ggtitle(title) +\n",
    "    xlab(\"Clusters\") +\n",
    "    ylab(\"Mean\") +\n",
    "    theme(text = element_text(size=16),\n",
    "          axis.text.x = element_text(size=18, face=\"bold\"),\n",
    "          legend.position = \"none\")\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz(as.numeric(frame_3_mean_sd$avg_avg_earnings_mean), as.numeric(frame_3_mean_sd$avg_avg_earnings_sd), \"The Average Quarterly Earnings Are Much Higher in Cluster 2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz(as.numeric(frame_3_mean_sd$avg_bottom_25_pctile_mean), as.numeric(frame_3_mean_sd$avg_bottom_25_pctile_sd), \"The Average 25th Percentile of Quarterly Earnings are Much Higher in Cluster 2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz(\n",
    "    as.numeric(frame_3_mean_sd$avg_top_75_pctile_mean), \n",
    "    as.numeric(frame_3_mean_sd$avg_top_75_pctile_sd), \n",
    "    \"The Average 75th Percentile of Quarterly Earnings are Much Higher in Cluster 2 \\n But Have More Spread\"\n",
    ")"
   ]
  },
  {
   "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 by matching the employers in the original data frame (`emp`) back to the UI wage records for the given year and `empr_nbr`. When examing industries by cluster, be aware that some employers may have multiple associated NAICS codes within a calendar year. If you run into this situation where there are multiple non-999999 industry codes, you will need to devise a set of consistent decision rules to apply across the entire table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting *k*\n",
    "\n",
    "How do we know if we chose an optimal number of clusters to describe our data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elbow method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the *Elbow method* as an input in selecting the optimal cluster number. Recall that *k*-means starts with k random cluster centers (centroids), assigns each data point to the closest centroid, and calculates the distances between each point and the centroid. Then it moves the positions of the centroids and repeats the previous steps until there is convergence. 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(1)\n",
    "\n",
    "# function to compute total within-cluster sum of square\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(k.values, wss_values, \n",
    "    type = \"b\", pch=19, frame=FALSE,\n",
    "    xlab = \"Number of clusters K\", \n",
    "    ylab = \"Total within-clusters sum of squares\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that SSE decreases as we increase k. Here, it decreases faster when k is small. As k increases, the reduction in SSE becomes smaller. 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`. 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": [
    "Let's run the model with 4 clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(1)\n",
    "k4 <- kmeans(emp_ml_scale, centers = 4, nstart = 20)\n",
    "k4$size"
   ]
  },
  {
   "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",
    "# remove empr_nbr, year, and num_quarters columns\n",
    "frame_4 <- frame_4 %>%\n",
    "    select(-c(empr_nbr, year, num_quarters))\n",
    "\n",
    "# summarize and add in sizes of each cluster\n",
    "frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarize_all(\"mean\") %>%\n",
    "    mutate(\n",
    "        size = k4$size\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The large employers are now even more separated into Cluster 3. These employers tend to pay reasonably well, with low turnover rates, but they do not pay as well as those in Cluster 2. Clusters 2 and 4 contain employers of similar size, but the difference in pay is quite glaring between the two. Cluster 1 consists primarily of much smaller employers with high turnover rates and low pay."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which clustering results - `frame_3` or `frame_4` - do you prefer? Do you think it could be optimal to choose more clusters?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, in clustering there is no single right answer - every time we run a different number of clusters, interesting patterns about our data can be exposed. However, what we do want to know is *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).\n",
    "\n",
    "Experiment with different numbers of clusters in the Checkpoint 1 below - given knowledge about Tennessee's labor market in 2017, which number of clusters makes most sense to you?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red><h3> Checkpoint 1: Run a K-Means clustering model </h3></font> \n",
    "\n",
    "1. Take a look again at the elbow curve, which number(s) do you think is (are) optimal?\n",
    "\n",
    "2. Choose a cluster number that you think is best (other than 3 or 4). Use `kmeans()` to run a k-means clustering model with the number you choose. Save your results and features in `frame_k`. \n",
    "\n",
    "3. Compare your results with the results we got previously. Do you find any differences? Are the results improved, in your opinion?\n",
    "\n",
    "Hint: in the Elbow method graph, it looks like 12 could be another optimal cluster - you can try with 12 clusters and see the differences, but remember: to achieve convergence, you’ll probably need to set nstart=1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run clustering algorithm using different k\n",
    "\n",
    "set.seed(1)\n",
    "__ <- kmeans(emp_ml_scale, centers = ___, nstart = __)\n",
    "__$size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# explore differences between clusters\n",
    "\n",
    "__ <- emp %>%                     \n",
    "    mutate(__.cluster = __$cluster)  \n",
    "\n",
    "# remove empr_nbr, year, and num_quarters columns\n",
    "___ <- ___ %>%\n",
    "    select(-c(empr_nbr, year, num_quarters))\n",
    "\n",
    "# summarize and add in sizes of each cluster\n",
    "___ %>%\n",
    "    group_by(___) %>%\n",
    "    summarize_all(\"mean\") %>%\n",
    "    mutate(\n",
    "        size = __$size\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cohort's Employers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we will take a look at our cohort's employers, and identify the clusters they belong to based on the `frame_4` clustering results. Additionally, we will look at our cohort's distribution across the four clusters and compare their employment outcomes.\n",
    "\n",
    "> We will need to subset `nb_cohort_wages_link` to just jobs in 2017 in order to match 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_tn_2021.dbo.nb_cohort_wages_link\n",
    "where year(job_date) = '2017'\n",
    "\"\n",
    "df_wages <- dbGetQuery(con, qry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will not subset to the dominant employer for each `hashed_ssn` in `df_wages`."
   ]
  },
  {
   "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_nbr')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's look at the number of employers that employed at least one individual in our original cohort in 2017 by cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Group by clusters and find number of unique employers in each cluster\n",
    "df_wages_clus %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_cohort = n_distinct(empr_nbr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare that to the percentage of all employers in our clusters that hired at least one individual in our cohort in the calendar year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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_nbr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of scrolling between these two outputs to compare the proportion of employers in each cluster that employed at least one individual in our cohort, we can combine them into one data frame using `inner_join()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save cohort and all employers dataframes\n",
    "\n",
    "cohort_emp <- df_wages_clus %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_cohort = n_distinct(empr_nbr))\n",
    "\n",
    "emp_all <- frame_4 %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(emp_all = n_distinct(empr_nbr))\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 3, the majority of employers in other clusters did not employ anyone in our cohort in 2017"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare average earnings of our cohort by cluster with average earnings of all employees in each cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate average earnings for TN graduates by cluster\n",
    "# also see number of individuals within each cluster\n",
    "df_wages_clus %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarise(\n",
    "        mean_earnings_cohort = mean(wge_amt),\n",
    "        num_individuals = n_distinct(SSN)\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-16 graduating cohort, those who worked for employers in clusters 2 earned more, 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. Additionally, We can see that out of those in the cohort who were employed in Tennessee in 2017, the majority were employed by employers in Clusters 3 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": "markdown",
   "metadata": {},
   "source": [
    "<font color=red> Note that you need to change the directory in read_csv() statements below. Replace \". .\" with your username.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# most common majors\n",
    "common_major <- read_csv(\"U:\\\\..\\\\TN Training\\\\Results\\\\common_major.csv\")\n",
    "\n",
    "common_major"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can limit employment information from `df_wages_clus` to just those receiving degrees in the most common majors with a `filter()` command. Afterwards, we can count the number of individuals within each  major/cluster combination then finding the proportion of individuals with each degree that are employed by employers in the four clusters."
   ]
  },
  {
   "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 %>% \n",
    "    filter(CIP_Family %in% common_major$CIP_Family) %>%\n",
    "    group_by(CIP_Family, k4.cluster) %>%\n",
    "    summarize(\n",
    "        n_inds = n_distinct(SSN)\n",
    "    ) %>%\n",
    "    ungroup() %>%\n",
    "    group_by(CIP_Family) %>%\n",
    "    mutate(\n",
    "        prop_by_cluster = n_inds / sum(n_inds)\n",
    "    ) %>%\n",
    "    ungroup() \n",
    "\n",
    "head(prop_by_cip_cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there are a lot of rows (one for each `CIP_Family`/`k4.cluster` combination), we can more easily discern the breakdown of employment by cluster within each major grouping with a simple bar plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prop_by_cip_cluster %>%\n",
    "    ggplot(aes(x=word(CIP_Family, 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 cluster 4, all degree recipients were most commonly employed by employers in cluster 4 in Tennessee in 2017. However, you can note that those receiving degrees in Engineering and Computer and Information Sciences were more likely to be employed by employers in Cluster 2, the cluster with the highest average quarterly pay, than those receiving degrees in the other three majors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=red><h3> Checkpoint 2: Cohort's Employers </h3></font> \n",
    "\n",
    "Repeat this exploration of employer clusters for the sample you created in Checkpoint 1. Which clusters hire the greatest proportion of our original cohort? Which pay the most?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "Feder, Benjamin, Barrett, Nathan, Simone, Sean, & Seith, David. (2022, March 30). Cluster Analysis using New Jersey Employment Data. Zenodo. https://doi.org/10.5281/zenodo.6399336"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}