K-means Clustering in R Programming with Example

What is Cluster analysis?

Cluster analysis belongs to unsupervised learning. A cluster is a group of data that share similar features. We can say, clustering analysis is more about discovery than a prediction. The machine searches for similarity in the data. For instance, you can use cluster analysis for the following application:

• Customer segmentation: Looks for similarity between groups of customers
• Stock Market clustering: Group stock based on performances
• Reduce dimensionality of a dataset by grouping observations with similar values

Clustering analysis is not too difficult to implement and is meaningful as well as actionable for business.

The most striking difference between supervised and unsupervised learning lies in the results. Unsupervised learning creates a new variable, the label, while supervised learning predicts an outcome. The machine helps the practitioner in the quest to label the data based on close relatedness. It is up to the analyst to make use of the groups and give a name to them.

Let’s make an example to understand the concept of clustering. For simplicity, we work in two dimensions. You have data on the total spend of customers and their ages. To improve advertising, the marketing team wants to send more targeted emails to their customers.

In the following graph, you plot the total spend and the age of the customers.

library(ggplot2)
df <- data.frame(age = c(18, 21, 22, 24, 26, 26, 27, 30, 31, 35, 39, 40, 41, 42, 44, 46, 47, 48, 49, 54),
    spend = c(10, 11, 22, 15, 12, 13, 14, 33, 39, 37, 44, 27, 29, 20, 28, 21, 30, 31, 23, 24)
)
ggplot(df, aes(x = age, y = spend)) +
    geom_point()

A pattern is visible at this point

1. At the bottom-left, you can see young people with a lower purchasing power
2. Upper-middle reflects people with a job that they can afford spend more
3. Finally, older people with a lower budget.

In the figure above, you cluster the observations by hand and define each of the three groups. This example is somewhat straightforward and highly visual. If new observations are appended to the data set, you can label them within the circles. You define the circle based on our judgment. Instead, you can use Machine Learning to group the data objectively.

In this tutorial, you will learn how to use the k-means algorithm.

K-Means Algorithm

K-means is the most widely used clustering method. The algorithm dates back to the 1950s and has been refined many times since.

The algorithm tries to find groups by minimizing the distance between the observations, called local optimal solutions. The distances are measured based on the coordinates of the observations. For instance, in a two-dimensional space the coordinates are simply x and y.

The algorithm works as follow:

• Step 1: Choose k initial centroids at random in the feature space
• Step 2: Assign every observation to its nearest cluster centre, the centroid. This produces k groups
• Step 3: Shift the initial centroid to the mean of the coordinates within a group.
• Step 4: Minimize the distance according to the new centroids. New boundaries are created. Thus, observations will move from one group to another
• Repeat until no observation changes groups

K-means normally measures the Euclidean distance between two observations x and y:

Other measures are available, such as the Manhattan or Minkowski distance. Note that, K-mean returns different groups each time you run the algorithm. Recall that the first initial guesses are random and compute the distances until the algorithm reaches a homogeneity within groups. That is, k-mean is very sensitive to the first choice, and unless the number of observations and groups are small, it is almost impossible to get the same clustering.

Selecting the number of clusters

Another difficulty found with k-mean is the choice of the number of clusters. You can set a high value of k, meaning a large number of groups, to improve within-group homogeneity, but you risk overfitting the data. Overfitting means performance drops sharply on new data, because the model has memorised the noise in this particular sample rather than the underlying pattern.

The number of clusters depends on the nature of the data set, the industry, business and so on. However, there is a rule of thumb to select the appropriate number of clusters:

where n is the number of observations in the dataset.

In practice it is worth spending time searching for the value of k that best fits the business need.

We will use the Prices of Personal Computers dataset to perform our clustering analysis. This dataset contains 6259 observations and 10 features. The dataset observes the price from 1993 to 1995 of 486 personal computers in the US. The variables are price, speed, ram, screen, cd among other.

You will proceed as follow:

• Import data
• Train the model
• Evaluate the model

Import data

K-means is not suitable for factor variables, because it works on distances and discrete labels have no meaningful distance between them. Drop the three categorical columns (cd, multi, premium) along with the row index X. This dataset contains no missing values.

library(dplyr)
PATH <-"https://raw.githubusercontent.com/guru99-edu/R-Programming/master/computers.csv"
df <- read.csv(PATH) %>%
select(-c(X, cd, multi, premium))
glimpse(df)

Output

