Cluster Analysis in R

1 Calculating distance between observations

1.1 What is cluster analysis?

1.1.1 When to cluster?

Example of scenarios where clustering methods would likely be appropriate:

• Using consumer behavior data to identify distinct segments within a market.
• Identifying distinct groups of stocks that follow similar trading patterns.

1.2 Distance between two observations

1.2.1 Calculate & plot the distance between two players

We will plot the positions of the 2 players and manually calculate the distance between them by using the Euclidean distance formula.

two_players <- as.data.frame(list(x = c(5, 15), y = c(4, 10)), .Names = c("x", "y"))

two_players

# Plot the positions of the players
ggplot(two_players, aes(x = x, y = y)) + 
  geom_point() +
  # Assuming a 40x60 field
  lims(x = c(-30,30), y = c(-20, 20))

# Split the players data frame into two observations
player1 <- two_players[1, ]
player2 <- two_players[2, ]

Calculate the distance between player1 and player2 by using the Euclidean distance formula

\(\sqrt{(x1−x2)^2+(y1−y2)^2}\)

# Calculate and print their distance using the Euclidean Distance formula
player_distance <- sqrt( (player1$x - player2$x)^2 + (player1$y - player2$y)^2 )
player_distance

[1] 11.6619

1.2.2 Using the dist() function

Using the Euclidean formula manually may be practical for 2 observations but can get more complicated rather quickly when measuring the distance between many observations.

The dist() function simplifies this process by calculating distances between our observations (rows) using their features (columns). In this case the observations are the player positions and the dimensions are their x and y coordinates.

The default distance calculation for the dist() function is Euclidean distance

three_players <- as.data.frame(list(x = c(5, 15, 0), y = c(4, 10, 20)), .Names = c("x", 
"y"))

# Calculate the Distance Between two_players
dist_two_players <- dist(two_players)
dist_two_players

        1
2 11.6619

# Calculate the Distance Between three_players
dist_three_players <- dist(three_players)
dist_three_players

         1        2
2 11.66190         
3 16.76305 18.02776

The dist() function makes life easier when working with many dimensions and observations.

1.2.3 Who are the closest players?

Which two players are closest to one another?

four_players <- as.data.frame(list(x = c(5, 15, 0, -5), y = c(4, 10, 20, 5)), .Names = c("x", 
"y"))

(dist(four_players))

         1        2        3
2 11.66190                  
3 16.76305 18.02776         
4 10.04988 20.61553 15.81139

Players 1 and 4 are the closest to one another.

1.3 The importance of scale

1.3.1 Effects of scale

When a variable is on a larger scale than other variables in your data it may disproportionately influence the resulting distance calculated between your observations. Lets see this in action by observing a sample of data from the trees data set.

trees <- as.data.frame(list(Girth = c(8.3, 8.6, 8.8, 10.5, 10.7, 10.8, 11, 
11, 11.1, 11.2, 11.3, 11.4, 11.4, 11.7, 12, 12.9, 12.9, 13.3, 
13.7, 13.8, 14, 14.2, 14.5, 16, 16.3, 17.3, 17.5, 17.9, 18, 18, 
20.6), Height = c(70, 65, 63, 72, 81, 83, 66, 75, 80, 75, 79, 
76, 76, 69, 75, 74, 85, 86, 71, 64, 78, 80, 74, 72, 77, 81, 82, 
80, 80, 80, 87), Volume = c(10.3, 10.3, 10.2, 16.4, 18.8, 19.7, 
15.6, 18.2, 22.6, 19.9, 24.2, 21, 21.4, 21.3, 19.1, 22.2, 33.8, 
27.4, 25.7, 24.9, 34.5, 31.7, 36.3, 38.3, 42.6, 55.4, 55.7, 58.3, 
51.5, 51, 77)), .Names = c("Girth", "Height", "Volume"))

We will leverage the scale() function which by default centers & scales our column features.

Our variables are the following:

• Girth - tree diameter in inches
• Height - tree height in inches

three_trees <- trees[1:3, 1:2]

# Calculate distance for three_trees 
dist_trees <- dist(three_trees)

# Scale three trees & calculate the distance  
scaled_three_trees <- scale(three_trees)
dist_scaled_trees <- dist(scaled_three_trees)

