---
title: "Comprehensive Guide to Clustering Validation in R"
author: "Clustering Validation Tutorial"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: show
    theme: cosmo
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, 
                      fig.width = 8, fig.height = 6)
```

# Introduction

Clustering validation is essential for evaluating the quality of clustering results. This tutorial demonstrates various clustering validation metrics in R using multiple packages and the famous `iris` dataset.

**Key concepts covered:**

- Internal validation measures (Silhouette, Dunn Index, Connectivity)
- Determining optimal number of clusters
- Stability measures
- Visual validation methods

# Load Required Packages

```{r load-packages}
# Install packages if needed
# install.packages(c("cluster", "factoextra", "clValid", "NbClust", "ggplot2"))

library(cluster)      # For silhouette, pam, and clustering functions
library(factoextra)   # For visualization functions
library(clValid)      # For comprehensive validation
library(NbClust)      # For multiple indices at once
library(ggplot2)      # For additional plotting
```

# Prepare Data

We'll use the iris dataset, removing the species column for unsupervised clustering.

```{r prepare-data}
# Load and prepare data
data("iris")
head(iris)

# Remove species column and scale the data
iris_data <- iris[, -5]
iris_scaled <- scale(iris_data)

# Check the scaled data
head(iris_scaled)
```

# Part 1: Determining Optimal Number of Clusters

## 1.1 Elbow Method (Within-Cluster Sum of Squares)

The elbow method plots the within-cluster sum of squares (WSS) against the number of clusters. The "elbow" point indicates a good number of clusters.

```{r elbow-method}
# Elbow method using factoextra
fviz_nbclust(iris_scaled, kmeans, method = "wss") +
  geom_vline(xintercept = 3, linetype = 2, color = "red") +
  labs(title = "Elbow Method for Optimal k",
       subtitle = "Look for the 'elbow' in the curve",
       x = "Number of clusters (k)",
       y = "Total within-cluster sum of squares") +
  theme_minimal()
```

**Interpretation:** The plot shows the total WSS decreases as k increases. Look for the point where the decrease slows down dramatically (the "elbow"). For iris data, this appears around k=3.

## 1.2 Silhouette Method

The silhouette method measures how similar an object is to its own cluster compared to other clusters. Higher average silhouette width indicates better clustering.

```{r silhouette-optimal}
# Silhouette method
fviz_nbclust(iris_scaled, kmeans, method = "silhouette") +
  labs(title = "Silhouette Method for Optimal k",
       subtitle = "Higher average silhouette width = better clustering",
       x = "Number of clusters (k)",
       y = "Average Silhouette Width") +
  theme_minimal()
```

**Interpretation:** The peak of the curve indicates the optimal number of clusters. Silhouette values range from -1 to 1:

- Close to 1: Well-clustered observations
- Close to 0: Borderline between two clusters
- Negative: Likely assigned to wrong cluster

## 1.3 Gap Statistic Method

The gap statistic compares the within-cluster variation to a null reference distribution of the data.

```{r gap-statistic}
# Compute gap statistic
set.seed(123)
gap_stat <- clusGap(iris_scaled, FUN = kmeans, nstart = 25, 
                    K.max = 10, B = 50)

# Print results
print(gap_stat, method = "firstmax")

# Visualize gap statistic
fviz_gap_stat(gap_stat) +
  labs(title = "Gap Statistic Method for Optimal k",
       subtitle = "Choose k where gap statistic is largest") +
  theme_minimal()
```

**Interpretation:** The optimal number of clusters is where the gap statistic is maximized. The method compares the observed within-cluster variation to that expected under a null reference distribution.

## 1.4 NbClust: Multiple Indices Simultaneously

NbClust provides 30 different indices for determining optimal cluster number in one function call.

```{r nbclust, results='hide', fig.show='hide'}
# Run NbClust with multiple indices
# Note: This may take some time as it computes 30 different indices
set.seed(123)
nb_result <- NbClust(data = iris_scaled, 
                     distance = "euclidean",
                     min.nc = 2, 
                     max.nc = 10, 
                     method = "kmeans", 
                     index = "all")
