Clustering.Rmd

---
title: "CLustering Unsupervised Machine Learning"
output:
  html_document: default
  html_notebook: default
  pdf_document: default
---

```{r}
#Load the data
data("USArrests")
mydata <- USArrests

#Remove mising values
mydata <- na.omit(mydata)

#Scale the variables
mydata <- scale(mydata)

head(mydata, n=5)
```
```{r}
#Data preparation example

set.seed(124)
ss <- sample(1:50,10)
df <- USArrests[ss, ]
df <- na.omit(df)
head(df,n=6)

df.scaled <- scale(df)
head(round(df.scaled, 2))
```


```{r}
#install.packages("cluster")

#For visualization
#if(!require(devtools))
#  install.packages("devtools")
#devtools::install_github("kassambara/factoextra")

#install.packages("factoextra")

#library(cluster)
#library(factoextra)

#install.packages("corrplot")
```

```{r}
#second argument take value where to apply the function 
#1. for function on rows
#2. for function on columns

desc_stats <- data.frame(
  Min = apply(USArrests, 2, min),
  Max = apply(USArrests, 2, max),
  Med = apply(USArrests, 2, median),
  SD = apply(USArrests, 2, sd),
  Mean = apply(USArrests, 2, mean)
)

desc_stats <- round(desc_stats,1)
head(desc_stats)
```

```{r}
#Functions for computing distances
#dist()#############################################################
library(stats)
eucl <- dist(df.scaled, method = "euclidean" )

#method can be euclidean, manhattan, correlation etc
round(as.matrix(eucl)[1:6,1:6],1)
#t() used for transposing data

#correlation based distance#########################################
#Cor computes coefficient between variables but we need observations 
#so using t() to transpose matrix

print("Correlation method for distance measure")
cor <- cor(t(df.scaled), method = "pearson")
dist_cor <- as.dist(1 - cor)
round(as.matrix(dist_cor)[1:6,1:6],1)

#daisy() to compute dissimilarity matrices between observations 

daisy(df.scaled, metric = c("euclidean", "manhattan", "gower"), stand = FALSE)

#stand: if TRUE, then the measurements in df.scaled are standardized before calculating the dissimilarities. Measurements are standardized for each variable (column), by subtracting the variable’s mean value and dividing by the variable’s mean absolute deviation

data("flower")
head(flower)

str(flower)

daisy_dist <- as.matrix(daisy(flower))
head(round(daisy_dist[1:6,1:6]),2)

#visualizing distance matrices

library(corrplot)
corrplot(as.matrix(eucl), is.corr = FALSE, method = "color")

corrplot(as.matrix(eucl), is.corr = FALSE, method = "color", order = "hclust", type = "upper")

#hierarchial clustering dendrogram 
plot(hclust(eucl, method = "ward.D2"))

#Heatmap
heatmap(as.matrix(eucl), symm = TRUE, distfun = function(x) as.dist(x))
```

```{r}
#Visualization ofdistance matrices using functions from factoextra
pearson.dist <- get_dist(USArrests, stand = TRUE, method = "pearson")
#for computing distance matrix between rows of a data mtrix 

fviz_dist(pearson.dist, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

```

