---
title: "Clustering"
author: "Vladan Devedzic"
output:
  html_document: default
  html_notebook: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## (Installing and) Loading the required R packages
```{r message=FALSE}
# install.packages("ggplot2")
library(ggplot2)
```

## Dataset

Read and the clean the dataset:  
`source("Prepare The Beatles song dataset for clustering.R")`  
`prepareTheBeatlesDatasetForClustering()`  
`# saveRDS(object = <dataframe or another R object>, file = "<filename>")  # save R object for the next session`  
`# <dataframe or another R object> <- readRDS(file = "<filename>")         # restore R object in the next session`  
```{r}
the.beatles.songs.num <- 
  readRDS("The Beatles songs dataset (numeric), v4.1.RData")
```  

## K-Means

### Outliers

Check if there are outliers in the data, using boxplots:  
`boxplot(<dataset>$<column name>, xlab = "<column name>")                    # basic boxplot for <column name>`  
`boxplot(<dataset>)                                                          # basic boxplots for all columns`  
`boxplot(<dataset>)$stats                                                    # basic boxplots for all columns, stats`  
`boxplot(<dataset>)$stats[c(1, 5), ]                                         # basic boxplots - whiskers`  
`<output var> <- ggplot(<dataset>,                       # ggplot2 boxplots`  
`+                      aes(x = "", `  
`+                          y = <column name>)) +        # show boxplot of <column name>`  
`+  geom_boxplot(width = 0.5, fill = "<color>") +        # boxplot width and color`  
`+  stat_boxplot(geom ='errorbar', width = 0.15) +       # show whiskers, control their width`  
`+  guides(fill = FALSE) +                               # no legend (it makes no sense here)`  
`+  xlab("")                                             # no x-axis label (it makes no sense here)`  

Compare outputs of boxplot(<dataset>)$stats and summary(dataset):  
```{r}
boxplot(the.beatles.songs.num)$stats
boxplot(the.beatles.songs.num$Duration, xlab = "Duration")
```  

Boxplotting all numeric features using the boxplotFeature() utility function:  
`source("Boxplotting.R")`  
`boxplotFeature(<dataset with numeric features>,   # dataframe with numeric features`  
`+              "<numeric feature>",               # numeric feature to boxplot (its name, passed as a string)`  
`+              "<color>")                         # boxplot fill color`  
```{r}
source("Boxplotting.R")
boxplotFeature(the.beatles.songs.num, "Duration", "red")
boxplotFeature(the.beatles.songs.num, "Other.releases", "chartreuse")
boxplotFeature(the.beatles.songs.num, "Covered.by", "orange")
boxplotFeature(the.beatles.songs.num, "Top.50.Billboard", "yellow")
```  

Fix the outlier values for each variable with outliers - replace each outlier value with a specific percentile value of the data, typically 90th or 95th:  
`boxplot(<dataset>$<column name>, xlab = "<column name>")    # basic boxplot for <column name>`  
`boxplot.stats(<dataset>$<column name>)                      # examine the boxplot more closely`  
`boxplot.stats(<dataset>$<column name>)$out                  # examine the outliers more closely`  
`boxplot.stats(<dataset>$<column name>)$stats[c(1, 5)]       # get the whiskers`  
`boxplot.stats(<dataset>$<column name>)$stats[1]             # get the lower whisker`  
`boxplot.stats(<dataset>$<column name>)$stats[5]             # get the upper whisker`  
`sort(boxplot.stats(<dataset>$<column name>)$out)            # get and sort the outliers`  
`<quantiles> <- quantile(<dataset>$<column name>,            # examine the 90th, 95th, ..., percentile`  
`+                       probs = seq(from = 0.9, to = 1, by = 0.025))`  
`<new max value> <-                                          # the value to replace the outliers`  
`+ as.numeric(quantile(<dataset>$<column name>,              # pick <percentile> closest to`  
`+                     probs = <percentile>))                # the upper whisker`  
`<dataset>$<column name>[<dataset>$<column name> >`  
`+                       <new max value>] <-                 # replace the outliers`  
`+ <new max value>`  
`<quantiles> <- quantile(<dataset>$<column name>,            # examine the 0th, 5th, ..., percentile`  
`+                       probs = seq(from = 0.0, to = 0.1, by = 0.025))`  
`<new min value> <-                                          # the value to replace the outliers`  
`+ as.numeric(quantile(<dataset>$<column name>,              # pick <percentile> closest to`  
`+                     probs = <percentile>))                # the lower whisker`  
`<dataset>$<column name>[<dataset>$<column name> < <new min value>] <- # replace the outliers`  
`+ <new min value>`  