```

```{r nbclust-results}
# View the optimal number suggested by different indices
cat("Optimal number of clusters based on majority rule:\n")
print(nb_result$Best.nc[1,])

# View how many indices proposed each k
table(nb_result$Best.nc[1,])
```

**Interpretation:** NbClust applies 30 different indices and suggests the optimal k based on majority rule. This provides a comprehensive assessment across multiple validation criteria.

# Part 2: Detailed Silhouette Analysis

After determining optimal k, we perform detailed silhouette analysis.

```{r silhouette-analysis}
# Perform k-means clustering with k=3
set.seed(123)
km_result <- kmeans(iris_scaled, centers = 3, nstart = 25)

# Compute silhouette information
sil <- silhouette(km_result$cluster, dist(iris_scaled))

# Summary of silhouette
summary(sil)

# Plot silhouette
fviz_silhouette(sil) +
  labs(title = "Silhouette Plot for k=3 Clusters",
       subtitle = "Width indicates clustering quality for each observation") +
  theme_minimal()
```

## Identify Poorly Clustered Observations

```{r negative-silhouette}
# Find observations with negative silhouette values
neg_sil_index <- which(sil[, "sil_width"] < 0)

if(length(neg_sil_index) > 0) {
  cat("Observations with negative silhouette (poorly clustered):\n")
  print(sil[neg_sil_index, , drop = FALSE])
  cat("\nThese", length(neg_sil_index), 
      "observations may be assigned to the wrong cluster.\n")
} else {
  cat("No observations with negative silhouette values.\n")
}

# Find observations with low silhouette (< 0.25)
low_sil_index <- which(sil[, "sil_width"] < 0.25 & sil[, "sil_width"] >= 0)
if(length(low_sil_index) > 0) {
  cat("\nObservations with low silhouette (< 0.25):\n")
  print(head(sil[low_sil_index, , drop = FALSE]))
  cat("\nThese observations are borderline between clusters.\n")
}
```

## Compare Different Numbers of Clusters

```{r compare-k}
# Compare silhouettes for k = 2 to 6
par(mfrow = c(2, 3))
for(k in 2:6) {
  set.seed(123)
  km_temp <- kmeans(iris_scaled, centers = k, nstart = 25)
  sil_temp <- silhouette(km_temp$cluster, dist(iris_scaled))
  plot(sil_temp, main = paste("k =", k), 
       col = 2:(k+1), border = NA)
}
par(mfrow = c(1, 1))
```

# Part 3: Comprehensive Validation with clValid

The clValid package provides internal, stability, and biological validation measures.

## 3.1 Internal Validation

Internal measures include Connectivity, Dunn Index, and Silhouette Width.

```{r clvalid-internal}
# Perform internal validation
# Compare multiple clustering methods and numbers of clusters
clmethods <- c("hierarchical", "kmeans", "pam")

intern <- clValid(iris_scaled, 
                  nClust = 2:6,
                  clMethods = clmethods,
                  validation = "internal")

# Summary of results
summary(intern)
```

```{r clvalid-internal-plot, fig.height=8}
# Visualize internal measures
par(mfrow = c(3, 1))
plot(intern, legend = FALSE)
par(mfrow = c(1, 1))

# Get optimal scores
optimalScores(intern)
```

**Interpretation of Internal Measures:**

- **Connectivity**: Measures how well observations are placed with their nearest neighbors. **Minimize** this value (range: 0 to ∞).
- **Dunn Index**: Ratio of smallest inter-cluster distance to largest intra-cluster distance. **Maximize** this value (range: 0 to ∞).
- **Silhouette Width**: Average silhouette value across all observations. **Maximize** this value (range: -1 to 1).

## 3.2 Stability Validation

Stability measures evaluate how consistent clustering is when data is perturbed (removing one sample at a time).

```{r clvalid-stability}
# Perform stability validation
# This takes longer as it requires resampling
stab <- clValid(iris_scaled, 
                nClust = 2:6,
                clMethods = clmethods,
                validation = "stability")

