3. Clustering discrete distributional data

Distances between groups

The distance between two rows is \(d(i, i') = \sqrt{\sum_j d_W^2(x_{ij}, x_{i'j})}\), where \(d_W\) is the L2 Wasserstein distance between the distributions of item \(j\). The distance matrix is computed in C++:

D <- WD_MAT_DIST(M)
round(D[1:6, 1:6], 3)
#>       F01_M F01_W F02_M F02_W F03_M F03_W
#> F01_M 0.000 1.857 2.848 2.334 1.525 1.398
#> F01_W 1.857 0.000 2.427 2.175 1.901 1.946
#> F02_M 2.848 2.427 0.000 0.949 3.103 2.895
#> F02_W 2.334 2.175 0.949 0.000 2.687 2.491
#> F03_M 1.525 1.901 3.103 2.687 0.000 0.994
#> F03_W 1.398 1.946 2.895 2.491 0.994 0.000

The variables can be weighted (argument weights), e.g. to give the same importance to items with different variability.

Hierarchical clustering

WD_hclust() performs an agglomerative clustering on this distance; Ward’s method ("ward.D2") is consistent with the Wasserstein decomposition of the total variance into a within and a between part.

hc <- WD_hclust(M, method = "ward.D2")
plot(hc, main = "Ward clustering of faculty x gender groups", xlab = "", sub = "")
rect.hclust(hc, k = 3)

cl_h <- cutree(hc, k = 3)
table(cluster = cl_h, gender = groups$Gender)
#>        gender
#> cluster M W
#>       1 5 5
#>       2 6 6
#>       3 2 2
same <- tapply(cl_h, groups$Faculty, function(z) length(unique(z)) == 1)
mean(same)   # share of faculties whose two groups are in the same cluster
#> [1] 0.8461538

Men and women are evenly spread across clusters and, for most faculties, the two groups fall in the same cluster: satisfaction with the orientation services depends on the faculty much more than on the gender of the students.

K-means (dynamic clustering)

WD_kmeans() partitions the rows into k clusters minimizing the sum of squared Wasserstein distances from the cluster prototypes, which are the Wasserstein barycenters of the clusters. Several random starts (rep) are run and the best solution is returned. The quality index is the ratio between the between-cluster and the total sum of squares.

set.seed(123)
qual <- sapply(2:7, function(k) WD_kmeans(M, k, rep = 10)$quality)
ggplot(data.frame(k = 2:7, quality = qual), aes(k, quality)) +
  geom_line() + geom_point() +
  labs(y = "between SS / total SS", title = "Quality of the partition") +
  theme_minimal()

The gain is large up to three clusters and smaller afterwards; we keep \(k = 3\).

set.seed(123)
km <- WD_kmeans(M, k = 3, rep = 20)
km$quality
#> [1] 0.5762092
km$solution$cardinality
#> IDX
#>  1  2  3 
#>  4 12 10
split(names(km$solution$IDX), km$solution$IDX)
#> $`1`
#> [1] "F06_M" "F06_W" "F13_M" "F13_W"
#> 
#> $`2`
#>  [1] "F01_W" "F02_M" "F02_W" "F04_M" "F05_M" "F05_W" "F07_M" "F07_W" "F08_M"
#> [10] "F08_W" "F12_M" "F12_W"
#> 
#> $`3`
#>  [1] "F01_M" "F03_M" "F03_W" "F04_W" "F09_M" "F09_W" "F10_M" "F10_W" "F11_M"
#> [10] "F11_W"

The partition agrees with the hierarchical one:

table(kmeans = km$solution$IDX, hclust = cl_h)
#>       hclust
#> kmeans  1  2  3
#>      1  0  0  4
#>      2  0 12  0
#>      3 10  0  0

Describing the clusters

The prototypes are a MatD (one row per cluster), so they can be plotted and summarized like the data:

Cluster labels are arbitrary; we sort the clusters by the mean global satisfaction of their prototypes (from the least to the most satisfied):

proto <- km$solution$centers
ord <- order(get.MatD.stats(proto, "mean")$mat[, "global"])
proto <- proto[ord, ]
rownames(proto@M) <- c("Low", "Medium", "High")
round(get.MatD.stats(proto, "mean")$mat, 2)
#>        informat willingn officeho compete global
#> Low        4.35     5.17     4.47    4.52   4.64
#> Medium     5.42     5.91     5.12    5.42   5.52
#> High       5.99     6.41     5.67    6.14   6.20
round(get.MatD.stats(proto, "quantile", prob = 0.1)$mat, 2)
#>        informat willingn officeho compete global
#> Low        2.25     2.75     2.25    2.75   2.50
#> Medium     3.50     4.40     3.00    3.70   4.00
#> High       4.50     5.17     3.83    4.83   5.08
plot(proto)

The three clusters correspond to low, medium and high satisfaction: in the least satisfied one the distributions are shifted to the left and are more dispersed, with a non-negligible share of low scores (look at the 10% quantiles), while the most satisfied cluster has distributions concentrated on 6 and 7.

How much of the separation is due to the means? The squared distances between the least and the most satisfied prototypes can be decomposed into position, size and shape components (see WassSqDistD()):

comp <- t(sapply(items, function(v)
  WassSqDistD(get.cell.MatD(proto, 1, v), get.cell.MatD(proto, 3, v), details = TRUE)))
round(comp, 3)
#>           SQ.W POSITION  SIZE SHAPE   rQQ
#> informat 3.081    2.688 0.264 0.129 0.955
#> willingn 2.329    1.541 0.698 0.091 0.963
#> officeho 1.618    1.459 0.060 0.098 0.972
#> compete  3.139    2.617 0.368 0.154 0.942
#> global   2.994    2.431 0.472 0.091 0.961

For all the items most of the distance is due to the different means (POSITION), but differences in dispersion and shape are not negligible, especially for willingn and global.

Standardized clustering

If the variables have very different variability, standardize = TRUE weights each variable by the inverse of its Wasserstein variance:

set.seed(123)
km_s <- WD_kmeans(M, k = 3, rep = 20, standardize = TRUE)
table(raw = km$solution$IDX, standardized = km_s$solution$IDX)
#>    standardized
#> raw  1  2  3
#>   1  0  4  0
#>   2 12  0  0
#>   3  0  0 10

Here the two partitions coincide, up to the numbering of the clusters.

Graphical tools for distributional data (Green Eyes plots, distributional heatmaps with dendrograms) are provided by the companion package DistrViz.

References

• Irpino A., Verde R. (2006). A new Wasserstein based distance for the hierarchical clustering of histogram symbolic data. In Data Science and Classification (IFCS 2006), Springer, 185-192.
• Irpino A., Verde R., Lechevallier Y. (2006). Dynamic clustering of histograms using Wasserstein metric. In COMPSTAT 2006 - Advances in Computational Statistics, Physica-Verlag, 869-876.
• Irpino A., Verde R. (2015). Basic statistics for distributional symbolic variables: a new metric-based approach. Advances in Data Analysis and Classification, 9(2), 143-175. doi:10.1007/s11634-014-0176-4
• Crocetta C., Irpino A., Antonucci L., Marin C. (2025). A performance indicator and its decomposition according to the impacts of different aspects based on distributional data. Quality & Quantity, 59(3), 2091-2110. doi:10.1007/s11135-024-01956-0