Exercise Solution for 5.1

The inner lapply function

The first lapply is generating a vector of n = 100 normally distributed random numbers, and creating two separate dataframes packed into a list that consist of the mean (my) and standard deviation, respectively. Each individual value is specifically assigned a 0 or an 8.

simdatmy = lapply(c(0,8), function(my) {
    
    tibble(y = rnorm(100, mean = my, sd = 2),
           class = paste(my, sep = ":"))
   })
summary(simdatmy)

##      Length Class  Mode
## [1,] 2      tbl_df list
## [2,] 2      tbl_df list

The outer lapply function

The second (outer) lapply uses the same idea to apply the same, random, 0 or 8 assignment to values in the mx function. The ouput is now four separate dataframes within a list that contain all of the mx data and all of the my data. Within the tibble function, they include the code class = to ensure that each row in each of the 4 the lists is assigned one of the four possible two-way combinations of 0 and 8. This is important to simulate a clustered dataset.

simdatmx = lapply(c(0, 8), function(mx) {
  lapply(c(0,8), function(my) {
    
    tibble(x = rnorm(100, mean = mx, sd = 2),
           y = rnorm(100, mean = my, sd = 2),
           class = paste(mx, my, sep = ":"))
    })})
summary(simdatmx)

##      Length Class  Mode
## [1,] 2      -none- list
## [2,] 2      -none- list

Putting it together

The last step is to bind the list of dataframes into one single dataframe. The final dataframe includes all of the x and y data, each with assigned classes, defined by a combination of 0 and 8.

set.seed(1)


simdat = lapply(c(0, 8), function(mx) {
  lapply(c(0,8), function(my) {
    
    tibble(x = rnorm(100, mean = mx, sd = 2),
           y = rnorm(100, mean = my, sd = 2),
           class = paste(mx, my, sep = ":"))
   }) %>% bind_rows
}) %>% bind_rows

head(simdat)

## # A tibble: 6 x 3
##        x       y class
##    <dbl>   <dbl> <chr>
## 1 -1.25  -1.24   0:0  
## 2  0.367  0.0842 0:0  
## 3 -1.67  -1.82   0:0  
## 4  3.19   0.316  0:0  
## 5  0.659 -1.31   0:0  
## 6 -1.64   3.53   0:0

unique(simdat$class)

## [1] "0:0" "0:8" "8:0" "8:8"

The final simdat dataframe includes 400 random points witih an assigned class to simulate clustering. We can look at the data using a simple ggplot scatterplot, color coded by the class of each point.

ggplot(simdat, aes(x = x, y = y, color = class)) + geom_point()

The next part of exploring the data is to compute the within-groups sum of squares (WSS) for the clusters that we just generated. The goal of this section is to observe how the WSS changes as the number of clusters is increased from 1 to 8 when using the k-means. Chapter 5 provides us with the following code and graph:

#1
wss = tibble(k = 1:8, value = NA_real_)

#2
wss$value[1] = sum(scale(simdatxy, scale = FALSE)^2)

#3
for (i in 2:nrow(wss)) {
  km  = kmeans(simdatxy, centers = wss$k[i])
  wss$value[i] = sum(km$withinss)
}

ggplot(wss, aes(x = k, y = value)) + geom_col()

What is really going on here?

This first chunk is setting up a one-column dataframe with blank NA values. The NAs will be filled in with values as the rest of the code processes.

#1
wss = tibble(k = 1:8, value = NA_real_)
wss

## # A tibble: 8 x 2
##       k value
##   <int> <dbl>
## 1     1    NA
## 2     2    NA
## 3     3    NA
## 4     4    NA
## 5     5    NA
## 6     6    NA
## 7     7    NA
## 8     8    NA

The second chunk of code is calculating the value for k = 1 individually. They use the scale function to scale down the value for k = 1 because it is so much larger than the rest of the k-values. Without scaling down the k = 1 value, it would be difficult to observe any sharp decreases that might indicate a “potential sweet spot” for the number of clusters.

#2
wss$value[1] = sum(scale(simdatxy, scale = FALSE)^2)
wss

## # A tibble: 8 x 2
##       k  value
##   <int>  <dbl>
## 1     1 15781.
## 2     2    NA 
## 3     3    NA 
## 4     4    NA 
## 5     5    NA 
## 6     6    NA 
## 7     7    NA 
## 8     8    NA

The last part of this chunk is running a k-means clustering on the remaining k 2 through 8 and then pulling out the withinss value for all of the observations, summing it, and assigning that value to each individual k-value.

#3
for (i in 2:nrow(wss)) {
  km  = kmeans(simdatxy, centers = wss$k[i])
  wss$value[i] = sum(km$withinss)
}

km$withinss

## [1] 275.2533 165.2368 199.8292 257.9090 285.8292 131.9047 239.8457 339.1308

wss

## # A tibble: 8 x 2
##       k  value
##   <int>  <dbl>
## 1     1 15781.
## 2     2  9055.
## 3     3  5683.
## 4     4  3088.
## 5     5  2755.
## 6     6  2441.
## 7     7  2152.
## 8     8  1895.

These corresponding values are then neatly displayed in a barchart of the WSS stastistic as a function of k. The sharp decrease between k = 3 and k = 4 (at the elbow) is indicative of the number of clusters present in the dataset.

ggplot(wss, aes(x = k, y = value)) + geom_col()

Computing the silhouette index for simdat

Next up is the code necessary to plot the silhouette index. The silhouette function comes from the cluster package, and the resulting graph provides an average silhouette width for k = 4 clusters.

pam4 = pam(simdatxy, 4)
sil = silhouette(pam4, 4, border = NA)

The pam (partitioning around medoids) function is doing the same thing as the kmeans call from the earlier chunk of code, but using the cluster package’s algorithm to calculate the k-means clustering. We use the pam function here because the we need the “pam” and “partition” output class to run the silhouette function. With this information, we can then compute the silhouette index and view the output summary and plot.

class(pam4)

## [1] "pam"       "partition"

sil = silhouette(pam4, 4, border = NA)
summary(sil)

## Silhouette of 400 units in 4 clusters from pam(x = simdatxy, k = 4) :
##  Cluster sizes and average silhouette widths:
##       103       100        99        98 
## 0.5279715 0.5047895 0.4815427 0.4785642 
## Individual silhouette widths:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.04754  0.41232  0.54916  0.49858  0.63554  0.71440

Before we explore the final output plot, it might be interesting to look at plots of the simulated values with their respective cluster assignments based on pam k-means clustering and the silhouette index. With some (a lot of) help from Brooke, we have the following code to view this.

For the most part, all of the points were assigned to the same cluster as the original, with the occational border point mis-assigned to the neighboring cluster. Interestingly, the silhouette index approaches zero when you near the border of of the cluster and is much higher near the center of the cluster. Although we would expect this, it can be helpful to view this graphically.

sil %>% 
  unclass() %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column(var = "orig_order") %>% 
  arrange(as.numeric(orig_order)) %>% 
  bind_cols(simdat) %>% 
  ggplot(aes(x = x, y = y, shape = as.factor(cluster), color = sil_width)) + 
  geom_point() + 
  facet_wrap(~ class)

And finally, a silhouette plot with the n = 400 data points. The average silhouette width is a metric that we can use to summarize everything at a level of the full clustering process. Essentially, the closer that this average is to 0.5, then the more accurate our number of clusters k is. This concept is further explored in the Part B.

plot(sil, col=c("red","green","blue","purple"), main="Silhouette", border = "NA")