# Summary
summary(stab)

# Optimal scores
optimalScores(stab)
```

**Interpretation of Stability Measures:**

- **APN (Average Proportion of Non-overlap)**: Average proportion of observations not placed in same cluster. **Minimize** (range: 0 to 1).
- **AD (Average Distance)**: Average distance between observations in same cluster. **Minimize**.
- **ADM (Average Distance between Means)**: Average distance between cluster centers. **Minimize**.
- **FOM (Figure of Merit)**: Average intra-cluster variance. **Minimize**.

# Part 4: Dunn Index Calculation

The Dunn Index is the ratio of the minimum inter-cluster distance to the maximum intra-cluster distance.

```{r dunn-index}
# Calculate Dunn Index for different k values
dunn_values <- sapply(2:10, function(k) {
  set.seed(123)
  km <- kmeans(iris_scaled, centers = k, nstart = 25)
  dunn(distance = dist(iris_scaled), clusters = km$cluster)
})

# Plot Dunn Index
dunn_df <- data.frame(k = 2:10, dunn = dunn_values)

ggplot(dunn_df, aes(x = k, y = dunn)) +
  geom_line(color = "steelblue", size = 1) +
  geom_point(color = "steelblue", size = 3) +
  geom_vline(xintercept = which.max(dunn_values) + 1, 
             linetype = "dashed", color = "red") +
  labs(title = "Dunn Index for Different Number of Clusters",
       subtitle = "Higher values indicate better clustering",
       x = "Number of clusters (k)",
       y = "Dunn Index") +
  theme_minimal()

cat("Optimal k based on Dunn Index:", which.max(dunn_values) + 1, "\n")
cat("Maximum Dunn Index value:", max(dunn_values), "\n")
```

# Part 5: PAM Clustering with Validation

Partitioning Around Medoids (PAM) is more robust to outliers than k-means.

```{r pam-clustering}
# Perform PAM clustering
set.seed(123)
pam_result <- pam(iris_scaled, k = 3)

# Summary
summary(pam_result)

# Visualize PAM clustering
fviz_cluster(pam_result, 
             data = iris_scaled,
             ellipse.type = "convex",
             palette = "jco",
             ggtheme = theme_minimal(),
             main = "PAM Clustering (k=3)")
```

```{r pam-silhouette}
# Silhouette plot for PAM
fviz_silhouette(pam_result) +
  labs(title = "Silhouette Plot for PAM Clustering (k=3)") +
  theme_minimal()
```

# Part 6: Hierarchical Clustering with Validation

```{r hierarchical}
# Perform hierarchical clustering
set.seed(123)
hc_result <- hcut(iris_scaled, k = 3, hc_method = "ward.D2")

# Visualize dendrogram
fviz_dend(hc_result, 
          rect = TRUE,
          cex = 0.5,
          main = "Hierarchical Clustering Dendrogram (k=3)",
          k_colors = c("#E7B800", "#00AFBB", "#FC4E07"))
```

```{r hierarchical-silhouette}
# Silhouette for hierarchical clustering
fviz_silhouette(hc_result) +
  labs(title = "Silhouette Plot for Hierarchical Clustering (k=3)") +
  theme_minimal()
```

# Part 7: Cluster Visualization and Comparison

```{r cluster-comparison, fig.height=10}
# Compare three clustering methods side by side
set.seed(123)

# K-means
km3 <- kmeans(iris_scaled, centers = 3, nstart = 25)
p1 <- fviz_cluster(list(data = iris_scaled, cluster = km3$cluster),
                   geom = "point",
                   ellipse.type = "convex",
                   palette = "jco",
                   main = "K-means Clustering") +
  theme_minimal()

# PAM
pam3 <- pam(iris_scaled, k = 3)
p2 <- fviz_cluster(pam3,
                   geom = "point",
                   ellipse.type = "convex",
                   palette = "jco",
                   main = "PAM Clustering") +
  theme_minimal()