Fix the outliers in the first variable manually:  
```{r}
boxplot(the.beatles.songs.num$Duration, xlab = "Duration")
boxplot.stats(the.beatles.songs.num$Duration)
boxplot.stats(the.beatles.songs.num$Duration)$out
boxplot.stats(the.beatles.songs.num$Duration)$stats[c(1, 5)]
sort(boxplot.stats(the.beatles.songs.num$Duration)$out)
quantile(the.beatles.songs.num$Duration, 
         probs = seq(from = 0.9, to = 1, by = 0.025))
new.max.duration <- 
  as.numeric(quantile(the.beatles.songs.num$Duration,         # the 92.5th percentile seems to be 
                      probs = 0.925))                         # a good cut-off point
the.beatles.songs.num$Duration[the.beatles.songs.num$Duration > new.max.duration] <- new.max.duration
quantile(the.beatles.songs.num$Duration, 
         probs = seq(from = 0, to = 0.1, by = 0.025))
new.min.duration <- 
  as.numeric(quantile(the.beatles.songs.num$Duration,         # the 2.5th percentile seems to be 
                      probs = 0.025))                         # a good cut-off point
the.beatles.songs.num$Duration[the.beatles.songs.num$Duration < new.min.duration] <- new.min.duration
boxplot(the.beatles.songs.num$Duration, xlab = "Duration")
# no more outliers in Duration
```  

Call fixOutliers() for the other columns:
```{r}
# source("Fix outliers.R")
# the.beatles.songs.num <- fixOutliers(the.beatles.songs.num, "<column name>")
source("Fix outliers.R")

boxplot(the.beatles.songs.num$Other.releases, xlab = "Other.releases")
boxplot.stats(the.beatles.songs.num$Other.releases)
boxplot.stats(the.beatles.songs.num$Other.releases)$out
boxplot.stats(the.beatles.songs.num$Other.releases)$stats[c(1, 5)]
the.beatles.songs.num <- fixOutliers(the.beatles.songs.num, "Other.releases")
boxplot(the.beatles.songs.num$Other.releases, xlab = "Other.releases")
# no more outliers in Other.releases
```  

```{r}
the.beatles.songs.num <- fixOutliers(the.beatles.songs.num, "Covered.by")
boxplot(the.beatles.songs.num$Covered.by, xlab = "Covered.by")
# no more outliers in Covered.by
```  

Demonstrate an attempt to fix outliers with highly skewed data:  
```{r}
boxplot(the.beatles.songs.num$Top.50.Billboard, xlab = "Top.50.Billboard")
boxplot.stats(the.beatles.songs.num$Top.50.Billboard)
boxplot.stats(the.beatles.songs.num$Top.50.Billboard)$out
boxplot.stats(the.beatles.songs.num$Top.50.Billboard)$stats[c(1, 5)]
temp <- the.beatles.songs.num$Top.50.Billboard # save current Top.50.Billboard, for restoring it later
the.beatles.songs.num <- fixOutliers(the.beatles.songs.num, "Top.50.Billboard")
boxplot(the.beatles.songs.num$Top.50.Billboard, xlab = "Top.50.Billboard")
the.beatles.songs.num$Top.50.Billboard <- temp # restore Top.50.Billboard
```  

Summarize the results so far:  
```{r}
summary(the.beatles.songs.num)
```  

### Patterns in the data

See if there are some patterns in the data, pairwise, to possibly indicate clusters:  
`pairs(~ <column 1 name> + <column 2 name> + ...,`  
`+     data = <dataframe>)`  
```{r}
pairs(~ Duration + Other.releases + Covered.by + Top.50.Billboard,  # no any striking pattern, i.e. 
      the.beatles.songs.num)                                        # no visual indication of clusters
```  

### K-Means with 2 variables

Try K-Means with 2 variables.  

Plot the data first:  
`<scatterplot> <-`
`+ ggplot(<dataset>, aes(x = <num.var.1>, y = <num.var.2>)) +`  
`+ geom_point(shape = <n>,         # <n> = 1: hollow circle, no fill;`  
`+                                 # <n> = 21: circle that can be filled`  
`+            fill = <color 1>,    # color of point fill (optional)`  
`+            color = <color 2>,   # color of point line (optional)`  
`+            size = <s>)          # size  of point line (optional)`  
`<scatterplot> <- <scatterplot> + xlab("<x label>")                    # label/caption on x-axis`  
`<scatterplot> <- <scatterplot> + ylab("<y label>")                    # label/caption on y-axis`  
`<scatterplot> <- <scatterplot> + ggtitle("<scatterplot title>")       # scatterplot title`  
`<scatterplot>                                                         # plot it`  