```{r}
#kmeans(x, centers, iter.max = 10, nstart = 1)

#x: numeric matrix, numeric data frame or a numeric vector
#centers: Possible values are the number of clusters (k) or a set of initial (distinct) cluster centers. If a number, a random set of (distinct) rows in x is chosen as the initial centers.
#iter.max: The maximum number of iterations allowed. Default value is 10.
#nstart: The number of random starting partitions when centers is a number. Trying nstart > 1 is often recommended.

#kmeans() function returns a list including:

#cluster: A vector of integers (from 1:k) indicating the cluster to which each point is allocated
#centers: A matrix of cluster centers (cluster means)
#totss: The total sum of squares (TSS), i.e ∑(xi−x¯)2∑(xi−x¯)2. TSS measures the total variance in the data.
#withinss: Vector of within-cluster sum of squares, one component per cluster
#tot.withinss: Total within-cluster sum of squares, i.e. sum(withinss)sum(withinss)
#betweenss: The between-cluster sum of squares, i.e. totss−tot.withinsstotss−tot.withinss
#size: The number of observations in each cluster

set.seed(123)
df2 <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2),
             matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))
colnames(df2) <- c("x","y")
head(df2)

set.seed(123)
km <- kmeans(df2, 2, nstart = 25)
km$cluster

km$size
km$centers

plot(df2, col = km$cluster, pch = 19, frame = FALSE)
points(km$cluster, col = 1:2, pch = 8, cex = 3)


set.seed(123)
km2 <- kmeans(USArrests, 4, nstart = 25)
plot(USArrests, col = km2$cluster, pch = 19)
points(km2$cluter, col = 1:4, pch = 8, cex = 3)
```
```{r}
#k means on USArrests data
library(factoextra)
library(ggplot2)

data("USArrests")
dataf <- na.omit(USArrests)
dataf <- scale(dataf)
set.seed(123)
#To find out how many cluster should be suggested 
fviz_nbclust(dataf, kmeans, method = "wss") + geom_vline(xintercept = 4, linetype = 2)

kmean <- kmeans(dataf, 4, nstart = 25)
print(kmean)

#compute mean of each variables
aggregate(USArrests, by = list(cluster=kmean$cluster), mean)

fviz_cluster(kmean , data = dataf)
```

```{r}
#Partitioning around medoids method(k-mediods)
#The pam algorithm is based on the search for k representative objects or medoids among the observations of the dataset. These observations should represent the structure of the data. After finding a set of k medoids, k clusters are constructed by assigning each observation to the nearest medoid. The goal is to find k representative objects which minimize the sum of the dissimilarities of the observations to their closest representative object.


library(cluster)
library(fpc)

data("USArrests")
pam.res <- pam(scale(USArrests),4)

pam.res$medoids

head(pam.res$cluster)
clusplot(pam.res, color= TRUE)
fviz_cluster(pam.res)

plot(silhouette(pam.res), col = 2:5)
fviz_silhouette(silhouette(pam.res))

#computing silhouette
sil <- silhouette(pam.res)[,1:3]
neg_sil <- which(sil[,'sil_width'] < 0)
sil[neg_sil,,drop = FALSE]
```


```{r}
#CLARA (Clustering LArge Applications)
#Algorithm : 
#Split randomly the data sets in multiple subsets with fixed size
#Compute PAM algorithm on each subset and choose the corresponding k representative objects (medoids). Assign each observation of the entire dataset to the nearest medoid.
#Calculate the mean (or the sum) of the dissimilarities of the observations to their closest medoid. This is used as a measure of the goodness of the clustering.
#Retain the sub-dataset for which the mean (or sum) is minimal. A further analysis is carried out on the final partition.

set.seed(1234)

#Generate 500 objects divided into 2 clusters
x <- rbind(cbind(rnorm(200,0,8), rnorm(200,0,8)),
           cbind(rnorm(300,50,8),rnorm(300,50,8)))

head(x)

#compute Clara

clarax <- clara(x,2,samples=50)
#fviz_cluster(clarax)
plot(silhouette(clarax), col = 2:3)
```

```{r}
#Visualization examples

set.seed(123)
# K-means clustering
km.res <- kmeans(scale(USArrests), 4, nstart = 25)
# Use clusplot function
library(cluster)
clusplot(scale(USArrests), km.res$cluster,  main = "Cluster plot",
         color=TRUE, labels = 2, lines = 0)

fviz_cluster(pam.res, frame.type = "t",
             frame.alpha = 0, frame.level = 0.7)

fviz_cluster(pam.res) + 
  scale_color_brewer(palette = "Set2")+
  scale_fill_brewer(palette = "Set2") +
  theme_minimal()

#fviz_cluster(object, data = NULL, stand = TRUE,
#             geom = c("point", "text"), 
#             frame = TRUE, frame.type = "convex")
```

```{r}
#Hierarchial Clustering

data("USArrests")
my_data <- scale(USArrests)

d <- dist(my_data, method = "euclidean")
res.hc <- hclust(d, method = "ward.D2")
grp <- cutree(res.hc,k=4)
plot(res.hc, cex = 0.6)
rect.hclust(res.hc, k = 4, border = 2:5)


library(factoextra)
res <- hcut(USArrests, k = 4, stand = TRUE)
fviz_dend(res, rect = TRUE, cex = 0.5, k_colors = c("Red", "Green", "Blue", "Yellow"))
```