# Output the results of both Matrices
print('Without Scaling')

[1] "Without Scaling"

dist_trees

         1        2
2 5.008992         
3 7.017834 2.009975

print('With Scaling')

[1] "With Scaling"

dist_scaled_trees

          1         2
2 1.8286961          
3 2.7778767 0.9691601

1.4 Measuring distance for categorical data

1.4.1 Calculating distance between categorical variables

We will explore how to calculate binary (Jaccard) distances. In order to calculate distances we will first have to dummify our categories using the dummy.data.frame() from the library dummies.

We will use a small collection of survey observations stored in the data frame job_survey with the following columns:

• job_satisfaction Possible options: “Hi”, “Mid”, “Low”
• is_happy Possible options: “Yes”, “No”

job_survey <- as.data.frame(list(job_satisfaction = structure(c(2L, 2L, 1L, 2L, 
3L), .Label = c("Hi", "Low", "Mid"), class = "factor"), is_happy = structure(c(1L, 
1L, 2L, 1L, 1L), .Label = c("No", "Yes"), class = "factor")), .Names = c("job_satisfaction", 
"is_happy"))

library(dummies)

# Dummify the Survey Data
dummy_survey <- dummy.data.frame(job_survey)

non-list contrasts argument ignorednon-list contrasts argument ignored

# Calculate the Distance
dist_survey <- dist(dummy_survey, method = "binary")

# Print the Original Data
job_survey

# Print the Distance Matrix
dist_survey

          1         2         3         4
2 0.0000000                              
3 1.0000000 1.0000000                    
4 0.0000000 0.0000000 1.0000000          
5 0.6666667 0.6666667 1.0000000 0.6666667

Notice that this distance metric successfully captured that observations 1 and 2 are identical (distance of 0)

2 Hierarchical clustering

2.1 Comparing more than two observations

2.1.1 Calculating linkage

Let us revisit the example with three players on a field. The distance matrix between these three players is shown below and is available as the variable dist_players.

(dist_players <- dist_three_players)

         1        2
2 11.66190         
3 16.76305 18.02776

From this we can tell that the first group that forms is between players 1 & 2, since they are the closest to one another with a Euclidean distance value of 11.

Now you want to apply the three linkage methods you have learned to determine what the distance of this group is to player 3.

# Extract the pair distances
distance_1_2 <- dist_players[1]
distance_1_3 <- dist_players[2]
distance_2_3 <- dist_players[3]

# Extract the pair distances
distance_1_2 <- dist_players[1]
distance_1_3 <- dist_players[2]
distance_2_3 <- dist_players[3]

# Calculate the complete distance between group 1-2 and 3
complete <- max(c(distance_1_3, distance_2_3))
complete

[1] 18.02776

# Calculate the single distance between group 1-2 and 3
single <- min(c(distance_1_3, distance_2_3))
single

[1] 16.76305

# Calculate the average distance between group 1-2 and 3
average <- mean(c(distance_1_3, distance_2_3))
average

[1] 17.39541

2.2 Capturing K clusters

2.2.1 Assign cluster membership

We will leverage the hclust() function to calculate the iterative linkage steps and we will use the cutree() function to extract the cluster assignments for the desired number (k) of clusters.

The positions of 12 players at the start of a 6v6 soccer match. This is stored in the lineup data frame.

lineup <- rbind(structure(list(x = c(-1, -2, 8, 7, -12, -15, -13, 15, 21, 12, 
-25, 26), y = c(1, -3, 6, -8, 8, 0, -10, 16, 2, -15, 1, 0)), .Names = c("x", 
"y"), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
-12L)))

We know that this match has two teams (k = 2), let’s use the clustering methods to assign which team each player belongs in based on their position.

Notes:

• The linkage method can be passed via the method parameter: hclust(distance_matrix, method = “complete”)
• Remember that in soccer opposing teams start on their half of the field.
• Because these positions are measured using the same scale we do not need to re-scale our data.

# Calculate the Distance
dist_players <- dist(lineup)

# Perform the hierarchical clustering using the complete linkage
hc_players <- hclust(dist_players, method = "complete")

# Calculate the assignment vector with a k of 2
clusters_k2 <- cutree(hc_players, k = 2)

