R ile Denetimsiz Öğrenme: K-Means, Hiyerarşik Kümeleme ve PCA

In the previous three posts we worked on problems with a known target: we predicted sales and classified stores as High or Low. Now we move to the side where there are no labels. All we have are indicators for 149 countries, and the question is: are there natural groups among them?

In this final post of the series I apply K-means clustering, hierarchical clustering and principal component analysis (PCA) to the World Happiness Report 2021 data.

“In unsupervised learning there is no right answer; there are only defensible interpretations.”

1. The Data and Scaling

The World Happiness Report 2021 holds a happiness score and six explanatory indicators for 149 countries. We go in with eight numeric variables, leaving out the confidence interval columns (upper and lower whisker and the standard error) because they are derived from the score itself and would count the same information three times.

hap <- read.csv("world-happiness-report-2021.csv", sep = ";", dec = ",")

numeric_vars <- hap[sapply(hap, is.numeric)]
numeric_vars <- numeric_vars[, !names(numeric_vars) %in% c(
  "Standard.error.of.ladder.score", "upperwhisker", "lowerwhisker")]

X <- scale(na.omit(numeric_vars))   # 149 x 8

2. K-Means

K-means divides observations into K clusters and iteratively updates each cluster's centroid. It has a single parameter, and that parameter decides the whole analysis: what should K be?

library(factoextra)
library(cluster)

set.seed(123)
km <- kmeans(X, centers = 2, nstart = 25)

km$size          # 87 62
km$betweenss / km$totss   # 0.3472

fviz_cluster(km, data = X)

3. How Many Clusters? Three Criteria

There is no single right way to choose K. Reading three common criteria together is the soundest approach.

The elbow method

The within-cluster sum of squares is computed for different values of K. It always falls as K grows; what you look for is the point where the fall clearly slows down.

The largest drop is from 1 to 2: 1184 down to 773, a fall of 35%. From two to three the drop is 19%, and from three to four 15%. The elbow sits between 2 and 3.

Silhouette score

The silhouette measures how close each observation is to its own cluster and how far from the others. It ranges from −1 to 1; higher means cleaner separation.

fviz_nbclust(X, kmeans, method = "silhouette")

# k=2 : 0.3095   <- highest
# k=3 : 0.2696
# k=4 : 0.2692
# k=5 : 0.2666
# k=6 : 0.2230

The gap statistic

set.seed(123)
gap <- clusGap(X, FUN = kmeans, nstart = 25, K.max = 10, B = 50)
fviz_gap_stat(gap)

The gap statistic compares the observed within-cluster dispersion against randomly generated reference data. All three criteria point to k = 2.

4. Interpreting the Clusters

The algorithm hands you labels "1" and "2"; giving them meaning is the analyst's job. The method is always the same: look at the cluster centres.

The table speaks for itself. The 87 countries in cluster 1 sit above average on happiness, income, social support, health and freedom, and below average on perceived corruption. The 62 countries in cluster 2 are the mirror image on every one of those indicators.

The interesting one is generosity: it is the only variable running the other way. Countries in cluster 2 score above average on generosity despite lower income and happiness. It is a small but meaningful detail showing that the clustering is not simply a rich–poor divide.

Only 34.7% of the total variation is explained between the clusters (betweenss / totss). The remaining 65.3% is within-cluster variation — both groups are quite heterogeneous internally.

Which countries are in each cluster?

Centroids give a cluster's average profile, but understanding a cluster means looking inside it. The most useful move is to list the observations closest to each centroid: those are the cluster's "typical" members.

for (k in 1:2) {
  members <- X[km$cluster == k, ]
  d <- apply(members, 1, function(row)
              sqrt(sum((row - km$centers[k, ])^2)))
  print(head(hap$Country.name[km$cluster == k][order(d)], 5))
}

The gap between the two mean happiness scores is 1.74 points, with a 95% confidence interval of [1.53, 1.95] (Welch t = 16.45, p < 0.001). Given that the clusters were built using the ladder score, that is no surprise — but the size of the gap shows how pronounced the split is.

The border countries are more interesting. An observation with a negative silhouette resembles its neighbouring cluster more than its own. There are four of them here: Albania, North Macedonia, Turkey and Ukraine.

Are these clusters stable?

Clustering returns a result on any dataset. The real question is whether the same clusters would appear if the data changed slightly. The standard check is bootstrap stability: resample the data repeatedly, rerun the clustering, and measure how well the original clusters survive using Jaccard similarity.

jaccard <- function(a, b)
  length(intersect(a, b)) / length(union(a, b))

set.seed(7)
for (b in 1:200) {
  i <- sample(nrow(X), replace = TRUE)
  new <- kmeans(X[unique(i), ], centers = 2, nstart = 10)
  # best-matching Jaccard for each original cluster
}
# Cluster 1: 0.9912   Cluster 2: 0.9871

The usual rule of thumb: above 0.85 is "highly stable", 0.75–0.85 is "stable", below 0.6 is unreliable. Both clusters sit above 0.98 and clear 0.75 in all 200 bootstrap samples.

And if we had chosen k = 3?

