Integrantes:

Caso aplicativo

El conjunto de datos USArrests contiene información sobre el número de delitos (agresión, asesinatos y violación) por cada 100,000 residentes para cada uno de los 50 estados de USA (1973). También se proporciona el porcentaje de la población que vive en zonas urbanas (Urbanpop). Se pretende estudiar si existe una agrupación subyacente de los estados mediante K-medoids clustering, ya que se sospecha de la presencia de outliers.

El siguiente código da inicio a las tareas de preprocesamiento de datos:

Lectura de Datos

data("USArrests")
DT::datatable(USArrests)
str(USArrests)
## 'data.frame':    50 obs. of  4 variables:
##  $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
##  $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
##  $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
##  $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...

Análisis descriptivo

library(funModeling)
df_status(USArrests)
##   variable q_zeros p_zeros q_na p_na q_inf p_inf    type unique
## 1   Murder       0       0    0    0     0     0 numeric     43
## 2  Assault       0       0    0    0     0     0 integer     45
## 3 UrbanPop       0       0    0    0     0     0 integer     36
## 4     Rape       0       0    0    0     0     0 numeric     48
summary(USArrests)
##      Murder          Assault         UrbanPop          Rape      
##  Min.   : 0.800   Min.   : 45.0   Min.   :32.00   Min.   : 7.30  
##  1st Qu.: 4.075   1st Qu.:109.0   1st Qu.:54.50   1st Qu.:15.07  
##  Median : 7.250   Median :159.0   Median :66.00   Median :20.10  
##  Mean   : 7.788   Mean   :170.8   Mean   :65.54   Mean   :21.23  
##  3rd Qu.:11.250   3rd Qu.:249.0   3rd Qu.:77.75   3rd Qu.:26.18  
##  Max.   :17.400   Max.   :337.0   Max.   :91.00   Max.   :46.00
library(ggplot2)
library(hrbrthemes)
library(viridis)

data_box %>%
  ggplot( aes(x=delito, y=valor, fill=delito)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6) + #paleta de colores
    geom_jitter(color="black", size=0.4, alpha=0.9) + # presentan las observaciones y las desfasa para que no se sobrepongan.
    theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Distribución de los delitos en USA") +
    xlab("")

Las variables a estudiar presentan heterogeneidad entre sus varianzas y desigualdad en sus rangos. Por ello, en este caso práctico es necesario estandarizar los datos.

Analizar la estandarización

datos<-as.data.frame(scale(USArrests))
summary(datos)
##      Murder           Assault           UrbanPop             Rape        
##  Min.   :-1.6044   Min.   :-1.5090   Min.   :-2.31714   Min.   :-1.4874  
##  1st Qu.:-0.8525   1st Qu.:-0.7411   1st Qu.:-0.76271   1st Qu.:-0.6574  
##  Median :-0.1235   Median :-0.1411   Median : 0.03178   Median :-0.1209  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.7949   3rd Qu.: 0.9388   3rd Qu.: 0.84354   3rd Qu.: 0.5277  
##  Max.   : 2.2069   Max.   : 1.9948   Max.   : 1.75892   Max.   : 2.6444
data_box_scale %>%
  ggplot( aes(x=delito, y=valor, fill=delito)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6) +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Distribución de los delitos en USA") +
    xlab("")

Hallando el valor óptimo de K o el número de cluster

Analizamos la matriz de distancias y el dendograma correspondiente para tener una noción de cuántos grupos significativamente distintos podríamos tener en el conjunto de datos. Por otro lado, se uso la distancia de Manhattan para el análisis de disimilaridad ya que esta medida se ve menos afectada por outliers (es más robusta).

  • Distancia de Manhattan:

\[\begin{gather*} d_{man}(p,q) = \sum_{i=1}^{n}\left | (p_{i}-q_{i}) \right | \end{gather*}\]

#distancias
library(factoextra)
res.dist <- get_dist(datos, stand = FALSE, method = "manhattan")
round(as.matrix(res.dist)[1:6, 1:6], 1)
##            Alabama Alaska Arizona Arkansas California Colorado
## Alabama        0.0    4.2     4.4      2.3        5.8      4.9
## Alaska         4.2    0.0     4.5      4.0        3.8      3.9
## Arizona        4.4    4.5     0.0      4.7        2.2      2.1
## Arkansas       2.3    4.0     4.7      0.0        6.2      4.4
## California     5.8    3.8     2.2      6.2        0.0      2.2
## Colorado       4.9    3.9     2.1      4.4        2.2      0.0
#matriz de distancias(grafico)
fviz_dist(res.dist)

