Análisis de Conglomerados

En la literatura estadística recibe el nombre de Cluster Analysis, o en español Análisis de Conglomerados. En la literatura de Inteligencia Artificial se utiliza la expresión clasificación no supervisada.

La idea general es que tenemos una muestra de observaciones multivariada, y deseamos agruparla en grupos. Muy breve la respuesta pero: ¿qué son grupos? Imaginemos una imagen aérea fija de un patio de un colegio. En esta imagen los datos son las posiciones de los niños. ¿Se agrupan los niños formando grupos o todos están jugando con todos y los grupos son una consecuencia pasajera del juego (delanteros y defensas que aparecen agrupados en un ataque)?. Aunque parece claro y simple el problema. Estadísticamente, es un tema de investigación complicado y muy relevante. Dado que definir por ejemplo, ¿Qué es un grupo? ¿Cómo defino un grupo? ¿Cuántos grupos distintos hay en los datos? Estamos viendo el efecto del ruido o realmente hay una estructura debajo que la vemos en un entorno con ruido.

¿Qué quiere decir encontrar grupos? Se trata de clasificar las observaciones en grupos de modo que las observaciones de un mismo grupo sean lo más similares que podamos y que los grupos entre sí sean muy distintos. El número de procedimientos que se han propuesto en la literatura es muy grande. La mayor parte de ellos no tienen un modelo probabilístico debajo, no son procedimientos basados en modelo. Son métodos que esencialmente utilizan el concepto de proximidad. Valoran de distintas formas lo cercanos que están los puntos, dentro de un mismo grupo y entre distintos grupos. Es pues, el primer punto a tratar: ¿cómo cuantificamos lo cerca o lejos que están los distintos puntos?. También sería necesario, valorar cuando dos conjuntos de puntos son más o menos parecidos, proximos o similares.

Supongamos que ya hemos clasificado en distintos grupos. ¿Cómo comprobar si hemos tenido éxito al hacerlo? Cuando tenemos un análisis discriminante, tenemos una muestra donde sabemos a qué grupo pertenece el individuo y dónde lo hemos clasificado. Esto nos permite valorar si nuestro procedimiento clasificación en bueno o malo. Pero cuando no se cuenta con dicha información, no vamos a tener esta referencia que nos da la muestra de entrenamiento. Entonces, ¿cómo valorarlo?. Rousseeuw en 1990, propone un concepto conocido como silueta, el cual nos va a servir para ello. No es ni tan simple ni tan satisfactorio como en análisis discriminante (como es de esperar, pues tenemos menos información para trabajar).

Entre los muchos procedimientos de obtener los grupos a partir de los datos, los más utilizados se pueden clasificar en dos tipos: procedimientos jerárquicos y métodos de particionamiento.

Ilustremos el problema con ejemplos

Ejemplo 1

Simulación: Generamos tres muestras correspondientes a distribuciones normales bivariadas. La matriz de covarianzas vamos a suponer que es la matriz identidad. Los vectores de medias son \(\mu_1 = c(1, 1), \, \mu_2 = c(3.5, 3)\) y \(\mu_3 = c(6, 5.5)\). En resumen, tenemos un vector de datos \(X_i \sim N_d(\mu_i; I_{2\times 2})\). Vamos a generar 100 datos de cada grupo, para ello es necesario el paquete mvtnor.

library(mvtnorm)
# Generamos los datos
x1 <- rmvnorm(n = 100, mean = c(1, 1))
x2 <- rmvnorm(n = 100, mean = c(3.5, 4))
x3 <- rmvnorm(n = 100, mean = c(6, 5.5))
# Gráfica de los datos
limx <- c(-1, 8)
limy <- limx
plot(x1, xlim = limx, ylim = limy,col="blue",xlab="Normales bivariadas",ylab=" ")
points(x2, pch = 2,col="red")
points(x3, pch = 3,col="purple")