# Create a new data frame storing these results
lineup_k2_complete <- mutate(lineup, cluster = clusters_k2)

2.2.2 Exploring the clusters

Because clustering analysis is always in part qualitative, it is incredibly important to have the necessary tools to explore the results of the clustering.

We will explore that data frame just created.

library(dplyr)

# Count the cluster assignments
count(lineup_k2_complete, cluster)

# Plot the positions of the players and color them using their cluster
ggplot(lineup_k2_complete, aes(x = x, y = y, color = factor(cluster))) +
  geom_point()

2.3 Vizualizing the dendrogram

2.3.1 Comparing average, single & complete linkage

We are now ready to analyze the clustering results of the lineup dataset using the dendrogram plot. This will give you a new perspective on the effect the decision of the linkage method has on your resulting cluster analysis.

# Prepare the Distance Matrix
dist_players <- dist(lineup)

# Generate hclust for complete, single & average linkage methods
hc_complete <- hclust(dist_players, method = "complete")
hc_single <- hclust(dist_players, method = "single")
hc_average <- hclust(dist_players, method = "average")

# Plot & Label the 3 Dendrograms Side-by-Side
# Hint: To see these Side-by-Side run the 4 lines together as one command
par(mfrow = c(1,3))
plot(hc_complete, main = 'Complete Linkage')
plot(hc_single, main = 'Single Linkage')

plot(hc_average, main = 'Average Linkage')

2.3.2 Height of the tree

An advantage of working with a clustering method like hierarchical clustering is that you can describe the relationships between your observations based on both the distance metric and the linkage metric selected (the combination of which defines the height of the tree).

2.4 Cutting the tree

2.4.1 Clusters based on height

We will leverage the visual representation of the dendrogram in order to group your observations into clusters using a maximum height (h), below which clusters form.

We will work the color_branches() function from the dendextend library in order to visually inspect the clusters that form at any height along the dendrogram.

library(dendextend)

library(dendextend)
dist_players <- dist(lineup, method = 'euclidean')
hc_players <- hclust(dist_players, method = "complete")

# Create a dendrogram object from the hclust variable
dend_players <- as.dendrogram(hc_players)

# Plot the dendrogram
plot(dend_players)

# Color branches by cluster formed from the cut at a height of 20 & plot
dend_20 <- color_branches(dend_players, h = 20)

# Plot the dendrogram with clusters colored below height 20
plot(dend_20)

# Color branches by cluster formed from the cut at a height of 40 & plot
dend_40 <- color_branches(dend_players, h = 40)

# Plot the dendrogram with clusters colored below height 40
plot(dend_40)

The height that we use to cut the tree greatly influences the number of clusters and their size

2.4.2 Exploring the branches cut from the tree

The cutree() function can also be used to cut a tree at a given height by using the h parameter.

dist_players <- dist(lineup, method = 'euclidean')
hc_players <- hclust(dist_players, method = "complete")

# Calculate the assignment vector with a h of 20
clusters_h20 <- cutree(hc_players, h = 20)

# Create a new data frame storing these results
lineup_h20_complete <- mutate(lineup, cluster = clusters_h20)

# Calculate the assignment vector with a h of 40
clusters_h40 <- cutree(hc_players, h = 40)

# Create a new data frame storing these results
lineup_h40_complete <- mutate(lineup, cluster = clusters_h40)

# Plot the positions of the players and color them using their cluster for height = 20
ggplot(lineup_h20_complete, aes(x = x, y = y, color = factor(cluster))) +
  geom_point()

# Plot the positions of the players and color them using their cluster for height = 40
ggplot(lineup_h40_complete, aes(x = x, y = y, color = factor(cluster))) +
  geom_point()

The height of any branch is determined by the linkage and distance decisions (in this case complete linkage and Euclidean distance). While the members of the clusters that form below a desired height have a maximum linkage+distance amongst themselves that is less than the desired height.

2.5 Making sense of the clusters

2.5.1 Segment wholesale customers

We’re now ready to use hierarchical clustering to perform market segmentation (i.e. use consumer characteristics to group them into subgroups).

customers_spend <- readRDS("ws_customers.rds")

We are provided with the amount spent by 45 different clients of a wholesale distributor for the food categories of Milk, Grocery & Frozen. This is stored in the data frame customers_spend. Assign these clients into meaningful clusters.

