Lecture 6 Demo: Gaussian Mixture Models on Penguins

library(palmerpenguins)
library(mclust)

1) Prepare the same penguins data used in Lecture 5

We reuse the same four numeric measurements and keep complete cases.

raw_df <- penguins[, c(
  "species", "island",
  "bill_length_mm", "bill_depth_mm",
  "flipper_length_mm", "body_mass_g"
)]

keep <- complete.cases(raw_df)
penguins_df <- raw_df[keep, ]

X <- as.matrix(penguins_df[, c(
  "bill_length_mm", "bill_depth_mm",
  "flipper_length_mm", "body_mass_g"
)])

X_scaled <- scale(X)

2) Quick baseline: K-means (hard assignments)

set.seed(32950)
km3 <- kmeans(X_scaled, centers = 3, nstart = 20)

table(km3$cluster)

  1   2   3 
132  87 123 

table(TrueSpecies = penguins_df$species, KmeansCluster = km3$cluster)

           KmeansCluster
TrueSpecies   1   2   3
  Adelie    127  24   0
  Chinstrap   5  63   0
  Gentoo      0   0 123

3) Compute GMM using MClust()

Model selection with BIC surface

Though GMM can model scale differences, the solution may not be stable if scales across variables differ a lot. So we still tend to scale the data first before running GMM in practice.

We use R function Mclust(). By default, it automatically select both number of components and covariance structure by BIC.

set.seed(32950)
gmm_auto <- Mclust(X_scaled, G = 1:12, verbose = F) ## G: the possible numbers of components

summary(gmm_auto)

---------------------------------------------------- 
Gaussian finite mixture model fitted by EM algorithm 
---------------------------------------------------- 

Mclust VEE (ellipsoidal, equal shape and orientation) model with 4 components: 

Clustering table:
  1   2   3   4 
154  66  57  65 

The BIC in MClust is defined as the reverse of regularly defined BIC, so higher BIC represents better model in MClust: \[\text{BIC}_{\text{mclust}} = -\text{BIC} =2l(\hat \Theta) - d\log n\]

We can actually visualize how BIC changes across different modeling choices.

plot(gmm_auto, what = "BIC")

You may also specify a specific convariance structure and number of components (for instance, fixed G = 3).

gmm_3 <- Mclust(X_scaled, G = 3, modelNames = "VEE", verbose = F)
gmm_3$bic

[1] -2510.485

The output of Mclust

We have the fitted model stored in gmm_auto$paramters, the estimated responsibilites stored in gmm_auto$z, and the corresponding hard clustering results in gmm_auto$classification and the uncertainties stored in gmm_auto$uncertainty. \[\gamma_{ik} = P(Z_i=k\mid x_i), \quad \text{uncertainty}_i=1-\max_k\gamma_{ik}\]

print("First 6 entries of the responsibilities:")

[1] "First 6 entries of the responsibilities:"

head(gmm_auto$z)

          [,1]         [,2]         [,3]         [,4]
[1,] 0.9999859 7.831331e-42 1.088997e-27 1.405042e-05
[2,] 0.9997910 1.289681e-26 3.344178e-19 2.089817e-04
[3,] 0.9917237 7.758721e-30 1.436243e-21 8.276320e-03
[4,] 0.9999983 2.581199e-43 4.760017e-29 1.671648e-06
[5,] 0.9999928 1.510400e-52 3.022682e-33 7.230942e-06
[6,] 0.9999509 1.237485e-35 1.283305e-24 4.905485e-05

print("First 6 entries of the uncertainty:")

[1] "First 6 entries of the uncertainty:"

head(gmm_auto$uncertainty)

[1] 1.405042e-05 2.089817e-04 8.276320e-03 1.671648e-06 7.230942e-06 4.905485e-05

print("First 6 entries of hard clustering labels:")

[1] "First 6 entries of hard clustering labels:"

head(gmm_auto$classification)

[1] 1 1 1 1 1 1

4) Visualize the clustering results

MClust provides visualzation of clustering results in the original data space

plot(gmm_auto, what = "classification")

plot(gmm_3, what = "classification")

Though BIC prefers \(K = 4\), compared to the results when \(K =3\), the additional cluster is not very well seperated in the data, so \(K=3\) seems a more reasonable choice.

plot(gmm_3, what = "uncertainty")

5) Visualize clustering results in PC space

pca <- prcomp(X_scaled)
score12 <- pca$x[, 1:2]

op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

plot(score12,
     col = penguins_df$species,
     pch = 19,
     xlab = "PC1", ylab = "PC2",
     main = "Ground truth (Species)")
legend("topright",
       legend = levels(as.factor(penguins_df$species)),
       col = 1:length(levels(as.factor(penguins_df$species))),
       pch = 19)

plot(score12,
     col = km3$cluster,
     pch = 19,
     xlab = "PC1", ylab = "PC2",
     main = "K-means clusters")

plot(score12,
     col = gmm_3$classification,
     pch = 19,
     xlab = "PC1", ylab = "PC2",
     main = paste("GMM:", gmm_3$modelName, "(G=", gmm_3$G, ")"))

plot(score12,
     col = "grey80",
     pch = 16,
     xlab = "PC1", ylab = "PC2",
     main = "GMM uncertainty (>0.1 highlighted)")

idx <- gmm_3$uncertainty > 0.1   # adjust threshold
points(score12[idx, ],
       col = "red",
       pch = 19)

par(op)

The GMM results are closer to the truth than K-means.

species <- penguins_df$species

tab_kmeans <- table(Cluster = km3$cluster, Species = species)
tab_gmm  <- table(Cluster = gmm_3$classification, Species = species)

tab_kmeans

       Species
Cluster Adelie Chinstrap Gentoo
      1    127         5      0
      2     24        63      0
      3      0         0    123

tab_gmm

       Species
Cluster Adelie Chinstrap Gentoo
      1    151         5      0
      2      0         0    123
      3      0        63      0

6) How covariance constraints change clustering (fixed G = 3)

To isolate covariance assumptions, hold G = 3 fixed and compare BIC values.

cov_models <- c(
  "EII", # spherical, equal volume
  "VII", # spherical, unequal volume
  "VEE", # ellipsoidal, equal shape and orientation
  "VVV"  # ellipsoidal, varying covariance (most flexible)
)

fits_cov <- lapply(cov_models, function(m) {
  Mclust(X_scaled, G = 3, modelNames = m, verbose = F)
})
names(fits_cov) <- cov_models

bic_cov <- sapply(fits_cov, function(f) f$bic)
comparison_table <- data.frame(
  model = cov_models,
  BIC = round(as.numeric(bic_cov), 2)
)
comparison_table <- comparison_table[order(-comparison_table$BIC), ]
comparison_table

Visualize how the covariance choices change the final clustering results

op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (m in cov_models) {
  plot(score12,
       col = fits_cov[[m]]$classification,
       pch = 19,
       xlab = "PC1", ylab = "PC2",
       main = paste0(m, " (G=3)"))
}

par(op)

9) Takeaways

• K-means gives hard assignments only.
• GMM gives soft assignments via posterior probabilities.
• Different covariance constraints imply different geometric assumptions, making GMM more flexible than K-means.
• BIC helps compare both number of components and covariance structure.