Se puede observar que hay tres grupos, pero estos no están claramente delimitados. Cabe resalta que nosotros no disponemos de esta información. Conocemos los valores que componen los vectores de datos, pero no conocemos el grupo al que podría pertenecer cada uno de ellos. Tampoco tenemos porqué tener prefijado el número de grupos. Aunque los datos son artificiales, ilustran bien el problema.

Ejemplo 2

Datos Ruspini: Cargamos los datos Ruspini que se encuentra dentro de la paquetería cluster

library(cluster)
data(ruspini)
plot(ruspini)

Son datos bivariados, donde se puede observar cómo se agrupan los puntos. Parece claro que podemos distinguir cuatro grupos. Para este conjunto de datos, no tenemos la clasificación real.

Disimilaridades

Disimilaridades entre observaciones

Empezamos tratando el problema de cuantificar el grado de proximidad o similaridad entre dos puntos. Tradicionalmente este tema en Matemáticas se ha formalizado a través del concepto de distancia o métrica. Una métrica es una función que a cada par de puntos \((x,y)\) le asocia un valor positivo de modo que cuando mayor es más distantes son, más alejados están. En la formalización matemática de un concepto intuitivo, es necesario pedir que se verifiquen ciertos axiomas que resulten razonables. En concreto la función de distancia \(d\) se dice que es una métrica si:

  • \(d(x,y) \geq 0\).
  • \(d(x,x) =0\).
  • \(d(x,y) = d(y,x)\).
  • \(d(x,z) \leq d(x,y) + d(y,z)\).

Las distancias más utilizadas en Análisis de Conglomerados son: la distancia euclídea y la distancia de Manhattan. Para dos vectores \(x\) e \(y\), la distancia euclídea se define como: \[d(x,y) = \sqrt{\sum_{k=1}^p (x_k-y_k)^2}\], con \(x,y \in \mathbb{R}^d\). La distancia de Manhattan esta dada por: \[d(x,y) = \sum_{k=1}^p \mid x_k-y_k \mid. \]

Las distancias euclídea y de Manhattan son adecuadas cuando trabajamos con variables continuas y que además estén en una misma escala. Notemos que cada una de las componentes del vector pesan igualmente. Si tenemos variables que no están igualmente escaladas, estas distancias pueden pesar más unas variables que otras y no por lo diferentes que sean entre los individuos sino simplemente por su escala.

Con mucha frecuencia nos encontramos trabajando con variables que aún siendo continuas están medidas en escalas muy diversas o bien tenemos variables que son continuas, otras que son binarias, otras categóricas con más de dos categorías o bien variable ordinales. En resumen, todos los posibles tipos de variables simultáneamente considerados. Es lo habitual. Una variable binaria la codícamos habitualmente como 1 y 0 indicando presencia o ausencia del atributo que estemos considerando. En una variable categórica la codificación es completamente arbitraria y por lo tanto no tiene sentido la aplicación de una de estas distancias. Todo esto plantea el hecho de que no es recomendable, ni tampoco fácil, en un banco de datos con distintos tipos de variables considerar una métrica o distancia, esto es, algo que verifique las propiedades anteriores. Lo que se hace es proponer medidas entre los vectores de características que tienen algunas de las propiedades y son, además, razonables. Por ello no hablaremos, en general, de una distancia o una métrica, sino de una medida de disimilaridad. Para eso, se han propuesto distintas medidas de disimilaridad entre variables cualitativas (binarias simétricas o asimétricas, cualitativas, ordinales) y cuantitativas. Hablaremos a continuación sobre la función daisy de la librería cluster. Es una opción muy genérica y bastante razonable.

Supongamos que el individuo o caso de estudio se describe mediante \(p\) variables binarias. Es natural construir la tabla de contingencia \(2\times 2\) que aparece en la siguiente tabla donde las filas corresponden con el individuo \(i\) y las columnas con el individuo \(j\).

