Unsupervised Learning in R

1 Unsupervised learning in R

1.1 Introduction

1.1.1 Identify clustering problems

Example of clustering problems:

• Determining the natural groupings of houses for sale based on size, number of bedrooms, etc.
• Determining if there are common patterns in the demographics of people at a commerce site

1.2 Introduction to k-means clustering

1.2.1 k-means clustering

We have created some two-dimensional data and stored it in a variable called x.

plot(x)

Our task is to create a k-means model of the x data using 3 clusters, then to look at the structure of the resulting model using the summary() function.

# Create the k-means model: km.out
km.out <- kmeans(x, centers = 3, nstart = 20 )

# Inspect the result
summary(km.out)

             Length Class  Mode   
cluster      300    -none- numeric
centers        6    -none- numeric
totss          1    -none- numeric
withinss       3    -none- numeric
tot.withinss   1    -none- numeric
betweenss      1    -none- numeric
size           3    -none- numeric
iter           1    -none- numeric
ifault         1    -none- numeric

1.2.2 Results of kmeans()

The kmeans() function produces several outputs.

We will access the cluster component directly. This is useful anytime we need the cluster membership for each observation of the data used to build the clustering model.

k-means models also have a print method to give a human friendly output of basic modeling results. This is available by using print() or simply typing the name of the model.

# Print the cluster membership component of the model
km.out$cluster

  [1] 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [35] 3 3 2 2 2 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3
 [69] 3 3 3 3 3 2 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 1 1
[103] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[137] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[171] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[205] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[239] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[273] 2 2 2 3 3 2 2 3 2 2 2 2 2 2 3 2 2 2 2 2 2 3 2 2 2 3 2 2

# Print the km.out object
print(km.out)

K-means clustering with 3 clusters of sizes 150, 52, 98

Cluster means:
        [,1]        [,2]
1 -5.0556758  1.96991743
2  0.6642455 -0.09132968
3  2.2171113  2.05110690

Clustering vector:
  [1] 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [35] 3 3 2 2 2 3 3 3 3 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3 3 3 3 3 3
 [69] 3 3 3 3 3 2 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 1 1
[103] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[137] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[171] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[205] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[239] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[273] 2 2 2 3 3 2 2 3 2 2 2 2 2 2 3 2 2 2 2 2 2 3 2 2 2 3 2 2