#análisis de agrupamiento jerárquico
library(cluster)
library(purrr)
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")
ac <- function(x) {
 + agnes(datos, method = x)$ac} # AC describe la fuerza de la estructura de agrupamiento que se ha obtenido mediante el enlace.

map_dbl(m, ac) #Aplica la funcion ac a cada elemento del vector "m" y devuelve variables tipo double. 
##   average    single  complete      ward 
## 0.7379371 0.6276128 0.8531583 0.9346210

Ward es el valor más grande aquí con 0.9346210, por lo que ward.D2 se determina como método para la agrupación jerárquica.

# metodo de enlace
library(stats)
res.hc <- hclust(res.dist, method = "ward.D2" )
head(res.hc$merge,10)
##       [,1] [,2]
##  [1,]  -15  -29
##  [2,]  -13  -32
##  [3,]  -23  -49
##  [4,]  -14  -16
##  [5,]  -20  -31
##  [6,]  -36  -46
##  [7,]  -17  -26
##  [8,]  -37  -47
##  [9,]  -19    1
## [10,]  -35    4
#alturas
alturas <- data.frame(etapa = 1:49, distancia = res.hc$height)
tail(alturas) #etapa = 1:n-1
##    etapa distancia
## 44    44  5.155297
## 45    45  5.519602
## 46    46  7.132458
## 47    47 11.102206
## 48    48 12.541788
## 49    49 25.477600
library(ggplot2)
library(gghighlight)
ggplot(alturas) + aes(x = etapa, y = distancia)  +
  geom_point() + geom_line()  + 
  scale_x_continuous(breaks = seq(1, 49,by=2)) + 
  geom_vline(xintercept = 46, col = "red", lty = 3,size = 1.2) + 
  geom_vline(xintercept = 47, col = "blue", lty = 3,size = 1.2) + 
  theme_bw()

#dendograma
library(factoextra)
fviz_dend(res.hc, cex = 0.5, k=2, rect = TRUE,show_labels=F)

Parsimonia (menos es mas)

Tener menos cluster es mejor, manteniendo las propiedades de compactos entre si y heterogéneos entre ellos. A partir de cambios más drásticos en el proceso de agrupamiento.

El criterio del Gráfico de Silueta

library(cluster)
library(factoextra)
fviz_nbclust(x = datos, FUNcluster = pam, method = "silhouette", k.max = 15,
             diss = dist(datos, method = "manhattan"))

Observaciones al medoid más cercano

set.seed(123)
pam_clusters <- pam(x= datos, k=2, metric = "manhattan");pam_clusters
## Medoids:
##            ID     Murder    Assault   UrbanPop       Rape
## New Mexico 31  0.8292944  1.3708088  0.3081225  1.1603196
## Nebraska   27 -0.8008247 -0.8250772 -0.2445636 -0.5052109
## Clustering vector:
##        Alabama         Alaska        Arizona       Arkansas     California 
##              1              1              1              2              1 
##       Colorado    Connecticut       Delaware        Florida        Georgia 
##              1              2              2              1              1 
##         Hawaii          Idaho       Illinois        Indiana           Iowa 
##              2              2              1              2              2 
##         Kansas       Kentucky      Louisiana          Maine       Maryland 
##              2              2              1              2              1 
##  Massachusetts       Michigan      Minnesota    Mississippi       Missouri 
##              2              1              2              1              1 
##        Montana       Nebraska         Nevada  New Hampshire     New Jersey 
##              2              2              1              2              2 
##     New Mexico       New York North Carolina   North Dakota           Ohio 
##              1              1              1              2              2 
##       Oklahoma         Oregon   Pennsylvania   Rhode Island South Carolina 
##              2              2              2              2              1 
##   South Dakota      Tennessee          Texas           Utah        Vermont 
##              2              1              1              2              2 
##       Virginia     Washington  West Virginia      Wisconsin        Wyoming 
##              2              2              2              2              2 
## Objective function:
##    build     swap 
## 2.563358 2.360113 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
datos.j = cbind(datos, cluster= pam_clusters$cluster)

ACP y visualización de grupos

fviz_cluster(object = pam_clusters, data = datos, ellipse.type = "t",
             repel = TRUE) +
  theme_bw() +
  labs(title = "Resultados clustering PAM") +
  theme(legend.position = "none")

Como hay más de 2 variables, se están representando las 2 primeras componentes de un PCA. Se tiene que calcular el PCA y extraer las proyecciones almacenadas en el elemento x.

medoids <- prcomp(datos)$x

Se seleccionan únicamente las proyecciones de las observaciones que son medoids.