1 0
1 \(A\) \(B\) \(A+B\)
0 \(C\) \(D\) \(C+D\)
\(A+B\) \(B+D\) \(m=A+B+C+D\)

Según la tabla los individuos \(i\) y \(j\) coincidirán en la presencia de \(A\) atributos y en la no presencia de \(D\) atributos. Tenemos \(B\) atributos en \(i\) que no están en \(j\) y \(C\) atributos que no están en \(i\) pero sí que están en \(j\).

El total de variables binarias es \(m= A+B+C+D\). Basándonos en esta tabla se pueden definir distintas medidas de disimilaridad. Vamos a considerar dos situaciones distintas, variables binarias simétricas y variables binarias no simétricas.

Si las variables son todas binarias simétricas, utilizamos como disimilaridad: \[d(i, j) = \dfrac{B + C}{m} = 1-\dfrac{A + D}{m}.\] La interpretación de esta medida de disimilaridad es simple. Dos individuos son tanto más disimilares cuantas más variables binarias tienen distintas.

Ahora, si las variables son todas binarias asimétricas, utilizamos como disimilaridad: \[d(i, j) = \dfrac{B + C}{A+B+C} = 1-\dfrac{A }{A+B+C}.\] La interpretación es que, simplemente no consideramos los acoplamientos negativos.

Consideremos ahora el caso de variables categóricas con más de dos categorías. Si todas las variables son de este tipo y tenemos un total de \(p\) variables, entonces los individuos \(i\) y \(j\) son tanto más disimilares cuanto más variables categóricas son distintas. Si denotamos por \(u\) el número de variables en las que coinciden los dos individuos, entonces la medida de disimilaridad sería: \[d(i, j) = \dfrac{p-u}{p}.\] Finalmente veamos cómo tratar las variables ordinales. Lo que se hace para variables de este tipo es, transformarlas al intervalo \([0, 1]\). Si \(x_{ij}\) denota la \(j\)-ésima variable del \(i\)-ésimo individuo, entonces consideramos la transformación: \[y_{ij} = \dfrac{x_{ij-1}}{M_j-1},\] donde \(1,2,\dots,M_j\) los valores que puede tomar la \(j\)-ésima variable ordinal. Lo que estamos haciendo con este procedimiento es transformar la variable ordinal en una variable numérica con una escala común. En la medida en que el número de categorías sea mayor esta transformación tendría más sentido.

Hemos visto cómo tratar cada tipo de variable aisladamente. El problema es combinar todas ellas en una sola medida de disimilaridad. La función daisy del paquete cluster utiliza la siguiente medida: \[ d(i,j)= \dfrac{\displaystyle\sum_{k=1}^p \delta_{ij}^{(k)}d_{ij}^{(k)}}{\displaystyle\sum_{k=1}^d \delta_{ij}^{(k)}}. \] Si todas las variables son categóricas, entonces nos da el número de acoplamientos del total de pares disponibles. Cuando todas las variables con las que trabajamos son numéricas, la medida de disimilaridad es la distancia de Manhattan donde cada variable está normalizada.

Dado un conjunto de datos \(x_i\) con \(i = 1,2,n\) tendremos, utilizando algunas de las medidas de disimilaridades comentadas, una matriz de dimensión \(n \times n\). Con esta matriz cuantificamos la disimilaridad que hay entre los elementos originales de la muestra.

Disimilaridades entre grupos de observaciones

En algunos procedimientos de agrupamiento (en particular, los jerárquicos) vamos a necesitar calcular disimilaridades entre conjuntos disjuntos de las observaciones originales. Estas disimilaridades las podemos calcular a partir de las disimilaridades comentadas en la sección anterior.

