Tienes razón en que la agrupación de k-means no debe hacerse con datos de tipos mixtos. Dado que k-means es esencialmente un simple algoritmo de búsqueda para encontrar una partición que minimice las distancias euclidianas al cuadrado dentro del clúster entre las observaciones agrupadas y el centroide del clúster, sólo debería utilizarse con datos en los que las distancias euclidianas al cuadrado fueran significativas.
Cuando sus datos consisten en variables de tipos mixtos, necesita utilizar la distancia de Gower. El usuario de CV @ttnphns tiene una gran descripción de la distancia de Gower ici . En esencia, se calcula una matriz de distancia para las filas de cada variable, utilizando un tipo de distancia apropiado para ese tipo de variable (por ejemplo, euclidiana para datos continuos, etc.); la distancia final de la fila $i$ a $i'$ es la media (posiblemente ponderada) de las distancias de cada variable. Una cosa que hay que tener en cuenta es que la distancia de Gower no es realmente una métrica . Sin embargo, con datos mixtos, la distancia de Gower es en gran medida el único juego en la ciudad.
En este punto, se puede utilizar cualquier método de clustering que pueda operar sobre una matriz de distancia en lugar de necesitar la matriz de datos original. (Tenga en cuenta que k-means necesita esta última). Las opciones más populares son partición en torno a los medoides (PAM, que es esencialmente lo mismo que k-means, pero utiliza la observación más central en lugar del centroide), varios enfoques de agrupación jerárquica (por ejemplo, mediana, enlace simple y enlace completo; con la agrupación jerárquica tendrá que decidir dónde ' cortar el árbol ' para obtener las asignaciones finales de los clusters), y DBSCAN que permite formas de cluster mucho más flexibles.
He aquí una sencilla R
demo (n.b., en realidad hay 3 clusters, pero los datos parecen ser en su mayoría 2 clusters apropiados):
library(cluster) # we'll use these packages
library(fpc)
# here we're generating 45 data in 3 clusters:
set.seed(3296) # this makes the example exactly reproducible
n = 15
cont = c(rnorm(n, mean=0, sd=1),
rnorm(n, mean=1, sd=1),
rnorm(n, mean=2, sd=1) )
bin = c(rbinom(n, size=1, prob=.2),
rbinom(n, size=1, prob=.5),
rbinom(n, size=1, prob=.8) )
ord = c(rbinom(n, size=5, prob=.2),
rbinom(n, size=5, prob=.5),
rbinom(n, size=5, prob=.8) )
data = data.frame(cont=cont, bin=bin, ord=factor(ord, ordered=TRUE))
# this returns the distance matrix with Gower's distance:
g.dist = daisy(data, metric="gower", type=list(symm=2))
Podemos empezar buscando sobre diferentes números de clusters con PAM:
# we can start by searching over different numbers of clusters with PAM:
pc = pamk(g.dist, krange=1:5, criterion="asw")
pc[2:3]
# $nc
# [1] 2 # 2 clusters maximize the average silhouette width
#
# $crit
# [1] 0.0000000 0.6227580 0.5593053 0.5011497 0.4294626
pc = pc$pamobject; pc # this is the optimal PAM clustering
# Medoids:
# ID
# [1,] "29" "29"
# [2,] "33" "33"
# Clustering vector:
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
# 1 1 1 1 1 2 1 1 1 1 1 2 1 2 1 2 2 1 1 1 2 1 2 1 2 2
# 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
# 1 2 1 2 2 1 2 2 2 2 1 2 1 2 2 2 2 2 2
# Objective function:
# build swap
# 0.1500934 0.1461762
#
# Available components:
# [1] "medoids" "id.med" "clustering" "objective" "isolation"
# [6] "clusinfo" "silinfo" "diss" "call"
Estos resultados pueden compararse con los de la agrupación jerárquica:
hc.m = hclust(g.dist, method="median")
hc.s = hclust(g.dist, method="single")
hc.c = hclust(g.dist, method="complete")
windows(height=3.5, width=9)
layout(matrix(1:3, nrow=1))
plot(hc.m)
plot(hc.s)
plot(hc.c)
El método de la mediana sugiere 2 (posiblemente 3) conglomerados, el simple sólo admite 2, pero el método completo podría sugerir 2, 3 o 4 a mi entender.
Por último, podemos probar el DBSCAN. Para ello, es necesario especificar dos parámetros: eps, la "distancia de alcance" (la proximidad que deben tener dos observaciones para estar vinculadas entre sí) y minPts (el número mínimo de puntos que deben estar conectados entre sí antes de estar dispuesto a llamarlos "cluster"). Una regla general para los minPts es utilizar uno más que el número de dimensiones (en nuestro caso 3+1=4), pero no se recomienda tener un número demasiado pequeño. El valor por defecto para dbscan
es 5; nos quedaremos con eso. Una forma de pensar en la distancia de alcance es ver qué porcentaje de las distancias es menor que un valor determinado. Podemos hacerlo examinando la distribución de las distancias:
windows()
layout(matrix(1:2, nrow=1))
plot(density(na.omit(g.dist[upper.tri(g.dist)])), main="kernel density")
plot(ecdf(g.dist[upper.tri(g.dist)]), main="ECDF")
Las propias distancias parecen agruparse en grupos visualmente discernibles de "más cerca" y "más lejos". Un valor de 0,3 parece distinguir más limpiamente entre los dos grupos de distancias. Para explorar la sensibilidad de la salida a diferentes opciones de eps, podemos probar también con .2 y .4:
dbc3 = dbscan(g.dist, eps=.3, MinPts=5, method="dist"); dbc3
# dbscan Pts=45 MinPts=5 eps=0.3
# 1 2
# seed 22 23
# total 22 23
dbc2 = dbscan(g.dist, eps=.2, MinPts=5, method="dist"); dbc2
# dbscan Pts=45 MinPts=5 eps=0.2
# 1 2
# border 2 1
# seed 20 22
# total 22 23
dbc4 = dbscan(g.dist, eps=.4, MinPts=5, method="dist"); dbc4
# dbscan Pts=45 MinPts=5 eps=0.4
# 1
# seed 45
# total 45
Utilizando eps=.3
da una solución muy limpia, que (al menos cualitativamente) coincide con lo que vimos en otros métodos anteriores.
Dado que no existe un clúster 1 , debemos tener cuidado de intentar hacer coincidir las observaciones que se denominan "cluster 1" de diferentes agrupaciones. En cambio, podemos formar tablas y si la mayoría de las observaciones denominadas "cluster 1" en un ajuste se denominan "cluster 2" en otro, veríamos que los resultados siguen siendo sustancialmente similares. En nuestro caso, las diferentes agrupaciones son en su mayoría muy estables y colocan las mismas observaciones en los mismos conglomerados cada vez; sólo difiere la agrupación jerárquica de vinculación completa:
# comparing the clusterings
table(cutree(hc.m, k=2), cutree(hc.s, k=2))
# 1 2
# 1 22 0
# 2 0 23
table(cutree(hc.m, k=2), pc$clustering)
# 1 2
# 1 22 0
# 2 0 23
table(pc$clustering, dbc3$cluster)
# 1 2
# 1 22 0
# 2 0 23
table(cutree(hc.m, k=2), cutree(hc.c, k=2))
# 1 2
# 1 14 8
# 2 7 16
Por supuesto, no hay garantía de que ningún análisis de clústeres recupere los verdaderos clústeres latentes en sus datos. La ausencia de las verdaderas etiquetas de los conglomerados (que estarían disponibles, por ejemplo, en una situación de regresión logística) significa que una enorme cantidad de información no está disponible. Incluso con conjuntos de datos muy grandes, los clusters pueden no estar lo suficientemente bien separados como para ser perfectamente recuperables. En nuestro caso, dado que conocemos la verdadera pertenencia a los clusters, podemos compararla con la salida para ver qué tan bien lo hizo. Como he señalado anteriormente, en realidad hay 3 conglomerados latentes, pero los datos dan la apariencia de 2 conglomerados:
pc$clustering[1:15] # these were actually cluster 1 in the data generating process
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# 1 1 1 1 1 2 1 1 1 1 1 2 1 2 1
pc$clustering[16:30] # these were actually cluster 2 in the data generating process
# 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
# 2 2 1 1 1 2 1 2 1 2 2 1 2 1 2
pc$clustering[31:45] # these were actually cluster 3 in the data generating process
# 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
# 2 1 2 2 2 2 1 2 1 2 2 2 2 2 2