## 11.1 K-Means

The k-means clustering algorithm randomly assigns all observations to one of $$k$$ clusters. K-means then iteratively calculates the cluster centroids and reassigns the observations to their nearest centroid. The centroid is the mean of the points in the cluster (Hence the name “k-means”). The iterations continue until either the centroids stabilize or the iterations reach a set maximum, iter.max (typically 50). The result is k clusters with the minimum total intra-cluster variation.

The centroid of cluster $$c_i \in C$$ is the mean of the cluster observations $$S_i$$: $$c_i = \frac{1}{|S_i|} \sum_{x_i \in S_i}{x_i}$$. The nearest centroid is the minimum squared euclidean distance, $$\underset{c_i \in C}{\operatorname{arg min}} dist(c_i, x)^2$$. A more robust version of k-means is partitioning around medoids (pam), which minimizes the sum of dissimilarities instead of a sum of squared euclidean distances. That’s what I’ll use.

The algorithm will converge to a result, but the result may only be a local optimum. Other random starting centroids may yield a different local optimum. Common practice is to run the k-means algorithm nstart times and select the lowest within-cluster sum of squared distances among the cluster members. A typical number of runs is nstart = 20.

#### Choosing K

What is the best number of clusters? You may have a preference in advance, but more likely you will use a scree plot or use the silhouette method. The scree plot is a plot of the total within-cluster sum of squared distances as a function of k. The sum of squares always decreases as k increases, but at a declining rate. The optimal k is at the “elbow” in the curve - the point at which the curve flattens. kmeans() returns an object of class kmeans, a list in which one of the components is the model sum of squares tot.withinss. In the scree plot below, the elbow may be k = 5.

wss <- map_dbl(2:10, ~ kmeans(dat_2_gwr, centers = .x)$tot.withinss) wss_tbl <- tibble(k = 2:10, wss) ggplot(wss_tbl, aes(x = k, y = wss)) + geom_point(size = 2) + geom_line() + scale_x_continuous(breaks = 2:10) + labs(title = "Scree Plot") The silhouette method calculates the within-cluster distance $$C(i)$$ for each observation, and its distance to the nearest cluster $$N(i)$$. The silhouette width is $$S = 1 - C(i) / N(i)$$ for $$C(i) < N(i)$$ and $$S = N(i) / C(i) - 1$$ for $$C(i) > N(i)$$. A value close to 1 means the observation is well-matched to its current cluster; A value near 0 means the observation is on the border between the two clusters; and a value near -1 means the observation is better-matched to the other cluster. The optimal number of clusters is the number that maximizes the total silhouette width. cluster::pam() returns a list in which one of the components is the average width silinfo$avg.width. In the silhouette plot below, the maximum width is at k = 6.

sil <- map_dbl(2:10, ~ pam(dat_2_gwr, k = .x)$silinfo$avg.width)
sil_tbl <- tibble(k = 2:10, sil)
ggplot(sil_tbl, aes(x = k, y = sil)) +
geom_point(size = 2) +
geom_line() +
scale_x_continuous(breaks = 2:10) +
labs(title = "Silhouette Plot") #### Summarize Results

Run pam() again and attach the results to the original table for visualization and summary statistics.

mdl <- pam(dat_2_gwr, k = 6)
dat_3 <- dat_2 %>% mutate(cluster = as.factor(mdl$clustering)) Here are the six medoids from our data set. dat_2[mdl$medoids, ] %>%
t() %>% as.data.frame() %>% rownames_to_column() %>%
flextable::flextable() %>% flextable::autofit()
 rowname V1 V2 V3 V4 V5 V6 EmployeeNumber 1171 35 65 221 747 1408 Attrition No No Yes No No No OverTime No No Yes No No No JobLevel 2 2 1 1 2 4 MonthlyIncome 5155 6825 3441 2713 5304 16799 YearsAtCompany 6 9 2 5 8 20 StockOptionLevel 0 1 0 1 1 1 YearsWithCurrManager 5 2 2 2 7 9 TotalWorkingYears 9 10 2 5 10 21 MaritalStatus Single Married Single Married Divorced Married Age 42 42 28 28 30 42 YearsInCurrentRole 4 7 2 2 7 7 JobRole Sales Executive Sales Executive Laboratory Technician Research Scientist Sales Executive Manager EnvironmentSatisfaction 2 3 3 3 3 3 JobInvolvement 3 3 3 3 3 3 BusinessTravel Travel_Rarely Travel_Rarely Travel_Rarely Travel_Rarely Travel_Rarely Travel_Rarely