Supongamos que tenemos un banco de datos con \(n\) individuos. Sean \(A\) y \(B\) dos subconjuntos disjuntos de la muestra. ¿Cómo podemos definir una disimilaridad entre \(A\) y \(B\) partiendo de las disimilaridades entre los datos individuales? Para resolver este problema, se han propuesto muchos procedimientos en la literarura. Si denotamos la disimilaridad entre \(A\) y \(B\) como \(d(A,B)\) entonces las disimilaridades más utilizadas son las siguientes:

  • Enlace simple: La disimilaridad entre los dos grupos es el mínimo de las disimilaridades entre todos los posibles pares. \[d(A,B) = \min_{\{a\in A, \, b \in B\}}d(a,b)\]

  • Enlace completo Ahora tomamos como disimilaridad entre los grupos como el máximo de las disimilaridadesentre entre todos los posibles pares. \[d(A,B) = \max_{\{a\in A, \, b \in B\}}d(a,b)\]

  • Promedio La disimilaridad entre los dos grupos es el promedio de las disimilaridades entre todos los posibles pares. \[d(A,B) = \dfrac{1}{|A| \times |B|} \sum_{a \in A,\, b \in B} d(a,b)\]

donde \(|A|\) es el cardinal del conjunto \(A\).

Es importante notar que solamente necesitamos conocer las disimilaridades entre los individuos para poder calcular las disimilaridades entre grupos de individuos. En la siguiente sección nos vamos a ocupar de los métodos jerárquicos en los cuales es fundamental el procedimiento que elijamos para calcular distintas entre grupos.

Cluster jerárquico

La idea de estos procedimientos es construir una jerarquía de particiones del conjunto de índices. Supongamos que \(\{C_1,C_2,\dots,Cr\}\) es una partición de un conjunto de observaciones.

Dada una partición de un conjunto de observaciones, podemos calcular la matriz \(r \times r\) que en la posición \((i,j)\) tiene la disimilaridad entre el conjunto \(C_i\) y \(C_j\), según alguno de los procedimientos antes indicados.

Veamos paso a paso el proceso que se realiza en un cluster jerárquicos:

Hay una representación gráfica muy utilizada para describir los resultados de un cluster jerárquico como el que acabamos de describir. Esta representación tiene el nombre de dendograma. En el dendograma se va mostrando a qué valor de la medida de disimilaridad se produce la unión de los grupos y simultáneamente qué grupos se están uniendo para esa disimilaridad. También nos permite una valoración rápida de cuántos grupos puede haber en el banco de datos. Simplemente trazando una linea horizontal a la altura en que tengamos el número de grupos que pensamos que puede haber.

Ejemplo paso a paso

Supongamos tenemos una base de datos con 7 observaciones. Los pasos que deberíamos seguir para conseguir un cluster jerárquico, son:

    1. (Paso 1) Supongamos que calculamos por alguno de los métodos anteriores la matriz de disimilaridades entre estas 7 observaciones y obtenemos:
\(x_1\) \(x_2\) \(x_3\) \(x_4\) \(x_5\) \(x_6\) \(x_7\)
\(x_1\) 0
\(x_2\) 2.15 0
\(x_3\) 0.7 1.53 0
\(x_4\) 1.07 1.14 0.43 0
\(x_5\) 0.85 1.38 0.21 0.29 0
\(x_6\) 1.16 1.01 0.55 0.22 0.41 0
\(x_7\) 1.56 2.83 1.86 2.04 2.02 2.05 0
    1. (Paso 2) Observamos que \(\min\{d(c_i,c_j)\}=d(x_3,x_5)=0.21\), por lo que el primer cluster que se forma es el cluster \((C,E)\).
    1. (Paso 1) La matriz de disimilaridades en este paso sería:
\(x_1\) \(x_2\) \((x_3,x_5)\) \(x_4\) \(x_6\) \(x_7\)
\(x_1\) 0
\(x_2\) 2.15 0
\((x_3,x_5)\) 0.7 1.38 0
\(x_4\) 1.07 1.14 0.29 0
\(x_6\) 1.16 1.01 0.41 0.22 0
\(x_7\) 1.56 2.83 1.86 2.04 2.05 0
    1. (Paso 2) Ahora, \(\min\{d(c_i,c_j)\}=d(x_4,x_6)=0.22\), formándose el cluster \((x_4,x_6)\).
    1. (Paso 1) La matriz de disimilaridades en este paso sería:
\(x_1\) \(x_2\) \((x_3,x_5)\) \((x_4,x_6)\) \(x_7\)
\(x_1\) 0
\(x_2\) 2.15 0
\((x_3,x_5)\) 0.7 1.38 0
\((x_4,x_6)\) 1.07 1.01 0.29 0
\(x_7\) 1.56 2.83 1.86 2.04 0
    1. (Paso 2) Para el nuevo paso, \(\min\{d(c_i,c_j)\}=d((x_3,x_5),(x_4,x_6))=0.29\), formándose el cluster \(((x_3,x_5),(x_4,x_6))\).
    1. (Paso 1) La matriz de disimilaridades en este paso sería:
\(x_1\) \(x_2\) \((x_3,x_5),(x_4,x_6)\) \(x_7\)
\(x_1\) 0
\(x_2\) 2.15 0
\((x_3,x_5),(x_4,x_6)\) 0.7 1.01 0
\(x_7\) 1.56 2.83 1.86 0
    1. (Paso 2) En este paso, \(\min\{d(c_i,c_j)\}=d(x_1,(x_3,x_5),(x_4,x_6))=0.7\), formándose el el cluster \((x_1,(x_3,x_5),(x_4,x_6))\).
    1. (Paso 1) La matriz de disimilaridades en este paso sería:
\((x_1,(x_3,x_5),(x_4,x_6))\) \(x_2\) \(x_7\)
\((x_1,(x_3,x_5),(x_4,x_6))\) 0
\(x_2\) 1.01 0
\(x_7\) 1.56 2.83 0
    1. (Paso 2) En este paso, \(\min\{d(c_i,c_j)\}=d(x_2,(x_1,(x_3,x_5),(x_4,x_6)))=1.01\), formándose el el cluster \((x_2,(x_1,(x_3,x_5),(x_4,x_6)))\).
    1. (Paso 1) La matriz de disimilaridades en este paso sería:
\((x_2,(x_1,(x_3,x_5),(x_4,x_6)))\) \(x_7\)
\((x_2,(x_1,(x_3,x_5),(x_4,x_6)))\) 0
\(x_7\) 1.56 0
    1. (Paso 3) Este seía el último paso, pues el próximo paso, evidentemente, se tendrá un único cluster formado por las 7 observaciones.

El dendograma asosiado de este cluster jerárquico será:

Ejemplo 1 en R

Para este primer ejemplo, generaremos los datos continuos,

# Datos simulados
library(mvtnorm)
x1 <- rmvnorm(n = 20, mean = c(1, 1))
x2 <- rmvnorm(n = 20, mean = c(3.5, 4))
x3 <- rmvnorm(n = 20, mean = c(6, 5.5))
x  <- c(x1,x2,x3)
plot(x,col="blue",xlab=" ",ylab=" ")

Utilizaremos variables numéricas para ete primer ejemplo, para simplemente demostrar cómo se realiza la agrupación en R. Las variables categóricas, como ya lo vimos, requieren un tratamiento especial, que ejemplificaremos más adelante.

Un primer paso necesario antes de iniciar el proceso es, escalar todas las variables numéricas. Escalar significa que cada variable ahora tendrá una media cero y una desviación estándar uno. Lo ideal es que quieras que una unidad en cada coordenada represente el mismo grado de diferencia. Escalar las variables hace que la desviación estándar sea la unidad de medida en cada coordenada. Esto se hace para evitar que el algoritmo de agrupamiento dependa de una unidad variable arbitraria. Puede escalar/estandarizar los datos usando la función scale de R:

dat <-scale(x)
plot(dat,col="red",xlab=" ",ylab=" ")