Alternatively:  
`<scatterplot> <-`  
`+  ggplot(<dataset>, aes(x = <num.var.1>, y = <num.var.2>)) +`  
`+  geom_point(shape = <n>,         # <n> = 1: hollow circle, no fill;`  
`+                                  # <n> = 21: circle that can be filled`  
`+             fill = <color 1>,    # color of point fill (optional)`  
`+             color = <color 2>,   # color of point line (optional)`  
`+             size = <s>)          # size  of point line (optional)`  
`<scatterplot> <- <scatterplot> +`  
`+ labs(x = "<x label>",                     # label/caption on x-axis`  
`+      y = "<y label>",                     # label/caption on y-axis`  
`+      title = "<scatterplot title>") +     # scatterplot title`  
`<scatterplot>                                                         # plot it`  
```{r}
scatterplot.Other.releases.vs.Covered.by <- 
  ggplot(the.beatles.songs.num, aes(x = Other.releases, y = Covered.by)) +
  geom_point(shape = 21, fill = "yellow", size = 2) + 
  labs(x = "Other.releases", y = "Covered.by", title = "Covered.by vs. Other.releases") +
  theme_bw()
scatterplot.Other.releases.vs.Covered.by
```  

Subset the original data to include only the variables to be used in K-Means:  
`<new dataframe> <- <dataframe>[, c("<col1 name>", "<col2 name>")]`  
`<new dataframe> <- <dataframe>[, <col1 index>:<col2 index>]`  
Alternatively:  
`<new dataframe> <- subset(<dataframe>, select = c("<col1 name>", "<col2 name>"))`  
`<new dataframe> <- subset(<dataframe>, select = c(<col1 index>:<col2 index>))`  
```{r}
the.beatles.songs.num.1 <- the.beatles.songs.num[, c("Other.releases", "Covered.by")]
# the.beatles.songs.num.1 <- subset(the.beatles.songs.num, select = c("Other.releases", "Covered.by"))
summary(the.beatles.songs.num.1)
head(the.beatles.songs.num.1)
```  

### Data normalization

Required by K-Means when the variables have different ranges.  
`range(<dataframe with numeric columns>$<numeric column 1> # check the range of <numeric column 1>`  
`range(<dataframe with numeric columns>$<numeric column 2> # check the range of <numeric column 2>`  
`...`  
`# install.packages("clusterSim")`  
`library(clusterSim)`  
`<dataframe with numeric columns> <-                       # works with vectors and matrices as well`  
`+ data.Normalization(<dataframe with numeric columns>,`  
`+                    type = "n4",                         # normalization: (x - min(x)) / (max(x) - min(x))`  
`+                    normalization = "column")            # normalization by columns`  
```{r}
range(the.beatles.songs.num.1$Other.releases)
range(the.beatles.songs.num.1$Covered.by)
library(clusterSim)
the.beatles.songs.num.2 <- 
  data.Normalization(the.beatles.songs.num.1, 
                     type = "n4", 
                     normalization = "column")
tail(the.beatles.songs.num.2)
```  

### K-Means for K = 3

Run K-Means for K = 3:  
`set.seed(<seed>)`  
`<clusters> <- kmeans(x = <normalized dataframe>,`  
`+                    centers = 3,                         # K = 3`  
`+                    iter.max = 20,                       # max number of iterations allowed`  
`+                    nstart = 1000)                       # no. of initial configurations`  
`+                                                         # (report generated based on the best one)`  
`<clusters>`  
```{r}
set.seed(888)
clusters.K3 <- kmeans(x = the.beatles.songs.num.2, centers = 3, iter.max = 20, nstart = 1000)
clusters.K3
```  

The meaning of parameters in the report:  

* withinss (within_SS) - within cluster sum of squares, i.e., sum of squared differences between individual data points in a cluster and the cluster center; it is computed for each cluster  
* totss (total_SS) - the sum of squared differences of each data point to the global sample mean  
* betweenss (between_SS) - the sum of squared differences of each cluster center to the global sample mean; the squared difference of each cluster center to the global sample mean is multiplied by the number of data points in that cluster  
* tot.withinss - the sum of squared differences between data points and cluster centers (the sum of within_SS for all the clusters)  
* between_SS / totat_SS - indicates how well the sample splits into clusters; the higher the ratio, the better clustering  