Within cluster sum of squares by cluster:
[1] 295.16925  95.50625 148.64781
 (between_SS / total_SS =  87.2 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"    
[5] "tot.withinss" "betweenss"    "size"         "iter"        
[9] "ifault"      

Because printing the whole model object to the console outputs many different things, we may wish to instead print a specific component of the model object using the $ operator.

1.2.3 Visualizing and interpreting results of kmeans()

One of the more intuitive ways to interpret the results of k-means models is by plotting the data as a scatter plot and using color to label the samples’ cluster membership. We will use the standard plot() function to accomplish this.

To create a scatter plot, we can pass data with two features (i.e. columns) to plot() with an extra argument col = km.out$cluster, which sets the color of each point in the scatter plot according to its cluster membership.

# Scatter plot of x
plot(x, col = km.out$cluster, main = "k-means with 3 clusters",
xlab = "", ylab = "" )

1.3 How k-means works and practical matters

1.3.1 Handling random algorithms

We saw how kmeans() randomly initializes the centers of clusters. This random initialization can result in assigning observations to different cluster labels. Also, the random initialization can result in finding different local minima for the k-means algorithm. This example will demonstrate both results.

At the top of each plot, the measure of model quality—total within cluster sum of squares error—will be plotted. Look for the model(s) with the lowest error to find models with the better model results.

Because kmeans() initializes observations to random clusters, it is important to set the random number generator seed for reproducibility.

# Set up 2 x 3 plotting grid
par(mfrow = c(2, 3))

# Set seed
set.seed(1)

for(i in 1:6) {
  # Run kmeans() on x with three clusters and one start
  km.out <- kmeans(x, centers = 3, nstart = 1)
  
  # Plot clusters
  plot(x, col = km.out$cluster, 
       main = km.out$tot.withinss, 
       xlab = "", ylab = "")
}

Because of the random initialization of the k-means algorithm, there’s quite some variation in cluster assignments among the six models.

1.3.2 Selecting number of clusters

The k-means algorithm assumes the number of clusters as part of the input. If we know the number of clusters in advance (e.g. due to certain business constraints) this makes setting the number of clusters easy. However, if we do not know the number of clusters and need to determine it, we will need to run the algorithm multiple times, each time with a different number of clusters. From this, we can observe how a measure of model quality changes with the number of clusters.

We will run kmeans() multiple times to see how model quality changes as the number of clusters changes. Plots displaying this information help to determine the number of clusters and are often referred to as scree plots.

The ideal plot will have an elbow where the quality measure improves more slowly as the number of clusters increases. This indicates that the quality of the model is no longer improving substantially as the model complexity (i.e. number of clusters) increases. In other words, the elbow indicates the number of clusters inherent in the data.

# Initialize total within sum of squares error: wss
wss <- 0

# For 1 to 15 cluster centers
for (i in 1:15) {
  km.out <- kmeans(x, centers = i, nstart = 20)
  # Save total within sum of squares to wss variable
  wss[i] <- km.out$tot.withinss
}

# Plot total within sum of squares vs. number of clusters
plot(1:15, wss, type = "b", 
     xlab = "Number of Clusters", 
     ylab = "Within groups sum of squares")

# Set k equal to the number of clusters corresponding to the elbow location
k <- 2

Looking at the scree plot, it looks like there are inherently 2 or 3 clusters in the data.

1.4 Introduction to the Pokemon data

1.4.1 Practical matters: working with real data

Dealing with real data is often more challenging than dealing with synthetic data. Synthetic data helps with learning new concepts and techniques, but we will deal with data that is closer to the type of real data you might find in your professional or academic pursuits.

The first challenge with the Pokemon data is that there is no pre-determined number of clusters. We will determine the appropriate number of clusters, keeping in mind that in real data the elbow in the scree plot might be less of a sharp elbow than in synthetic data. Use your judgement on making the determination of the number of clusters.

This example includes plotting the outcomes of the clustering on two dimensions, or features, of the data. These features were chosen somewhat arbitrarily for this exercise. Think about how you would use plotting and clustering to communicate interesting groups of Pokemon to other people.

An additional note: this example utilizes the iter.max argument to kmeans(). As we’ve seen, kmeans() is an iterative algorithm, repeating over and over until some stopping criterion is reached. The default number of iterations for kmeans() is 10, which is not enough for the algorithm to converge and reach its stopping criterion, so we’ll set the number of iterations to 50 to overcome this issue. To see what happens when kmeans() does not converge, try running the example with a lower number of iterations (e.g. 3). This is another example of what might happen when we encounter real data and use real cases.

library(readr)
pokemon <- read_csv("Pokemon.csv") %>% 
  select(HitPoints, Attack, Defense, SpecialAttack, SpecialDefense, Speed)

Parsed with column specification:
cols(
  Number = col_double(),
  Name = col_character(),
  Type1 = col_character(),
  Type2 = col_character(),
  Total = col_double(),
  HitPoints = col_double(),
  Attack = col_double(),
  Defense = col_double(),
  SpecialAttack = col_double(),
  SpecialDefense = col_double(),
  Speed = col_double(),
  Generation = col_double(),
  Legendary = col_logical()
)

# Initialize total within sum of squares error: wss
wss <- 0

# Look over 1 to 15 possible clusters
for (i in 1:15) {
  # Fit the model: km.out
  km.out <- kmeans(pokemon, centers = i, nstart = 20, iter.max = 50)
  # Save the within cluster sum of squares
  wss[i] <- km.out$tot.withinss
}

# Produce a scree plot
plot(1:15, wss, type = "b", 
     xlab = "Number of Clusters", 
     ylab = "Within groups sum of squares")

# Select number of clusters
k <- 2

# Build model with k clusters: km.out
km.out <- kmeans(pokemon, centers = k, nstart = 20, iter.max = 50)

# View the resulting model
km.out

K-means clustering with 2 clusters of sizes 374, 426

Cluster means:
  HitPoints   Attack  Defense SpecialAttack SpecialDefense    Speed
1  54.87968 57.94118 55.49733      52.33155       53.44118 54.04278
2  81.88263 97.49061 89.94836      90.80751       88.11033 80.77465

Clustering vector:
  [1] 1 1 2 2 1 1 2 2 2 1 1 2 2 1 1 1 1 1 1 2 1 1 2 2 1 1 1 2 1 2 1 2 1 2
 [35] 1 1 2 1 1 2 1 2 1 2 1 1 1 2 1 1 2 1 1 1 2 1 1 1 1 1 2 1 2 1 2 1 1 2
 [69] 1 1 2 2 1 1 2 1 1 2 1 2 1 1 2 1 2 1 2 2 1 2 1 1 2 1 2 1 2 1 2 1 1 2
[103] 2 1 1 2 1 2 1 2 1 2 1 2 2 2 1 1 2 1 2 1 2 2 2 1 2 1 2 1 2 2 2 2 2 2
[137] 2 2 2 1 2 2 2 1 1 2 2 2 1 1 2 1 2 2 2 2 2 2 2 1 1 2 2 2 2 2 1 1 2 1
[171] 1 2 1 1 2 1 1 1 2 1 1 1 1 2 1 2 1 1 1 1 1 1 2 1 1 2 2 2 1 1 1 2 1 1
[205] 2 1 1 2 1 1 2 2 2 1 2 2 1 1 2 1 2 1 1 2 2 1 2 1 2 2 2 2 2 1 1 2 1 1
[239] 1 2 1 1 2 1 2 2 1 2 2 2 1 2 2 2 1 1 2 1 1 1 2 2 2 2 2 1 1 2 2 2 2 2
[273] 1 1 2 2 1 1 2 2 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 2 1 1 2
[307] 2 1 1 1 2 1 2 2 1 2 1 1 1 2 1 2 1 1 1 1 1 2 1 2 1 2 2 2 1 1 2 1 2 2
[341] 1 1 1 1 1 1 2 1 2 2 1 2 1 2 2 2 1 2 1 1 1 2 1 2 1 2 2 2 2 2 2 1 2 1
[375] 2 1 2 1 2 1 2 1 2 1 2 1 2 2 1 2 2 1 2 2 1 1 2 2 1 1 2 1 2 2 2 1 1 1
[409] 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 1 1 2 1 1 2 1
[443] 1 2 1 1 1 1 1 1 2 1 2 1 2 1 2 1 2 2 2 2 1 2 1 1 2 1 2 1 2 2 1 2 1 2
[477] 2 2 2 1 2 1 1 2 1 2 1 1 1 1 2 1 1 2 2 1 1 2 2 1 2 1 2 1 2 2 1 2 1 1
[511] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[545] 2 2 2 2 2 2 2 2 2 2 1 1 2 1 1 2 1 1 2 1 1 1 1 2 1 2 1 2 1 2 1 2 1 2
[579] 1 1 2 1 2 1 1 2 1 1 1 2 2 2 1 1 2 1 1 2 2 2 1 1 2 1 1 2 1 2 1 2 2 1
[613] 1 2 1 2 2 2 1 2 1 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 2 1 1 2 1 2 1 1 2
[647] 1 2 1 1 2 1 2 1 2 2 1 2 1 2 1 2 2 1 1 2 1 2 1 1 2 1 1 2 1 2 2 1 2 2
[681] 1 2 2 1 2 1 2 2 1 2 1 2 2 2 1 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[715] 2 2 2 2 1 1 2 1 1 2 1 1 2 1 1 1 1 2 1 1 1 1 2 1 1 2 1 2 1 2 2 1 2 2
[749] 1 2 2 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 2 2 1 2 1 2 2 2 1 2 1 1
[783] 1 1 2 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2

Within cluster sum of squares by cluster:
[1]  911966.2 2007145.9
 (between_SS / total_SS =  31.9 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"    
[5] "tot.withinss" "betweenss"    "size"         "iter"        
[9] "ifault"      

# Plot of Defense vs. Speed by cluster membership
plot(pokemon[, c("Defense", "Speed")],
     col = km.out$cluster,
     main = paste("k-means clustering of Pokemon with", k, "clusters"),
     xlab = "Defense", ylab = "Speed")

# Select number of clusters
k <- 3

# Build model with k clusters: km.out
km.out <- kmeans(pokemon, centers = k, nstart = 20, iter.max = 50)
km.pokemon <- km.out
# View the resulting model
km.out

K-means clustering with 3 clusters of sizes 175, 270, 355

Cluster means:
  HitPoints   Attack   Defense SpecialAttack SpecialDefense    Speed
1  79.30857 97.29714 108.93143      66.71429       87.04571 57.29143
2  81.90370 96.15926  77.65556     104.12222       86.87778 94.71111
3  54.68732 56.93239  53.64507      52.02254       53.04789 53.58873

Clustering vector:
  [1] 3 3 2 2 3 3 2 2 2 3 3 1 2 3 3 3 3 3 3 2 3 3 2 2 3 3 3 2 3 2 3 2 3 1
 [35] 3 3 1 3 3 2 3 2 3 2 3 3 3 2 3 3 2 3 1 3 2 3 3 3 2 3 2 3 2 3 2 3 3 1
 [69] 3 2 2 2 3 1 1 3 3 2 3 2 3 1 1 3 2 3 1 1 3 2 3 3 2 3 1 3 1 3 1 3 2 2
[103] 2 1 3 1 3 1 3 2 3 2 3 1 1 1 3 3 1 3 1 3 1 1 1 3 2 3 1 3 2 2 2 2 2 2
[137] 1 1 1 3 1 1 1 3 3 2 2 2 3 3 1 3 1 2 2 1 2 2 2 3 3 2 2 2 2 2 3 3 1 3
[171] 3 2 3 3 1 3 3 3 2 3 3 3 3 2 3 2 3 3 3 3 3 3 2 3 3 2 2 1 3 3 1 2 3 3
[205] 2 3 3 3 3 3 1 2 1 3 1 2 3 3 2 3 1 3 1 1 1 3 1 3 1 1 1 1 1 3 3 1 3 1
[239] 3 1 3 3 2 3 2 1 3 2 2 2 3 1 2 2 3 3 1 3 3 3 1 2 2 2 1 3 3 1 1 2 2 2
[273] 3 3 2 2 3 3 2 2 3 3 1 1 3 3 3 3 3 3 3 3 3 3 3 2 3 3 2 3 3 3 1 3 3 2
[307] 2 3 3 3 1 3 3 2 3 2 3 3 3 2 3 1 3 1 3 3 3 1 3 1 3 1 1 1 3 3 2 3 2 2
[341] 3 3 3 3 3 3 1 3 2 2 3 2 3 2 1 1 3 2 3 3 3 2 3 2 3 1 2 2 2 2 1 3 1 3
[375] 1 3 1 3 1 3 1 3 2 3 1 3 2 2 3 1 1 3 2 2 3 3 2 2 3 3 2 3 1 1 1 3 3 1
[409] 2 2 3 1 1 2 1 1 1 2 2 2 2 2 2 1 2 2 2 2 2 2 1 2 3 1 1 3 3 2 3 3 2 3
[443] 3 2 3 3 3 3 3 3 2 3 2 3 1 3 1 3 1 1 1 2 3 1 3 3 2 3 2 3 1 2 3 2 3 2
[477] 2 2 2 3 2 3 3 2 3 1 3 3 3 3 1 3 3 2 2 3 3 2 2 3 1 3 1 3 2 1 3 2 3 3
[511] 2 1 2 2 1 1 1 2 2 2 2 1 2 1 2 2 2 2 1 1 2 2 2 2 2 2 2 1 2 2 2 2 2 2
[545] 2 2 1 2 2 2 2 2 2 2 3 3 2 3 3 2 3 3 2 3 3 3 3 1 3 2 3 2 3 2 3 2 3 1
[579] 3 3 2 3 2 3 1 1 3 2 3 2 1 1 3 1 1 3 3 2 1 1 3 3 2 3 3 2 3 2 3 2 2 3
[613] 3 2 3 1 2 2 3 1 3 1 2 3 1 3 1 3 2 3 1 3 2 3 2 3 3 1 3 3 2 3 2 3 3 2
[647] 3 2 2 3 1 3 1 3 2 1 3 2 3 1 3 1 1 3 3 2 3 2 3 3 2 3 1 1 3 1 2 3 2 1
[681] 3 2 1 3 1 3 1 1 3 1 3 1 2 1 3 3 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[715] 2 2 2 2 3 1 1 3 3 2 3 3 2 3 3 3 3 2 3 3 3 3 2 3 3 2 3 2 3 1 2 3 2 2
[749] 3 1 2 1 3 1 3 2 3 1 3 1 3 1 3 2 3 2 3 1 3 2 2 2 2 1 3 2 2 1 3 1 3 3
[783] 3 3 1 1 1 1 3 1 3 2 2 2 1 1 2 2 2 2

Within cluster sum of squares by cluster:
[1]  709020.5 1018348.0  812079.9
 (between_SS / total_SS =  40.8 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"    
[5] "tot.withinss" "betweenss"    "size"         "iter"        
[9] "ifault"      

# Plot of Defense vs. Speed by cluster membership
plot(pokemon[, c("Defense", "Speed")],
     col = km.out$cluster,
     main = paste("k-means clustering of Pokemon with", k, "clusters"),
     xlab = "Defense", ylab = "Speed")

2 Hierarchical clustering

2.1 Introduction

2.1.1 Hierarchical clustering with results

We will create your first hierarchical clustering model using the hclust() function.

We will create a hierarchical clustering model of x. The first step to hierarchical clustering is determining the similarity between observations, which you will do with the dist() function.

# Create hierarchical clustering model: hclust.out
hclust.out <- hclust(dist(x))

# Inspect the result
summary(hclust.out)

            Length Class  Mode     
merge       598    -none- numeric  
height      299    -none- numeric  
order       300    -none- numeric  
labels        0    -none- NULL     
method        1    -none- character
call          2    -none- call     
dist.method   1    -none- character

Now that we’ve made your first hierarchical clustering model, let’s learn how to use it to solve problems.

2.2 Selecting number of clusters

plot(hclust.out)
abline(h = c(3.5,4.5, 6.9, 9.0), col = "red")

If we cut the tree at a height of 6.9, we’re left with 3 branches representing 3 distinct clusters.

2.2.1 Cutting the tree

cutree() is the R function that cuts a hierarchical model. The h and k arguments to cutree() allow you to cut the tree based on a certain height h or a certain number of clusters k.

# Cut by height
cutree(hclust.out, h = 7)

  [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 2 1 1 2 1 1 1 2 1 2 1 2 1 1 1 2
 [35] 1 2 2 2 2 1 1 1 1 2 2 1 1 1 2 1 1 2 1 1 1 1 1 1 2 1 2 1 1 1 2 1 1 1
 [69] 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 2 1 2 2 1 1 3 3
[103] 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
[137] 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
[171] 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
[205] 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
[239] 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 1 2 2 2 2 2 1 1 2 2 2 2 1 2 2 2 2 2
[273] 2 2 2 1 2 2 1 1 1 2 2 2 2 2 1 2 2 2 2 2 2 1 2 2 2 1 2 2

# Cut by number of clusters
cutree(hclust.out, k = 3)

  [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 2 1 1 2 1 1 1 2 1 2 1 2 1 1 1 2
 [35] 1 2 2 2 2 1 1 1 1 2 2 1 1 1 2 1 1 2 1 1 1 1 1 1 2 1 2 1 1 1 2 1 1 1
 [69] 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 2 1 2 2 1 1 3 3
[103] 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
[137] 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
[171] 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
[205] 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
[239] 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 1 2 2 2 2 2 1 1 2 2 2 2 1 2 2 2 2 2
[273] 2 2 2 1 2 2 1 1 1 2 2 2 2 2 1 2 2 2 2 2 2 1 2 2 2 1 2 2

The output of each cutree() call represents the cluster assignments for each observation in the original dataset.

2.3 Clustering linkage and practical matters

2.3.1 Linkage methods

We will produce hierarchical clustering models using different linkages and plot the dendrogram for each, observing the overall structure of the trees.

# Cluster using complete linkage: hclust.complete
hclust.complete <- hclust(dist(x), method = "complete")

# Cluster using average linkage: hclust.average
hclust.average <- hclust(dist(x), method = "average")

# Cluster using single linkage: hclust.single
hclust.single <- hclust(dist(x), method = "single")

# Plot dendrogram of hclust.complete
plot(hclust.complete, main = "Complete")

# Plot dendrogram of hclust.average
plot(hclust.average, main = "Average")

# Plot dendrogram of hclust.single
plot(hclust.single, main = "Single")

Whether we want balanced or unbalanced trees for our hierarchical clustering model depends on the context of the problem we’re trying to solve. Balanced trees are essential if we want an even number of observations assigned to each cluster. On the other hand, if we want to detect outliers, for example, an unbalanced tree is more desirable because pruning an unbalanced tree can result in most observations assigned to one cluster and only a few observations assigned to other clusters.

2.3.2 Practical matters: scaling

Clustering real data may require scaling the features if they have different distributions.

We will go back to working with “real” data, the pokemon dataset introduced in the first chapter. We will observe the distribution (mean and standard deviation) of each feature, scale the data accordingly, then produce a hierarchical clustering model using the complete linkage method.

# View column means
colMeans(pokemon)

     HitPoints         Attack        Defense  SpecialAttack 
      69.25875       79.00125       73.84250       72.82000 
SpecialDefense          Speed 
      71.90250       68.27750 

# View column standard deviations
apply(pokemon, 2, sd)

     HitPoints         Attack        Defense  SpecialAttack 
      25.53467       32.45737       31.18350       32.72229 
SpecialDefense          Speed 
      27.82892       29.06047 

# Scale the data
pokemon.scaled = scale(pokemon)

# Create hierarchical clustering model: hclust.pokemon
hclust.pokemon = hclust(dist(pokemon.scaled), method = "complete")

We first checked to see if the column means and standard deviations vary. Because they do, we scaled the data, converted the scaled data to a similarity matrix and passed it into the hclust() function.

2.3.3 Comparing kmeans() and hclust()

Comparing k-means and hierarchical clustering, we’ll see the two methods produce different cluster memberships. This is because the two algorithms make different assumptions about how the data is generated.

# Apply cutree() to hclust.pokemon: cut.pokemon
cut.pokemon <- cutree(hclust.pokemon, k = 3)

# Compare methods
table(km.pokemon$cluster, cut.pokemon)

   cut.pokemon
      1   2   3
  1 171   3   1
  2 267   3   0
  3 350   5   0

Looking at the table, it looks like the hierarchical clustering model assigns most of the observations to cluster 1, while the k-means algorithm distributes the observations relatively evenly among all clusters. It’s important to note that there’s no consensus on which method produces better clusters. The job of the analyst in unsupervised clustering is to observe the cluster assignments and make a judgment call as to which method provides more insights into the data.

3 Dimensionality reduction with PCA

3.1 Introduction

3.1.1 PCA using prcomp()

We will create a s principal component analysis (PCA) model and observe the diagnostic results.

pokemon_pca is a 4 dimensions version of our pokemon dataset.

head(pokemon_pca)

           HitPoints Attack Defense Speed
Quilava           58     64      58    80
Goodra            90    100      70    80
Mothim            70     94      50    66
Marowak           60     80     110    45
Chandelure        60     55      90    80
Helioptile        44     38      33    70

# Perform scaled PCA: pr.out
pr.out <- prcomp(x = pokemon_pca,
scale = TRUE)

# Inspect model output
summary(pr.out)

Importance of components:
                          PC1    PC2    PC3     PC4
Standard deviation     1.4420 1.0013 0.7941 0.53595
Proportion of Variance 0.5199 0.2507 0.1577 0.07181
Cumulative Proportion  0.5199 0.7705 0.9282 1.00000

What is the minimum number of principal components that are required to describe at least 75% of the cumulative variance in this dataset? 2. The first two principal components describe around 77% of the variance.

3.1.2 Additional results of PCA

PCA models in R produce additional diagnostic and output components:

• center: the column means used to center to the data, or FALSE if the data weren’t centered
• scale: the column standard deviations used to scale the data, or FALSE if the data weren’t scaled
• rotation: the directions of the principal component vectors in terms of the original features/variables. This information allows you to define new data in terms of the original principal components
• x: the value of each observation in the original dataset projected to the principal components

We can access these the same as other model components. For example, use pr.out$rotation to access the rotation component.

The x component is a table with the same dimensions as the original data.

dim(pr.out$x)

[1] 50  4

dim(pokemon_pca)

[1] 50  4

The data were centered prior to performing PCA.

pr.out$center

HitPoints    Attack   Defense     Speed 
    71.08     81.22     78.44     66.58 

The data were scaled prior to performing PCA.

pr.out$scale

HitPoints    Attack   Defense     Speed 
 25.62193  33.03078  32.05809  27.51036 

The directions of the principal component vectors are presented in a table with different dimensions as the original data.

dim(pr.out$rotation)

[1] 4 4

dim(pokemon_pca)

[1] 50  4

3.2 Visualizing and interpreting PCA results

3.2.1 Interpreting biplots

The biplot() function plots both the principal components loadings and the mapping of the observations to their first two principal component values.

biplot(pr.out)

which two original variables have approximately the same loadings in the first two principal components? Attack and HitPoints

which two Pokemon are the least similar in terms of the second principal component? Kadabra and Torkoal

3.2.2 Variance explained

The second common plot type for understanding PCA models is a scree plot. A scree plot shows the variance explained as the number of principal components increases. Sometimes the cumulative variance explained is plotted as well.

# Variability of each principal component: pr.var
pr.var <- pr.out$sdev^2

# Variance explained by each principal component: pve
pve <- pr.var / sum(pr.var)

3.2.3 Visualize variance explained

Now we will create a scree plot showing the proportion of variance explained by each principal component, as well as the cumulative proportion of variance explained.

These plots can help to determine the number of principal components to retain. One way to determine the number of principal components to retain is by looking for an elbow in the scree plot showing that as the number of principal components increases, the rate at which variance is explained decreases substantially. In the absence of a clear elbow, you can use the scree plot as a guide for setting a threshold.

# Plot variance explained for each principal component
plot(pve, xlab = "Principal Component",
     ylab = "Proportion of Variance Explained",
     ylim = c(0, 1), type = "b")

# Plot cumulative proportion of variance explained
plot(cumsum(pve), xlab = "Principal Component",
     ylab = "Cumulative Proportion of Variance Explained",
     ylim = c(0, 1), type = "b")

Notice that when the number of principal components is equal to the number of original features in the data, the cumulative proportion of variance explained is 1.

3.3 Practical issues

3.3.1 Scaling

Scaling your data before doing PCA changes the results of the PCA modeling. Here, we will perform PCA with and without scaling, then visualize the results using biplots.

Sometimes scaling is appropriate when the variances of the variables are substantially different. This is commonly the case when variables have different units of measurement, for example, degrees Fahrenheit (temperature) and miles (distance). Making the decision to use scaling is an important step in performing a principal component analysis.

# Mean of each variable
colMeans(pokemon)

     HitPoints         Attack        Defense  SpecialAttack 
      69.25875       79.00125       73.84250       72.82000 
SpecialDefense          Speed 
      71.90250       68.27750 

# Standard deviation of each variable
apply(pokemon, 2, sd)

     HitPoints         Attack        Defense  SpecialAttack 
      25.53467       32.45737       31.18350       32.72229 
SpecialDefense          Speed 
      27.82892       29.06047 

# PCA model with scaling: pr.with.scaling
pr.with.scaling <- prcomp(x = pokemon_pca, scale = TRUE)

# PCA model without scaling: pr.without.scaling
pr.without.scaling <- prcomp(x = pokemon_pca, scale = FALSE)

# Create biplots of both for comparison
biplot(pr.with.scaling)

biplot(pr.without.scaling)

The new Total column contains much more variation, on average, than the other four columns, so it has a disproportionate effect on the PCA model when scaling is not performed. After scaling the data, there’s a much more even distribution of the loading vectors.