Clustering of very skewed, count data: any suggestions to go about (transform etc)?

Question

Basic problem

Here is my basic problem: I am trying to cluster a dataset containing some very skewed variables with counts. The variables contain many zeros and are therefore not very informative for my clustering procedure - which is likely to be k-means algorithm.

Fine, you say, just transform the variables using square root, box cox, or logarithm. But since my variables are based on categorical variables, I fear that I might introduce a bias by handling a variable (based on one value of the categorical variable), while leaving others (based on others values of the categorical variable) the way they are.

Let's go into some more detail.

The dataset

My dataset represents purchases of items. The items have different categories, for example color: blue, red, and green. The purchases are then grouped together, e.g. by customers. Each of these customers is represented by one row of my dataset, so I somehow have to aggregate purchases over customers.

The way I do this is by counting the number of purchases, where the item is a certain color. So instead of a single variable color, I end up with three variables count_red, count_blue, and count_green.

Here is an example for illustration:

-----------------------------------------------------------
customer | count_red  |    count_blue   | count_green     |
-----------------------------------------------------------
c0       |    12      |        5        |       0         |
-----------------------------------------------------------
c1       |     3      |        4        |       0         |
-----------------------------------------------------------
c2       |     2      |       21        |       0         |
-----------------------------------------------------------
c3       |     4      |        8        |       1         |
-----------------------------------------------------------

Actually, I do not use absolute counts in the end, I use ratios (fraction of green items of all purchased items per customer).

-----------------------------------------------------------
customer | count_red  |    count_blue   | count_green     |
-----------------------------------------------------------
c0       |    0.71    |        0.29     |       0.00      |
-----------------------------------------------------------
c1       |    0.43    |        0.57     |       0.00      |
-----------------------------------------------------------
c2       |    0.09    |        0.91     |       0.00      |
-----------------------------------------------------------
c3       |    0.31    |        0.62     |       0.08      |
-----------------------------------------------------------

The result is the same: For one of my colors, e.g. green (nobody likes green), I get a left-skewed variable containing many zeros. Consequently, k-means fails to find a good partitioning for this variable.

On the other hand, if I standardize my variables (subtract mean, divide by standard deviation), the green variable "blows up" due to its small variance and takes values from a much larger range than the other variables, which makes it look more important to k-means than it actually is.

The next idea is to transform the sk(r)ewed green variable.

Transforming the skewed variable

If I transform the green variable by applying the square root it looks a little less skewed. (Here the green variable is plotted in red and green to ensure confusion.)

Red: original variable; blue: transformed by square root.

Let's say I am satisfied with the result of this transformation (which I am not, since the zeros still strongly skew the distribution). Should I now also scale the red and blue variables, although their distributions look fine?

Bottom line

In other words, do I distort the clustering results by handling the color green on one way, but not handling red and blue at all? In the end, all three variables belong together, so shouldn't they be handled in the same way?

EDIT

To clarify: I am aware that k-means is probably not the way to go for count-based data. My question however really is about the treatment of dependent variables. Choosing the correct method is a separate matter.

The inherent constraint in my variables is that

count_red(i) + count_blue(i) + count_green(i) = n(i), where n(i) is the total number of purchases of customer i.

(Or, equivalently, count_red(i) + count_blue(i) + count_green(i) = 1 when using relative counts.)

If I transform my variables differently, this corresponds to giving different weights to the three terms in the constraint. If my goal is to optimally separate groups of customers, do I have to care about violating this constraint? Or does "the end justify the means"?

Answer 1

It is not wise to transform the variables individually because they belong together (as you noticed) and to do k-means because the data are counts (you might, but k-means is better to do on continuous attributes such as length for example).