We can assume, because the data is all of the same type (amount spent) and you will not need to scale it.

# Calculate Euclidean distance between customers
dist_customers <- dist(customers_spend)

# Generate a complete linkage analysis 
hc_customers <- hclust(dist_customers, method = "complete")

# Plot the dendrogram
plot(hc_customers)

# Create a cluster assignment vector at h = 15000
clust_customers <- cutree(hc_customers, h = 15000)

# Generate the segmented customers data frame
segment_customers <- mutate(customers_spend, cluster = clust_customers)

2.5.2 Explore wholesale customer clusters

We are now ready to analyze the characteristics of these clusters.

Since we are working with more than 2 dimensions it would be challenging to visualize a scatter plot of the clusters, instead we will rely on summary statistics to explore these clusters. We will analyze the mean amount spent in each cluster for all three categories.

dist_customers <- dist(customers_spend)
hc_customers <- hclust(dist_customers)
clust_customers <- cutree(hc_customers, h = 15000)
segment_customers <- mutate(customers_spend, cluster = clust_customers)

# Count the number of customers that fall into each cluster
count(segment_customers, cluster)

# Color the dendrogram based on the height cutoff
dend_customers <- as.dendrogram(hc_customers)
dend_colored <- color_branches(dend_customers, h = 15000)

# Plot the colored dendrogram
plot(dend_colored)

# Calculate the mean for each category
segment_customers %>% 
  group_by(cluster) %>% 
  summarise_all(list(mean))

• Customers in cluster 1 spent more money on Milk than any other cluster.
• Customers in cluster 3 spent more money on Grocery than any other cluster.
• Customers in cluster 4 spent more money on Frozen goods than any other cluster.
• The majority of customers fell into cluster 2 and did not show any excessive spending in any category. press

3 K-means clustering

3.1 Introduction

3.1.1 K-means on a soccer field

We will use the lineup dataframe.

Note that in the kmeans() function k is specified using the centers parameter.

# Build a kmeans model
model_km2 <- kmeans(lineup, centers = 2)

# Extract the cluster assignment vector from the kmeans model
clust_km2 <- model_km2$cluster

# Create a new data frame appending the cluster assignment
lineup_km2 <- mutate(lineup, cluster = clust_km2)

# Plot the positions of the players and color them using their cluster
ggplot(lineup_km2, aes(x = x, y = y, color = factor(cluster))) +
  geom_point()

Knowing the desired number of clusters ahead of time can be very helpful when performing a k-means analysis.

3.1.2 K-means on a soccer field (part 2)

We successfully used the k-means algorithm to cluster the two teams from the lineup data frame. This time, let’s explore what happens when you use a k of 3.

# Build a kmeans model
model_km3 <- kmeans(lineup, centers = 3)

# Extract the cluster assignment vector from the kmeans model
clust_km3 <- model_km3$cluster

# Create a new data frame appending the cluster assignment
lineup_km3 <- mutate(lineup, cluster = clust_km3)

# Plot the positions of the players and color them using their cluster
ggplot(lineup_km3, aes(x = x, y = y, color = factor(cluster))) +
  geom_point()

Does this result make sense? Remember we only have 2 teams on the field. It’s very important to remember that k-means will run with any k that is more than 2 and less than your total observations, but it doesn’t always mean the results will be meaningful.

3.2 Evaluating different values of K by eye

3.2.1 Many K’s many models

While the lineup dataset clearly has a known value of k, often times the optimal number of clusters isn’t known and must be estimated.

We will leverage map_dbl() from the purrr library to run k-means using values of k ranging from 1 to 10 and extract the total within-cluster sum of squares metric from each one. This will be the first step towards visualizing the elbow plot.

library(purrr)

# Use map_dbl to run many models with varying value of k (centers)
tot_withinss <- map_dbl(1:10,  function(k){
  model <- kmeans(x = lineup, centers = k)
  model$tot.withinss
})

# Generate a data frame containing both k and tot_withinss
elbow_df <- data.frame(
  k = 1:10 ,
  tot_withinss = tot_withinss
)

3.2.2 Elbow (Scree) plot

