January 22, 2015

The power of randomness: Pham vs Gap Part II

This is part two of two looking at the performance of methods for determining the number of clusters in k-means clustering. For part one, see Optimal “k” when there’s no cluster? Pham vs Gap Part I.

Code

All of the examples, as well as plots can be found at this Github repo under the directory part2.

Overview

In the previous post, we looked at how the gap statistic and Pham et. al. method behave in several conditions under the reference distribution.

In this post we’ll look at how both methods behave under simulations performed in the gap statistic paper. Unlike Part I where we looked at how those methods behaved when there is no clustering, here we will look at how methods behave when there is clustering.

For references on the gap statistic or the Pham method, please see Part I. Unlike Part I, in this post I’ve also included the “prediction strength” method, also by Tibshirani et. al..

Prediction strength algorithm

We’ll review the prediction strength algorithm briefly here. The paper is pretty good, and the technique is pretty simple and clever.

First, let us define the term co-membership. Co-membership is defined by Tibshirani et. al. as an indicator variable which is 1 if the two observations belong to the same cluster, as the clustering is defined, and 0 otherwise.

The prediction strength idea depends heavily on cross validation and this idea of co-membership. The overview of the algorithm is below:

Randomly partition all of your data into a training set and a testing set
Cluster the training set, then cluster the test set
Compute the co-membership matrix of the test set using the clustering from the training data (this is the “prediction”).
Compute the co-membership matrix of the test set using the testing data (this is the “truth”)
Evaluate how well the training clustering predicts the co-membership of the testing clustering

This procedure is performed using cross-validation several times for each candidate k. Then, for each k, the mean cross-validation prediction score is computed. The largest k such that the mean prediction strength is greater than some value (in the paper 0.8 or 0.9) is chosen as the optimal k.

Simulations

I recreated the simulations from the gap statistic paper. For the most part, they are the same, or similar, but vary the sample size. I’ll quote the paper, then give any differences below.

For each configuration, I did 100 realizations. I did not include simulation (a) as Part I does a more in depth analysis of it.

All simulations were implemented in R. I used the following implementations:

clusGap() from the the cluster package for the gap statistic (see my parallelized fork of the cluster package)
prediction.strength() from the fpc package for the prediction strength method
kselection() from the kselection package for the Pham method

Simulation (b)

three clusters in two dimensions – the clusters are standard normal variable with (25, 25, 30) observations, centered at (0, 0), (0, 5), and (5, -3).

The implementation, (sim_b()), follows this description exactly. There are three additional sample sizes tested. In total, these are the ones simulated: (25, 50, 25), (50, 100, 50), (125, 250, 125), (2500, 5000, 2500).

Simulation (c)

four clusters in three dimensions – each cluster was randomly chosen to have 25 or 50 standard normal observations, with centres randomly chosen as $N(0, 5I)$ (any simulation with clusters having a minimum distance less than 1.0 units between them was discarded;

The implementation (sim_c()) follows the description pretty similarly, except I simulated in four dimensions instead of three due to a typo. The sample sizes tested were: (25, 50), (50, 100), (125, 250), (2500, 5000).

Simulation (d)

four clusters in 10 dimensions – each cluster was randomly chosen to have 25 or 50 standard normal observations, with centres randomly chosen as $N(0, 1.9I)$ (any simulation with clusters having a minimum distance less than 1.0 units between them was discarded; in this and the previous scenario, the settings are such that about half of the random realizations were discarded);

The implementation (sim_d()) is exactly as described. The sample sizes tested were: (25, 50), (50, 100), (125, 250), (2500, 5000).

Simulation (e)

two elongated clusters in three dimensions – each cluster is generated as follows Set $x_1 = x_2 = x_3 = t$ with $t$ taking 100 equally spaced values from -0.5 to 0.5 and then Gaussian noise with standard deviation 0.1 added to each feature. Cluster 2 is generated the same way, except that the value 10 is added to each feature at the end. The result is two elongated clusters, stretching out along the main diagonal of a three dimensional cube.

sim_e() is implemented exactly as described. The sample sizes tested were: 100, 250, 500, 7500. The sample size 15,000 was also tested, but isn’t displayed.

Results

Correct classification of k

Perhaps the most important factor of all of these methods is accuracy. The figure below shows how each method did in the different scenarios:

The header of each subfigure shows the sample size for each simulation.

The table of counts can be seen below in the Github repository by loading part2/results.RData and looking at the variable X_table, where X is the simulation identifier.