Ejercicio Monte Carlo con conglomerados

Basado en la propuesta de Halkidi, Batistakis y Varzirgiannis¹.

library(mclust)
library(clue)
library(MASS)
library(fpc)
library(clusterCrit)

Suponga un conjunto de datos dado, X, que contiene \(100\) vectores tridimensionales.

Los puntos de \(X\) forman cuatro conglomerados de \(25\) puntos cada uno. Cada conglomerado está generado por una distribución normal. Los vectores medios de las cuatro distribuciones son \([0.2, 0.2, 0.2]^t\), \([0.5, 0.2, 0.8]^t\), \([0.5, 0.8, 0.2]^t\) y \([0.8, 0.8, 0.8]^t\). Las matrices de covarianza de estas distribuciones son todas iguales a \(0,2I\) , donde \(I\) es una matriz de identidad de \(3 \times 3\).

Agrúpese independientemente el conjunto de datos \(X\) en cuatro grupos según una partición \(P\) de tal modo que los primeros \(25\) vectores pertenezcan al primer grupo \(P1\), los \(25\) siguientes pertenezcan al segundo grupo \(P2\), los \(25\) siguientes al tercer grupo \(P3\) y los últimos \(25\) vectores al grupo \(P4\).

set.seed(123)
sigma <- sqrt(0.2)
means <- list(c(0.2, 0.2, 0.2), c(0.5, 0.2, 0.8), c(0.5, 0.8, 0.2), c(0.8, 0.8, 0.8))
x1 <- matrix(data = rnorm(75, mean = 0.2, sd = sigma), nrow = 25, byrow = F)
x2 <- matrix(data = c(rnorm(25, mean = 0.5, sd = sigma),
                      rnorm(25, mean = 0.2, sd = sigma),
                      rnorm(25, mean = 0.8, sd = sigma)), 
             nrow = 25, byrow = F)
x3 <- matrix(data = c(rnorm(25, mean = 0.5, sd = sigma),
                      rnorm(25, mean = 0.8, sd = sigma),
                      rnorm(25, mean = 0.2, sd = sigma)), 
             nrow = 25, byrow = F)
x4 <- matrix(data = rnorm(75, mean = 0.8, sd = sigma), nrow = 25, byrow = F)

X <- rbind(x1, x2, x3, x4)
P <- rep(1:4, each = 25)

Sea una visualización de \(X\):

Xgraph <- as.data.frame(X)
g1 <- ggplot() +
  geom_point(mapping = aes(V1, V2, color = factor(P)), data = Xgraph)
g2 <- ggplot() +
  geom_point(mapping = aes(V1, V3, color = factor(P)), data = Xgraph)
g3 <- ggplot() +
  geom_point(mapping = aes(V2, V3, color = factor(P)), data = Xgraph)
g1 / g2 / g3

Ejecute un algoritmo de agrupación k-means para \(k = 4\) grupos y suponga que \(C\) es la estructura de agrupación resultante.

set.seed(123)
C <- kmeans(X, centers = 4)

Sea una visualización de \(X\) respecto a \(C\):

Xgraph <- as.data.frame(X)
g1 <- ggplot() +
  geom_point(mapping = aes(V1, C$cluster, color = factor(P)), data = Xgraph)
g2 <- ggplot() +
  geom_point(mapping = aes(V2, C$cluster, color = factor(P)), data = Xgraph)
g3 <- ggplot() +
  geom_point(mapping = aes(V3, C$cluster, color = factor(P)), data = Xgraph)
g1 / g2 / g3

Obsérvese que aún coincidiendo, k-means coloca el número de cada algoritmo de manera aleatoria:

table(P, C$cluster)

##    
## P    1  2  3  4
##   1  5  1  4 15
##   2 12  0 10  3
##   3  2  3  4 16
##   4  4 16  3  2

Considere que \(C = {C_1, \dots, C_m}\) es una estructura de agrupamiento de un conjunto de datos \(X\) y \(P = {P_1, \dots, P_s}\) es una partición definida de los datos.

Refiérase a un par de puntos \((x_v, x_u)\) del conjunto de datos utilizando los siguientes términos:

\(II\): si ambos puntos pertenecen al mismo conglomerado de la estructura C y al mismo grupo de partición P.
\(ID\): si los puntos pertenecen al mismo conglomerado de C y a distintos grupos de P.
\(DI\): si los puntos pertenecen a distintos conglomerados de C y al mismo grupo de P.
\(DD\): si ambos puntos pertenecen a distintos conglomerados de C y a distintos grupos de P.

Suponiendo ahora que a, b, c y d son el número de pares \(II, ID, DI\) y \(DD\) respectivamente, entonces \(a + b + c + d = M\), que es el número máximo de todos los pares del conjunto de datos, es decir, \(M = n(n - 1)/2\) donde \(n\) es el número total de puntos del conjunto de datos.