We have calculated the total within-cluster sum of squares for values of k ranging from 1 to 10. We can visualize this relationship using a line plot to create what is known as an elbow plot (or scree plot).

When looking at an elbow plot we want to see a sharp decline from one k to another followed by a more gradual decrease in slope. The last value of k before the slope of the plot levels off suggests a “good” value of k.

# Use map_dbl to run many models with varying value of k (centers)
tot_withinss <- map_dbl(1:10,  function(k){
  model <- kmeans(x = lineup, centers = k)
  model$tot.withinss
})

# Generate a data frame containing both k and tot_withinss
elbow_df <- data.frame(
  k = 1:10,
  tot_withinss = tot_withinss
)

# Plot the elbow plot
ggplot(elbow_df, aes(x = k, y = tot_withinss)) +
  geom_line() +
  scale_x_continuous(breaks = 1:10)

We have learned how to create and visualize elbow plots as a tool for finding a “good” value of k.

We can see that there is a sharp change in the slope of this line that makes an "elbow" shape. Furthermore, this is supported by the prior knowledge that there are two teams in this data and a k of 2 is desired.

3.3 Silhouette analysis: observation level performance

3.3.1 Silhouette analysis

Silhouette analysis allows you to calculate how similar each observations is with the cluster it is assigned relative to other clusters. This metric (silhouette width) ranges from -1 to 1 for each observation in your data and can be interpreted as follows:

• Values close to 1 suggest that the observation is well matched to the assigned cluster
• Values close to 0 suggest that the observation is borderline matched between two clusters
• Values close to -1 suggest that the observations may be assigned to the wrong cluster

We will leverage the pam() and the silhouette() functions from the cluster library to perform silhouette analysis to compare the results of models with a k of 2 and a k of 3.

library(cluster)

# Generate a k-means model using the pam() function with a k = 2
pam_k2 <- pam(lineup, k = 2)

# Plot the silhouette visual for the pam_k2 model
plot(silhouette(pam_k2))

# Generate a k-means model using the pam() function with a k = 3
pam_k3 <- pam(lineup, k = 3)

# Plot the silhouette visual for the pam_k3 model
plot(silhouette(pam_k3))

Did you notice that for k = 2, no observation has a silhouette width close to 0? What about the fact that for k = 3, observation 3 is close to 0 and is negative? This suggests that k = 3 is not the right number of clusters.

3.4 Making sense of the K-means clusters

3.4.1 Revisiting wholesale data: “Best” k

We explored wholesale distributor data customers_spend using hierarchical clustering. This time we will analyze this data using the k-means clustering tools. The first step will be to determine the “best” value of k using average silhouette width.

# Use map_dbl to run many models with varying value of k
sil_width <- map_dbl(2:10,  function(k){
  model <- pam(x = customers_spend, k = k)
  model$silinfo$avg.width
})

# Generate a data frame containing both k and sil_width
sil_df <- data.frame(
  k = 2:10,
  sil_width = sil_width
)

# Plot the relationship between k and sil_width
ggplot(sil_df, aes(x = k, y = sil_width)) +
  geom_line() +
  scale_x_continuous(breaks = 2:10)

From the plot notice that k = 2 has the highest average sillhouette width and is the “best” value of k we will move forward with.

3.4.2 Exploration

We have found that k = 2 has the highest average silhouette width. We will continue to analyze the wholesale customer data by building and exploring a kmeans model with 2 clusters.

set.seed(42)

# Build a k-means model for the customers_spend with a k of 2
model_customers <- kmeans(customers_spend, centers = 2)

# Extract the vector of cluster assignments from the model
clust_customers <- model_customers$cluster

# Build the segment_customers data frame
segment_customers <- mutate(customers_spend, cluster = clust_customers)

# Calculate the size of each cluster
count(segment_customers, cluster)

# Calculate the mean for each category
segment_customers %>% 
  group_by(cluster) %>% 
  summarise_all(list(mean))

It seems that in this case cluster 1 consists of individuals who proportionally spend more on Frozen food while cluster 2 customers spent more on Milk and Grocery.

When we explored this data using hierarchical clustering, the method resulted in 4 clusters while using k-means we got 2. Both of these results are valid, but which one is appropriate for this would require more subject matter expertise.

“Generating clusters is a science, but interpreting them is an art.”