Add the vector of clusters to the dataframe:  
`<normalized dataframe>$<new column> <- factor(<clusters>$cluster)   # <clusters>: from the previous step`  
`head(<normalized dataframe>)`  
```{r}
the.beatles.songs.num.2$Cluster <- factor(clusters.K3$cluster)
head(the.beatles.songs.num.2)
```  

Plot the clusters in a new scatterplot, using plotClusters() utility function:  
`source("Plot clusters.R")`  
`plotClusters() <- function(dataset,           # dataset with the cluster column`  
`+                          xcol,              # dataset column for the x-axis, passed as a string`  
`+                          ycol,              # dataset column for the y-axis, passed as a string`  
`+                          clustercol,        # dataset column showing the clusters, passed as a string`  
`+                          title,             # plot title`  
`+                          x.label,           # x-axis label`  
`+                          y.label,           # y-axis label`  
`+                          legend.label,      # plot legend label`  
`+                          show.centers.flag, # plot cluster centers if TRUE`  
`+                          clusters) {        # clusters computed by kmeans() in a previous step`  
```{r}
source("Plot clusters.R")
plotClusters(the.beatles.songs.num.2,
             "Other.releases",
             "Covered.by",
             "Cluster",
             "Clusters: (Other.releases, Covered.by), normalized",
             "Other.releases",
             "Covered.by",
             "Cluster",
             TRUE,
             clusters.K3)
```  

### Optimal value for K
Find the optimal value for K, using the Elbow method (a call to the appropriate utility function); an appropriate dataframe should be passed as the parameter (just numeric variables, no clusters):  
`source("Elbow method.R")`  
`<elbow parameters> <- getElbowMethodParameters(<dataframe>[, c(<n1>, <n2>, ...)])   # leave out the cluster column`  
`<elbow parameters>`  
`plotElbow(<elbow parameters>)`  
```{r}
source("Elbow method.R")
elbow.2 <- getElbowMethodParameters(the.beatles.songs.num.2[, c(1,2)])   # remove the Cluster column
elbow.2
plotElbow(elbow.2)
```  

Show differences in tot.withinss for different values of K more precisely, using the getDifferences() utility function:  
`source("Elbow method.R")`  
`<diff dataframe> <-`  
`+ data.frame(K = <n1>:<n2>,`  
`+            diff.tot.withinss =`  
`+              getDifferences(<elbow stats>$<tot.withinss>),                 # from the previous step`  
`+            diff.ratio =`  
`+              getDifferences(<elbow stats>$<ratio between_SS / total_SS>))  # from the previous step`  
`names(<diff dataframe>) <- c("K", "Difference in tot.withinss", "Difference in between_SS / total_SS")`  
`<diff dataframe>`  
```{r}
df.differences <- data.frame(K = 2:8, 
                             getDifferences(elbow.2[, 2]), 
                             getDifferences(elbow.2[, 3]))
names(df.differences) <- c("K", "Difference in tot.withinss", "Difference in between_SS / total_SS")
df.differences
```  

### K-Means for K = 4

Run K-Means also for K = 4:  
```{r}
set.seed(818)
the.beatles.songs.num.2 <- the.beatles.songs.num.2[, -3]    # remove the Cluster column for the new run
clusters.K4 <- kmeans(x = the.beatles.songs.num.2, centers = 4, iter.max = 20, nstart = 1000)
clusters.K4
the.beatles.songs.num.2$Cluster <- factor(clusters.K4$cluster)
head(the.beatles.songs.num.2)
source("Plot clusters.R")
plotClusters(the.beatles.songs.num.2,
             "Other.releases",
             "Covered.by",
             "Cluster",
             "Clusters: (Other.releases, Covered.by), normalized",
             "Other.releases",
             "Covered.by",
             "Cluster",
             TRUE,
             clusters.K4)
```  

Examine clusters more closely by looking into the cluster centers (means) and standard deviations from the centers. In doing so, use 'regular' (not normalized) features and the summarizeClusterStats() utility function:  
`source("Summary statistics about clusters.R")`  
`<stats dataframe> <- `  
`+ summarizeClusterStats(feature.set = <dataframe with 'regular' (not normalized) features>, `  
`+                       clusters = <clusters>$cluster,   # <clusters>: result of kmeans()`  
`+                       cl.num = <K>)                    # <K>: number of clusters`  
Compare dispersion (sd) in this representation of the data and in the corresponding plots.`  
```{r}
source("Summary statistics about clusters.R")
clusters.K3.stats <- summarizeClusterStats(feature.set = the.beatles.songs.num.1, 
                                           clusters = clusters.K3$cluster, cl.num = 3)
