3. Clustering histogram data: k-means and fuzzy c-means

The clustering methods of HistDAWass group the individuals (rows) of a MatH object using the L2 Wasserstein distance between histograms. The prototype of each cluster is itself a histogram for each variable (the Wasserstein mean of the histograms of the cluster).

We use the OzoneFull data: 78 stations in the USA described by the histograms of four variables.

The algorithms start from random partitions, so we fix the seed and repeat each algorithm rep times, keeping the best solution.

K-means

set.seed(1)
km <- WH_kmeans(OzoneFull, k = 3, rep = 5,
                simplify = TRUE, qua = 20, standardize = TRUE)
names(km$solution)
#> [1] "IDX"         "cardinality" "centers"     "Crit"
km$solution$cardinality
#> IDX
#> Cl.1 Cl.2 Cl.3 
#>   31   13   34
km$quality   # proportion of the total sum of squares explained by the partition
#> [1] 0.4845037
head(km$solution$IDX)
#> I1 I2 I3 I4 I5 I6 
#>  1  1  3  3  1  3

The prototypes are a MatH object, so they can be plotted:

plot(km$solution$centers)

Choosing the number of clusters

The quality (between / total sum of squares) increases with k; a common choice is where the gain begins to flatten (elbow):

set.seed(1)
qual <- sapply(2:7, function(k)
  WH_kmeans(OzoneFull, k = k, rep = 5, simplify = TRUE, qua = 20,
            standardize = TRUE)$quality)

round(qual, 3)
#> [1] 0.345 0.485 0.571 0.639 0.682 0.713
plot(2:7, qual, type = "b", xlab = "k", ylab = "quality")

K-means with adaptive distances

WH_adaptive.kmeans() also estimates a weight for each variable, so that the variables that better separate the clusters weigh more. With schema = 2 there is a weight for the position (mean) and one for the variability (centred distribution) of each variable; with schema = 3 or 4 the weights also change with the cluster.

set.seed(1)
akm <- WH_adaptive.kmeans(OzoneFull, k = 3, schema = 2, rep = 5,
                          simplify = TRUE, qua = 20, standardize = TRUE)

round(akm$weights[, 1], 3)   # the same for all the clusters with schema = 2
#>         Ozone.Conc.ppb(Cen.)         Ozone.Conc.ppb(Var.) 
#>                        0.990                        1.100 
#>          Temperature.C(Cen.)          Temperature.C(Var.) 
#>                        0.502                        1.699 
#> Solar.Radiation.WattM2(Cen.) Solar.Radiation.WattM2(Var.) 
#>                        2.172                        0.834 
#>        Wind.Speed.mSec(Cen.)        Wind.Speed.mSec(Var.) 
#>                        0.925                        0.642
akm$cardinality
#> IDX
#> Cl.1 Cl.2 Cl.3 
#>   28   13   37
akm$quality
#> [1] 0.5701985

By default the product of the weights is one (weight.sys = "PROD"); with weight.sys = "SUM" they sum to one.

Fuzzy c-means

WH_fcmeans() assigns to each individual a degree of membership to each cluster; m (greater than 1) controls the fuzziness.

set.seed(1)
fcm <- WH_fcmeans(OzoneFull, k = 3, m = 1.6, rep = 5,
                  simplify = TRUE, qua = 20, standardize = TRUE)

round(head(fcm$solution$membership), 3)
#>     [,1]  [,2]  [,3]
#> I1 0.955 0.002 0.043
#> I2 0.599 0.017 0.383
#> I3 0.185 0.119 0.696
#> I4 0.010 0.004 0.986
#> I5 0.920 0.002 0.077
#> I6 0.074 0.018 0.908
fcm$solution$cardinality   # crisp assignment to the cluster of maximum membership
#> IDX
#> Cl.1 Cl.2 Cl.3 
#>   33   13   32
fcm$quality
#> [1] 0.4470096

Individuals with no clearly dominant membership lie between clusters:

maxmemb <- apply(fcm$solution$membership, 1, max)
sum(maxmemb < 0.6)
#> [1] 12

Fuzzy c-means with adaptive distances

set.seed(1)
afcm <- WH_adaptive_fcmeans(OzoneFull, k = 3, schema = 2, m = 1.6, rep = 3,
                            simplify = TRUE, qua = 20, standardize = TRUE)

round(afcm$solution$var.weights[, 1], 3)
#>         Ozone.Conc.ppb(Cen.)         Ozone.Conc.ppb(Var.) 
#>                        0.944                        1.110 
#>          Temperature.C(Cen.)          Temperature.C(Var.) 
#>                        0.599                        1.776 
#> Solar.Radiation.WattM2(Cen.) Solar.Radiation.WattM2(Var.) 
#>                        1.776                        0.794 
#>        Wind.Speed.mSec(Cen.)        Wind.Speed.mSec(Var.) 
#>                        0.995                        0.638
afcm$quality
#> [1] 0.5286997

Comparing partitions

Hierarchical clustering (WH_hclust(), a wrapper of stats::hclust()) gives a dendrogram that can be cut at any number of clusters:

hc <- WH_hclust(OzoneFull, simplify = TRUE, qua = 20, standardize = TRUE,
                method = "ward.D2")
plot(hc, labels = FALSE, main = "OzoneFull, Ward method")

table(hclust = cutree(hc, k = 3), kmeans = km$solution$IDX)
#>       kmeans
#> hclust  1  2  3
#>      1 12  0 34
#>      2 19  0  0
#>      3  0 13  0

Self-organising maps (WH_2d_Kohonen_maps() and WH_2d_Adaptive_Kohonen_maps()) are also available.

References

Irpino, A., Verde, R., Lechevallier, Y. (2006) Dynamic clustering of histograms using Wasserstein metric. In: COMPSTAT 2006, Physica-Verlag, 869-876.

Irpino, A., Verde, R., De Carvalho, F.A.T. (2014) Dynamic clustering of histogram data based on adaptive squared Wasserstein distances. Expert Systems with Applications, 41, 3351-3366. https://doi.org/10.1016/j.eswa.2013.12.001