Cluster Analysis
Department of Mathematics
The Chinese University of Hong Kong
HSMMC Math Modelling Advanced Training Workshop
June 27, 2026
1 Cluster Analysis (CA)
Cluster analysis is a data analytical technique designed to identify and aggregate homogeneous observations or data points—entities that exhibit comparable traits and inherent properties. Classified as a data mining method, it constitutes an unsupervised learning algorithm that uncovers latent patterns and intrinsic structural relationships embedded within datasets.
Centroid-based clustering and hierarchical clustering constitute two distinct approaches widely adopted for cluster analysis.
2 Python Code
2.1 Step 1: Load Data
# Python equivalent of the provided R code
import pandas as pd
import numpy as np
df = pd.read_csv("df_coffee.csv")Random sample of 20 observations, set random seed for reproducibility
Descriptive statistics
## === Descriptive Statistics ===## Aroma Flavor ... Quakers Category.Two.Defects
## count 20.000000 20.000000 ... 20.0 20.000000
## mean 7.576000 7.605000 ... 0.0 4.050000
## std 0.266959 0.367058 ... 0.0 6.021671
## min 7.170000 6.920000 ... 0.0 0.000000
## 25% 7.420000 7.420000 ... 0.0 0.000000
## 50% 7.540000 7.625000 ... 0.0 2.500000
## 75% 7.750000 7.750000 ... 0.0 4.000000
## max 8.080000 8.670000 ... 0.0 20.000000
##
## [8 rows x 15 columns]Drop columns with zero standard deviation (Uniformity, Clean.Cup, Sweetness, Category.Two.Defects)
drop_cols = ["Uniformity", "Clean.Cup", "Sweetness", "Category.Two.Defects"]
df = df.drop(columns=drop_cols)Remove last column (character categorical variable)
2.2 Step 2: Complete Cases
Drop rows with any missing values (complete cases only)
2.3 Step 3: Scaling
Scaling the dataset is critical because features often have vastly different measurement ranges. Without scaling, variables with larger numerical magnitudes will disproportionately skew clustering outputs.
Two widely used scaling methods for cluster analysis are standardization, which transforms data by subtracting each feature’s mean and dividing by its standard deviation, and normalization, which rescales all values to a fixed bounded range.
\[ X_{stand} = \frac{X - \mu_x}{\sigma_x} \\ X_{norm} = \frac{X - X_{min}}{X_{max}-X_{min}} \]
import pandas as pd
from sklearn.preprocessing import StandardScaler, MinMaxScaler
# ----------------------
# Standardization (Z-score)
# ----------------------
scaler_standard = StandardScaler()
df_scaled_standard = pd.DataFrame(
scaler_standard.fit_transform(df_com),
columns=df_com.columns
)
# ----------------------
# Min-Max Normalization (0~1 range)
# ----------------------
scaler_minmax = MinMaxScaler()
df_scaled_norm = pd.DataFrame(
scaler_minmax.fit_transform(df_com),
columns=df_com.columns
)
# If you only need standardized data:
df_scaled = df_scaled_standard 2.4 Step 4: Optimal Number of Clusters
The optimal number of clusters varies depending on the dataset and its intended application. That said, several methods exist to help identify the ideal cluster count for centroid-based clustering.
2.4.1 Elbow Method
This method involves plotting the within-cluster sum of squares (WSS) against the number of clusters. WSS quantifies how tightly grouped data points are within their respective clusters. The elbow method identifies the optimal cluster count by locating the point on the WSS curve where the rate of decline begins to flatten out; this inflection point is typically chosen as the ideal number of clusters.
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
import numpy as np
# Calculate within-cluster sum of squares (WSS) for elbow method
wss = []
max_k = 10 # You can adjust the maximum number of clusters to test
for k in range(1, max_k + 1):
kmeans = KMeans(n_clusters=k, random_state=100)
kmeans.fit(df_scaled)
wss.append(kmeans.inertia_) # inertia = WSSKMeans(n_clusters=10, random_state=100)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Parameters
Fitted attributes
10 features
| kmeans0 |
| kmeans1 |
| kmeans2 |
| kmeans3 |
| kmeans4 |
| kmeans5 |
| kmeans6 |
| kmeans7 |
| kmeans8 |
| kmeans9 |
# Plot elbow curve
plt.figure(figsize=(8, 5))
plt.plot(range(1, max_k + 1), wss, marker='o')
plt.xlabel("Number of clusters k")
plt.ylabel("Within-cluster sum of squares (WSS)")
plt.title("Elbow method")
plt.axvline(x=2, linestyle='--', color='red') # Vertical dashed line at k=2
plt.show()
2.4.2 Silhouette Score/Width
The silhouette score quantifies how closely each data point aligns with its assigned cluster relative to neighboring clusters. Its values span from −1 to 1: a score of 1 means the point fits perfectly within its cluster, a score of 0 means the point lies directly on the boundary between two clusters, and a score of −1 means the point would be a better fit for an alternative cluster. The optimal number of clusters is typically the value that yields the highest average silhouette score.
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples
import numpy as np
# ----------------------
# Equivalent to fviz_nbclust(..., method = "silhouette")
# Plot average silhouette score for different k values
# ----------------------
max_k = 10
sil_avg_scores = []
k_range = range(2, max_k + 1) # silhouette undefined for k=1
for k in k_range:
km = KMeans(n_clusters=k, random_state=100)
cluster_labels = km.fit_predict(df_scaled)
avg_sil = silhouette_score(df_scaled, cluster_labels)
sil_avg_scores.append(avg_sil)
plt.figure(figsize=(8, 5))
plt.plot(k_range, sil_avg_scores, marker="o")
plt.xlabel("Number of clusters k")
plt.ylabel("Average Silhouette Score")
plt.title("Silhouette method")
plt.show()
# ----------------------
# Equivalent to silhouette() + plot + fviz_silhouette for k=2
# ----------------------
# Fit KMeans with 2 clusters
km_res = KMeans(n_clusters=2, random_state=100)
labels = km_res.fit_predict(df_scaled)
# Compute silhouette values for every single sample
sil_vals = silhouette_samples(df_scaled, labels)
# Simple silhouette plot (matches base R plot(sil))
fig, ax = plt.subplots(figsize=(8, 5))
y_lower = 10
for cluster_id in np.unique(labels):
# Extract silhouette scores for current cluster and sort
cluster_sil = sil_vals[labels == cluster_id]
cluster_sil.sort()
y_upper = y_lower + len(cluster_sil)
ax.fill_betweenx(np.arange(y_lower, y_upper), 0, cluster_sil, alpha=0.7)
y_lower = y_upper + 10
ax.axvline(x=np.mean(sil_vals), color="red", linestyle="--")
ax.set_xlabel("Silhouette Coefficient")
ax.set_ylabel("Sample Index")
ax.set_title("Silhouette Plot (k=2)")
plt.show()
A positive silhouette coefficient means an observation fits well within its assigned cluster.
2.4.3 Gap Statistic Method
The gap statistic helps determine the optimal number of clusters by comparing the observed within-cluster dispersion with a reference distribution. The idea is to select the number of clusters that maximizes the gap statistic. A larger gap indicates a more significant deviation from randomness, suggesting a better-defined clustering structure. By comparing the gaps for different numbers of clusters, we can identify the number of clusters that provides the most distinct and meaningful clustering pattern.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.metrics import pairwise_distances_argmin_min
def calculate_gap_statistic(data, max_k, n_references=10):
"""
Replicate factoextra::fviz_nbclust(..., method = "gap_stat")
Compute gap statistic and standard error for each k
"""
# Convert dataframe to numpy array
X = data.values if hasattr(data, "values") else np.array(data)
n_samples, n_features = X.shape
k_range = np.arange(1, max_k + 1)
gap_scores = []
gap_standard_errors = []
for k in k_range:
# 1. Calculate within-cluster dispersion Wk on real data
kmeans = KMeans(n_clusters=k, random_state=100)
labels = kmeans.fit_predict(X)
dists, _ = pairwise_distances_argmin_min(X, kmeans.cluster_centers_)
wk = np.sum(dists ** 2)
# 2. Generate random reference datasets & compute reference Wkb
ref_wk_log = []
for _ in range(n_references):
ref_data = np.random.uniform(
low=X.min(axis=0),
high=X.max(axis=0),
size=X.shape
)
ref_kmeans = KMeans(n_clusters=k, random_state=100)
ref_labels = ref_kmeans.fit_predict(ref_data)
ref_dists, _ = pairwise_distances_argmin_min(ref_data, ref_kmeans.cluster_centers_)
ref_wk = np.sum(ref_dists ** 2)
ref_wk_log.append(np.log(ref_wk))
# Gap statistic formula
log_wk = np.log(wk)
mean_log_ref = np.mean(ref_wk_log)
gap = mean_log_ref - log_wk
gap_scores.append(gap)
# Standard error for error bars
sd_ref = np.std(ref_wk_log)
se = sd_ref * np.sqrt(1 + 1 / n_references)
gap_standard_errors.append(se)
return k_range, np.array(gap_scores), np.array(gap_standard_errors)
# ---------------------- Run Gap Statistic Plot ----------------------
# Replace df_scaled with your actual standardized dataframe variable
max_k_test = 10
k_values, gap_values, gap_error = calculate_gap_statistic(df_scaled, max_k=max_k_test)
# Plot equivalent to R fviz_nbclust gap_stat chart
plt.figure(figsize=(8, 5))
plt.errorbar(k_values, gap_values, yerr=gap_error, marker="o", capsize=4)
plt.xlabel("Number of clusters k")
plt.ylabel("Gap statistic")
plt.title("Gap statistic method")
plt.xticks(k_values)## ([<matplotlib.axis.XTick object at 0x00000208975922D0>, <matplotlib.axis.XTick object at 0x00000208954B81A0>, <matplotlib.axis.XTick object at 0x00000208954B65D0>, <matplotlib.axis.XTick object at 0x00000208975DACF0>, <matplotlib.axis.XTick object at 0x00000208975D9F10>, <matplotlib.axis.XTick object at 0x00000208975D93A0>, <matplotlib.axis.XTick object at 0x00000208975D8560>, <matplotlib.axis.XTick object at 0x0000020897591A30>, <matplotlib.axis.XTick object at 0x00000208975DB470>, <matplotlib.axis.XTick object at 0x00000208975DBE60>], [Text(1, 0, '1'), Text(2, 0, '2'), Text(3, 0, '3'), Text(4, 0, '4'), Text(5, 0, '5'), Text(6, 0, '6'), Text(7, 0, '7'), Text(8, 0, '8'), Text(9, 0, '9'), Text(10, 0, '10')])
2.4.4 Combination of All
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, pairwise_distances_argmin_min
# --------------------------
# 1. Helper Function: Elbow Method (WSS)
# --------------------------
def elbow_curve(data, max_k=10, seed=100):
X = data.values
wss = []
for k in range(1, max_k+1):
km = KMeans(n_clusters=k, random_state=seed)
km.fit(X)
wss.append(km.inertia_)
return np.arange(1, max_k+1), np.array(wss)
# --------------------------
# 2. Helper Function: Silhouette Score
# --------------------------
def silhouette_curve(data, max_k=10, seed=100):
X = data.values
k_range = np.arange(2, max_k+1)
sil_scores = []
for k in k_range:
km = KMeans(n_clusters=k, random_state=seed)
labels = km.fit_predict(X)
sil_scores.append(silhouette_score(X, labels))
return k_range, np.array(sil_scores)
# --------------------------
# 3. Helper Function: Gap Statistic
# --------------------------
def gap_stat_curve(data, max_k=10, n_ref=10, seed=100):
X = data.values
n, p = X.shape
k_list = np.arange(1, max_k+1)
gaps = []
np.random.seed(seed)
for k in k_list:
# Real data Wk
km = KMeans(n_clusters=k, random_state=seed)
labels = km.fit_predict(X)
dists, _ = pairwise_distances_argmin_min(X, km.cluster_centers_)
wk = np.sum(dists ** 2)
log_wk = np.log(wk)
# Reference uniform datasets
ref_log = []
for _ in range(n_ref):
ref = np.random.uniform(X.min(axis=0), X.max(axis=0), X.shape)
km_r = KMeans(n_clusters=k, random_state=seed)
r_labels = km_r.fit_predict(ref)
r_dists, _ = pairwise_distances_argmin_min(ref, km_r.cluster_centers_)
ref_log.append(np.log(np.sum(r_dists ** 2)))
gap = np.mean(ref_log) - log_wk
gaps.append(gap)
return k_list, np.array(gaps)
# --------------------------
# 4. Run all three methods on scaled data
# --------------------------
max_cluster = 10
k_elbow, wss_vals = elbow_curve(df_scaled, max_k=max_cluster)
k_sil, sil_vals = silhouette_curve(df_scaled, max_k=max_cluster)
k_gap, gap_vals = gap_stat_curve(df_scaled, max_k=max_cluster)
# Find optimal k for each method
# Elbow: pick elbow inflection point (simple heuristic: max drop rate)
wss_drop = np.diff(wss_vals)
opt_elbow = np.argmin(wss_drop) + 2
# Silhouette: k with maximum silhouette score
opt_sil = k_sil[np.argmax(sil_vals)]
# Gap Statistic: k with maximum gap value
opt_gap = k_gap[np.argmax(gap_vals)]
# Aggregate votes (match R n_clust majority logic)
votes = [opt_elbow, opt_sil, opt_gap]
from collections import Counter
vote_count = Counter(votes)
majority_k = vote_count.most_common(1)[0][0]
# Print summary (equivalent to print(n_clust) in R)
print("=== Optimal Cluster Count Votes from Multiple Methods ===")## === Optimal Cluster Count Votes from Multiple Methods ===## Elbow Method optimal k: 2## Silhouette Method optimal k: 2## Gap Statistic optimal k: 1##
## Vote tally: {np.int64(2): 2, np.int64(1): 1}## Majority-supported optimal number of clusters: 2# --------------------------
# Plot combined results (equivalent to plot(n_clust) in R)
# --------------------------
fig, axes = plt.subplots(3, 1, figsize=(8, 10))
# Subplot 1: Elbow Method
axes[0].plot(k_elbow, wss_vals, marker="o")
axes[0].axvline(opt_elbow, color="red", ls="--", label=f"Opt k={opt_elbow}")
axes[0].set_title("Elbow Method (WSS)")
axes[0].set_xlabel("Number of clusters k")
axes[0].set_ylabel("Within-cluster Sum of Squares")
axes[0].legend()
axes[0].grid(alpha=0.3)
# Subplot 2: Silhouette Score
axes[1].plot(k_sil, sil_vals, marker="o", color="green")
axes[1].axvline(opt_sil, color="red", ls="--", label=f"Opt k={opt_sil}")
axes[1].set_title("Silhouette Score Method")
axes[1].set_xlabel("Number of clusters k")
axes[1].set_ylabel("Average Silhouette Score")
axes[1].legend()
axes[1].grid(alpha=0.3)
# Subplot3: Gap Statistic
axes[2].plot(k_gap, gap_vals, marker="o", color="orange")
axes[2].axvline(opt_gap, color="red", ls="--", label=f"Opt k={opt_gap}")
axes[2].set_title("Gap Statistic Method")
axes[2].set_xlabel("Number of clusters k")
axes[2].set_ylabel("Gap Statistic Value")
axes[2].legend()
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.show()
2.5 Step 5: Hierarchical Clustering
2.5.1 Step 5.1: Distance
2.5.1.1 Numeric Data
Use the function dist():
import pandas as pd
import numpy as np
from scipy.spatial.distance import pdist
# Convert scaled dataframe to numpy array
X = df_scaled.values
# ----------------------
# Numeric Data Distance Matrices (equivalent to R dist())
# pdist returns condensed distance matrix same as R dist object
# ----------------------
# Euclidean distance (R method = 'euclidian' / euclidean)
df_dist_euc = pdist(X, metric="euclidean")
# Maximum / Chebyshev distance
# df_dist_max = pdist(X, metric="chebyshev")
# Manhattan / L1 distance
# df_dist_man = pdist(X, metric="cityblock")
# Canberra distance
# df_dist_can = pdist(X, metric="canberra")
# Minkowski distance (p=2 = Euclidean by default)
# df_dist_mink = pdist(X, metric="minkowski", p=3)
# ----------------------
# Categorical binary distance (R method = "binary")
# For dummified binary categorical data matrix dummified_data
# ----------------------
# dummified_data_np = dummified_data.values
# dist_binary = pdist(dummified_data_np, metric="jaccard")2.5.2 Step 5.2: Linkage Criteria
Linkage defines the rule used to merge clusters during hierarchical clustering. Multiple linkage criteria exist, and the selected criterion can drastically alter clustering outputs:
- Complete linkage: Uses the maximum distance between two clusters
- Single linkage: Uses the minimum distance between two clusters
- Average linkage: Uses the average pairwise distance between two clusters
- Ward’s linkage: Minimizes the total within-cluster variance
Below is code to generate hierarchical clusters from numeric data using Euclidean distance and each of the above linkage methods.
from scipy.cluster.hierarchy import linkage
# df_dist_euc is the condensed Euclidean distance matrix from pdist()
# linkage() takes condensed distance matrix as input, equivalent to R hclust()
hclust_single = linkage(df_dist_euc, method="single")
hclust_complete = linkage(df_dist_euc, method="complete")
hclust_average = linkage(df_dist_euc, method="average")
hclust_centroid = linkage(df_dist_euc, method="centroid")
hclust_ward = linkage(df_dist_euc, method="ward")Get the cluster information:
from scipy.cluster.hierarchy import fcluster
# Cut tree into exactly 2 clusters (equivalent to cutree(..., k = 2))
clusters_k2 = fcluster(hclust_ward, t=2, criterion="maxclust")
# Cut tree at height h = 0.7 (equivalent to cutree(..., h = 0.7))
clusters_h07 = fcluster(hclust_ward, t=0.7, criterion="distance")import pandas as pd
from scipy.cluster.hierarchy import fcluster
# Assign cluster labels by cutting dendrogram into 2 clusters
cluster_labels = fcluster(hclust_ward, t=2, criterion="maxclust")
# Attach cluster column and compute mean summary per cluster
df_num["cluster"] = cluster_labels
cluster_summary = df_num.groupby("cluster").mean().add_suffix("_Average")
print(cluster_summary)## Aroma_Average ... Category.One.Defects_Average
## cluster ...
## 1 7.425455 ... 0.636364
## 2 7.760000 ... 0.222222
##
## [2 rows x 10 columns]Dendrogram will help to visualize as well as to determine the number of clusters to use or the height where to cut.
import numpy as np
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram, fcluster
# Total sample count
n_samples = hclust_ward.shape[0] + 1
# ====================== 1. Basic dendrogram + safe blue cutoff line ======================
plt.figure(figsize=(10, 6))
dendrogram(hclust_ward)## {'icoord': [[35.0, 35.0, 45.0, 45.0], [25.0, 25.0, 40.0, 40.0], [65.0, 65.0, 75.0, 75.0], [95.0, 95.0, 105.0, 105.0], [85.0, 85.0, 100.0, 100.0], [70.0, 70.0, 92.5, 92.5], [55.0, 55.0, 81.25, 81.25], [32.5, 32.5, 68.125, 68.125], [15.0, 15.0, 50.3125, 50.3125], [5.0, 5.0, 32.65625, 32.65625], [115.0, 115.0, 125.0, 125.0], [135.0, 135.0, 145.0, 145.0], [185.0, 185.0, 195.0, 195.0], [175.0, 175.0, 190.0, 190.0], [165.0, 165.0, 182.5, 182.5], [155.0, 155.0, 173.75, 173.75], [140.0, 140.0, 164.375, 164.375], [120.0, 120.0, 152.1875, 152.1875], [18.828125, 18.828125, 136.09375, 136.09375]], 'dcoord': [[0.0, np.float64(2.8211994490582724), np.float64(2.8211994490582724), 0.0], [0.0, np.float64(3.3416529959187184), np.float64(3.3416529959187184), np.float64(2.8211994490582724)], [0.0, np.float64(1.3541797122858186), np.float64(1.3541797122858186), 0.0], [0.0, np.float64(1.2640743886952885), np.float64(1.2640743886952885), 0.0], [0.0, np.float64(2.428994380285572), np.float64(2.428994380285572), np.float64(1.2640743886952885)], [np.float64(1.3541797122858186), np.float64(2.7691404691045887), np.float64(2.7691404691045887), np.float64(2.428994380285572)], [0.0, np.float64(3.559154299785624), np.float64(3.559154299785624), np.float64(2.7691404691045887)], [np.float64(3.3416529959187184), np.float64(4.403923017090526), np.float64(4.403923017090526), np.float64(3.559154299785624)], [0.0, np.float64(6.027631053839452), np.float64(6.027631053839452), np.float64(4.403923017090526)], [0.0, np.float64(6.590497475713074), np.float64(6.590497475713074), np.float64(6.027631053839452)], [0.0, np.float64(2.2713867810103223), np.float64(2.2713867810103223), 0.0], [0.0, np.float64(1.41825729301792), np.float64(1.41825729301792), 0.0], [0.0, np.float64(1.3122498495396393), np.float64(1.3122498495396393), 0.0], [0.0, np.float64(1.6038764997622879), np.float64(1.6038764997622879), np.float64(1.3122498495396393)], [0.0, np.float64(2.2724840629843372), np.float64(2.2724840629843372), np.float64(1.6038764997622879)], [0.0, np.float64(3.137028030598917), np.float64(3.137028030598917), np.float64(2.2724840629843372)], [np.float64(1.41825729301792), np.float64(4.843050954713168), np.float64(4.843050954713168), np.float64(3.137028030598917)], [np.float64(2.2713867810103223), np.float64(7.30855846442715), np.float64(7.30855846442715), np.float64(4.843050954713168)], [np.float64(6.590497475713074), np.float64(12.194460625748167), np.float64(12.194460625748167), np.float64(7.30855846442715)]], 'ivl': ['8', '16', '19', '2', '3', '6', '5', '15', '1', '11', '17', '4', '13', '0', '9', '14', '12', '7', '10', '18'], 'leaves': [8, 16, 19, 2, 3, 6, 5, 15, 1, 11, 17, 4, 13, 0, 9, 14, 12, 7, 10, 18], 'color_list': ['C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C0'], 'leaves_color_list': ['C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2']}# Safe calculation for cut height k=2
clusters_k2 = fcluster(hclust_ward, t=2, criterion="maxclust")
dist_vals = np.unique(hclust_ward[:, 2])
dist_vals_sorted = np.sort(dist_vals)
if len(dist_vals_sorted) >= 2:
cut_height = dist_vals_sorted[-2]
else:
cut_height = dist_vals_sorted[-1]
plt.axhline(y=cut_height, c="blue", lw=2, linestyle="--")
plt.title("Ward Dendrogram with Cut Line k=2")
plt.xlabel("Samples")
plt.ylabel("Distance")
plt.xticks(rotation=45, ha="right")## (array([ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125,
## 135, 145, 155, 165, 175, 185, 195]), [Text(5, 0, '8'), Text(15, 0, '16'), Text(25, 0, '19'), Text(35, 0, '2'), Text(45, 0, '3'), Text(55, 0, '6'), Text(65, 0, '5'), Text(75, 0, '15'), Text(85, 0, '1'), Text(95, 0, '11'), Text(105, 0, '17'), Text(115, 0, '4'), Text(125, 0, '13'), Text(135, 0, '0'), Text(145, 0, '9'), Text(155, 0, '14'), Text(165, 0, '12'), Text(175, 0, '7'), Text(185, 0, '10'), Text(195, 0, '18')])
# ====================== 2. Simple dendrogram (REMOVED broken branch color function) ======================
plt.figure(figsize=(10, 6))
dendrogram(hclust_ward)## {'icoord': [[35.0, 35.0, 45.0, 45.0], [25.0, 25.0, 40.0, 40.0], [65.0, 65.0, 75.0, 75.0], [95.0, 95.0, 105.0, 105.0], [85.0, 85.0, 100.0, 100.0], [70.0, 70.0, 92.5, 92.5], [55.0, 55.0, 81.25, 81.25], [32.5, 32.5, 68.125, 68.125], [15.0, 15.0, 50.3125, 50.3125], [5.0, 5.0, 32.65625, 32.65625], [115.0, 115.0, 125.0, 125.0], [135.0, 135.0, 145.0, 145.0], [185.0, 185.0, 195.0, 195.0], [175.0, 175.0, 190.0, 190.0], [165.0, 165.0, 182.5, 182.5], [155.0, 155.0, 173.75, 173.75], [140.0, 140.0, 164.375, 164.375], [120.0, 120.0, 152.1875, 152.1875], [18.828125, 18.828125, 136.09375, 136.09375]], 'dcoord': [[0.0, np.float64(2.8211994490582724), np.float64(2.8211994490582724), 0.0], [0.0, np.float64(3.3416529959187184), np.float64(3.3416529959187184), np.float64(2.8211994490582724)], [0.0, np.float64(1.3541797122858186), np.float64(1.3541797122858186), 0.0], [0.0, np.float64(1.2640743886952885), np.float64(1.2640743886952885), 0.0], [0.0, np.float64(2.428994380285572), np.float64(2.428994380285572), np.float64(1.2640743886952885)], [np.float64(1.3541797122858186), np.float64(2.7691404691045887), np.float64(2.7691404691045887), np.float64(2.428994380285572)], [0.0, np.float64(3.559154299785624), np.float64(3.559154299785624), np.float64(2.7691404691045887)], [np.float64(3.3416529959187184), np.float64(4.403923017090526), np.float64(4.403923017090526), np.float64(3.559154299785624)], [0.0, np.float64(6.027631053839452), np.float64(6.027631053839452), np.float64(4.403923017090526)], [0.0, np.float64(6.590497475713074), np.float64(6.590497475713074), np.float64(6.027631053839452)], [0.0, np.float64(2.2713867810103223), np.float64(2.2713867810103223), 0.0], [0.0, np.float64(1.41825729301792), np.float64(1.41825729301792), 0.0], [0.0, np.float64(1.3122498495396393), np.float64(1.3122498495396393), 0.0], [0.0, np.float64(1.6038764997622879), np.float64(1.6038764997622879), np.float64(1.3122498495396393)], [0.0, np.float64(2.2724840629843372), np.float64(2.2724840629843372), np.float64(1.6038764997622879)], [0.0, np.float64(3.137028030598917), np.float64(3.137028030598917), np.float64(2.2724840629843372)], [np.float64(1.41825729301792), np.float64(4.843050954713168), np.float64(4.843050954713168), np.float64(3.137028030598917)], [np.float64(2.2713867810103223), np.float64(7.30855846442715), np.float64(7.30855846442715), np.float64(4.843050954713168)], [np.float64(6.590497475713074), np.float64(12.194460625748167), np.float64(12.194460625748167), np.float64(7.30855846442715)]], 'ivl': ['8', '16', '19', '2', '3', '6', '5', '15', '1', '11', '17', '4', '13', '0', '9', '14', '12', '7', '10', '18'], 'leaves': [8, 16, 19, 2, 3, 6, 5, 15, 1, 11, 17, 4, 13, 0, 9, 14, 12, 7, 10, 18], 'color_list': ['C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C0'], 'leaves_color_list': ['C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2']}plt.title("Ward Dendrogram")
plt.xlabel("Samples")
plt.ylabel("Distance")
plt.xticks(rotation=45, ha="right")## (array([ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125,
## 135, 145, 155, 165, 175, 185, 195]), [Text(5, 0, '8'), Text(15, 0, '16'), Text(25, 0, '19'), Text(35, 0, '2'), Text(45, 0, '3'), Text(55, 0, '6'), Text(65, 0, '5'), Text(75, 0, '15'), Text(85, 0, '1'), Text(95, 0, '11'), Text(105, 0, '17'), Text(115, 0, '4'), Text(125, 0, '13'), Text(135, 0, '0'), Text(145, 0, '9'), Text(155, 0, '14'), Text(165, 0, '12'), Text(175, 0, '7'), Text(185, 0, '10'), Text(195, 0, '18')])
2.6 Step 5: K-means Clustering
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
# 1. Run K-Means clustering (equivalent to kmeans(df_scaled, centers = 2))
km_model = KMeans(n_clusters=2, random_state=100)
cluster_assign = km_model.fit_predict(df_scaled)
# Print model information, same as printing km object in R
print("K-Means Model Output")## K-Means Model Output## Cluster centers:
## [[ 0.55342061 0.65939883 0.68248618 0.61146386 0.62287806 0.42464149
## 0.64082913 0.67537713 -0.15838297 -0.23102817]
## [-0.67640297 -0.80593191 -0.83414977 -0.74734472 -0.76129541 -0.51900626
## -0.7832356 -0.82546094 0.19357919 0.28236777]]## Inertia (total within-cluster sum of squares): 124.3613## Cluster labels for each sample:
## [0 1 1 0 0 0 1 0 1 0 0 1 0 0 0 1 1 1 0 1]# 2. Equivalent to fviz_cluster: PCA 2D projection for cluster visualization
# fviz_cluster automatically uses PCA to reduce dimensions to 2
pca = PCA(n_components=2)
pca_result = pca.fit_transform(df_scaled)
pca_df = pd.DataFrame(pca_result, columns=["PC1", "PC2"])
pca_df["cluster"] = cluster_assign
plt.figure(figsize=(8, 5))
sns.scatterplot(data=pca_df, x="PC1", y="PC2", hue="cluster", palette="Set1", s=60)
# Plot cluster centers projected to PCA space
centers_pca = pca.transform(km_model.cluster_centers_)
plt.scatter(centers_pca[:, 0], centers_pca[:, 1], marker="X", c="black", s=200, label="Centers")
plt.legend(title="Clusters")
plt.title("K-Means Clusters (PCA Projection)")
sns.set_style("white") # match theme_minimal()
plt.tight_layout()
plt.show()
# 3. Scatter plot of raw variables Aroma vs Aftertaste colored by cluster
df_km = df_num.copy()
df_km["cluster"] = cluster_assign
plt.figure(figsize=(8, 5))
sns.scatterplot(data=df_km, x="Aroma", y="Aftertaste", hue="cluster", palette="Set1", s=60)
plt.legend(title="Clusters")
plt.title("Clusters: Aroma vs Aftertaste")
sns.set_style("white")
plt.tight_layout()
plt.show()
# 4. Cluster mean summary table (dplyr group_by + summarise all mean values)
cluster_summary = df_km.groupby("cluster").mean().add_suffix("_Average")
print("\nCluster Variable Average Summary:")##
## Cluster Variable Average Summary:## Aroma_Average ... Category.One.Defects_Average
## cluster ...
## 0 7.72 ... 0.181818
## 1 7.40 ... 0.777778
##
## [2 rows x 10 columns]