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

d  <- dist(Xs, method = "euclidean")  # default distance
hc <- hclust(d, method = "complete")  # default linkage

k  <- 4
clusters <- cutree(hc, k = k)         # cluster membership per observation

# dendrogram fit
coph <- cophenetic(hc)
cor(d, coph)                          # cophenetic correlation

library(factoextra)

fviz_dend(hc, k = k, cex = 0.7,
          k_colors = "jco",           # cluster colour palette
          rect = TRUE, rect_fill = TRUE,
          color_labels_by_k = TRUE,
          main = "Cluster Dendrogram", xlab = "Objects", ylab = "Height")

library(factoextra)
library(cluster)

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

# Silhouette method - hcut() uses its own internal Ward/Euclidean default
sil <- fviz_nbclust(Xs, FUN = hcut, method = "silhouette")
suggested_k_silhouette <- sil$data$clusters[which.max(sil$data$y)]

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

library(factoextra)
library(ggplot2)

# Elbow method
fviz_nbclust(Xs, FUNcluster = kmeans, method = "wss") +
  geom_vline(xintercept = suggested_k_elbow, linetype = "dashed", color = "steelblue", linewidth = 0.8) +
  labs(title = "Elbow Method")

# Silhouette method
fviz_nbclust(Xs, FUN = hcut, method = "silhouette") +
  labs(title = "Average Silhouette Method")

# Gap statistic
fviz_gap_stat(gap) +
  labs(title = "Gap Statistic Method")

library(FactoMineR)

res.pca  <- PCA(Xs, graph = FALSE)          # PCA on the scaled data
res.hcpc <- HCPC(res.pca, nb.clust = k, graph = FALSE)  # Ward + k-means consolidation

res.hcpc$data.clust$clust                   # cluster membership
res.hcpc$desc.var$quanti                    # top variables per cluster - v.test, Mean/sd in category,
                                             # Overall mean/sd, p.value, and n (observations per cluster)

# cluster quality: Total / Within / Between sum-of-squares (PCA space)
total_ss    <- sum(scale(res.pca$ind$coord, scale = FALSE)^2)
between_ss  <- sum(sapply(split(as.data.frame(res.pca$ind$coord), res.hcpc$data.clust$clust),
                          function(g) nrow(g) * sum((colMeans(g))^2)))
within_ss   <- total_ss - between_ss
cluster_seperation_ss <- between_ss / total_ss                       # Between/Total ratio (separation indicator)

# PCA eigenvalues, variable contributions, and variable correlations with components
res.pca$eig                                  # eigenvalue / % variance / cumulative %
res.pca$var$contrib                          # % contribution of each variable per component
res.pca$var$cor                              # correlation of each variable with each component

# Cluster-wise means - raw units (uses the cutree()-based clusters, not HCPC)
cluster_means <- aggregate(X, by = list(cluster = clusters), FUN = mean)

# Cluster-wise means - standardized scale ("Clustering Summary", uses HCPC clusters)
hcpc_cluster <- res.hcpc$data.clust$clust
tab_means <- aggregate(as.data.frame(Xs), by = list(cluster = hcpc_cluster), FUN = mean)

# Intra-/inter-cluster distances (uses the cutree()-based clusters + distance matrix)
dmat        <- as.matrix(dist(Xs, method = "euclidean"))
cluster_ids <- sort(unique(clusters))

intra <- sapply(cluster_ids, function(cl) {
  idx <- which(clusters == cl)
  m   <- dmat[idx, idx]
  mean(m[upper.tri(m)])                      # mean pairwise distance within the cluster
})
names(intra) <- cluster_ids

inter <- combn(cluster_ids, 2, function(pair) {
  mean(dmat[clusters == pair[1], clusters == pair[2]])   # mean pairwise distance between the pair
})
names(inter) <- combn(cluster_ids, 2, paste, collapse = "-")

library(dendextend)

hc1 <- hclust(dist(Xs, method = "euclidean"), method = "average")
hc2 <- hclust(dist(Xs, method = "euclidean"), method = "ward.D2")

dl <- dendlist(as.dendrogram(hc1), as.dendrogram(hc2))
names(dl) <- c("average", "ward.D2")
dl <- untangle(dl, method = "step1side")
entanglement(dl)                             # 0 = identical, 1 = maximal tangling

# Plot the tanglegram (k = number of clusters highlighted on each side)
k_highlight <- 4
tanglegram(dl, k_labels = k_highlight, k_branches = k_highlight,
           main_left = names(dl)[1], main_right = names(dl)[2],
           margin_inner = 10, lwd = 2)

# Linkage method correlation matrix (cophenetic agreement across 6 linkage methods)
methods <- c("single", "complete", "average", "centroid", "median", "ward.D2")
dends   <- lapply(methods, function(m) as.dendrogram(hclust(dist(Xs), method = m)))
cor.dendlist(do.call(dendlist, dends), method = "cophenetic")