data(USArrests)
X  <- USArrests                      # Murder, Assault, UrbanPop, Rape (all numeric)
Xs <- scale(X)                       # z-score standardisation (RAISINS default)

set.seed(123)
km <- kmeans(Xs, centers = 2, iter.max = 50, nstart = 20)

km$centers                           # cluster centroids (scaled units)
km$cluster                           # cluster membership of each observation
km$totss                             # total sum of squares (TSS)
km$tot.withinss                      # total within-cluster sum of squares (TWSS)
km$betweenss                         # between-cluster sum of squares (BSS)
100 * km$betweenss / km$totss        # explained variance (%)
km$withinss                          # within-cluster sum of squares, per cluster
km$size                              # number of observations, per cluster

means_by_cluster <- aggregate(X, by = list(cluster = km$cluster), mean)  # raw (unscaled) cluster means of variables
dd <- cbind(X, cluster = km$cluster)                                     # cluster membership with variables and observations

library(factoextra)
library(cluster)
library(ggplot2)

#Elbow method - WCSS by k
elbow <- fviz_nbclust(Xs, FUNcluster = kmeans, method = "wss")
wss   <- elbow$data$y
suggested_k_elbow <- which.min(diff(wss)) + 1

elbow +
  geom_vline(xintercept = suggested_k_elbow, linetype = "dashed", color = "skyblue", linewidth = 0.8) +
  theme_minimal() +
  labs(title = "Elbow Method")

#Silhouette method - average silhouette width by k
sil <- fviz_nbclust(Xs, FUN = kmeans, method = "silhouette",
                     iter.max = 50, nstart = 20)
suggested_k_silhouette <- sil$data$clusters[which.max(sil$data$y)]

sil +
  ggtitle("Average Silhouette Method") +
  theme(plot.title = element_text(hjust = 0.5))

#Gap statistic
gap <- clusGap(Xs, FUN = kmeans, nstart = 20, K.max = 10, B = 50)
suggested_k_gap <- maxSE(gap$Tab[, "gap"], gap$Tab[, "SE.sim"], method = "firstSEmax")

fviz_gap_stat(gap) +
  ggtitle("Gap Statistic Method") +
  theme(plot.title = element_text(hjust = 0.5))

library(factoextra)
library(clValid)

distance   <- get_dist(Xs, method = "euclidean")   # pairwise Euclidean distance matrix
pca        <- prcomp(Xs, scale. = FALSE)           # principal component analysis
validation <- clValid(Xs, nClust = 2:6, clMethods = "kmeans", validation = "internal")

as.matrix(distance)                  # full pairwise distance matrix
validation@measures[, , 1]           # Connectivity, Dunn Index, Silhouette Width by k

# best k under each internal validation measure
val_df <- as.data.frame(validation@measures[, , 1])
opt_scores <- data.frame(
  Measure  = rownames(val_df),
  Score    = c(min(val_df["Connectivity", ]), max(val_df["Dunn", ]), max(val_df["Silhouette", ])),
  Clusters = c(
    as.numeric(colnames(val_df)[which.min(val_df["Connectivity", ])]),
    as.numeric(colnames(val_df)[which.max(val_df["Dunn", ])]),
    as.numeric(colnames(val_df)[which.max(val_df["Silhouette", ])])
  )
)