Ahora se pueden definir los siguientes índices para medir el grado de similitud entre \(C\) y \(P\):

Estadístico de Rand: \(R = (a + d)/M\),
Coeficiente de Jaccard: \(J = a/(a + b + c)\),

Los dos índices anteriores toman valores entre \(0\) y \(1\), y se maximizan cuando \(m = s\), como es el caso actual en que el número de conglomerados buscado es igual al número de particiones.

Índice de Folkes y Mallows: \(FM = \frac{a}{\sqrt{m1m2}} = \sqrt{\frac{a}{a + b}\cdot \frac{a}{a + c}}\), donde \(m1 = (a + b)\), \(m2 = (a + c)\).

Para los tres índices anteriores se ha comprobado que los valores altos de los índices indican una gran similitud entre C y P. Cuanto más altos son los valores de estos índices, más similares son C y P. Siendo que se definieron los cuatro grupos con medias diferentes e igualñ dispersión, se espera que se parezcan.

Otro índice es:

Estadístico \(\hat{\Gamma}\) de Huberts: \(huberts = corr(X, Y) = \frac{\sum_{i<j}{(X_{i,j}-\mu_{x})(Y_{i,j}-\mu_{y}})}{M\sigma_{x}\sigma_{y}}\)

Los valores altos de este índice indican una gran similitud entre X e Y.

Obsérvese que \(X_{i,j}, Y_{i,j}\) son matrices binarias (indicadoras) de dimensión \(n \times n\) acerca de si, para la posición i menor a j de algún vector, los valores son iguales.

Todas estas estadísticas presentan distribuciones con una alta probabilidad de observar valores extremos en su cola derecha, y este comportamiento se analiza bajo el supuesto de que los eventos generadores de los datos son aleatorios, es decir, ocurren sin patrones o tendencias subyacentes.

Calcule los valores de los índices para la agrupación \(C\) y la partición \(P\).

calcule_pares <- function(congl1, congl2) {
  n <- length(congl1)
  a <- b <- c <- d <- 0
  
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      if (congl1[i] == congl1[j] && congl2[i] == congl2[j]) {
        a <- a + 1  # SS
      } else if (congl1[i] == congl1[j] && congl2[i] != congl2[j]) {
        b <- b + 1  # SD
      } else if (congl1[i] != congl1[j] && congl2[i] == congl2[j]) {
        c <- c + 1  # DS
      } else {
        d <- d + 1  # DD
      }
    }
  }
  list(a = a, b = b, c = c, d = d, M = a + b + c + d)
}

rand_statistic <- function(a, d, M) {
  (a + d) / M
}

jaccard_coefficient <- function(a, b, c) {
  a / (a + b + c)
}

fowlkes_mallows_index <- function(congl1, congl2, a, b, c) {
  table <- table(congl1, congl2)
  sum_squares <- sum(table^2)
  sum_congl1 <- sum(rowSums(table)^2)
  sum_congl2 <- sum(colSums(table)^2)
  return(sum_squares / sqrt(sum_congl1 * sum_congl2))
}

huberts_gamma <- function(congl1, congl2, a, b, c) {
  n <- length(congl1)
  M = n*(n-1)/2
  X <- matrix(0, n, n)
  Y <- matrix(0, n, n)
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      X[i, j] <- ifelse(congl1[i] == congl1[j], 1, 0)
      Y[i, j] <- ifelse(congl2[i] == congl2[j], 1, 0)
    }
  }
  mu_congl1 <- sum(X)/M
  mu_congl2 <- sum(Y)/M
  sigma_congl1 <- sqrt((a + b)/M - ((a + b)/M)^2)
  sigma_congl2 <- sqrt((a + c)/M - ((a + c)/M)^2)
  
  sum_val <- 0
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      sum_val <- sum_val + ((X[i, j] - mu_congl1)*(Y[i, j] - mu_congl2))
    }
  }
  gamma <- sum_val / (M * sigma_congl1 * sigma_congl2)
  return(gamma)
}

res <- calcule_pares(C$cluster, P)
R <- rand_statistic(a = res$a, d = res$d, M = res$M)
J <- jaccard_coefficient(a = res$a, b = res$b, c = res$c)
FM <- fowlkes_mallows_index(congl1 = C$cluster, congl2 = P)
hubert <- huberts_gamma(congl1 = C$cluster, congl2 = P, a = res$a, b = res$b, c = res$c)

Se obtiene \(R = 0.698\), \(J = 0.249\), \(FM = 0.422\) y \(hubert = 0.198\).

Se ha dicho que valores altos implican semejanza entre los dos pares. Ninguno es alto, debido a la alta dispersión de los datos globulares.

Si bien se han realizado todos los cálculos de manera manual, mediante el paquete clusterCrit se obtiene el calculo de manera rápida:

extCriteria(part1 = C$cluster, part2 = P, crit = list("rand", "jaccard", "folkes_mallows", "hubert"))

