Supplemental_Machine_Learning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> <img style=\"float: center;\" src=\"images/CI_horizontal.png\" width=\"450\">\n",
    "<center>\n",
    "    <span style=\"font-size: 1.5em;\">\n",
    "        <a href='https://www.coleridgeinitiative.org'>Website</a>\n",
    "    </span> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Unsupervised Machine Learning</center>\n",
    "<br>\n",
    "<center>Benjamin Feder, Brian Kim, Ekaterina Levitskaya, Allison Nunez. </center>\n",
    "<br>\n",
    "<a href=\"https://doi.org/10.5281/zenodo.6463919\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.6463919.svg\" alt=\"DOI\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are problems where the goal is not necessarily prediction, but instead it is to discover any inherent groupings or patterns in the data. Unsupervised machine learning methods can help tackle these situations. Clustering is the most common unsupervised machine learning technique, but other methods such as principal components analysis (PCA) or neural networks implementations (e.g. self-organizing maps) are also widely-used in the field. In fact, you have actually already done a form of unsupervised learning already in the text analysis notebook. Topic modeling with Latent Dirichlet Allocation is a form of unsupervised learning. This notebook will introduce unsupervised machine learning through the lens of 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 the distance metric. The most common distance metric 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 another 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 ones are task-dependent and thus must be evaluated by the user. There may be different clusterings that can be generated with the same data. One can imagine clustering similar news stories based on the topic content, based 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 re-running 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: one can look at the centroid of a cluster, frequency distributions of different features (and compare them to the prior distribution of each feature), or other aspects.\n",
    "\n",
    "This notebook focuses 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 to generate).\n",
    "2. Initialize the algorithm 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 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 to obtain different clusters. 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 one's needs. When the mean of the data points cannot be computed (i.e. for categorical variables), a related method called *K-medoids* can be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates how *k*-means clustering can be employed to better understand the types of PhD students based on funding history. This notebook covers a few different values of *k* to see how to best characterize the PhD funding experiences 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 primary R package used in this notebook for clustering is called `cluster`. The usual packages for database connection and data manipulation/visualization will also be imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(warn=-1)\n",
    "\n",
    "#database interaction imports\n",
    "suppressMessages(library(odbc))\n",
    "\n",
    "# for data manipulation/visualization\n",
    "suppressMessages(library(tidyverse))\n",
    "suppressMessages(library(ggplot2))\n",
    "\n",
    "# clustering\n",
    "suppressMessages(library(cluster))\n",
    "suppressMessages(library(reshape2))\n",
    "\n",
    "options(warn=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connect to the database\n",
    "con <- DBI::dbConnect(odbc::odbc(),\n",
    "                     Driver = \"SQL Server\",\n",
    "                     Server = \"msssql01.c7bdq4o2yhxo.us-gov-west-1.rds.amazonaws.com\",\n",
    "                     Trusted_Connection = \"True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Read in the Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To start the clustering example, the proper dataset must be crafted. In this example, the dataset will contain funding history from UMETRICS, as well as time measures regarding degree completion from SED. The individuals in this dataset must individually appear in both the SED and UMETRICS tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read into R\n",
    "qry <- \"\n",
    "select sed.drf_id, team_size, any_non_federal\n",
    "from ds_nsf_ncses.dbo.nsf_sed sed\n",
    "join tr_ncses_2021.dbo.sed_umetrics_xwalk xwalk\n",
    "on sed.drf_id = xwalk.drf_id\n",
    "join ds_iris_umetrics.dbo.semester sem\n",
    "on xwalk.emp_number = sem.emp_number\n",
    "where sed.phdfy = '2015' \n",
    "\"\n",
    "students <- dbGetQuery(con, qry)\n",
    "\n",
    "# see the dataframe\n",
    "head(students)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add total number of semesters of funding:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add total number of semesters of funding\n",
    "students <- students %>%\n",
    "    group_by(drf_id) %>%\n",
    "    mutate(sem_funding = n())  # count the number of rows per each unique drf_id\n",
    "\n",
    "students %>% head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add a flag for federal and non-federal funding:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If team_size > 0, then flag as federal (add 1), if team_size is missing, then add 1 in the \"any_non_federal\" variable\n",
    "students <- students %>%\n",
    "    mutate(federal = ifelse(team_size > 0, 1, 0),\n",
    "           non_federal = ifelse(is.na(team_size) | any_non_federal == 1, 1, 0))\n",
    "\n",
    "head(students)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count total number of semesters with federal and non-federal funding:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "students <- students %>%\n",
    "            group_by(drf_id) %>%                                 # for each unique drf_id\n",
    "            mutate(sem_federal = sum(federal, na.rm = TRUE),     # sum values in the federal flag\n",
    "                   sem_non_federal = sum(non_federal)) %>%       # sum values in the any_non_federal flag\n",
    "            select(-c(team_size, any_non_federal, federal, non_federal)) %>%  # remove these columns\n",
    "            slice(1) %>%                                         # deduplicate rows, choose just 1 row per unique drf_id\n",
    "            replace(is.na(.), 0)                                 # fill NA with 0\n",
    "\n",
    "head(students)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Clean the Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Missing Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running k-means clustering on `students`, the data frame must be further cleaned. Clustering algorithms will not work if there are any missing values in any of the features. Here, the `na.omit()` function removes all rows with any NA values. (If a student has missing information in any of the columns, a row will be dropped).\n",
    "\n",
    "> Note: **never remove data** if possible - in a real world setting any missing data would be likely filled with an imputation or baseline assumption. Imputation of missing data will be discussed during the Inference session."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check number of rows (where each row is a unique student)\n",
    "nrow(students)\n",
    "head(students)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove all missing data points before running clustering\n",
    "# na.omit will remove any rows with any NA values\n",
    "students_ml <- na.omit(students)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check number of rows after dropping rows with any NA values\n",
    "nrow(students_ml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explanatory power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The **drf_id** variable should be removed from the data frame since the feature does not provide any explanatory power for the k-means algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove drf_id\n",
    "students_features <- students_ml %>%\n",
    "    ungroup() %>%\n",
    "    select(-drf_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(students_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variable Types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Recall that 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, such as k-medoids as mentioned in the Introduction to Clustering [section](#Introduction-to-Clustering). R packages such as `klaR` and `cba` can be used to implement algorithms to handle categorical variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check data type of all variables - make sure all of them are numeric\n",
    "str(students_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Scaling Numerical Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to consider scaling these features before computing *k*-means clustering, especially if the metrics are on a variety of numerical scales. The scales of each of the variables can be verified using `summary`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get descriptions of each variable using \"summary\" function\n",
    "summary(students_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the variables are on different numerical scales, the `scale()` function can be used to set all of these variables onto the same scale.\n",
    "\n",
    "> Note: It is not always the case that the variables will be on different scales. Exercise caution before implementing any variable scaling.\n",
    "\n",
    "In this case the variables are comparable, as all of them represent semesters of funding, therefore, the code below is commented out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is the code for scaling of the features, if needed\n",
    "#students_features <- scale(students_features)\n",
    "\n",
    "# View first rows after scaling\n",
    "#head(students_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data frame is now ready for clustering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Choose the Number of Clusters, *K*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running a *k*-means model is simple once the underlying dataset is ready: the `kmeans()` function will implement the algorithm as long as the number of clusters (called `centers`) is specified. What k value is optimal when starting? With multiple features, it is hard to visualize the data and decide the proper number by using the eyes. K will be set to a low number, such as 3, to see these results.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because *k*-means clustering will generate different results (due to different starting points), set a seed so that the work in this notebook can be reproducible using `set.seed`. To get the same results, the same seed must be set before running the clustering algorithm every time. Luckily, if the same seed is set and all collaborators are running the same *k*-means algorithm, the results will be consistent, even if collaborators are working in different environments, i.e. Jupyter notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### k = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model and run on students_features\n",
    "set.seed(1)\n",
    "k3 <- kmeans(students_features, centers=3, nstart=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The `nstart` argument 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 components of k3\n",
    "str(k3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`kmeans` function returns the following components - among the most useful are:\n",
    "- `cluster` - integers indicating the cluster assignment for each row\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": [
    "As indicated above, the size of each cluster can be found by analyzing the `size` component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see cluster size\n",
    "k3$size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the students are concentrated in clusters 1 and 3. In the perfect world, the students would be distributed more evenly across clusters, but oftentimes, it may make sense that they would not. Most importantly, the goal is for high intra-cluster similarity and low inter-cluster similarity.\n",
    "\n",
    "To determine if `k=3` may accomplish the goal, basic descriptive statistics of the students within these clusters can be analyzed. To start the process, the clustering results can be added to the original data frame, `students` before finding the means of the various features within each cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " # add cluster number to the original dataframe and call frame_3\n",
    "frame_3 <- data.frame(students_ml, k3$cluster)\n",
    "\n",
    "# find within-cluster means\n",
    "frame_3 %>% select(-drf_id) %>%\n",
    "    group_by(k3.cluster) %>%\n",
    "    summarize_all(\"mean\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although cluster means should not be the only measure analyzed (e.g. standard deviation, median, etc.), the means can allow for a quick characterization of each of the clusters. How can these three clusters be described?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 2 and 3 are characterized by a higher number of semesters of funding in comparison with cluster 1. Cluster 2 is characterized by students who had a higher number of non-federal funding and cluster 3 is characterized by students who had, on the contrary, a higher number of federal funding. Cluster 1 is characterized by students who in general didn't have many semesters of any funding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate clusters\n",
    "\n",
    "One simple way to evaluate resulting clusters is to compare the summary statistics between the key variables of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The differences between the clusters can also be visualized in more detail by finding the mean and standard deviation of the variables. In this example, the variables of interest are primary source of support and time to degree completion from the SED data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read into R\n",
    "qry <- \"\n",
    "select drf_id, srceprim, ttddoc\n",
    "from ds_nsf_ncses.dbo.nsf_sed\n",
    "where phdfy = '2015'\n",
    "\"\n",
    "sed <- dbGetQuery(con, qry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add in more SED data to frame_3\n",
    "frame_3 <- frame_3 %>% \n",
    "    inner_join(sed, by = c('drf_id'))\n",
    "\n",
    "# head(frame_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recode the primary source of support into the following categories, using case_when\n",
    "frame_3 <- frame_3 %>%\n",
    "            mutate(srceprim = case_when((srceprim == 'A' | srceprim == 'B') ~ 'Fellowship, scholarship or dissertation grant',\n",
    "                                        srceprim == 'C' ~ 'Teaching assistantship',\n",
    "                                      (srceprim == 'D' | srceprim == 'E' | srceprim == 'F' | srceprim == 'G') ~ 'Research assistantship or traineeship',\n",
    "                                       (srceprim == 'H' | srceprim == 'I' | srceprim == 'J' | srceprim == 'K') ~ 'Own resources',\n",
    "                                       (srceprim == 'L' | srceprim == 'M' | srceprim == 'N' | srceprim == '') ~ 'Other sources'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The breakdown by primary source of support within each cluster can now be analyzed both numerically and visually.\n",
    "\n",
    "> The proportions will be used since there are different counts within each of the clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Proportions by sex within each cluster\n",
    "cluster_counts_by_srceprim <- frame_3 %>%\n",
    "    count(k3.cluster, srceprim) %>%\n",
    "    group_by(k3.cluster) %>%\n",
    "    mutate(prop = n/sum(n)) %>%\n",
    "    select(-n) %>%\n",
    "    ungroup()\n",
    "\n",
    "\n",
    "cluster_counts_by_srceprim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(cluster_counts_by_srceprim, aes(x = k3.cluster, y=prop, fill=srceprim)) +\n",
    "    geom_bar(stat = 'identity', position = 'stack') +\n",
    "    labs(\n",
    "        x = \"Cluster\",\n",
    "        y = \"Proportion\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same process can be followed for the time to degree completion, with the addition of the standard deviation to account for the spread within each of the clusters for the numerical variable of interest (**ttddoc**)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# average time to degree by cluster\n",
    "time_to_degree <- frame_3 %>%\n",
    "    group_by(k3.cluster) %>% \n",
    "    summarize(\n",
    "        avg_time_to_degree = mean(ttddoc, na.rm = TRUE),\n",
    "        sd = sd(ttddoc, na.rm = TRUE)\n",
    "    )\n",
    "\n",
    "time_to_degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(time_to_degree, aes(x=k3.cluster, y=avg_time_to_degree, fill=k3.cluster)) +\n",
    "    geom_bar(stat=\"identity\", position = position_dodge()) +   # plot bars for the mean values\n",
    "    geom_errorbar(aes(ymax= avg_time_to_degree + sd, ymin = avg_time_to_degree-sd),            # add standard deviation bars\n",
    "                  width=.2,\n",
    "                  position = position_dodge(.9)) +\n",
    "    labs(\n",
    "        title = \"No Significant Difference in TTD by Cluster\",  # add title\n",
    "        x = \"Clusters\",                                       \n",
    "        y = \"Mean\"\n",
    "    ) +                                              \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": [
    "It looks like federal vs non-federal funding doesn't matter for time to degree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting *k*\n",
    "\n",
    "Does it seem as though three clusters may be sufficient after quickly taking a look at the differentiations amongst these key variables? How can one be confident?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elbow method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *Elbow method* can be used as one input in determining 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*, the sum of squared errors (`SSE`) is calculated after the model converges for different values of *k*. All the resulting `SSE` values are plotted by increasing value of *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(students_features, k)$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 the SSE decreases as k increases. Here, the SSE decreases faster when k is small. As k increases, the reduction in SSE becomes smaller. 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 this graph, the elbow curve begins to flatten around k = 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try running the model with 4 clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(1)\n",
    "k4 <- kmeans(students_features, centers = 4)\n",
    "k4$size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save these results to a dataframe called `frame_4`, and check the characteristics of students in each cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame_4 <- data.frame(students_ml, k4$cluster)  # add cluster number to the original dataframe\n",
    "\n",
    "frame_4 %>% select(-drf_id) %>%\n",
    "    group_by(k4.cluster) %>%\n",
    "    summarize_all(\"mean\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which clustering results - `frame_3` or `frame_4` - seem to be more useful for the data interpretation? Is there any additional information that should be considered?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In clustering, there is often no single right answer - everytime a different number of clusters is used, interesting patterns about the data can be exposed. However, it is crucial to think about and decide if whether the clusters represent true subgroups in the 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": [
    "### Adding Categorical Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clustering can also be performed using categorical variables. There are many ways to approach using a mix of categorical and numerical variables. These methods are all centered around the idea that there must be a distance measure that can incorporate both numerical and categorical variables. This can be done through either alternative distance metrics or through using categorical variables as numerical variables. \n",
    "\n",
    "In this example, we will go with the latter and the categorical variables will need to be converted to binary (0, 1) values and then scaled with the numerical variables. Note that this makes a pretty big assumption about how we should think about differences in categorical variables and differences in numerical variables. That is, we are assuming the maximum difference in any given numerical variable is the same as a difference on any individual indicator variable. Depending on the topic, this may or may not make sense. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Modal funder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read into R\n",
    "qry <- \"\n",
    "select sed.drf_id, team_size, any_non_federal, sem.modal_funder\n",
    "from ds_nsf_ncses.dbo.nsf_sed sed\n",
    "join tr_ncses_2021.dbo.sed_umetrics_xwalk xwalk\n",
    "on sed.drf_id = xwalk.drf_id\n",
    "join ds_iris_umetrics.dbo.semester sem\n",
    "on xwalk.emp_number = sem.emp_number\n",
    "where sed.phdfy = '2015'\n",
    "\"\n",
    "students <- dbGetQuery(con, qry)\n",
    "\n",
    "# see the dataframe\n",
    "head(students)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the variables with the number of semesters of funding, federal and non-federal funding, like described at the beginning of the notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "students <- students %>%\n",
    "    group_by(drf_id) %>%\n",
    "    mutate(sem_funding = n()) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If team_size > 0, then flag as federal (add 1), if team_size is missing, then add 1 in the \"any_non_federal\" variable\n",
    "students <- students %>%\n",
    "    mutate(federal = ifelse(team_size > 0, 1, 0),\n",
    "           non_federal = ifelse(is.na(team_size) | any_non_federal == 1, 1, 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "students <- students %>%\n",
    "            group_by(drf_id) %>%                                 # for each unique drf_id\n",
    "            mutate(sem_federal = sum(federal, na.rm = TRUE),     # sum values in the federal flag\n",
    "                   sem_non_federal = sum(non_federal)) %>%       # sum values in the any_non_federal flag\n",
    "            select(-c(team_size, any_non_federal, federal, non_federal)) %>%  # remove these columns\n",
    "            slice(1) %>% ungroup()                                         # deduplicate rows, choose just 1 row per unique drf_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Segment the data frame by 5 major agencies (NIH, NSF, DOE, DOD, USDA), otherwise overwriting **modal_funder** as **\"OTHER\"**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "students <- students %>%\n",
    "    mutate(modal_funder = ifelse(modal_funder %in% c(\"NIH\", \"NSF\", \"DOE\", \"DOD\", \"USDA\"), modal_funder, \"OTHER\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "students <- students %>%\n",
    "            replace(is.na(.), 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modal_funder_df <- students %>%\n",
    "                    select(c(drf_id, modal_funder))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use dcast function to get create 0 and 1 values for the primary source of support variable\n",
    "modal_funder_df <- dcast(data = modal_funder_df, drf_id ~ modal_funder, length)\n",
    "\n",
    "merged <- inner_join(students, modal_funder_df, by=c('drf_id'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(merged)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged <- merged %>%\n",
    "            select(-c(drf_id, modal_funder))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the head of the table with modal funder coded as 0 and 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(merged)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scale the features\n",
    "merged_scale <- scale(merged)\n",
    "\n",
    "# View first rows after scaling\n",
    "head(merged_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove all missing data points before running clustering\n",
    "# na.omit will remove any rows with any NA values\n",
    "merged_scale <- na.omit(merged_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_scale <- as.data.frame(merged_scale)"
   ]
  },
  {
   "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(merged_scale, k)$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": [
    "Based on the Elbow method, it looks like 7 clusters are optimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the model \n",
    "set.seed(1)\n",
    "k7 <- kmeans(merged_scale, centers=7, nstart=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged <- na.omit(merged)                     # remove missing values\n",
    "frame_7 <- data.frame(merged, k7$cluster)  # add cluster number to the original dataframe \n",
    "\n",
    "frame_7 %>% \n",
    "    group_by(k7.cluster) %>%\n",
    "    summarize_all(\"mean\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How might you interpret this? Can you identify clusters of interest? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Close the database connection\n",
    "dbDisconnect(con)"
   ]
  },
  {
   "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"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}