Hay varias funciones disponibles en R para realizar cluster jerárquicos. Los más utilizados comunmente son: hclust y agnes.

Para la función hclust, se requieren los valores de distancia que se pueden calcular en R utilizando la función dist. La medida predeterminada para la función dist es la Euclidiana, sin embargo, puede cambiarla con el argumento del método. Con esto, también necesitamos especificar el método que queremos usar (es decir, “simple”,“completo”,“promedio”).

d   <- dist(dat, method = "euclidean")
hc1 <- hclust(d, method = "complete" )
plot(hc1,cex=.5,hang=-1)

Otra alternativa es la función agnes. Ambas funciones son bastante similares; sin embargo, con la esta función también puede obtener el coeficiente de aglomeración, que mide la cantidad de estructura de agrupamiento encontrada (los valores más cercanos a 1 sugieren una estructura de agrupación fuerte).

hc2 <- agnes(dat, method = "complete")
pltree(hc2,cex=0.5)

Además, puedes explorar con que método consigues mejor coeficiente de aglomeración

library(purrr)
m <- c("average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")
ac <- function(x) {
  agnes(dat, method = x)$ac
}
map_dbl(m, ac) 
##   average    single  complete      ward 
## 0.9878299 0.8564875 0.9935320 0.9981470

Si suponemos que tenemos cuatro grupos, la división, sería:

pltree(hc2, hang=-1, cex = 0.6)
rect.hclust(hc2, k = 4, border = 1:4)

c1 <-cutree(hc2, 4)
plot(dat, col = c1)

Ejemplo 2 en R

Para este segundo ejemplo, generaremos los datos categoricos,

x1 <- rbinom(n = 15,10,0.5)
x2 <- rbinom(n = 15,15,0.5)

data  <- matrix(c(x1,x2),ncol=1)
plot(data,col="green",xlab=" ",ylab=" ")

hc1 <-  agnes(daisy(data), diss = TRUE, method = "complete")
pltree(hc1,cex=0.5,hang=-1,)

Si suponemos que tenemos dos grupos, la división, sería:

pltree(hc1, hang=-1, cex = 0.6)
rect.hclust(hc1, k = 2, border = 2:3)

cc <-cutree(hc1, 2)
plot(data, col = cc)

Ejemplo 3 en R

Consideremos los datos ruspini. Apliquemos un cluster de conglomerados jerárquico utilizando como disimilaridad entre grupos el promedio de las disimilaridades y como medida de disimilaridad la distancia euclídea.

agj <- agnes(ruspini, metric = "euclidean", method = "average")
plot(agj, which = 2)

Supongamos que decidimos quedarnos con 4 grupos. Ahora podemos representar los datos de modo que ponemos en un mismo color los datos que hemos clasificado en un grupo.

c1 <-cutree(agj, 4)
plot(ruspini, col = c1)

Ejemplo 4 en R

Veamos un análisis de Conglomerados jerárquico. Los datos son los porcentajes de votos que recibió un candidato republicano en las elecciones americanas entre los años 1856 y 1976. Cada observación corresponde a un estado y las variables son las distintas elecciones.

library(cluster)
data(votes.repub)
agj1 <- agnes(votes.repub, metric = "manhattan", stand = TRUE)
plot(agj1,which=2,xlab="",ylab=" ")

Ahora si cambiamos la matriz de distancias, el procedimiento para el cálculo de las disimilaridades entre grupos es:

agj2 <- agnes(daisy(votes.repub), diss = TRUE, method = "complete")
plot(agj2, which = 2)

Un procedimiento jerárquico tiene interés cuando la base de datos es pequeña. En otros casos es difícil de interpretar un dendograma o de seleccionar el número de grupos que hay. En todo caso, hasta el momento lo vemos como parte de un análisis exploratorio para decidir el número de grupos y posteriormente utilizar algún método de particionamiento como los que vamos a ver en las próximas clases.