So far throughout the course, most of the models we have discussed have been examples of supervised learning. This means that we have an outcome variable that we are trying to predict, and we use the observed outcome values to help train the model. For example, with linear regression we may know the price of a house and try to predict that price using characteristics of the house, and with KNN we may know the species of a flower and use nearby flowers to classify a new observation. But, what happens if there is no outcome variable and we just want to see if there are natural groups hidden within the data?
This type of problem falls under the umbrella of unsupervised learning. Instead of trying to predict a known outcome, we are interested in finding structure within the predictor variables themselves. One of the most common examples of unsupervised learning is clustering, where we attempt to group observations that are similar to each other. There are many possible clustering algorithms, but we will begin with one of the most approachable methods: K-Means Clustering.
A common application might be customer segmentation. Imagine we have information about how much customers spend, how often they shop, and how long they have been customers. We may not have a variable telling us which “type” of customer each person is, but perhaps there are natural groups in the data. We may find one group of frequent low-spending customers, another group of occasional high-spending customers, and another group somewhere in between. K-Means will try to identify these groups for us.
Before we discuss exactly how K-Means works, we should distinguish clustering from the classification methods we have already seen. With classification, we are given the correct classes while training the model. With clustering, those classes do not exist ahead of time.
| Method | Outcome Known? | Main Goal |
|---|---|---|
| Linear Regression | Yes | Predict a quantitative value |
| Logistic Regression | Yes | Predict a binary outcome |
| KNN | Yes | Predict a categorical outcome using nearby observations |
| K-Means | No | Group similar observations together |
In KNN, \(K\) refers to the number of nearby observations we use to classify a new observation. In K-Means, \(K\) refers to the number of clusters we want the algorithm to find. KNN already knows the classifications of the training observations, while K-Means is trying to discover groups without being given any classifications.
For instance, suppose we have the following data:
Without being told anything else about these points, we might look at the plot and decide that there appear to be three groups. There is a group toward the bottom left, another near the center, and another toward the upper right. The goal of K-Means is essentially to formalize this process so the computer can decide which observations should be placed together, as high dimensional data will be impossible for us humans to do by just looking at it.
The main question then becomes: what do we mean when we say two observations are “similar”? Much like with KNN, we will define similarity using distance. For two-dimensional data, the Euclidean distance between an observation \((x_i,y_i)\) and the center of a cluster \((\bar{x}_c,\bar{y}_c)\) can be written as
\[d_i = \sqrt{(x_i-\bar{x}_c)^2 + (y_i-\bar{y}_c)^2}.\]
More generally, if we have \(p\) features, we can calculate the distance between observation \(i\) and the center of cluster \(c\) as
\[d(x_i,\mu_c)=\sqrt{\sum_{j=1}^{p}(x_{ij}-\mu_{cj})^2}.\]
The notation may look a little intimidating, but all we are really doing is finding the straight-line distance between the observation and the center of the cluster. We already encountered this same general idea when discussing KNN.
The K-Means algorithm has a fairly simple objective: place observations into \(K\) clusters so that observations within the same cluster are as close to each other as possible. Each cluster will have a centroid, which is the mean location of all observations assigned to that cluster. The algorithm then repeatedly moves observations and centroids until the groups stop changing.
The general process is:
Let us use the small dataset from above and select three initial centers. These are intentionally not perfect centers because the whole point of the algorithm is that they will be updated.
x <- c(1, 1, 2, 2, 7, 8, 8, 9, 4, 5, 5, 6)
y <- c(1, 2, 1, 2, 8, 8, 9, 8, 5, 4, 5, 4)
toy_data <- data.frame(x, y)
initial_centers <- matrix(c(1, 4,
4, 7,
8, 3),
nrow=3, byrow=TRUE)
plot(toy_data$x, toy_data$y, pch=19,
xlab="x", ylab="y", main="Starting Centroids")
points(initial_centers[,1], initial_centers[,2], pch=8, cex=2, lwd=3)
For every observation, K-Means will calculate the distance to each of these three centers and assign the observation to whichever center is closest. For instance, we could manually calculate the distances from the first observation \((1,1)\) to each starting centroid:
point <- c(1,1)
sqrt(sum((point - initial_centers[1,])^2))[1] 3sqrt(sum((point - initial_centers[2,])^2))[1] 6.708204sqrt(sum((point - initial_centers[3,])^2))[1] 7.28011The observation would be assigned to the cluster associated with the smallest of these three distances. After every observation has been assigned, we calculate a new mean \(x\) and mean \(y\) for each cluster. Those means become the new centroids, and then we repeat the process until the clusters do not change (this may take multiple iterations to settle down).
Luckily, we do not need to do all of these calculations manually. The kmeans() function in R will repeat the process for us until it converges.
toy_kmeans <- kmeans(toy_data, centers=initial_centers)
toy_kmeansK-means clustering with 3 clusters of sizes 4, 4, 4
Cluster means:
x y
1 1.5 1.50
2 8.0 8.25
3 5.0 4.50
Clustering vector:
[1] 1 1 1 1 2 2 2 2 3 3 3 3
Within cluster sum of squares by cluster:
[1] 2.00 2.75 3.00
(between_SS / total_SS = 95.8 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" Notice that the cluster centers are not necessarily actual observations in the dataset. They are simply the means of the observations within each cluster. This is why we call the method K-Means (…the naming people did us a favor on this one).
The algorithm is trying to minimize the within-cluster sum of squares, which we will abbreviate as WCSS. This measures how far observations are from their own cluster center:
\[\text{WCSS}=\sum_{c=1}^{K}\sum_{i\in C_c}\|x_i-\mu_c\|^2.\]
This idea should feel somewhat familiar. In regression we tried to find a line that minimized the squared residuals. Here, instead of measuring the vertical distance from an observation to a regression line, we are measuring the distance from an observation to its cluster centroid. We square these distances and try to make the total as small as possible.
We can access this information from the K-Means model:
toy_kmeans$withinss[1] 2.00 2.75 3.00toy_kmeans$tot.withinss[1] 7.75toy_kmeans$size[1] 4 4 4toy_kmeans$centers x y
1 1.5 1.50
2 8.0 8.25
3 5.0 4.50The withinss value tells us the within-cluster sum of squares for each individual cluster, while tot.withinss adds all of them together. The size output tells us how many observations were assigned to each cluster, and centers gives us the coordinates of the final centroids.
Let us now move to a dataset we have encountered before. The iris dataset contains measurements of 150 flowers from three different species. We have previously used the species as an outcome variable for classification, but this time we are going to pretend that we do not know the species. We will only give the algorithm the four quantitative measurements and see whether it can discover natural groups on its own.
head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosairis_data <- iris[,1:4]We happen to know that there are three species in the data, so we will begin by asking K-Means to create three clusters. The basic syntax is very simple: pass the quantitative data into kmeans() and specify the number of centers.
set.seed(101)
iris_kmeans <- kmeans(iris_data, centers=3)
iris_kmeansK-means clustering with 3 clusters of sizes 38, 62, 50
Cluster means:
Sepal.Length Sepal.Width Petal.Length Petal.Width
1 6.850000 3.073684 5.742105 2.071053
2 5.901613 2.748387 4.393548 1.433871
3 5.006000 3.428000 1.462000 0.246000
Clustering vector:
[1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[38] 3 3 3 3 3 3 3 3 3 3 3 3 3 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[75] 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 1 1 2 1 1 1 1
[112] 1 1 2 2 1 1 1 1 2 1 2 1 2 1 1 2 2 1 1 1 1 1 2 1 1 1 1 2 1 1 1 2 1 1 1 2 1
[149] 1 2
Within cluster sum of squares by cluster:
[1] 23.87947 39.82097 15.15100
(between_SS / total_SS = 88.4 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" There are several important pieces of information in the output. The cluster means are stored in centers, the cluster assignment for each observation is stored in cluster, and the number of observations in each cluster is stored in size.
iris_kmeans$centers Sepal.Length Sepal.Width Petal.Length Petal.Width
1 6.850000 3.073684 5.742105 2.071053
2 5.901613 2.748387 4.393548 1.433871
3 5.006000 3.428000 1.462000 0.246000iris_kmeans$size[1] 38 62 50head(iris_kmeans$cluster)[1] 3 3 3 3 3 3We can visualize the clusters using Petal Length and Petal Width:
plot(iris$Petal.Length, iris$Petal.Width,
col=iris_kmeans$cluster, pch=19,
xlab="Petal Length", ylab="Petal Width",
main="K-Means Clustering of Iris")
points(iris_kmeans$centers[,"Petal.Length"],
iris_kmeans$centers[,"Petal.Width"],
pch=8, cex=2, lwd=3)
Even though we are only displaying two variables in the graph, remember that the model above used all four quantitative variables when deciding which observations belonged together, so that is why the observations might not appear to belong to the closest centroid.
Since we actually know the flower species, we can compare the clusters K-Means discovered with the true species after fitting the model. The species variable was not used by K-Means at all.
table(Actual_Species=iris$Species, Cluster=iris_kmeans$cluster) Cluster
Actual_Species 1 2 3
setosa 0 0 50
versicolor 2 48 0
virginica 36 14 0We should be careful when interpreting this table. Cluster 1 does not automatically mean virginica, Cluster 2 does not automatically mean versicolor, and so on. The numbers assigned to the clusters are arbitrary labels. If we ran the model again, the exact cluster numbers could switch even if the underlying groupings remained essentially the same.
What we should look for is whether the clusters line up reasonably well with the known species. We will usually see that setosa forms a very distinct group, while there is more overlap between versicolor and virginica. This makes sense if we look at the flower measurements, as setosa tends to be separated more clearly from the other species.
This example is useful because it lets us check how well the clusters correspond to a known grouping. In most real clustering problems though, we will not have a “correct” outcome variable sitting off to the side. If we already had that outcome variable, we would probably be doing classification instead!
There is one additional complication we should discuss before moving on. The K-Means algorithm begins with starting centroids, and different starting locations can sometimes produce different final clusters. This means that K-Means can settle on a solution that is good, but not necessarily the absolute best possible solution.
To reduce this issue, we can use the nstart argument. If we specify nstart=25, R will run K-Means 25 times using different starting centroids and keep the solution with the smallest total within-cluster sum of squares.
set.seed(101)
iris_kmeans <- kmeans(iris_data, centers=3, nstart=25)
iris_kmeans$tot.withinss[1] 78.85144As a general rule, I would recommend using multiple starts rather than relying on a single random initialization. There is usually very little reason to trust one random starting point when R is perfectly willing to try several for us.
K-Means is a distance-based algorithm, which means the scale of the variables matters a great deal. We ran into the same issue with KNN. Suppose one feature ranges from 0 to 1 while another ranges from 0 to 10,000. A difference of 1 unit in the second variable is numerically much larger than a difference of 1 unit in the first variable, even though the first difference may be much more meaningful. The variable with the larger scale will tend to dominate the distance calculation.
A good example of this issue can be seen using the USArrests dataset, which contains violent crime statistics for the 50 states in 1973.
head(USArrests) Murder Assault UrbanPop Rape
Alabama 13.2 236 58 21.2
Alaska 10.0 263 48 44.5
Arizona 8.1 294 80 31.0
Arkansas 8.8 190 50 19.5
California 9.0 276 91 40.6
Colorado 7.9 204 78 38.7summary(USArrests) Murder Assault UrbanPop Rape
Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30
1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07
Median : 7.250 Median :159.0 Median :66.00 Median :20.10
Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23
3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18
Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00 apply(USArrests, 2, sd) Murder Assault UrbanPop Rape
4.355510 83.337661 14.474763 9.366385 Notice that the variables are measured on very different scales. Assault, for example, is much larger numerically than Murder. If we perform K-Means directly on the raw data, the variables with the largest numerical scales will have the greatest influence on the clusters.
set.seed(101)
us_raw_kmeans <- kmeans(USArrests, centers=3, nstart=25)
us_raw_kmeans$centers Murder Assault UrbanPop Rape
1 11.812500 272.5625 68.31250 28.37500
2 4.270000 87.5500 59.75000 14.39000
3 8.214286 173.2857 70.64286 22.84286To avoid this problem, we can standardize each feature before clustering. A standardized value, or z-score, can be calculated as
\[z=\frac{x-\bar{x}}{s}.\]
This transforms every variable so that it has a mean of approximately 0 and a standard deviation of 1. Therefore, each variable begins on a comparable scale.
USArrests_scaled <- scale(USArrests)
apply(USArrests_scaled, 2, mean) Murder Assault UrbanPop Rape
1.543210e-16 1.143530e-16 -3.996803e-16 8.526513e-16 apply(USArrests_scaled, 2, sd) Murder Assault UrbanPop Rape
1 1 1 1 We can now fit the model using the standardized data.
set.seed(101)
us_kmeans <- kmeans(USArrests_scaled, centers=3, nstart=25)
us_kmeansK-means clustering with 3 clusters of sizes 13, 20, 17
Cluster means:
Murder Assault UrbanPop Rape
1 -0.9615407 -1.1066010 -0.9301069 -0.9667633
2 1.0049340 1.0138274 0.1975853 0.8469650
3 -0.4469795 -0.3465138 0.4788049 -0.2571398
Clustering vector:
Alabama Alaska Arizona Arkansas California
2 2 2 3 2
Colorado Connecticut Delaware Florida Georgia
2 3 3 2 2
Hawaii Idaho Illinois Indiana Iowa
3 1 2 3 1
Kansas Kentucky Louisiana Maine Maryland
3 1 2 1 2
Massachusetts Michigan Minnesota Mississippi Missouri
3 2 1 2 2
Montana Nebraska Nevada New Hampshire New Jersey
1 1 2 1 3
New Mexico New York North Carolina North Dakota Ohio
2 2 2 1 3
Oklahoma Oregon Pennsylvania Rhode Island South Carolina
3 3 3 3 2
South Dakota Tennessee Texas Utah Vermont
1 2 2 3 1
Virginia Washington West Virginia Wisconsin Wyoming
3 3 1 1 3
Within cluster sum of squares by cluster:
[1] 11.95246 46.74796 19.62285
(between_SS / total_SS = 60.0 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" When we interpret the final clusters, though, standardized values are not always the easiest numbers to explain. It may be more meaningful to return to the original variables and calculate the average crime statistics within each cluster.
aggregate(USArrests, by=list(Cluster=us_kmeans$cluster), mean) Cluster Murder Assault UrbanPop Rape
1 1 3.600000 78.53846 52.07692 12.17692
2 2 12.165000 255.25000 68.40000 29.16500
3 3 5.841176 141.88235 72.47059 18.82353This output allows us to describe what makes each cluster different. For instance, one group may contain states with relatively low values across the crime variables while another may contain states with higher violent-crime measurements. The important part is that we are interpreting the clusters after the algorithm creates them. K-Means does not tell us that a group should be called “low crime” or “high crime”; those descriptions are ones we create by examining the characteristics of each cluster.
Since K-Means relies on Euclidean distance, I would typically recommend scaling quantitative variables unless there is a meaningful reason not to. If all variables are already measured on the same scale and a one-unit difference means roughly the same thing across variables, then scaling may not be necessary. But if the scales are substantially different, we should at least stop and think about it before running the model.
So far we have conveniently told K-Means how many clusters to create. With iris, we knew there were three species, so using \(K=3\) was an obvious place to start. In a real unsupervised learning problem, however, the whole point may be that we do not know how many groups exist. This creates one of the most important questions in K-Means clustering: how do we choose K?
One possibility is to compare the total within-cluster sum of squares for several different values of \(K\). Remember that we want this value to be small because it means observations are close to their own cluster centers. Let us calculate the total within-cluster sum of squares for values of \(K\) from 1 through 10.
wss <- c()
for(k in 1:10){
set.seed(101)
model <- kmeans(USArrests_scaled, centers=k, nstart=25)
wss[k] <- model$tot.withinss
}
wss [1] 196.00000 102.86240 78.32327 56.40317 48.94420 42.83303 38.38847
[8] 34.27969 29.94611 26.18348We can then plot these values:
plot(1:10, wss, type="b", pch=19,
xlab="Number of Clusters (K)",
ylab="Total Within-Cluster Sum of Squares",
main="Elbow Plot")
As \(K\) increases, the within-cluster sum of squares will always decrease. This makes sense because it becomes easier to keep observations close to a centroid when we are allowed to create more clusters. In the ridiculous extreme case where every observation had its own cluster, the within-cluster distance would be 0. That would technically minimize our WCSS, but it would not be a very useful clustering model.
Instead, we look for an elbow in the graph, which is where we see a big change in slope (this might be 2 or 4 for this graph). This is the point where adding another cluster stops producing a dramatic improvement in the within-cluster sum of squares. The graph may drop sharply at first and then begin to level off. The bend between those two patterns is a reasonable candidate for \(K\).
The elbow method is helpful, but it is not magic. Sometimes there is a very obvious elbow and other times the plot looks more like an arm without an elbow (…which is not especially useful to us). Therefore, we may want another method to help support our decision.
One option is the silhouette score. The silhouette score compares how close an observation is to the other observations in its own cluster with how close it is to observations in neighboring clusters. The value ranges from \(-1\) to \(1\).
We can calculate the average silhouette score for several possible values of \(K\) using the cluster library.
library(cluster)
silhouette_mean <- c()
for(k in 2:10){
set.seed(101)
model <- kmeans(USArrests_scaled, centers=k, nstart=25)
sil <- silhouette(model$cluster, dist(USArrests_scaled))
silhouette_mean[k-1] <- mean(sil[,3])
}
silhouette_mean[1] 0.4084890 0.3094312 0.3396889 0.3030781 0.2859821 0.2916528 0.2579320
[8] 0.2576100 0.2623532Then we can plot the average silhouette score:
plot(2:10, silhouette_mean, type="b", pch=19,
xlab="Number of Clusters (K)",
ylab="Average Silhouette Score",
main="Silhouette Method")
With the silhouette method, larger values are preferred. Therefore, we might choose the value of \(K\) that produces the highest average silhouette score.
It is completely possible for the elbow method, silhouette score, and our practical interpretation of the clusters to suggest slightly different values of \(K\). This is not necessarily a problem. Unlike supervised learning, we often do not have a true outcome that tells us whether we are right or wrong. Choosing the number of clusters is partly a statistical question and partly an interpretation question.
Let us finish with one more example using mtcars. We will cluster the cars based on fuel efficiency (mpg), horsepower (hp), weight (wt), and quarter-mile time (qsec).
car_data <- mtcars[,c("mpg", "hp", "wt", "qsec")]
head(car_data) mpg hp wt qsec
Mazda RX4 21.0 110 2.620 16.46
Mazda RX4 Wag 21.0 110 2.875 17.02
Datsun 710 22.8 93 2.320 18.61
Hornet 4 Drive 21.4 110 3.215 19.44
Hornet Sportabout 18.7 175 3.440 17.02
Valiant 18.1 105 3.460 20.22Since horsepower, weight, fuel efficiency, and quarter-mile time are all measured on different scales, we should standardize the data before fitting the model.
car_scaled <- scale(car_data)
set.seed(101)
car_kmeans <- kmeans(car_scaled, centers=3, nstart=50)We can visualize the clusters using weight and fuel efficiency:
plot(mtcars$wt, mtcars$mpg, col=car_kmeans$cluster, pch=19,
xlab="Weight (1000 lbs)", ylab="Miles per Gallon",
main="K-Means Clusters for mtcars")
text(mtcars$wt, mtcars$mpg, labels=rownames(mtcars), pos=3, cex=.6)
Just like we did with USArrests, we should interpret the cluster characteristics using the variables in their original units.
aggregate(car_data, by=list(Cluster=car_kmeans$cluster), mean) Cluster mpg hp wt qsec
1 1 21.03636 103.3636 2.995909 19.13273
2 2 15.40667 206.9333 3.917267 16.68733
3 3 30.06667 75.5000 1.873000 18.39833The cluster numbers themselves are meaningless. There is nothing inherently “Cluster 1-ish” about the cars assigned to Cluster 1. Instead, we might examine the group averages and find a cluster of lighter, more fuel-efficient cars, another cluster of heavier and more powerful cars, and a third group somewhere between those extremes. We give the clusters meaning by examining the variables that characterize them.
K-Means is appealing because the underlying idea is fairly simple: group observations around the nearest mean. But, much like the other methods we have seen, most of the important decisions happen around the model rather than inside the function itself. We still need to decide which variables to include, whether they need to be scaled, how many clusters make sense, and whether the clusters tell us anything useful. The function can find the groups, but it is still our job as Data Scientists to decide what those groups mean.