clusters.K4.stats <- summarizeClusterStats(feature.set = the.beatles.songs.num.1, 
                                           clusters = clusters.K4$cluster, cl.num = 4)
```  

### K-Means with all variables 

Data normalization:  
```{r}
library(clusterSim)
the.beatles.songs.num.4 <- 
  data.Normalization(the.beatles.songs.num, 
                     type = "n4", 
                     normalization = "column")
tail(the.beatles.songs.num.4)
```  

Find the optimal value for K, using the Elbow method (a call to the appropriate utility function):  
```{r}
source("Elbow method.R")
elbow.4 <- getElbowMethodParameters(the.beatles.songs.num.4)
elbow.4
plotElbow(elbow.4)
```  

Show differences in tot.withinss for different values of K more precisely, using the getDifferences() utility function:  
```{r}
df.differences <- data.frame(K = 2:8, 
                             getDifferences(elbow.4[, 2]), 
                             getDifferences(elbow.4[, 3]))
names(df.differences) <- c("K", "Difference in tot.withinss", "Difference in between_SS / total_SS")
df.differences
```  

Run K-Means for K = 3, since K = 3 seems to be the best value for K:  
```{r}
set.seed(888)
clusters.K3.all.vars <- kmeans(x = the.beatles.songs.num.4, centers = 3, iter.max = 20, nstart = 1000)
clusters.K3.all.vars
```  

Run K-Means also for K = 4, since K = 4 seems to be the next best value for K:  
```{r}
set.seed(888)
clusters.K4.all.vars <- kmeans(x = the.beatles.songs.num.4, centers = 4, iter.max = 20, nstart = 1000)
clusters.K4.all.vars
```  

Examine and compare the cluster centers for K = 3 and K = 4:  
```{r}
clusters.K3.all.vars.stats <- summarizeClusterStats(feature.set = the.beatles.songs.num, 
                                                    clusters = clusters.K3.all.vars$cluster, cl.num = 3)
clusters.K3.all.vars.stats

clusters.K4.all.vars.stats <- summarizeClusterStats(feature.set = the.beatles.songs.num, 
                                                    clusters = clusters.K4.all.vars$cluster, cl.num = 4)
clusters.K4.all.vars.stats
```  

### Compare multiple clustering results/schemes

`# install.packages("fpc")`  
`library(fpc)`  
`?cluster.stats`  
`<comparison criteria> <-              # specify criteria (from cluster.stats()) for comparing`  
`+ c("<criterion 1>",                  # different clusterings (e.g., "max.diameter", "min.separation",`  
`+   "<criterion 2>", ...)             # "average.between", "average.within", "within.cluster.ss", ...)`  
`<distance matrix> <-`  
`+ dist(x = <normalized dataset>)`  
`<comparison> <- sapply(list(<clustering 1 name> = <clustering 1>$cluster,    # <clustering 1> computed by kmeans()`  
`+                           <clustering 2 name> = <clustering 2>$cluster,    # <clustering 2> computed by kmeans()`  
`+                           ...)`  
`+                      FUN = function(x)`  
`+                        cluster.stats(<distance matrix>, x))[<comparison criteria>, ]`  
Alternative 1 - show output as a table in the console:  
`# install.packages("knitr")`  
`library(knitr)`  
`kable(x = <comparison>, format = "rst")`  
Alternative 2 - show output as a dataframe, using an appropriate utility function:  
`source("Summary statistics about clusters.R")`  
`<comparison.df> <- `  
`+ compareMultipleClusterings(<comparison>)      # show comparison as a dataframe, using a corresponding utility function`  
`<comparison.df>`  
```{r}
library(fpc)
comparison.criteria <- c("max.diameter", "min.separation", "average.between", 
                         "average.within", "within.cluster.ss")  
d <- dist(x = the.beatles.songs.num.4)
comparison <- sapply(list(c.K3.var2 = clusters.K3$cluster, # clustering: 3 clusters, 2 variables
                          c.K4.var2 = clusters.K4$cluster, # clustering: 4 clusters, 2 variables
                          c.K3.var4 = clusters.K3.all.vars$cluster,  # 3 clusters, 4 variables
                          c.K4.var4 = clusters.K4.all.vars$cluster), # 4 clusters, 4 variables
                     FUN = function(x) cluster.stats(d, x))[comparison.criteria, ]

# Alternative 1 - show output as a table in the console: 
# library(knitr)
# kable(x = comparison, format = "rst")

# Alternative 2 - show output as a dataframe, using an appropriate utility function: 
comparison.df <- compareMultipleClusterings(comparison)
comparison.df
```  