medoids <- medoids[rownames(pam_clusters$medoids), c("PC1", "PC2")]
medoids <- as.data.frame(medoids)
medoids
##                  PC1        PC2
## New Mexico -1.960124  0.1414131
## Nebraska    1.252916 -0.1920044

Se emplean los mismos nombres que en el objeto ggplot.

colnames(medoids) <- c("x", "y")
medoids
##                    x          y
## New Mexico -1.960124  0.1414131
## Nebraska    1.252916 -0.1920044

Creación del gráfico

fviz_cluster(object = pam_clusters, data = datos, ellipse.type = "t",
             repel = TRUE,palette = c("#a4c64f", "#F99D28")) +
  theme_bw() +
  # Se resaltan las observaciones que actúan como medoids
  geom_point(data = medoids, color = "firebrick", size = 2) +
  labs(title = "Resultados clustering PAM") +
  theme(legend.position = "none")

Verificar que el cluster esté bien

datos.j2 = datos.j %>% dplyr::select(cluster,Murder,Assault, UrbanPop, Rape)
colnames(datos.j2) = c('<span style="color:red">Cluster</span>', "Asesinato", "Agresión", "PropUrbano", "Violación")
DT::datatable(round(datos.j2,4), escape = FALSE)
datos.j$cluster = as.factor(datos.j$cluster) 
str(datos.j)
## 'data.frame':    50 obs. of  5 variables:
##  $ Murder  : num  1.2426 0.5079 0.0716 0.2323 0.2783 ...
##  $ Assault : num  0.783 1.107 1.479 0.231 1.263 ...
##  $ UrbanPop: num  -0.521 -1.212 0.999 -1.074 1.759 ...
##  $ Rape    : num  -0.00342 2.4842 1.04288 -0.18492 2.06782 ...
##  $ cluster : Factor w/ 2 levels "1","2": 1 1 1 2 1 1 2 2 1 1 ...

Diagrama de líneas de promedio por cluster

library(dplyr)
datos.j %>%   
  group_by(cluster) %>%
  summarise_all(list(mean)) -> medias
medias
## # A tibble: 2 × 5
##   cluster Murder Assault UrbanPop   Rape
##   <fct>    <dbl>   <dbl>    <dbl>  <dbl>
## 1 1        1.00    1.01     0.198  0.847
## 2 2       -0.670  -0.676   -0.132 -0.565
datos.j %>%  summarise_if(is.numeric,mean) %>%
  round(4) -> general
general
##   Murder Assault UrbanPop Rape
## 1      0       0        0    0
general <- cbind(cluster="general",general) #agregamos una nueva columna
general
##   cluster Murder Assault UrbanPop Rape
## 1 general      0       0        0    0
medias  <- as.data.frame(rbind(medias,general)) #juntamos las filas de ambas tablas
medias
##   cluster    Murder    Assault   UrbanPop       Rape
## 1       1  1.004934  1.0138274  0.1975853  0.8469650
## 2       2 -0.669956 -0.6758849 -0.1317235 -0.5646433
## 3 general  0.000000  0.0000000  0.0000000  0.0000000

Convirtiendo la data a formato tidy

Pasamos a filas todas las columnas menos la de “cluster”, las etiquetas las almacenamos en “variable” y los valores de cada etiqueta en “valor”.

library(tidyr)
gathered_datos.j <- pivot_longer(data=medias,
                                 -cluster,
                                 names_to="variable",
                                 values_to = "valor")
head(gathered_datos.j)
## # A tibble: 6 × 3
##   cluster variable  valor
##   <fct>   <chr>     <dbl>
## 1 1       Murder    1.00 
## 2 1       Assault   1.01 
## 3 1       UrbanPop  0.198
## 4 1       Rape      0.847
## 5 2       Murder   -0.670
## 6 2       Assault  -0.676

Manera de hallar patrones usando las VARIABLES ACTIVAS

ggplot(gathered_datos.j) + aes(x=variable,y=valor,color=cluster) + 
  geom_point() + 
  geom_line(aes(group = cluster)) +
  theme_bw() +
  theme(legend.position = "bottom",legend.title=element_blank()) +
  labs(title="Diagrama de líneas de Cluster por Variable",
       x="Variable",y="")

Conclusión: En EE.UU se presentan dos grupos de estados notablemente diferenciados en cuanto a índices de delincuencia y tasa de población, siendo uno de los grupos el que posee mayores valores que superan a la media de las variables que miden la taza de delincuencia (asalto, agresión y violación). Por lo tanto, en caso de que una persona este pensando mudarse a EE.UU se le pondría indicar los estados que conforman el grupo que tiene menor taza de delincuencia.