# Hierarchical
hc3 <- hcut(iris_scaled, k = 3, hc_method = "ward.D2")
p3 <- fviz_cluster(list(data = iris_scaled, cluster = hc3$cluster),
                   geom = "point",
                   ellipse.type = "convex",
                   palette = "jco",
                   main = "Hierarchical Clustering") +
  theme_minimal()

# Arrange plots
gridExtra::grid.arrange(p1, p2, p3, ncol = 1)
```

# Part 8: Comparing Results to True Labels

Since iris has known species labels, we can compare clustering results.

```{r compare-true-labels}
# Compare k-means results to true species
table(Predicted = km3$cluster, Actual = iris$Species)

# Calculate adjusted Rand index
library(mclust)
adjustedRandIndex(km3$cluster, iris$Species)

# For PAM
cat("\nAdjusted Rand Index for PAM:", 
    adjustedRandIndex(pam3$clustering, iris$Species), "\n")

# For Hierarchical
cat("Adjusted Rand Index for Hierarchical:", 
    adjustedRandIndex(hc3$cluster, iris$Species), "\n")
```

**Interpretation:** Adjusted Rand Index (ARI) measures agreement between two clusterings, corrected for chance:

- ARI = 1: Perfect agreement
- ARI = 0: Random clustering
- ARI < 0: Worse than random

# Summary and Best Practices

## Key Findings for Iris Dataset

```{r summary-table}
# Create summary of all methods
summary_df <- data.frame(
  Method = c("Elbow", "Silhouette", "Gap Statistic", 
             "NbClust Majority", "Dunn Index", "clValid Internal"),
  Optimal_k = c(3, 2, 3, 3, which.max(dunn_values) + 1, 2),
  Notes = c("Visual inspection of elbow",
            "Maximum average silhouette",
            "First maximum of gap statistic",
            "Majority vote across 30 indices",
            "Maximum Dunn index",
            "Best connectivity & silhouette")
)

knitr::kable(summary_df, 
             caption = "Summary of Optimal k by Different Methods")
```

## Best Practices for Clustering Validation

1. **Use Multiple Validation Measures**: Different measures may suggest different optimal k. Consider the consensus.

2. **Visual Inspection**: Always visualize your clusters and silhouette plots to understand the clustering structure.

3. **Domain Knowledge**: Statistical measures should be combined with domain expertise. The "optimal" k should make sense in your application context.

4. **Data Scaling**: Always scale/standardize your data before clustering, especially when variables have different units.

5. **Stability Check**: Use stability measures to ensure your clustering is robust to small perturbations.

6. **Interpretation Guidelines**:
   - Average Silhouette > 0.7: Strong clustering
   - Average Silhouette > 0.5: Reasonable clustering
   - Average Silhouette > 0.25: Weak clustering
   - Average Silhouette < 0.25: No substantial structure

7. **Computational Considerations**: For large datasets, use sampling or efficient algorithms (CLARA for PAM, sampling for validation).

## Session Information

```{r session-info}
sessionInfo()
```

# References

1. Rousseeuw, P.J. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. *Journal of Computational and Applied Mathematics*, 20, 53-65.

2. Dunn, J.C. (1974). Well separated clusters and fuzzy partitions. *Journal on Cybernetics*, 4, 95-104.

3. Brock, G., Pihur, V., Datta, S., & Datta, S. (2008). clValid: An R package for cluster validation. *Journal of Statistical Software*, 25(4), 1-22.

4. Charrad, M., Ghazzali, N., Boiteau, V., & Niknafs, A. (2014). NbClust: An R package for determining the relevant number of clusters in a data set. *Journal of Statistical Software*, 61(6), 1-36.

5. Tibshirani, R., Walther, G., & Hastie, T. (2001). Estimating the number of clusters in a data set via the gap statistic. *Journal of the Royal Statistical Society: Series B*, 63(2), 411-423.