## Observations: 6,259
## Variables: 7
## $ price <int> 1499, 1795, 1595, 1849, 3295, 3695, 1720, 1995, 2225, 2...
## $ speed <int> 25, 33, 25, 25, 33, 66, 25, 50, 50, 50, 33, 66, 50, 25, ...
## $ hd <int> 80, 85, 170, 170, 340, 340, 170, 85, 210, 210, 170, 210...
## $ ram <int> 4, 2, 4, 8, 16, 16, 4, 2, 8, 4, 8, 8, 4, 8, 8, 4, 2, 4, ...
## $ screen <int> 14, 14, 15, 14, 14, 14, 14, 14, 14, 15, 15, 14, 14, 14, ...
## $ ads <int> 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, ...
## $ trend  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...

From the summary statistics, you can see the data has large values. A good practice before any distance-based method is to standardise the data so that every variable has a mean of zero and a standard deviation of one. Without it, price, whose values run into the thousands, would dominate the distance calculation and screen size would count for almost nothing.

summary(df)

Output:

##      price          speed              hd              ram        
##  Min.   : 949   Min.   : 25.00   Min.   :  80.0   Min.   : 2.000  
##  1st Qu.:1794   1st Qu.: 33.00   1st Qu.: 214.0   1st Qu.: 4.000
##  Median :2144   Median : 50.00   Median : 340.0   Median : 8.000  
##  Mean   :2220   Mean   : 52.01   Mean   : 416.6   Mean   : 8.287  
##  3rd Qu.:2595   3rd Qu.: 66.00   3rd Qu.: 528.0   3rd Qu.: 8.000  
##  Max.   :5399   Max.   :100.00   Max.   :2100.0   Max.   :32.000  
##      screen           ads            trend      
##  Min.   :14.00   Min.   : 39.0   Min.   : 1.00  
##  1st Qu.:14.00   1st Qu.:162.5   1st Qu.:10.00  
##  Median :14.00   Median :246.0   Median :16.00  
##  Mean   :14.61   Mean   :221.3   Mean   :15.93  
##  3rd Qu.:15.00   3rd Qu.:275.0   3rd Qu.:21.50  
##  Max.   :17.00   Max.   :339.0   Max.   :35.00

You rescale the variables with the scale() function of the dplyr library. The transformation reduces the impact of outliers and allows to compare a sole observation against the mean. If a standardized value (or z-score) is high, you can be confident that this observation is indeed above the mean (a large z-score implies that this point is far away from the mean in term of standard deviation. A z-score of two indicates the value is 2 standard deviations away from the mean. Note, the z-score follows a Gaussian distribution and is symmetrical around the mean.

# Note: speed is deliberately left out of the scaled data frame
rescale_df <- df %>%
    mutate(price_scal = scale(price),
    hd_scal = scale(hd),
    ram_scal = scale(ram),
    screen_scal = scale(screen),
    ads_scal = scale(ads),
    trend_scal = scale(trend)) %>%
select(-c(price, speed, hd, ram, screen, ads, trend))

R base has a function to run the k mean algorithm. The basic function of k mean is:

kmeans(df, k)
arguments:
-df: dataset used to run the algorithm
-k: Number of clusters

Train the model

In figure three, you detailed how the algorithm works. You can watch each step graphically with the animation package written by Yihui Xie, who also created knitr for R Markdown. The package is not in the conda channels, so install it from CRAN:

install.packages("animation")

After you load the library, you add .ani after kmeans and R will plot all the steps. For illustration purpose, you only run the algorithm with the rescaled variables hd and ram with three clusters.

set.seed(2345)
library(animation)
kmeans.ani(rescale_df[2:3], 3)

Code Explanation

• kmeans.ani(rescale_df[2:3], 3): Select the columns 2 and 3 of rescale_df data set and run the algorithm with k sets to 3. Plot the animation.

You can interpret the animation as follow:

• Step 1: R randomly chooses three points
• Step 2: Compute the Euclidean distance and draw the clusters. You have one cluster in green at the bottom left, one large cluster colored in black at the right and a red one between them.
• Step 3: Compute the centroid, i.e. the mean of the clusters
• Repeat until no data changes cluster

The algorithm converged after seven iterations. You can run the k-mean algorithm in our dataset with five clusters and call it pc_cluster.

pc_cluster <-kmeans(rescale_df, 5)

The list pc_cluster contains seven useful elements:

• pc_cluster$cluster: The cluster assigned to each observation
• pc_cluster$centers: The cluster centres
• pc_cluster$totss: The total sum of squares
• pc_cluster$withinss: Within sum of squares, one value per cluster
• pc_cluster$tot.withinss: The sum of withinss
• pc_cluster$betweenss: Total sum of squares minus the within sum of squares
• pc_cluster$size: Number of observation within each cluster

You will use the sum of the within sum of square (i.e. tot.withinss) to compute the optimal number of clusters k. Finding k is indeed a substantial task.

How to Find the Optimal k with the Elbow Method

One technique to choose the best k is called the elbow method. This method uses within-group homogeneity or within-group heterogeneity to evaluate the variability. In other words, you are interested in the percentage of the variance explained by each cluster. You can expect the variability to increase with the number of clusters, alternatively, heterogeneity decreases. Our challenge is to find the k that is beyond the diminishing returns. Adding a new cluster does not improve the variability in the data because very few information is left to explain.

In this tutorial, we find this point using the heterogeneity measure. The Total within clusters sum of squares is the tot.withinss in the list return by kmean().

You can construct the elbow graph and find the optimal k as follow:

• Step 1: Construct a function to compute the total within clusters sum of squares
• Step 2: Run the algorithm over a range of k values
• Step 3: Create a data frame with the results of the algorithm
• Step 4: Plot the results

Step 1) Construct a function to compute the total within clusters sum of squares