## $rand
## [1] 0.6983839
## 
## $jaccard
## [1] 0.248994
## 
## $folkes_mallows
## [1] 0.3989342
## 
## $hubert
## [1] 0.197914

Para conocer los cuatro datos, \(II, ID, DI\) y \(DD\) se usa:

clusterCrit::concordance(C$cluster, P)

##     y    n
## y 495  788
## n 705 2962

En segundo lugar, genere \(100\) conjuntos de datos \(X_i, i = 1, \dots, 100\), cada uno de ellos constando de \(100\) vectores aleatorios, en \(3\) dimensiones, utilizando la distribución uniforme. De acuerdo con la partición \(P\) definida anteriormente para cada \(X_i\), asigne los primeros \(25\) de sus vectores a \(P1\) y el segundo, tercer y cuarto grupo de \(25\) vectores a \(P2\), \(P3\) y \(P4\) respectivamente.

# Genera 100 conjuntos de datos nuevos con distribución uniforme. 
set.seed(123)
numero_conjuntos <- 100
numero_vectores <- 100
X_multiple <- lapply(1:numero_conjuntos, 
              function(i) matrix(runif(numero_vectores * 3), ncol = 3))

A continuación, ejecute k-means i veces, una vez por cada \(X_i\), para definir las respectivas estructuras de agrupación de los conjuntos de datos. Denomínense \(C_i\). Para cada uno de ellos calcule los valores de los índices \(R_i, J_i, FM_i, hubert_i, i = 1, \dots, 100\).

R_values <- numeric(numero_conjuntos)
J_values <- numeric(numero_conjuntos)
FM_values <- numeric(numero_conjuntos)
hubert_values <- numeric(numero_conjuntos)

set.seed(1515)
for (i in 1:numero_conjuntos) {
  nuevo_X <- X_multiple[[i]]
  kmeans_nuevo <- kmeans(nuevo_X, centers = 4)
  res <- calcule_pares(kmeans_nuevo$cluster, P)
  R_values[i] <- rand_statistic(a = res$a, d = res$d, M = res$M)
  J_values[i] <- jaccard_coefficient(a = res$a, b = res$b, c = res$c)
  FM_values[i] <- fowlkes_mallows_index(congl1 = kmeans_nuevo$cluster, congl2 = P)
  hubert_values[i] <- huberts_gamma(congl1 = kmeans_nuevo$cluster, congl2 = P,
                                    a = res$a, b = res$b, c = res$c)
}

Se obtiene, para el primer subconjunto de datos, \(R = 0.621\), \(J = 0.137\), \(FM = 0.271\) y \(hubert = -0.011\).

Sea una visualización de \(X_1\) respecto a \(C_1\):

set.seed(1515)
nuevo_X <- X_multiple[[1]]
kmeans_nuevo <- kmeans(nuevo_X, centers = 4)
Xgraph <- as.data.frame(nuevo_X)
g1 <- ggplot() +
  geom_point(mapping = aes(V1, kmeans_nuevo$cluster, color = factor(P)), data = Xgraph)
g2 <- ggplot() +
  geom_point(mapping = aes(V2, kmeans_nuevo$cluster, color = factor(P)), data = Xgraph)
g3 <- ggplot() +
  geom_point(mapping = aes(V3, kmeans_nuevo$cluster, color = factor(P)), data = Xgraph)
g1 / g2 / g3

Se espera que no tengan los índices ningún valor alto.

Fije el nivel de significación \(\alpha = 0.05\) y compare estos valores con los valores \(R, J, FM \text{ y } hubert\) correspondientes a X.

La hipótesis nula que se trabaja es que los datos ocurren sin patrones o tendencias subyacentes.

Rechace la hipótesis nula si \((1- \alpha)\cdot n = (1 - 0.05)100 = 95\) o menos valores de \(R_i, J_i, FM_i, hubert_i\) son mayores que los valores correspondientes de \(R, J, FM, hubert\). No la rechace de lo contrario.

nivel_significancia <- 0.05
umbral <- (1 - nivel_significancia) * numero_conjuntos

reject_R <- sum(R_values < R) > umbral
reject_J <- sum(J_values < J) > umbral
reject_FM <- sum(FM_values < FM) > umbral
reject_hubert <- sum(hubert_values < hubert) > umbral

¿Se rechaza la hipótesis nula para …

R? Cierto

J? Cierto

FM? Cierto

hubert? Cierto

Preguntas

¿Qué se hace cuando se pregunta sum(R_values < R) > umbral?

¿Qué distribución parecen tener los valores de \(R_i, J_i, FM_i, hubert_i\)?

On Clustering Validation Techniques. Journal of Intelligent Information Systems, 17:2/3, 107–145, 2001↩︎

Ejercicio Monte Carlo con conglomerados

Fernando López-Torrijos F.