In your place, I would compute chi-square distance (perfect for counts) between every pair of customers, based on the variables containing counts. Then do hierarchical clustering (for example, average linkage method or complete linkage method - they do not compute centroids and threfore don't require euclidean distance) or some other clustering working with arbitrary distance matrices.

Copying example data from the question:

-----------------------------------------------------------
customer | count_red  |    count_blue   | count_green     |
-----------------------------------------------------------
c0       |    12      |        5        |       0         |
-----------------------------------------------------------
c1       |     3      |        4        |       0         |
-----------------------------------------------------------
c2       |     2      |       21        |       0         |
-----------------------------------------------------------
c3       |     4      |        8        |       1         |
-----------------------------------------------------------

Consider pair c0 and c1 and compute Chi-square statistic for their 2x3 frequency table. Take the square root of it (like you take it when you compute usual euclidean distance). That is your distance. If the distance is close to 0 the two customers are similar.

It may bother you that sums in rows in your table differ and so it affects the chi-square distance when you compare c0 with c1 vs c0 with c2. Then compute the (root of) the Phi-square distance: Phi-sq = Chi-sq/N where N is the combined total count in the two rows (customers) currently considered. It is thus normalized distance wrt to overall counts.

Here is the matrix of sqrt(Chi-sq) distance between your four customers
 .000   1.275   4.057   2.292
1.275    .000   2.124    .862
4.057   2.124    .000   2.261
2.292    .862   2.261    .000

And here is the matrix of sqrt(Phi-sq) distance 
.000    .260    .641    .418
.260    .000    .388    .193
.641    .388    .000    .377
.418    .193    .377    .000

So, the distance between any two rows of the data is the (sq. root of) the chi-square or phi-square statistic of the 2 x p frequency table (p is the number of columns in the data). If any column(s) in the current 2 x p table is complete zero, cut off that column and compute the distance based on the remaining nonzero columns (it is OK and this is how, for example, SPSS does when it computes the distance). Chi-square distance is actually a weighted euclidean distance.

Answer 2

@ttnphns has provided a good answer.

Doing clustering well is often about thinking very hard about your data, so let's do some of that. To my mind, the most fundamental aspect of your data is that they are compositional.

On the other hand, your primary concern seems to be that you have a lot of 0s for green products and specifically wonder if you can transform only the green values to make it more similar to the rest. But because these are compositional data, you cannot think about one set of counts independently from the rest. Moreover, it appears that what you are really interested in are customers' probabilities of purchasing different colored products, but because many have not purchased any green ones, you worry that you cannot estimate those probabilities. One way to address this is to use a somewhat Bayesian approach in which we nudge customers' estimated proportions towards a mean proportion, with the amount of the shift influenced by how far they are from the mean and how much data you have to estimate their true probabilities.

Below I use your example dataset to illustrate (in R) one way to approach your situation. I read in the data and convert them into rowwise proportions, and then compute mean proportions by column. I add the means back to each count to get adjusted counts and new rowwise proportions. This nudges each customer's estimated proportion towards the mean proportion for each product. If you wanted a stronger nudge, you could use a multiple of the means (such as, 15*mean.props) instead.

d = read.table(text="id  red    blue    green
...
c3  4   8   1", header=TRUE)
tab = as.table(as.matrix(d[,-1]))
rownames(tab) = paste0("c", 0:3)
tab
#    red blue green
# c0  12    5     0
# c1   3    4     0
# c2   2   21     0
# c3   4    8     1
props = prop.table(tab, 1)
props
#           red       blue      green
# c0 0.70588235 0.29411765 0.00000000
# c1 0.42857143 0.57142857 0.00000000
# c2 0.08695652 0.91304348 0.00000000
# c3 0.30769231 0.61538462 0.07692308
mean.props = apply(props, 2, FUN=function(x){ weighted.mean(x, rowSums(tab)) })
mean.props
#        red       blue      green 
# 0.35000000 0.63333333 0.01666667 
adj.counts = sweep(tab, 2, mean.props, FUN="+");  adj.counts
#            red        blue       green
# c0 12.35000000  5.63333333  0.01666667
# c1  3.35000000  4.63333333  0.01666667
# c2  2.35000000 21.63333333  0.01666667
# c3  4.35000000  8.63333333  1.01666667
adj.props = prop.table(adj.counts, 1);  adj.props
#             red         blue        green
# c0 0.6861111111 0.3129629630 0.0009259259
# c1 0.4187500000 0.5791666667 0.0020833333
# c2 0.0979166667 0.9013888889 0.0006944444
# c3 0.3107142857 0.6166666667 0.0726190476

There are several results of this. One of which is that you now have non-zero estimates of the underlying probabilities of purchasing green products, even when a customer doesn't actually have any record of having purchased any green products yet. Another consequence is that you now have somewhat continuous values, whereas the original proportions were more discrete; that is, the set of possible estimates is less constricted, so a distance measure like the squared Euclidean distance might make more sense now.

We can visualize the data to see what happened. Because these are compositional data, we only actually have two pieces of information, and we can plot these in a single scatterplot. With most of the information in the red and blue categories, it makes sense to use those as the axes. You can see that the adjusted proportions (the red numbers) are shifted a little from their original positions.

windows()
  plot(props[,1], props[,2], pch=as.character(0:3),
       xlab="Proportion Red", ylab="Proportion Blue", xlim=c(0,1), ylim=c(0,1))
  points(adj.props[,1], adj.props[,2], pch=as.character(0:3), col="red")

At this point, you have data and a lot of people would begin by standardizing them. Again, because these are compositional data, I would run cluster analyses without doing any standardization—these values are already commensurate and standardization would destroy some of the relational information. In fact, from looking at the plot I think you really have only one dimension of information here. (At least in the sample dataset; your real dataset may well be different.) Unless, from a business point of view, you think it's important to recognize people who have any substantial probability of purchasing green products as a distinct cluster of customers, I would extract scores on the first principal component (which accounts for 99.5% of the variance in this dataset) and just cluster that.

pc.a.props = prcomp(adj.props[,1:2], center=T, scale=T)
cumsum(pc.a.props$sdev^2)/sum(pc.a.props$sdev^2)
# [1] 0.9946557 1.000000
pc.a.props$x
#           PC1         PC2
# c0 -1.7398975 -0.03897251
# c1 -0.1853614 -0.04803648
# c2  1.6882400 -0.06707115
# c3  0.2370189  0.15408015
library(mclust)
mc = Mclust(pc.a.props$x[,1])
summary(mc)
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm 
# ----------------------------------------------------
# 
# Mclust E (univariate, equal variance) model with 3 components:
# 
#  log.likelihood n df       BIC       ICL
#       -2.228357 4  6 -12.77448 -12.77448
# 
# Clustering table:
# 1 2 3 
# 1 2 1