You create the function that runs the k-mean algorithm and store the total within clusters sum of squares

kmean_withinss <- function(k) {
    cluster <- kmeans(rescale_df, k)
    return (cluster$tot.withinss)
}

Code Explanation

• function(k): Set the number of arguments in the function
• kmeans(rescale_df, k): Run the algorithm for this value of k
• return(cluster$tot.withinss): Store the total within clusters sum of squares

Test the function with k equal to 2.

Output:

## Try with 2 cluster

kmean_withinss(2)

Output:

## [1] 27087.07

Step 2) Run the algorithm n times

You will use the sapply() function to run the algorithm over a range of k. This technique is faster than creating a loop and store the value.

# Set maximum cluster 
max_k <-20 
# Run algorithm over a range of k 
wss <- sapply(2:max_k, kmean_withinss)

Code Explanation

• max_k <- 20: Set the maximum value of k to 20
• sapply(2:max_k, kmean_withinss): Run the function kmean_withinss() over a range 2:max_k, i.e. 2 to 20.

Step 3) Create a data frame with the results of the algorithm

With the function written and tested, run it across the range 2 to 20 and store every tot.withinss value.

# Create a data frame to plot the graph
elbow <-data.frame(2:max_k, wss)

Code Explanation

• data.frame(2:max_k, wss): Create a data frame with the output of the algorithm store in wss

Step 4) Plot the results

You plot the graph to visualize where is the elbow point

# Plot the graph with gglop
ggplot(elbow, aes(x = X2.max_k, y = wss)) +
    geom_point() +
    geom_line() +
    scale_x_continuous(breaks = seq(1, 20, by = 1))

From the graph, you can see the optimal k is seven, where the curve is starting to have a diminishing return.

Once you have our optimal k, you re-run the algorithm with k equals to 7 and evaluate the clusters.

Examining the clusters

pc_cluster_2 <-kmeans(rescale_df, 7)

As mention before, you can access the remaining interesting information in the list returned by kmean().

pc_cluster_2$cluster
pc_cluster_2$centers
pc_cluster_2$size

Evaluation is subjective and depends on what the clusters are for. The goal here is to group computers with similar specifications. A domain expert could do this by hand, but the process would be slow and error prone. K-means does the grouping objectively and leaves the expert to interpret and name the result.

As a prior evaluation, you can examine the size of the clusters.

pc_cluster_2$size

Output:

## [1] 608 1596 1231  580 1003  699  542

The largest cluster, number 2, holds 1,596 observations, while the smallest, number 7, holds only 542 computers. It might be good to have homogeneity between clusters, if not, a thinner data preparation might be required.

You get a deeper look at the data with the center component. The rows refer to the numeration of the cluster and the columns the variables used by the algorithm. The values are the average score by each cluster for the interested column. Standardization makes the interpretation easier. Positive values indicate the z-score for a given cluster is above the overall mean. For instance, cluster 4 has the highest average price (price_scal = 1.09), while cluster 5 has the lowest (-0.82).

center <-pc_cluster_2$centers
center

Output:

##   price_scal    hd_scal     ram_scal screen_scal   ads_scal trend_scal
## 1 -0.6372457 -0.7097995 -0.691520682  -0.4401632  0.6780366 -0.3379751
## 2 -0.1323863  0.6299541  0.004786730   2.6419582 -0.8894946  1.2673184
## 3  0.8745816  0.2574164  0.513105797  -0.2003237  0.6734261 -0.3300536
## 4  1.0912296 -0.2401936  0.006526723   2.6419582  0.4704301 -0.4132057
## 5 -0.8155183  0.2814882 -0.307621003  -0.3205176 -0.9052979  1.2177279
## 6  0.8830191  2.1019454  2.168706085   0.4492922 -0.9035248  1.2069855
## 7  0.2215678 -0.7132577 -0.318050275  -0.3878782 -1.3206229 -1.5490909

You can create a heat map with ggplot to help us highlight the difference between categories.

The default colors of ggplot need to be changed with the RColorBrewer library. You can use the conda library and the code to launch in the terminal:

conda install -c r r-rcolorbrewer

To create a heat map, you proceed in three steps:

• Build a data frame with the values of the center and create a variable with the number of the cluster
• Reshape the data with the gather() function of the tidyr library. You want to transform data from wide to long.
• Create the palette of colors with colorRampPalette() function

Step 1) Build a data frame

Let’s create the reshape dataset

library(tidyr)

# create dataset with the cluster number

cluster <- c(1: 7)
center_df <- data.frame(cluster, center)

# Reshape the data

center_reshape <- gather(center_df, features, values, price_scal: trend_scal)
head(center_reshape)

Output:

##   cluster   features     values
## 1       1 price_scal -0.6372457
## 2       2 price_scal -0.1323863
## 3       3 price_scal  0.8745816
## 4       4 price_scal  1.0912296
## 5       5 price_scal -0.8155183
## 6       6 price_scal  0.8830191		

Step 2) Create the colour palette

The code below builds the palette of colours used by the heat map.

library(RColorBrewer)
# Create the palette
hm.palette <-colorRampPalette(rev(brewer.pal(10, 'RdYlGn')),space='Lab')

Step 3) Visualize

You can plot the graph and see what the clusters look like.

# Plot the heat map
ggplot(data = center_reshape, aes(x = features, y = cluster, fill = values)) +
    scale_y_continuous(breaks = seq(1, 7, by = 1)) +
    geom_tile() +
    coord_equal() +
    scale_fill_gradientn(colours = hm.palette(90)) +
    theme_classic()

How to Get Reproducible K-Means Results with set.seed() and nstart

The tutorial notes that k-means returns different groups on every run. That is not a quirk to live with, it is a problem with two standard fixes, and neither is used in the examples above.

1. Fix the starting point with set.seed(). The initial centroids are drawn at random, so the same call produces different clusters each time. Setting a seed makes the whole analysis reproducible:

set.seed(123)
pc_cluster_2 <- kmeans(rescale_df, 7)

2. Run several starts and keep the best with nstart. A single random start can converge to a poor local optimum. The nstart argument runs the algorithm that many times from different random starts and returns the one with the lowest total within sum of squares:

set.seed(123)
pc_cluster_2 <- kmeans(rescale_df, centers = 7, nstart = 25)

Twenty-five starts is the usual recommendation and costs very little on a dataset of this size. Without it, the elbow curve itself becomes unstable, because each point on the plot comes from a single unlucky or lucky start. The kmean_withinss() function defined earlier should therefore be written as:

kmean_withinss <- function(k) {
    cluster <- kmeans(rescale_df, centers = k, nstart = 25)
    return (cluster$tot.withinss)
}

3. Raise iter.max if the algorithm warns. If R reports “did not converge in 10 iterations”, pass iter.max = 50 to give it more room.

K-Means vs Hierarchical Clustering in R

K-means is not the only clustering method available, and the choice usually comes down to the size of the data and whether you already know how many groups you want.

# Hierarchical alternative on the same scaled data
hc <- hclust(dist(rescale_df), method = "ward.D2")
plot(hc)
groups <- cutree(hc, k = 7)

With 6,259 observations, hierarchical clustering has to build a distance matrix of roughly 19.6 million pairs, which is why k-means is the practical choice for this dataset. Note too that k-means assumes clusters of similar size and roughly spherical shape; where that assumption fails, DBSCAN handles irregular shapes and identifies outliers instead of forcing every point into a group.

K-Means Clustering in R: Function Reference

Every function used in this tutorial is listed below:

FAQs