We’re most concerned about attrition. Do high-attrition employees fall into a particular cluster? Yes! 79.7% of cluster 3 left the company - that’s 59.5% of all turnover in the company.

dat_3_smry <- dat_3 %>%
count(cluster, Attrition) %>%
group_by(cluster) %>%
mutate(cluster_n = sum(n),
turnover_rate = scales::percent(n / sum(n), accuracy = 0.1)) %>%
ungroup() %>%
filter(Attrition == "Yes") %>%
mutate(pct_of_turnover = scales::percent(n / sum(n), accuracy = 0.1)) %>%
select(cluster, cluster_n, turnover_n = n, turnover_rate, pct_of_turnover)
dat_3_smry %>% flextable::flextable() %>% flextable::autofit()
 cluster cluster_n turnover_n turnover_rate pct_of_turnover 1 268 23 8.6% 9.7% 2 280 24 8.6% 10.1% 3 177 141 79.7% 59.5% 4 364 29 8.0% 12.2% 5 202 8 4.0% 3.4% 6 179 12 6.7% 5.1%

You can get some sense of the quality of clustering by constructing the Barnes-Hut t-Distributed Stochastic Neighbor Embedding (t-SNE).

dat_4 <- dat_3 %>%
left_join(dat_3_smry, by = "cluster") %>%
rename(Cluster = cluster) %>%
mutate(
MonthlyIncome = MonthlyIncome %>% scales::dollar(),
description = str_glue("Turnover = {Attrition}
MaritalDesc = {MaritalStatus}
Age = {Age}
Job Role = {JobRole}
Job Level {JobLevel}
Overtime = {OverTime}
Current Role Tenure = {YearsInCurrentRole}
Professional Tenure = {TotalWorkingYears}
Monthly Income = {MonthlyIncome}

Cluster: {Cluster}
Cluster Size: {cluster_n}
Cluster Turnover Rate: {turnover_rate}
Cluster Turnover Count: {turnover_n}
"))

tsne_obj <- Rtsne(dat_2_gwr, is_distance = TRUE)

tsne_tbl <- tsne_obj$Y %>% as_tibble() %>% setNames(c("X", "Y")) %>% bind_cols(dat_4) %>% mutate(Cluster = as_factor(Cluster)) g <- tsne_tbl %>% ggplot(aes(x = X, y = Y, colour = Cluster, text = description)) + geom_point() ggplotly(g) Another common approach is to take summary statistics and draw your own conclusions. You might start by asking which attributes differ among the clusters. The box plots below show the distribution of the numeric variables. All of the numeric variable distributions appear to vary among the clusters. my_boxplot <- function(y_var){ dat_3 %>% ggplot(aes(x = cluster, y = !!sym(y_var))) + geom_boxplot() + geom_jitter(aes(color = Attrition), alpha = 0.2, height = 0.10) + theme_minimal() + theme(legend.position = "none") + labs(x = "", y = "", title = y_var) } vars_numeric <- dat_3 %>% select(-EmployeeNumber) %>% select_if(is.numeric) %>% colnames() g <- map(vars_numeric, my_boxplot) gridExtra::marrangeGrob(g, nrow=1, ncol = 2)   You can perform an analysis of variance to confirm. The table below collects the ANOVA results for each of the numeric variables. The results indicate that there are significant differences among clusters at the .01 level for all of the numeric variables. km_aov <- vars_numeric %>% map(~ aov(rlang::eval_tidy(expr(!!sym(.x) ~ cluster)), data = dat_3)) km_aov %>% map(anova) %>% map(~ data.frame(F = .x$F value[], p = .x$Pr(>F)[])) %>% bind_rows() %>% bind_cols(Attribute = vars_numeric) %>% select(Attribute, everything()) %>% flextable::flextable() %>% flextable::autofit() %>% flextable::colformat_num(j = 2, digits = 2) %>% flextable::colformat_num(j = 3, digits = 4)  Attribute F p MonthlyIncome 718.88 0.0000 YearsAtCompany 134.58 0.0000 YearsWithCurrManager 65.53 0.0000 TotalWorkingYears 342.88 0.0000 Age 98.76 0.0000 YearsInCurrentRole 69.53 0.0000 We’re particularly interested in cluster 3, the high attrition cluster, and what sets it apart from the others. Right away, we can see that clusters 3 and 4 similar in almost all of theses attributes. They tend to have the smaller incomes, company and job tenure, years with their current manager, total working experience, and age. I.e., they tend to be lower on the career ladder. To drill into cluster differences to determine which clusters differ from others, use the Tukey HSD post hoc test with Bonferroni method applied to control the experiment-wise error rate. That is, only reject the null hypothesis of equal means among clusters if the p-value is less than $$\alpha / p$$, or $$.05 / 6 = 0.0083$$. The significantly different cluster combinations are shown in bold. Clusters 3 and 4 differ from the others on all six measures. However, there are no significant differences between c3 and c4 (highlighted green) for the numeric variables. km_hsd <- map(km_aov, TukeyHSD) map(km_hsd, ~ .x$cluster %>%
data.frame() %>%
rownames_to_column() %>%
filter(str_detect(rowname, "-"))) %>%
map2(vars_numeric, bind_cols) %>%
bind_rows() %>%
select(predictor = ...6, everything()) %>%
mutate(cluster_a = str_sub(rowname, start = 1, end = 1),
cluster = paste0("c", str_sub(rowname, start = 3, end = 3))) %>%
names_from = cluster_a, values_from = p.adj,
names_prefix = "c") %>%
flextable::flextable() %>%
flextable::colformat_num(j = c(3:7), digits = 4) %>%
flextable::bold(i = ~ c2 < .05 / length(vars_numeric), j = ~ c2, bold = TRUE) %>%
flextable::bold(i = ~ c3 < .05 / length(vars_numeric), j = ~ c3, bold = TRUE) %>%
flextable::bold(i = ~ c4 < .05 / length(vars_numeric), j = ~ c4, bold = TRUE) %>%
flextable::bold(i = ~ c5 < .05 / length(vars_numeric), j = ~ c5, bold = TRUE) %>%
flextable::bold(i = ~ c6 < .05 / length(vars_numeric), j = ~ c6, bold = TRUE) %>%
flextable::bg(i = ~ cluster == "c3", j = ~ c4, bg = "#B6E2D3") %>%
flextable::autofit() 
 predictor cluster c2 c3 c4 c5 c6 MonthlyIncome c1 0.6587 0.0000 0.0000 0.9720 0.0000 MonthlyIncome c2 0.0000 0.0000 0.9899 0.0000 MonthlyIncome c3 0.3975 0.0000 0.0000 MonthlyIncome c4 0.0000 0.0000 MonthlyIncome c5 0.0000 YearsAtCompany c1 0.9999 0.0000 0.0000 0.9922 0.0000 YearsAtCompany c2 0.0000 0.0000 0.9989 0.0000 YearsAtCompany c3 0.9935 0.0000 0.0000 YearsAtCompany c4 0.0000 0.0000 YearsAtCompany c5 0.0000 YearsWithCurrManager c1 1.0000 0.0000 0.0000 0.9763 0.0000 YearsWithCurrManager c2 0.0000 0.0000 0.9595 0.0000 YearsWithCurrManager c3 0.9965 0.0000 0.0000 YearsWithCurrManager c4 0.0000 0.0000 YearsWithCurrManager c5 0.0000 TotalWorkingYears c1 0.9840 0.0000 0.0000 0.9980 0.0000 TotalWorkingYears c2 0.0000 0.0000 0.8934 0.0000 TotalWorkingYears c3 0.9998 0.0000 0.0000 TotalWorkingYears c4 0.0000 0.0000 TotalWorkingYears c5 0.0000 Age c1 0.9980 0.0000 0.0000 0.9989 0.0000 Age c2 0.0000 0.0000 0.9691 0.0000 Age c3 0.0532 0.0000 0.0000 Age c4 0.0000 0.0000 Age c5 0.0000 YearsInCurrentRole c1 0.1013 0.0000 0.0000 0.6471 0.0000 YearsInCurrentRole c2 0.0000 0.0000 0.9573 0.0000 YearsInCurrentRole c3 0.9854 0.0000 0.0000 YearsInCurrentRole c4 0.0000 0.0000 YearsInCurrentRole c5 0.0000

So at this point we have kind of an incomplete picture. We know the high-attrition employees are low on the career ladder, but cluster 4 is also low on the career ladder and they are not high-attrition.

How about the factor variables? The tile plots below show that clusters 3 and 4 are lab technicians and research scientists. They are both at the lowest job level. But three factors distinguish these clusters from each other: cluster 3 is far more likely to work overtime, have no stock options, and be single.

my_tileplot <- function(y_var){
dat_3 %>%
count(cluster, !!sym(y_var)) %>%
ungroup() %>%
group_by(cluster) %>%
mutate(pct = n / sum(n)) %>%
ggplot(aes(y = !!sym(y_var), x = cluster, fill = pct)) +
geom_tile() +
geom_text(aes(label = scales::