```{r}
#Clustering tendecy

library("factoextra")
my_data <- scale(iris[, -5])
get_clust_tendency(my_data, n = 50,
                   gradient = list(low = "steelblue",  high = "white"))
```

```{r}
#Finding optimal number of clusters

library(NbClust)
data("iris")
iris.scaled <- scale(iris[,-5])

## K means method
set.seed(123)
km_res <- kmeans(iris.scaled,3,nstart=25)
km_res$cluster

fviz_cluster(km_res,data=iris.scaled,geom="point", stand = FALSE, frame.type = "norm")

library(cluster)

##PAM cluster method

pam_res <- pam(iris.scaled, 3)
fviz_cluster(pam_res, stand = FALSE, geom="point", frame.type = "norm")

##Hierarchial Cluster
dist_res <- dist(iris.scaled, method = "euclidean")
h_c <- hclust(dist_res, method = "ward.D2")
plot(h_c, labels = FALSE, hang = -1)
rect.hclust(h_c, k = 3, border = 2:4)

h_c.cut <- cutree(h_c, k = 3)
head(h_c.cut, 20)

############### Elbow Method for k means clustering #####################
set.seed(123)
k.max <- 15
data <- iris.scaled
wss <- sapply(1:k.max, 
              function(k){kmeans(data,k,nstart=10)$tot.withinss})

plot(1:k.max, wss,
     type = "b", pch = 19, frame = FALSE)
abline(v = 3, lty = 2)

fviz_nbclust(iris.scaled, kmeans, method = "wss") + geom_vline(xintercept = 3, linetype = 2)


############### Elbow Method for PAM clustering #####################
fviz_nbclust(iris.scaled, pam, method = "wss") + geom_vline(xintercept = 3, linetype = 2)

############### Elbow Method for hierarchial clustering #####################
fviz_nbclust(iris.scaled, hcut, method = "wss") + geom_vline(xintercept = 3, linetype = 2)

############### Average Silhouette Method for k means clustering#################

sil <- rep(0, k.max)
# Compute the average silhouette width for 
# k = 2 to k = 15

for(i in 2:k.max){
  km.res <- kmeans(data, centers = i, nstart = 25)
  ss <- silhouette(km.res$cluster, dist(data))
  sil[i] <- mean(ss[, 3])}

# Plot the  average silhouette width
plot(1:k.max, sil, type = "b", pch = 19, 
     frame = FALSE, xlab = "Number of clusters k")
abline(v = which.max(sil), lty = 2)

fviz_nbclust(iris.scaled, kmeans, method = "silhouette")
fviz_nbclust(iris.scaled, pam, method = "silhouette")
fviz_nbclust(iris.scaled, hcut, method = "silhouette", hc_method = "complete")


##############Gap Statistic method ##################################

gap_stat <- clusGap(iris.scaled, FUN = kmeans, nstart = 25,
                    K.max = 10, B = 50)
# Print the result
print(gap_stat, method = "firstmax")

plot(gap_stat, frame = FALSE, xlab = "Number of clusters k")
abline(v = 3, lty = 2)

fviz_gap_stat(gap_stat, 
          maxSE = list(method = "Tibs2001SEmax", SE.factor = 2))

gap_stat <- clusGap(iris.scaled, FUN = pam, K.max = 10, B = 50)
# Plot gap statistic
fviz_gap_stat(gap_stat)

gap_stat <- clusGap(iris.scaled, FUN = hcut, K.max = 10, B = 50)
# Plot gap statistic
fviz_gap_stat(gap_stat)

```

```{r}
#To choose appropriate clustering algo for your data

library(clValid)

intern <- clValid(my_data, nClust = 2:6, 
              clMethods = c("hierarchical","kmeans","pam"),
              validation = "internal")

summary(intern)

library(pvclust)
set.seed(123)
data("lung")
ss <- sample(1:73,30)
my_data2 <- lung[,ss]

res.pv <- pvclust(my_data2, method.dist="cor", 
                  method.hclust="average", nboot = 10)

plot(res.pv, hang = -1, cex = 0.5)
pvrect(res.pv)
```