All three criteria pointed at k = 2, but it is worth seeing what that decision costs.

The k = 3 solution lifts the explained share from 34.7% to 47.1% and produces a meaningful middle group. But the variable doing the new separating is perceived corruption: the second cluster is set apart from the third by institutional trust rather than by wealth.

5. Hierarchical Clustering

Unlike K-means, hierarchical clustering does not ask for the number of clusters up front. It treats every observation as its own cluster, merges the closest ones, and produces a tree (a dendrogram); you decide where to cut afterwards.

The critical choice is the linkage method: how is the distance between two clusters defined?

d <- dist(X)

for (m in c("single", "complete", "average", "ward.D2")) {
  hc <- hclust(d, method = m)
  cat(m, cor(d, cophenetic(hc)), "\n")
}

# Highest cophenetic correlation: average
hc <- hclust(d, method = "average")
fviz_dend(hc, k = 2, cex = 0.5)

There is an important distinction here: the highest cophenetic correlation does not mean "the best clustering". It only says the tree represents the distances most faithfully. Ward's method produces more interpretable clusters in many applications; the choice should follow the purpose.

But does this tree say the same thing as k-means?

That sentence reads like a general warning. On this dataset it can be tested directly: cut each tree into two and compare the result with the split k-means found.

for (m in c("average", "ward.D2")) {
  hc <- hclust(d, method = m)
  print(table(kmeans = km$cluster, hierarchical = cutree(hc, 2)))
}

The result is striking. Cut in two, the average method — the one with the highest cophenetic correlation — puts 148 countries in one cluster and a single country in the other. It represents the distances faithfully while producing no usable split. Its adjusted Rand index against k-means is 0.005, no better than chance agreement.

Ward.D2, despite having the lowest cophenetic correlation, yields a balanced 82 + 67 split that overlaps 92.6% with the k-means solution (ARI = 0.725). Only 11 of the 149 countries land in a different cluster.

The adjusted Rand index (ARI) measures agreement between two clusterings: 1 is a perfect match, 0 is the level of similarity you would get by chance. It is invariant to the cluster labels — swapping "Cluster 1" and "Cluster 2" does not change it.

6. Principal Component Analysis

Visualizing eight-variable data in two dimensions needs a method. PCA builds new axes from linear combinations of the variables, carrying as much variation as possible.

pca <- prcomp(numeric_vars, center = TRUE, scale. = TRUE)
summary(pca)

fviz_eig(pca)                      # scree plot
fviz_pca_var(pca, repel = TRUE)   # variable loadings
fviz_pca_biplot(pca, habillage = factor(km$cluster))

The first two components gather 65.5% of the information in eight variables onto two axes. Adding the third takes us to 80%. That makes a two-dimensional plot a reasonable way to show eight-dimensional data — while remembering that 34.5% of the information is not visible in it.

PC1 is a prosperity axis. Happiness, income, health, social support and freedom all load in the same direction, with perceived corruption against them. This single axis explains half the variation and reads in practice as "level of development".

PC2 is a generosity axis. With a loading of 0.702, generosity dominates this component, accompanied by perceived corruption (−0.395) and freedom of choice (+0.348). In other words, there is a second structure that prosperity does not explain, and it concerns social behaviour.

How many components should we keep?

"The first two components explain 65.5%" is an observation, not a decision. Several criteria exist for choosing how many to keep, and on this data two of them disagree.

Kaiser says three components, the broken stick says one. The clash is not accidental: the "eigenvalue > 1" rule is known to be generous and tends to retain more components than it should. The broken-stick criterion compares eigenvalues with the lengths expected from randomly breaking a stick, and is markedly more conservative.

Does clustering after PCA give the same answer?

A common workflow is to reduce dimensions with PCA first and cluster on the components afterwards, on the grounds that it removes noise and makes distances more meaningful. That claim is testable here.

Z <- pca$x[, 1:3]                 # first three components
set.seed(123)
km_pca <- kmeans(scale(Z), centers = 2, nstart = 25)

table(raw = km$cluster, pca = km_pca$cluster)
# agreement: 84.6%

Clustering on the first three components overlaps 84.6% with clustering on the raw scaled data. Twenty-three of the 149 countries land in a different cluster — not a trivial difference.

The reason is straightforward: the first three components carry 80.1% of the total variation and discard the remaining 19.9%. What gets discarded may be noise, but it may equally be information the clustering needed. On this data, part of what was dropped is exactly what separated the second group on institutional trust.

7. Common Mistakes

Conclusion

In unsupervised learning there is no test set to measure accuracy against. In its place you have internal consistency measures — silhouette, cophenetic correlation, explained variance — and, above all, whether the result makes sense.

All three methods pointed to the same structure here: 149 countries divide roughly in two along prosperity indicators, and that split explains about a third of the total variation. PCA then revealed a second dimension that prosperity does not account for, organized around generosity — a finding the clustering alone does not tell you.

This four-part series ends here. From data preparation through regression, classification and clustering, one theme ran through all of it: fitting a model is easy; deciding whether to trust it is hard.

“Clustering produces a result on any data. The real question is whether that result was actually there.”